Reputation: 2150
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?
Upvotes: 0
Views: 3309
Reputation: 12689
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.
Upvotes: 1