Stability Analysis for different Iterative Solvers on a System of Linear equations in MATLAB

Stability Analysis for different Iterative Solvers on a System of Linear equations in MATLAB

Problem Specification:-

Solving the system of linear equations in MATLAB using different iterative solvers and showing:

1. The impact of the spectral radius on the convergence of the solution.
2. The effect of diagonal magnification on the spectral radius for different iterative solvers.

Given:-

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

X=`[x_1,x_2,x_3]`

B= [10; -14; 33]

Increase the spectral radius to 1.1 by changing the diagonal terms suitably.

Determine:-

The eigenvalue of the iteration matrix.

The spectral radius of the iteration matrix.

Solve the system using Jacobi, Gauss-Seidel, and Successive Over-Relaxation (SOR) iterative solvers.

------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

THEORY:-

To solve a system of linear equations, we can use different iterative solvers like Jacobi, Gauss-Seidel, and Successive Over-Relaxation (SOR). Some important factors render the solution unstable. One factor is the spectral radius of the iteration matrix.

Note: We can use the inbuilt function to solve the system of linear equations, but for this problem, we will use different iterative solvers to solve the system of linear equations because we need to assess the consistency of iterative solvers.

Given the system of the linear equation can be written in matrix form as:

Splitting the matrix [A] into diagonal matrix [D], upper triangular matrix [U], and a lower triangular matrix [L]:

The iterative matrix for the different iterative solvers is given by:

Jacobi: `[D]^(-1) [L+U]`

Gauss-Seidel: `[D+L]^(-1) [U]`

SOR: `[D+[ω*L]]^(-1)*[ω*U+(ω-1)*D]`

Eigen Value of the matrix:

The eigenvalue (`λ`) of the matrix are the roots of the equation `|A-λI|=0`; where `I` is the identity matrix.

Calculating the eigenvalues for matrix A:

The characteristic equation for calculation of eigenvalue is given by:

After solving the above determinant, we get the equation:

`λ^3-7λ^2-60λ+402=0`

After solving the above polynomial equation; we get three eigenvalues:

`λ_1=-7.6663, λ_2=8.4900, and λ_3=6.1763`

Spectral Radius:

The spectral radius of a matrix is the maximum of the absolute values of the eigenvalues of the matrix, i.e. `ρ(A)=max⁡|λ|.`

Let A be an `(n X n)` matrix with complex or real elements with eigenvalues `λ_1,λ_2,…,λ_n.`

Then, the spectral radius `ρ(A)` is defined as:

`ρ(A)=max⁡{|λ_1 |,|λ_2 |,…,|λ_n |}`

For iterative solvers like Jacobi, Gauss-Seidel, and SOR; the characteristic equation is given by:

`X=TX+C`

`X-TX=C`

`X[I-T]=C` …… (I)

And, `AX=B` …… (II)

On comparing equation (I) and (II); we get:

`A=[I-T]`

The above criterion is only applicable when the spectral radius of the iteration matrix is less than 1 (`i.e. ρ(T)<1`).

However, it doesn’t mean that the equation has no solution; it only shows that our iterative solver will not be going to converge solution. We can still calculate the solution by the direct method, but it would be a time-consuming process and require a lot of labor.

The generalized equation for different iterative solvers:

Jacobi: `x^(k+1)=D^(-1) B+D^(-1) [-L-U] .x^k`

Gauss-Seidel: `x^(k+1)=[D+L]^(-1) B+[D+L]^(-1) [-U] .x^k`

SOR: `x^(k+1)=(1-ω) .x^k+ω[[D+L]^(-1) B+[D+L]^(-1) [-U] .x^k ]`

MATLAB Code:-

Creating a user-defined function for Jacobi iterative solver:

%Creating a user-defined function for Jacobi iterative solver
function [iteration,x]=func_Jacobi(L,U,D,B)
 
error=9e9;
tol=1e-4; %tolerance
x_old=[0;0;0];
iteration=1;
 
% convergence loop
while (error>tol)
    x=(inv(D)*B)+((inv(D)*(-L-U))*x_old);
    
    error=max(abs(x_old-x));
    x_old=x;
    iteration=iteration+1;
end
end

Creating a user-defined function for Gauss-Seidel iterative solver:

%Creating a user-defined function for Gauss-Seidel iterative solver
function [iteration,x]=func_Gauss_Seidel(L,U,D,B)
 
error=9e9;
tol=1e-4; %tolerance
x_old=[0;0;0];
iteration=1;
 
% convergence loop
while (error>tol)
    x=(inv(D+L)*B)+(((inv(D+L)*(-U))*x_old));
    
    error=max(abs(x_old-x));
    x_old=x;
    iteration=iteration+1;
end
end

Creating a user-defined function for SOR iterative solver:

%Creating a user-defined function for SOR iterative solver
function [iteration,x]=func_SOR(L,U,D,B,w)
 
error=9e9;
tol=1e-4; %tolerance
x_old=[0;0;0];
iteration=1;
 
