Understanding Linear Systems

Code:

clear all 
close all 
clc

A = [5 1 2 ; -3 9 4; 1 2 -7];

B=[10;-14;33];
%magnification factor m
m=[0.5 1 2 3 4 5 6];


%calculating the eigen value A-L*I=0;roots of this eqn will give the igen values
l=[0 0 0 ;-3 0 0;1 2 0];  %lower diag matrix
u=[0 1 2;0 0 4;0 0 0];    %upper diad marix
d=[5 0 0 ;0 9 0;0 0 -7];  %diag matrix
I=[1 0 0;0 1 0;0 0 1];    %identity matrix

syms f(L);
solver=input('enter the solver 1)jacobi 2)gauss sidel 3)sor');
if solver==1
   print(45);
   for i=1:length(m) 
    d=m(i)*d;
    sum=(l+u);
    t_jac=(inv(d))*(sum);
    f(L)=det(t_jac-L*I);
    c=double(coeffs(f(L),'all'))
    eign_val=roots(c);
    spectral_rad(i)=max(abs(eign_val));
      if spectral_rad(i)<1 
 [jac_iter(i),xnew]=methord(t_jac,d,B);
      else
          jac_iter(i)=nan;
          txt=sprintf('the jacobi solver could not converge for the spectral radius equal to: %d',spectral_rad(i))
          disp(txt);
      end
      end
     subplot(2,1,1);
     plot(m,spectral_rad);
     xlabel('diagonal magnification');
     ylabel('spectral radius ');
    
     subplot(2,1,2);
    loglog(spectral_rad,jac_iter);  
     xlabel('spectral radius');
     ylabel('number of iterations');
end



if solver==2
   
   for i=1:length(m) 
    d=m(i)*d;
    sum=(u);
    t_gs=(inv(d+l))*(sum);
    f(L)=det(t_gs-L*I)
    c=double(coeffs(f(L),'all'))
    eign_val=roots(c)
    spectral_rad(i)=max(abs(eign_val))
    if spectral_rad(i)<1  
 [gs_iter(i),xnew]=methord(t_gs,d,B);
    else 
        gs_iter(i)=nan;
         txt=sprintf('the gauss sidel solver could not converge for the spectral radius equal to: %d',spectral_rad(i))
          disp(txt);
    end
    
   end
     subplot(2,1,1);
     plot(m,spectral_rad);
      xlabel('diagonal magnification');
     ylabel('spectral radius ');
     
     subplot(2,1,2);
    loglog(spectral_rad,gs_iter);  
     xlabel('spectral radius');
     ylabel('number of iterations');
end

w=1.2

if solver==3
   
   for i=1:length(m) 
    d=m(i)*d;
    sum=(1-w)*d-w*u;
    t_sor=(inv(d+l*w))*(sum);
    f(L)=det(t_sor-L*I)
    c=double(coeffs(f(L),'all'))
    eign_val=roots(c)
    spectral_rad(i)=max(abs(eign_val))
 if spectral_rad(i)<1  
 [sor_iter(i),xnew]=methord(t_sor,d,B*w);
 else
     sor_iter(i)=nan;
     txt=sprintf('the sor solver could not converge for the spectral radius equal to: %d',spectral_rad(i))
     disp(txt);
 end
   end
     subplot(2,1,1);
     plot(m,spectral_rad);
     xlabel('diagonal magnification');
     ylabel('spectral radius ');
     subplot(2,1,2);
    loglog(spectral_rad,sor_iter);  
     xlabel('spectral radius');
     ylabel('number of iterations');
end

Function Used:

The function used in the above-mentioned program is as follows

function  [iter,xnew]=methord(t,d,B)
   tol=10e-3;
   x=[0;0;0]
   xold=x;
   error=9e9;
   iter=1
   
   while (error>tol)
    
    x=(t)*xold+(inv(d))*B ;  
    error=max(abs(xold-x))
    iter=iter+1
    xold=x;
   end
xnew=x;
end

GRAPHS:

The following graphs were obtained

1. Jacobi Method

2. Gauss-Seidel Method

The Gauss-Seidel solver could not converge for the spectral radius equal to : 1.381472e+00

3. SOR Method