%% MATLAB Iteration Demo
% This file demonstrates Jacobi iteration and convergence.

%% Problem setup

A=[4,-1,0,1,0;-1,4,-1,0,1;0,-1,4,-1,0;1,0,-1,4,-1;0,1,0,-1,4];
b=[100 100 100 100 100]';

toler=0.001;          %Iterate until relative error within tolerance
maxiter=100;        %Max Iterations (before giving up)
n=size(b);

%% Jacobi Iteration

display('Jacobi Iteration');

xk=zeros(n);
x0=ones(n);
k=0;

%Loop until converged, or failed
while ((max(abs((xk-x0)./x0))>toler) & (k<maxiter)) 
    x0=xk;              %Retain prior solution
    k=k+1;              %Count iterations
    R=b-A*x0;           %Note that we have used Matrix Multiply!
    xk=x0+R./diag(A);   %Note that only the diagonal is used here.
    display(['k=',int2str(k),'|',num2str(xk','  %1.16e')]);
end;

R

%% Gauss-Seidel Iteration

display('Gauss-Seidel Iteration');

xk=zeros(size(b));
x0=ones(size(b));
k=0;

%Loop until converged, or failed
while ((max(abs((xk-x0)./x0))>toler) & (k<maxiter)) 
    x0=xk;              %Retain prior solution
    k=k+1;              %Count iterations
    for i=1:n
        R=b(i)-A(i,:)*xk;
        xk(i)=xk(i)+R/A(i,i);
    end;
    display(['k=',int2str(k),'|',num2str(xk','  %1.16e')]);
end;

R

%% SOR 

display('Successive-Over-Relaxation');

xk=zeros(size(b));
x0=ones(size(b));
k=0;
omega=1.12;

%Loop until converged, or failed
while ((max(abs((xk-x0)./x0))>toler) & (k<maxiter)) 
    x0=xk;              %Retain prior solution
    k=k+1;              %Count iterations
    for i=1:n
        R=b(i)-A(i,:)*xk;
        xk(i)=xk(i)+omega*R/A(i,i);
    end;
    display(['k=',int2str(k),'|',num2str(xk','  %1.16e')]);
end;

R