% convergence loop
while (error>tol)
    x=((1-w)*x_old)+(w*((inv(D+L)*B)+(inv(D+L)*(-U)*x_old)));
    
    error=max(abs(x_old-x));
    x_old=x;
    iteration=iteration+1;
end
end

Main code:

clear all
close all
clc

% Given coefficient matrix
A=[5 1 2;-3 9 4;1 2 -7];
% Given RHS matrix
B=[-10;-14;33];

% relaxation factor used in SOR solver
w=1.1;

% values of unkowns [x1;x2;x3] is given by
x=A\B;
x_old=x;

% identity matrix
I=[1 0 0;0 1 0;0 0 1];

% Decomposition of matrix A:-

% Diagonal matrix
D1=[5 0 0;0 9 0;0 0 -7];
% Upper triangular matrix
U=[0 1 2;0 0 4;0 0 0];
% Lower triangular matrix
L=[0 0 0;-3 0 0;1 2 0];

% Assigning value of magnification factor to see the variation in spectral radius
magnification_factor=[0.5 1 2 3 4 5];
for i=1:length(magnification_factor)
    D=D1*magnification_factor(i);
    
% Operator for calling different iterative solver
% for Jacobi: iterative_solver=1, for Gauss-Seidel: iterative_solver=2, and for SOR: iterative_solver=3
iterative_solver=1                                                                                                                         ;

% JACOBI iterative solver
if iterative_solver==1
    
    % iteration matrix for Jacobi
    T_Jacobi=inv(D)*(L+U);
    
    % Calculating eigenvalues for T_Jacobi
    syms m
    eqn_Jacobi=det(T_Jacobi-I*m)==0;
    eigenvalue_T_Jacobi=double(solve(eqn_Jacobi,m));
    
    % Calculating spectral radius
    spectral_radius(i)=max(abs(eigenvalue_T_Jacobi));
    
    % Calling the function for Jacobi solver
    [iteration(i),x]=func_Jacobi(L,U,D,B);
end %end of Jacobi iterative solver

% GAUSS-SEIDEL iterative solver
if iterative_solver==2
    
    % iteration matrix for Gauss-Seidel
    T_Gauss=inv(D+L)*(U);
    
    % Calculating eigenvalues for T_Gauss
    syms n
    eqn_Gauss=det(T_Gauss-I*n)==0;
    eigenvalue_T_Gauss=double(solve(eqn_Gauss,n));
    
    % Calculating spectral radius
    spectral_radius(i)=max(abs(eigenvalue_T_Gauss));
    
    % Calling the function for Gauss-Seidel solver
    [iteration(i),x]=func_Gauss_Seidel(L,U,D,B);
end %end of Gauss-Seidel iterative solver
    
    
% Successive Over-Relaxation (SOR) iterative solver
if iterative_solver==3
    
    % iteration matrix for SOR
    T_SOR=inv(D+(w*L))*(((1-w)*D)-(w*U));
    
    % Calculating eigenvalue for T_SOR
    syms o
    eqn_SOR=det(T_SOR-I*o)==0;
    eigenvalue_T_SOR=double(solve(eqn_SOR,o));
    
    % Calculating spectral radius
    spectral_radius(i)=max(abs(eigenvalue_T_SOR));
    
    % Calling the function for SOR
    [iteration(i),x]=func_SOR(L,U,D,B,w);

end %end of SOR iterative solver
end
% Plots
if iterative_solver==1
    title_text=sprintf('JACOBI ITERATIVE SOLVER');
end
if iterative_solver==2
    title_text=sprintf('GAUSS-SEIDEL ITERATIVE SOLVER');
end
if iterative_solver==3
    title_text=sprintf('SUCCESSIVE OVER-RELAXATION (SOR) ITERATIVE SOLVER');
end
figure (1)
subplot (2,1,1)
plot (magnification_factor, spectral_radius,'-o','linewidth',2)
grid on
xlabel ('Magnification Factor', 'FontSize', 15)
ylabel ('Spectral Radius','Fontsize',15)
title (title_text,'FontSize',15)
subplot (2,1,2)
plot (spectral_radius, iteration,'-o','linewidth',2)
grid on
xlabel ('Spectral Radius', 'FontSize', 15)
ylabel ('No. of Iterations','FontSize',15)

       

RESULTS:-

Data showing the impact of different magnification factor on the spectral radius, number of iterations, and on the values of (`x_1,x_2, and x_3`):

CONCLUSION:-

From the above plots and data, we can conclude that:

1. Spectral radius decreases as we increase the magnification factor.
2. For the Gauss-Seidel and SOR solvers, the declivity in the values of the spectral radius is sharper than compared to the Jacobi solver.
3. As the spectral radius increases; the number of iterations to get a converged solution also increases.
4. When the spectral radius is greater than 1, then the iterative solvers take many iterations to run and even after that produces the ‘Nan’ values as the solution. Hence, we can say that the system becomes unstable for all the iterative solvers when the spectral radius for the iteration matrix is greater than 1.