%
%Feb 4: Examples for solving Ax=b by Jacobi, Gauss-Seidel and SOR methods
%
%Initialize the matrix A, x0 and b
%
>> A=[4 -1 1;-1 5 -1;1 -1 7];
>> b=[1;2;3];
>> x0=zeros(3,1);
%
% Partitioning A as: A=D-L-U
%
>> D=diag(diag(A));
>> L=[0 0 0;1 0 0;-1 1 0];
>> U=[0 1 -1;0 0 1;0 0 0];
%
% Check if D-L-U=A:
>> D-L-U
%
%Jacobi Method:
%
>> x1=inv(D)*(L+U)*x0+inv(D)*b;
>> x2=inv(D)*(L+U)*x1+inv(D)*b;
%
%Gauss-Seidel Method:
%
>> x1=inv(D-L)*U*x0+inv(D-L)*b;
>> x2=inv(D-L)*U*x1+inv(D-L)*b;
%
%SOR Method with w=1.2
%
>> w=1.2;
>> x1=inv(D-w*L)*((1-w)*D+w*U)*x0+inv(D-w*L)*w*b;
>> x2=inv(D-w*L)*((1-w)*D+w*U)*x1+inv(D-w*L)*w*b;
%
%Compute the spectral radius for each T
%
>> Tjac=inv(D)*(L+U);
>> abs(eig(Tjac))
%ans =
%
%    0.3886
%    0.2260
%    0.1627
%
%Spectral Radius of Tjac is 0.386<1 so Jacobi Method converges for this example.
%
>> Tgs=inv(D-L)*U;
>> abs(eig(Tgs))
%ans =
%
%         0
%    0.0845
%    0.0845
%
%Spectral Radius of Tgs is 0.0845<1 so Gauss-Seidel Method converges for this example.
%
>> Tsor=inv(D-1.2*L)*((1-1.2)*D+1.2*U);
>> abs(eig(Tsor))
%ans =
%
%    0.2362
%    0.2362
%    0.1434
%
%Spectral Radius of Tsor is 0.2353<1 so SOR Method with w=1.2 converges for this example.
%