Significance of Spectral radius on the converge and stability of a linear system

Objective: To find the significance of Spectral radius on the converge and stability of a linear system

Theory: 

The solution of a system of linear equations can be found by implementing the iterative solvers namely, Jacobi, Gauss-Seidal and SOR. However the before starting to solve a system, we can predict the converge and stability of a system. This can only be stated if we know about the a parameter called "Spectral Radius".  Before coming to the definition of spectral radius, let's define what are eigenvalues are. 

So Mathematically, eigenvalues (`lamda`) of a system of linear equations represented in the form of matrix (say matrix A), are the roots of the equation $\det (A-\lambda I)=0$`det(A-lamdaI)=0` and spectral radius is the maximum of the absolute values of these eigenvalues `rho(A) = max(|lamda|)`

The matrix representation of the system of linear equations are as follows: 

`ubrace{[[5,1,2],[-3,9,4],[1,2,-7]]}_{[A]} ubrace{[[x_1],[x_2],[x_3]]}_{[X]} = ubrace{[[10],[-14],[33]]}_{[B]}`

The matrix [A] can be split and represented in as follows

[A] = [D] - [L] - [U]

Where, 

Diagonal Matrix,[D] = `[[5, 0, 0],[0, -9, 0], [0, 0,-7]]`

Lower triangular matrix, [L] = 

Upper triangular matrix,[U] =

Here we have split the matrix A into 3 parts so that we can actually implement the iterative methods in this case utilizing the matrix format. The 3 iterative methods in their matrix formats have been shown below: 

1. Jacobi method: x = T*xold + c 

Where iterative matrix, T =inv(D)*(L+U)

and coefficient matrix c = inv(D)*B

2. Gauss-Seidel Method: x = T*xold + c 

Where iterative matrix, T = inv(D-L)*U

and coefficient matrix c = inv(D-L)*B

3. SOR Method: x = T*xold + c 

Where iterative matrix, T = inv((1/w)*D-L)*[((1/w)-1)*D + U]

and coefficient matrix c = inv((1/w)*D-L)*B

Then we would manipulate the diagonal matrix, D with a scalar magnifier to observe the effect of spectral radius on the converge and performance of the iterative methods. 

Code: 

1. Main Program

clear all
close all
clc
    
%% Part-1: INPUTS

    %Co-efficient Matrix 
    A=[5 1 2;-3 9 4;1 2 -7];
    %RHS or Residual Matrix
    B=[10;-14; 33];
 
    % Matrix Decomposition of Coefficient Matrix, A
    % Here, A=D-U-L
    U = -triu(A,1); 
    D = diag(diag(A));
    L = -(tril(A)-D);
    Dinitial =D;
    multiplier = linspace(0.3,5,20);
    tol =1e-4;
    
for k =1:length(multiplier)
    D = multiplier(k)*Dinitial; 
    
    %Modified Co-efficient Matrix 
    A=D-U-L;
    %% Part-2: Iterative Matrix
    %Jacobi Iterative Matrix
    Tjac = inv(D)*(L+U);
    %Gauss Seidal Iterative Matrix
    Tgs = inv(D-L)*U;
    %SOR Iterative Matrix
    omega =0.92;
    Tsor = inv(D-(omega*L))*[((1-omega)*D) + (omega*U)];

    %% Part-3: Spectral radius
    rho_jac(k) = spectral_rad(Tjac);
    rho_gs(k) = spectral_rad(Tgs);
    rho_sor(k)= spectral_rad(Tsor);

    %% Part-4: Computing x
    [iter_jac(k) x_j] = jacobi(A,B,U,D,L, multiplier, tol);
    x_jac(1,k) = x_j(1);
    x_jac(2,k) = x_j(2);
    x_jac(3,k) = x_j(3);

    [iter_gs(k) x_G] = gs(A,B,U,D,L, multiplier, tol);
    x_gs(1,k) = x_G(1);
    x_gs(2,k) = x_G(2);
    x_gs(3,k) = x_G(3);

    [iter_sor(k) x_S] = sor(A,B,U,D,L, multiplier, tol, omega);
    x_sor(1,k) = x_S(1);
    x_sor(2,k) = x_S(2);
    x_sor(3,k) = x_S(3);

    
