Jacobi to Gauss-Seidel

Question

I need to restructure my Jacobi code to Gauss-Seidel. I have tried many things but nothing has worked so far. Here is the current code that I have. What do I need to do to change this code into Gauss-Seidel?

function z=jacobi(A,B)
C=A; D=B; [N N]=size(A);
for i=1:N
    C(i,:)=-C(i,:)/A(i,i);
    C(i,i)=0;
    D(i)=D(i)/A(i,i);
end
x=zeros(N,1);
y=C*x+D;
m=0;
while(norm(x-y)>0.00000001 && m<100)
    x=y;
    y=C*x+D;
    m=m+1;
end
z=y;
end

score 0 · Accepted Answer · edited May 23 '17 at 12:20

The difference between a Jacobi solver and a Gauss-Seidel solver is that when you're solving for the solution of a variable x_i at the current iteration, you need to use the information from the previous variables (x_1, x_2, ..., x_{i-1}) as part of the solution for the current variable x_i. For Jacobi, you are simply using the previous iteration's solution to formulate the current solution. For Gauss-Seidel, for each variable that you solve for, you must use the solutions of the previous variables calculated from the current iteration as part of the solution for the variable you are focusing on.

As such, for your particular version of the code (though not optimal...), you simply need to add in a for loop where we solve for each variable one at a time, then keep feeding this information into the other variables. This has to be done in the while loop, because that's where you're doing the iterations.

Therefore:

function z=gaussseidel(A,B) %// Change the function name
C=A; D=B; [N N]=size(A);
for i=1:N
    C(i,:)=-C(i,:)/A(i,i);
    C(i,i)=0;
    D(i)=D(i)/A(i,i);
end
x=zeros(N,1);
y=C*x+D;
m=0;
while(norm(x-y)>0.00000001 && m<100)
    x=y;
    %// Change here
    for idx = 1 : N
        if idx == 1 %// Case where we are solving for the first variable
            y(idx) = C(idx,:)*x + D(idx);
        else %// Everything else
            y(idx) = C(idx,1:idx-1)*y(1:idx-1) + C(idx,idx:end)*x(idx:end) + D(idx);
        end
    end
    m=m+1;
end
z=y;
end

Look very closely at the code. If we are working on calculating the first variable, then we don't have any new information yet, and so we simply just update the solution using all of the previous information from the previous iteration. Once we pass this point, then we use the variables calculated from the current iteration from the first variable up to the point where we are calculating, then use the rest of the variables for calculating the solution.

If you want, there is another way to calculate Jacobi and Gauss-Seidel, and you can take a look at that post here: Jacobi iteration doesn't end. This code was originally written for Jacobi, and I provided the modification for making it into Gauss-Seidel.

Good luck!

Hey Ray, the code works pretty well but it takes too much time for make the calculus in comparisson of my jacobi code. It's not suppose to be that Seidel its faster tan Jacobi? — Pether Zavala, Nov 23 '14 at 18:55
@PetherZavala - No problem at all. Yes, it's supposed to **converge** much faster than Jacobi, but in terms of computation speed, what I have there isn't faster. The reason why is due to the `for` loop. That was the only way that I could think of given your code where I could implement the Gauss-Seidel changes. You need to use the variables from the current iteration to calculate the solutions for the other variables remaining. As such, there wasn't a way that I could think of doing it unless you used a loop. — rayryeng, Nov 23 '14 at 19:18

score 0 · Answer 2 · answered Nov 24 '14 at 06:35

Another version of the restructure code. By the way this code isn't quicker than the jacobi code but finds the solution in less iterations.

function [z,m]=seidel(A,B)
C=A; D=B; [N N]=size(A);
for i=1:N
 C(i,:)=-C(i,:)/A(i,i);
 C(i,i)=0;
 D(i)=D(i)/A(i,i);
end
x=zeros(N,1);
m=0;
while (m<100)
 t=x;
 m=m+1;
 for j=1:N
   y(j)=C(j,:)*x+D(j);
   x(j)=y(j);
 end
 if(norm(x-t)<0.00000001)
    break;
 end
end
z=y';
end

Jacobi to Gauss-Seidel

2 Answers2