Jacobi Method example

I wrote this code a while back. Apparently, this post is getting a few hits. I refined it a bit. Also wanted to make sure it was looked right.

What is the Jacobi method? It's an iterative way of solving a system of equations. Ax = b. This is learned in Numerical Analysis classes.

Generally, you decompose a matrix into its diagonal and lower and upper triangles. Such that,

A = D + L + U

x0 is the initial guess.

Now you want to solve for x. x_approx = inverse of D * (b - (L+U)x0)

Get x_approx from this formula, then plug it into x0 and do it all over again until you reach the desired error.

% Code is for matlab / octave
format long e;
A = [ 2 -1; -1 2];
b = [ 2; 1];
k = 20;
xi = [0;0]; % initial guess

%we want
%x' = Dinv(b - (L+U)x)

% just take the diagonal of A
D=[A(1,1) 0;0 A(2,2)]; 

Di=inv(D);

%lower and upper triangles of A
L=[0 0; -A(2,1) 0];
U=[0 -A(1, 2); 0 0];

% this part doesn't change, so calc now
C=Di*b;
% this is the main term of the algo
M=Di*(L+U);

x = [5/3; 4/3]; % solution

abs_err = [ 100 ; 100];
abs_err = abs_err - x;
i=0;

%stop until error is <= 1e-6
while abs_err(1) >= 1e-6 && abs_err(2) >= 1e-6
 xi=M*xi+C
 abs_err = abs(xi - x);
 fprintf('i = %f abs-err = %f %f \n',i,abs_err(1),abs_err(2));
 i=i+1;
end