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;
