Gauss-Seidel Method in MATLAB

Question

I am trying to implement the Gauss-Seidel method in MATLAB. But there are two major mistakes in my code, and I could not fix them:

1. My code converges very well on small matrices, but it never converges on large matrices.

2. The code makes redundant iterations. How can I prevent from redundant iterations?

Gauss-Seidel Method on wikipedia.

N=5;
A=rand(N,N);
b=rand(N,1);
x = zeros(N,1);
sum = 0;
xold = x;
tic
for n_iter=1:1000
    for i = 1:N
        for j = 1:N
            if (j ~= i)
                sum = sum + (A(i,j)/A(i,i)) * xold(j);
            else
                continue;
            end
        end
        x(i) = -sum + b(i)/A(i,i);
        sum = 0;
    end
    if(abs(x(i)-xold(j))<0.001)
        break;
    end
    xold = x;
end
gs_time=toc;
prompt1='Gauss-Seidel Method Time';
prompt2='x Matrix';
disp(prompt2);
disp(x);
disp(prompt1);
disp(gs_time);

Answer 1

First off, a generality. The Gauß-Seidel and Jacobi methods only apply to diagonally dominant matrices, not generic random ones. So to get correct test examples, you need to actually constructively ensure that condition, for instance via

A = rand(N,N)+N*eye(N)

or similar.

Else the method will diverge towards infinity in some or all components.

---

Now to some other strangeness in your implementation. What does

if(abs(x(i)-xold(j))<0.001)

mean? Note that this instruction is outside the loops where i and j are the iteration variables, so potentially, the index values are undefined. By inertia they will accidentally both have the value N, so this criterion makes at least a little sense.

What you want to test is some norm of the difference of the vectors as a whole, thus using sum(abs(x-xold))/N or max(abs(x-xold)). On the right side you might want to multiply with the same norm construction applied to x so that the test is for the relative error, taking the scale of the problem into account.

---

By the instructions in the given code, you are implementing the Jacobi iteration, computing all the updates first and then advancing the iteration vector. For the Gauß-Seidel variant you would need to replace the single components in-place, so that newly computed values are immediately used.

---

Also, you could shorten/simplify the inner loop

xold = x;
for i = 1:N
    sum = b(i);
    for j = 1:N
        if (j ~= i)
            sum = sum - A(i,j) * x(j);
        end
    end
    x(i) = sum/A(i,i);
end
err = norm(x-xold)

or even shorter using the language features of matlab

xold = x
for i = 1:N
    J = [1:(i-1) (i+1):N];
    x(i) = ( b(i) - A(i,J)*x(J) )/A(i,i);
end   
err = norm(x-xold)