Jacobi Iteration with Shift

98 Views Asked by At

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