0

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
  • 1
    It would help to introduce the shifted jacobi method for us with formulas instead of your code. – Thomas Dec 08 '14 at 07:23
  • Well all I know from class is that the shift is applied to the original matrix, and the normal Jacobi process is followed after that. There was not really any discussion about specifically using the shift, so I might be missing some steps there. Basically I'm solving the linear equation: (A+muI)x=b – jb68321 Dec 08 '14 at 18:34
  • I don't know if this solves the problem, but shouldn't x_ML=A\b since that is your original problem – Thomas Dec 08 '14 at 18:53
  • Sonystarmap: that line is just to compare a MATLAB answer/timing to the other methods. I could compare to either A or Ashift, depending on what my reasoning is. – jb68321 Dec 09 '14 at 20:41

0 Answers0