Chebyshev smoother, I cannot get this to run, I dont understand why.
2 views (last 30 days)
Show older comments
dtVx = rho*resolx.^2./(4*eta_nbv*(1.+eta_b)*4.1); %dtVx is a very large vector
A = (4*eta_nbv*(1.+eta_b)*4.1); %Here I am using Ax=b for the above dtVx vector
A = diag(A);
[V,D] = eigs(A);
lam1 = min(diag(D))
lam2 = max(diag(D))
%END OF ARNOLDI
S = sparse(A);
S = S \ A;
b = rho*resolx.^2; %This b is from the dtVx vector above
r = S.*(b - A.*dtVx); %redidual
thet = 0.5*(lam2+lam1); %Theta value
delt = 0.5*(lam2-lam1); %Delta value
sig = thet/delt; %sigma value
rho1 = 1/sig;
dk = r/thet;
xk = 0;
nbv = 44229;
for k = 1:nbv
xk = xk + dk;
r = r - S*A*dk;
rhok = 1/(2*sig-rhok);
dk = (rhok*rho1*dk + 2*rhok*r)/(delt+1);
rho1 = rhok
end
xk = xk + dk;
0 Comments
Accepted Answer
Paul
on 22 Apr 2023
Hi Nathan,
Can't run the code because not all of the information needed to run it is provided.
dtVx = rho*resolx.^2./(4*eta_nbv*(1.+eta_b)*4.1); %dtVx is a very large vector
Comparing the code to the Algorithm, it looks like there may be a few issues.
A = (4*eta_nbv*(1.+eta_b)*4.1); %Here I am using Ax=b for the above dtVx vector
A = diag(A);
[V,D] = eigs(A);
lam1 = min(diag(D))
lam2 = max(diag(D))
%END OF ARNOLDI
S = sparse(A);
This line should result in S being the identity matrix, to within numerical precision.
S = S \ A;
b = rho*resolx.^2; %This b is from the dtVx vector above
This line using element-wise multiplication and not matrix multiplication. Is thar correct? Also, this value of r is used and updated in the loop, even though the Algorithm doesn't initialize it as r_1. Is that a typo in the Algorithm?
r = S.*(b - A.*dtVx); %redidual
thet = 0.5*(lam2+lam1); %Theta value
delt = 0.5*(lam2-lam1); %Delta value
sig = thet/delt; %sigma value
rho1 = 1/sig;
dk = r/thet;
xk = 0;
nbv = 44229;
for k = 1:nbv
xk = xk + dk;
r = r - S*A*dk;
rhok = 1/(2*sig-rhok);
The next line doesn't look correct. Accoding to the algorithm it should be
dk = (rhok*rho1*dk + 2*rhok*r)/(delt+1);
% dk = rhok*rho1*dk + 2*rhok*r/delt % should be this?
rho1 = rhok;
end
xk = xk + dk;
0 Comments
More Answers (0)
See Also
Categories
Find more on Sparse Matrices in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!