Jacobi/Gauss Seidel Methods in Matlab

Question

I need to implement the Jacobi and Guass Seidel, methods in Matlab.

I found this link which has code that produces correct results (on the one sample I tried) for each.

See post 3.

My problem is that the implementation is different to that described here, and here

I'm happy enough using someone else's code (indeed I prefer something tried and tested) but I want to understand exactly how it works.

Can someone point me to a description of the implementations used in the post?

Alternately are there other implementations of these algortihms that might be considered benchmarks in Matlab?

Answer 1

I would like to show you how the code on Seidel works, hopefully you can do the same analysis on Jacobi yourself.

Q=tril(A); % Q == L
r=b-A*x;
dx=Q\r; 

This part mathematically means x(:,k+1) = inv(L) * (b - A*x(:,k)) = inv(L) * (b - L*x(:,k) - U*x(:,k));

while in the wikipedia page you provides, it requires inv(L) * (b - U*x(:,k));

but they are equivalent since inv(L) * (b - L*x(:,k) - U*x(:,k)) = inv(L) * (b - U*x(:,k)) - x(:,k); so if you follow the formula in wikipedia, the iteration update should be: x(:,k+1)=(dx + x(:,k)); while it is the same in the code you provided :x(:,k+1) = x(:,k) + lambda * dx;

Please note that lambda is a relaxation coefficient, mainly functions on the convergence speed adjuestment. You can set to 1 in the code, which makes it exactly the same as the formula in wikipeida.