The question is to solve a linear system using Jacobi iterations with a shift of mu = 5. My code converges very quickly, but it does not yield the results that MATLAB gives with the backslash operator. Perhaps I'm implementing the shift part incorrectly. Do you apply it to the "b" side also or take it out later? Thanks! Note that the A matrix is created using some specific criteria, and I believe that part is all correct.
% Centered Difference Scheme - Create "A" Matrix
% Boundary Conditions
x_0 = 0;
x_n = 1;
u_0 = 2;
u_n = -1;
n = 200; % Number of discrete points calculated
h = (x_n-x_0)/(n-1); % Distance between points
X = (x_0:h:x_n)'; % Initializing x vector
% Create the tri-diagonal coefficient matrix "A"
A = 1/h^2*(diag(-2*ones(1,n))+diag(ones(1,n-1),1)+diag(ones(1,n-1),-1));
% Create "b" matrix
f = 1/2-X(1:n);
b = f;
b(1) = f(1)-u_0/(h^2);
b(n) = f(n)-u_n/(h^2);
% Incorporating Given Shift
mu=5;
I=eye(size(A));
A_shift=A+mu*I;
% Jacobi Iteration Method
[m,n2]=size(A_shift);
S=zeros(m,n2);
for i=1:n2;
S(i,i)=A_shift(i,i);
end
T=S-A_shift;
x_J=zeros(size(b));
xnew=zeros(size(b));
for k=1:1000
for j=1:n2
xnew(j)=b(j);
for i=1:m
xnew(j)=xnew(j)+T(j,i)*x_J(i);
end
xnew(j)=xnew(j)/S(j,j);
end
x_J(:,k)=xnew;
end;
x_J % Jacobi Solution
% Compare to Matlab
x_ML = A_shift\b;
Update: I found that my actual algorithm was the problem, not the shifting. I re-wrote the entire thing and had better results! Here is the end result:
i_max = 1000000;
D=diag(diag(A_shift));
D_inv=inv(D);
E=A_shift-D;
x_J=rand((length(A_shift)),1);
T=-D_inv*E;
C=D_inv*b;
tol=.100000;
k_J=0;
d=1;
while d>tol
if k_J<i_max
k_J=1+k_J;
x_J=T*x_J+C;
d = max(abs(A_shift*x_J-b));
else
Time_J=toc;
tic
x_ML = A_shift\b;
Time_ML=toc;
return
end
end