Understanding linear systems- Program to determine the Eigen Value and Spectral Radius of the given system of equations. Also solve them using iterative techniques available.

The Gauss-Seidel method is a technique for solving the 'n' equations of the linear system of equation Ax=B  one at a time in sequence, and uses previously computed results as soon as they are available,

There are two important characteristics of the Gauss-Seidel method should be noted.

Firstly, the computations appear to be serial. Since each component of the new iterate depends upon all previously computed components, the updates cannot be done simultaneously as in the Jacobi Method.

Secondly, the new iterate `x^((k+1))` depends upon the order in which the equations are examined. If this ordering is changed, the components of the new iterates (and not just their order) will also change.

In terms of matrices, the definition of the Gauss-Seidel method can be expressed as

where the matrices D, -L, and -U represent the diagonal, strictly lower triangular, and strictly upper triangular parts of A, respectively.

The Gauss-Seidel method is applicable to strictly diagonally dominant, or symmetric positive definite matrices  A.

3. Successive Overrelaxation Method:

The successive overrelaxation method (SOR) is a method of solving a linear system of equation Ax=B  derived by extrapolating the Gauss-seidel method. This extrapolation takes the form of a weighted average between the previous iterate and the computed Gauss-Seidel iterate successively for each component,

In matrix terms, the SOR algorithm can be written as

where the matrices D, -L and -U represent the diagonal, strictly lower-triangular, and strictly upper-triangular parts of  A, respectively.

If  , the SOR method simplifies to the gauss-seidel method A theorem due to Kahan shows that SOR fails to converge if  is outside the interval (0,2).

So the iteration matrix for the above methods we have:

_______For Jacobi Method

_______For Gauss-seidel Method

_______For SOR-GS Method

Here, new guess values will be calculated as,

where, T- Iteration matrix of different methods

Main Program file:

%% By Gaurav Kharwade%%

clear all
close all
clc

% Understanding of linear equation-
% 5x1+x2+2x3=10
% -3x1+9x2+4x3=-14
% x1+2x2-7x3=33

% RHS   
B=[10;-14;33];

% Decomposition of matrix 
L=[0 0 0;-3 0 0;1 2 0];
U=[0 1 2;0 0 4;0 0 0];
D=[5 0 0;0 9 0;0 0 -7];

% co-efficient matrix 
A=L+U+D;

% Calculating for x
x_exact=inv(A)*B

disp('Type 1 to select Jacobi method')
disp('Type 2 to select GS method')
disp('Type 3 to select SOR method')
solving_method= input('Select solving method= ');

if solving_method==1

% For magnifying diagonal terms
 m=linspace(0.4602,5,10);
 
 Dold=D;

 [jac_iter,x_jac,error_jac,Computational_time,SpectralR_jac]=jac(A,B,D,Dold,L,U,m)

 
figure(1)
subplot(2,1,1)
plot(m,SpectralR_jac,'b')
xlabel('Diagonal Magnification')
ylabel('Spectral Radius')
title('Effect of Diagonal Magnification on Spectral Radius')
hold on
subplot(2,1,2)
plot(SpectralR_jac,jac_iter,'r','linewidth',2)
ylabel('Nos. of iteration')
xlabel('Spectral Radius')
title(sprintf('No. of iterations vs Spectral Radius graph'))

endif


if solving_method==2

% For magnifying diagonal terms
m=linspace(0.55,5,10);

Dold=D; 
 
[gs_iter,x_gs,SpectralR_gs,error_gs]=gs(A,B,Dold,D,L,U,m)

figure(1)
subplot(2,1,1)
plot(m,SpectralR_gs,'b')
xlabel('Diagonal Magnification')
ylabel('Spectral Radius')
title('Effect of Diagonal Magnification on Spectral Radius')
hold on
subplot(2,1,2)
plot(SpectralR_gs,gs_iter,'r','linewidth',2)
ylabel('Nos. of iteration')
xlabel('Spectral Radius')
title(sprintf('No. of iterations vs Spectral Radius graph'))

endif

if solving_method==3

% For magnifying diagonal terms
m=linspace(0.45,5,10);

disp('Enter value for relaxation factor suitably to get spectral radius to 1.1, i.e. 0.9')
omega=input('Enter value for relaxation factor=')

Dold=D;   