end   
    %% plots
    figure(1)
    plot(multiplier, rho_jac,'o-',multiplier, rho_gs,'o-',multiplier, rho_sor,'o-')
    xlabel("Diagonal Multiplier")
    ylabel("Spectral Radius")
    legend('Jacobi', 'Gauss Seidal','SOR')
    title("Effect on Spectral Radius due to magnification of diagonal matrix")
    grid on

    figure(2)
    plot(rho_jac, iter_jac,'o-')
    xlabel("Spectral Radius")
    ylabel("No of Iteration")
    grid on
    title("Jacobi")
    
    figure(3)
    plot(rho_gs, iter_gs,'o-')
    xlabel("Spectral Radius")
    ylabel("No of Iteration")
    grid on
    title("Gauss Seidal")
    
    figure(4)
    plot(rho_sor, iter_sor,'o-')
    xlabel("Spectral Radius")
    ylabel("No of Iteration")
    grid on
    title("SOR")
    
   
    

 2. User-Defined Functions

a. spectral_rad(T)

function output = spectral_rad(T)
%% spectral_radius
output = max(abs(roots(charpoly(T))));
end

b.  jacobi()

function [iter x] = jacobi(A,B,U,D,L, multiplier, tol)

n=length(B); 
xold=ones(n,1);
err = 9e9;
iter=0;
while(err>tol)
    T =inv(D)*(L+U);
    c = inv(D)*B;
    x = T*xold + c;
    err = max(abs(x-xold));
    xold=x;
    iter =iter +1;
end

end

c. gs()

function [iter x] = gs(A,B,U,D,L, multiplier, tol)

n=length(B); 
xold=ones(n,1);
err = 9e9;
iter=0;

while(err>tol)
    T = inv(D-L)*U;
    c = inv(D-L)*B;
    x = T*xold +c;
    err = max(abs(x-xold));
    xold=x;
    iter =iter +1;
end

end

d.  sor()

function [iter x] = sor(A,B,U,D,L, multiplier, tol, w)

n=length(B); 
xold=ones(n,1);
err = 9e9;
iter=0;
while(err>tol)
    %T = inv(D-w*L)*[(1-w)*D + w*U];
    T = inv((1/w)*D-L)*[((1/w)-1)*D + U];
    %c = inv(D-w*L)*w*B;
    c = inv((1/w)*D-L)*B;
    x = T*xold  + c;
    err = max(abs(x-xold));
    xold=x;
    iter =iter +1;
end

 Algorithm: 

1. At first we have defined the matrix A and calculated the upper, lower and diagonal matrix out of it and also defined a range of 20 diagonal multiplier from 0.3 to 5. 

2.  Then for each of the diagonal multiplier, we calculated the iterative matrix for each of the iterative methods. Calculated the spectral radius using the function called spectral_rad() and then computed the values of the solution matrix, X = [x1;x2;x3]. 

3. At the end, we have plotted the graphs to see the relation between the spectral radius and the scalar magnifier and also the effect of spectral radius on the converge on the solvers.

Graphs:

Results and Discussion:

1. So from the plots, it can be easily concluded that with the increase in the diagonal magnification factor the spectral radius decreases and with the increase in spectral radius the linear system of equations diverges.  

2. If we carefully observe the last 3 plots then we can state after divergence increase to a great extent after spectral radius 1 irrespective of the iterative methods we are using. 

3.  We have checked the solutions of the system with spectral radius,1 and it has been observed that the system doesn't converge at such scenarios. So we can conclude that the convergence of a system depends in spectral radius and it is ensured only when the spectral radius is less than 1.

For jacobi method we can see the system doesn't converge for spectral radius greater than 1

Similar is the case for other methods also.

and at spectral radius =1 the divergence is the maximum