Writing A MATLAB/Octave Code To Solve A Linear System Using Different Implicit Solvers To Compute The Spectral Radius

                                                             #User Defined Function (UDF) for Jacobi Method
-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------

function [jac_iter jac_error x] = Jacobi_function (L,d,U,B)
  
  x = zeros(3,1);
  x_old = x;
  jac_error = 9e9;
  tol = 1e-04;
  jac_iter = 1;
  
  T_jac = inv(d)*(L+U);
  
  while (jac_error > tol)
    
    
    x = (inv(d)*B - T_jac*x_old);
    jac_error(jac_iter) = max(abs(x - x_old));
    x_old = x;
    
    
    jac_iter = jac_iter + 1;
    
    
  endwhile
  
endfunction

                                                          #User Defined Function (UDF) for Gauss seidel Implicit method
-----------------------------------------------------------------------------------------------------------------------------------------------------------------------


function [gs_iter gs_error x] = GS_function (L,d,U,B)
  
  x = zeros(3,1);
  x_old = x;
  gs_error = 9e9;
  tol = 1e-04;
  gs_iter = 1;
  
  T_gs = (inv(d+L))*U;
  
  while (gs_error > tol)
    
    
    x = (inv(d+L)*B - T_gs*x_old);
    gs_error(gs_iter) = max(abs(x - x_old));
    x_old = x;
    
    
    gs_iter = gs_iter + 1;
    
    
  endwhile
  
endfunction

                                                        #User Defined Function(UDF)for Successive Over Relaxation (SOR) Implicit Method  
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------

function [sor_iter sor_error x] = SOR_function (L,d,U,B,omega)
  
  x = zeros(3,1);
  x_old = x;
  sor_error = 9e9;
  tol = 1e-04;
  sor_iter = 1;
  
  T_sor = (inv(d+(omega*L)))*(omega*U-((1-omega)*d));
  
  while (sor_error > tol)
    
    
    x = T_sor*x_old + (inv(d+(omega*L))*B);
    sor_error(sor_iter) = max(abs(x - x_old));
    x_old = x;
    
    
    sor_iter = sor_iter + 1;
    
    
  endwhile
  
endfunction

                                                        #Main program code to solve the Linear System using Jacobi, Gauss-seidel & SOR methods
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

clear all
close all
clc

A = [ 5 1 2; -3 9 4; 1 2 -7];
B = [10; -14; -33];

D = diag(diag(A)); %Diagonal Matrix
U = triu(A)-D; %Upper Triangular Matrix
L = tril(A)-D; %Lower Triangular MAtrix


%Calculate Eigen Value manually
%characteristic equation of 3*3 matrix
% -lamba^3 + lamba*trace(A) + (trace(A)^2 - trace(A^2))*lambda + det(A)
% we get [-1*lambda + 7*lambda + 60*lambda - 402*lambda]

co_eff = [-1 7 60 -402];

Eigen_value = roots(co_eff)
ev_check = eig(A);
spectral_radius = max(abs(Eigen_value));
sprintf('Eigen value of matrixc A = %f, %f, %f', Eigen_value(1), Eigen_value(2), Eigen_value(3))


e = 0.2:0.2:1.2; %changinf the diagonal elements (Diagonal Magnification)


%Solve using Different Iterative Methods


     for i = 1:length (e)
         d = e(i)*D;
         sprintf('nn Diagonal Magnification Factor = %f', e(i));
         
         
         %Jacobi  Method
         
         T_jac = inv(d)*(L+U);
         jac_spectral(i) = max(abs(eig(T_jac)));
         [jac_iter(i) jac_error x] = Jacobi_function (L,d,U,B);
         sprintf('Jacobi Method') 
         sprintf('Spectral Radius = %f', jac_spectral(i))
         sprintf('Number of iterations = %f', jac_iter(i))
         sprintf('Solution x1 = %f, x2 = %f, x3 = %f', x(1), x(2), x(3))
         
         
         
         %Gauss-Sidel (GS) MEthod
         
         T_gs = (inv(d+L))*U;
         gs_spectral(i) = max(abs(eig(T_gs)));
         [gs_iter(i) gs_error x] = GS_function (L,d,U,B);
         sprintf('Gauss Sidel Method') 
         sprintf('Spectral Radius = %f', gs_spectral(i))
         sprintf('Number of iterations = %f', gs_iter(i))
         sprintf('Solution x1 = %f, x2 = %f, x3 = %f', x(1), x(2), x(3))
       
         
         
         %Successive Over Relaxation (SOR) Method
         
         omega = 2-(2/(1+sqrt(1-(jac_spectral(i)^2))));
         T_SOR = (inv(d+(omega*L)))*(omega*U-((1-omega)*d));
         SOR_spectral(i) = max(abs(eig(T_SOR)));
         [SOR_iter(i) SOR_error x] = SOR_function (L,d,U,B,omega);
         sprintf('Successive Over Relaxation (SOR) Method') 
         sprintf('Spectral Radius = %f', SOR_spectral(i))
         sprintf('Number of iterations = %f', SOR_iter(i))
         sprintf('Solution x1 = %f, x2 = %f, x3 = %f', x(1), x(2), x(3))
         
         
     endfor

     
figure(1)
subplot(2,1,1)
plot(e, jac_spectral,'r')
xlabel('Diagonal Magnification')
ylabel('Spectral Radius')
title('Jacobi Method')
subplot(2,1,2)
plot(e, jac_iter, 'b')
xlabel('Diagonal Magnification')
ylabel('Number of Iteration')


figure(2)
subplot(2,1,1)
plot(e, gs_spectral,'r')
xlabel('Diagonal Magnification')
ylabel('Spectral Radius')
title('Gauss Seidel (GS) Method')
subplot(2,1,2)
plot(e, gs_iter, 'b')
xlabel('Diagonal Magnification')
ylabel('Number of Iteration')


figure(3)
subplot(2,1,1)
plot(e, SOR_spectral,'r')
xlabel('Diagonal Magnification')
ylabel('Spectral Radius')
title('Successive Over Relaxation (SOR) Method')
subplot(2,1,2)
plot(e, SOR_iter, 'b')
xlabel('Diagonal Magnification')
ylabel('Number of Iteration')

Results & Observations:-   

• From the first plots we see a similar trend for Jacobi, Gauss siedel & Successive Over Relaxation i.e. the spectral radius decreases as the diagonal magnification factor is increasing. However, for SOR method the spectral radius drop was highest.
• Similar observation can be made for the second plot, as we increased the diagonal maginification factor the number of iterations started decreasing.