[sor_iter,x_sor,SpectralR_sor,error_sor]=sor(A,B,Dold,D,L,U,m,omega)
figure(1)
subplot(2,1,1)
plot(m,SpectralR_sor,'b')
xlabel('Diagonal Magnification')
ylabel('Spectral Radius')
title('Effect of Diagonal Magnification on Spectral Radius')
hold on
subplot(2,1,2)
plot(SpectralR_sor,sor_iter,'r','linewidth',2)
xlabel('Spectral Radius')
ylabel('Nos. of iteration')
title(sprintf('No. of iterations vs Spectral Radius graph'))

endif

User_defined Function Program file- JACOBI METHOD

function [iter,x_jac,error_jac,Computational_time,SpectralR_jac]=jac(A,B,Dold,D,L,U,m)

tic
for j=1:length(m)
  
  D=m(j)*Dold; 
  
  n=length(B); 
  x_old=zeros(n,1);
  %xnew_jac=zeros(n,1);
  error_jac=9e9;
  tol=1e-3;

  Tjac=inv(D)*(L+U); % for Jacobi  
 
 iter(j)=1;
while(error_jac>tol) 
   
  xnew_jac=(-(inv(D)*(L+U))*x_old)+(inv(D)*B);
  error_jac=max(abs(xnew_jac-x_old));
  x_old=xnew_jac;
  iter(j)=iter(j)+1;
endwhile
 x_jac=xnew_jac;
 
eigen_value=eig(Tjac);
SpectralR_jac(j)=max(real(abs(eigen_value)));
 

endfor
Computational_time=toc;


endfunction

User_defined Function Program file- GAUSS-SEIDEL METHOD

unction [gs_iter,x_gs,SpectralR_gs,error_gs]=gs(A,B,Dold,D,L,U,m)
 
  for j=1:length(m)
   D=m(j)*Dold;
 
  n=length(B); 
  x_old=zeros(n,1);
  %xnew_gs=zeros(n,1);
  error_gs=9e9;
  tol=1e-5;
  
  
 Tgs=(L+D)^-1*(-U); % for GS

 gs_iter(j)=1;
while(error_gs>tol) 
   
  xnew_gs= ((-inv(D+L)*(U))*x_old)+(inv(D+L)*B);
  error_gs=max(abs(xnew_gs-x_old));
  x_old=xnew_gs;
  gs_iter(j)=gs_iter(j)+1;
  
  
endwhile
  x_gs=xnew_gs;
 Computational_time=toc;
 
 
eigen_value=eig(Tgs);
SpectralR_gs(j)=max(real(abs(eigen_value)));
 
 endfor
endfunction

User_defined Function Program file- SOR METHOD

function [sor_iter,x_sor,SpectralR_sor,error_sor]=sor(A,B,Dold,D,L,U,m,omega)

for j=1:length(m)
  D=m(j)*Dold;
  
  n=length(B); 
  x_old=zeros(n,1);
  %xnew_sor=zeros(n,1);
  error_sor=1e9;
  tol=1e-5;
 

Tsor=(L+(1/omega)*D)^-1*(((1/omega)*D)-D-U); %for sor 
  
 sor_iter(j)=1;
while(error_sor>tol) 
   
  xnew_sor= -(inv(D+omega*L)*((1-omega)*L+U)*x_old)+(inv(D+omega*L)*B);
  error_sor=max(abs(xnew_sor-x_old));
  x_old=xnew_sor;
  sor_iter(j)=sor_iter(j)+1;
  
  
endwhile
  x_sor=xnew_sor;
 Computational_time=toc;
 
 eigen_value=eig(Tsor);
 SpectralR_sor(j)=max(abs(eigen_value));
 
 endfor
endfunction

Following are results obtained:

PLOTS:

1. JACOBI METHOD

Values of Nos od iteration and Spectral radius is as below:

2. GAUSS-SEIDEL METHOD

Values of Nos od iteration and Spectral radius is as below:

3. SOR METHOD

The idea of over-relaxation is to choose the factor ω in such a way that the spectral radius ρ (P) is made as small as possible, thus giving the optimum rate of convergence.

Values of Nos od iteration and Spectral radius is as below:

From the graph, we can easily say that if the value for the spectral radius is more than 1(`rho>1`) numbers of iteration to converge will be much higher as it takes more updated guess values for convergence. 

As spectral radius coming down from 1.1, there is a drastic decrease in a number of iterations for convergence.

Conclusion:

The iterative scheme   is convergent if and only if every
the eigenvalue of [T] satisfies  |λ|<1, that is, the spectral radius .