Iterative solution of a system of linear equations and an analysis of spectral radius of a matrix

UNDERSTANDING LINEAR SYSTEMS
(ANALYSIS OF VARIOUS ITERATIVE SCHEMES TO SOLVE A SYSTEM OF LINEAR EQUATIONS TO FIND THE EIGEN VALUES AND SPECTRAL RADIUS)

(A) PROBLEM STATEMENT:

Given coefficient matrix:

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

Given Solution Matrix:

`X = [[x1],[x2],[x3]]`

Given RHS Matrix:

`B = [[10],[-14],[33]]`

We will be solving the system of linear equations using various iterative methods such as:

i. Gauss Jacobi method

ii. Gauss Siedel method

iii. Successive Over Relaxation (SOR) method

Theory:

For a general system of 3 linear equations in 3 variables we have:

`[[a11,a12,a13],[a21,a22,a23],[a31,a32,a33]][[x1],[x2],[x3]]=[[b1],[b2],[b3]]`

or `[A][X]=[B]` The solution is given by `[x]=[A]^-1[B]`

The Eigen Values of the matrix is given by the equation: `|A-lambdaI|=0`

`|([a11,a12,a13],[a21,a22,a23],[a31,a32,a33])-lambda([1,0,0],[0,1,0],[0,0,1])|=0`

 This equation is known as the the characteristic polynomial of the matrix and the roots are the eigen values `lambda` of the matrix:

`alambda^3+blambda^2+clambda+d=0`

The maximum value of the Eigen Values (roots of the eqn.) is known as the Spectral Radius `rho` of the matrix.

(B) MATLAB PROGRAMS:

1. Main Syntax of the program:

%% MATLAB Program to compute the Eigen Values and Spectral Radius of a 3x3 Matrix
clear all;
close all;
clc;
%Given Matrix:
A = [5 1 2; -3 9 4; 1 2 -7];
% Defining the Coefficient Matrix:
D = [ A(1,1) 0 0; 0 A(2,2) 0; 0 0 A(3,3) ]; %Diagonal elements
L = [ 0 0 0; A(2,1) 0 0; A(3,1) A(3,2) 0 ]; % Lower triangular elements
U = [ 0 A(1,2) A(1,3); 0 0 A(2,3); 0 0 0 ]; % Upper triangular elements

B = [10 ; -14; 33];
m = 1; %Diagonal Magnification Parameter
% Computation of Eigen Values:
% Solving the equation |A-xI|=0
% Characteristic equation for the system: -x^3+7x^2+60x-402=0

coeff = [-1 7 60 -402];
Eigen_value = roots(coeff);
% Checking the Eigen Values of the Matrix
Eigen_actual = eig(A);
%Checking the spectral radius of A:
Spect_rad = max(Eigen_value);

%% Solving the matrix by Iterative methods:
% Analytical Solution as reference:
exact_sol = AB;
% Solution by Gauss Jacobi
[iter_jac,x_jac,eig_jacobi,rho_jacobi] = jacobi_iter(D,L,U,B,m);
% Solution by Gauss Siedel
[iter_gs,x_gs,eig_gs,rho_gs] = siedel_iter(D,L,U,B,m);
% Solution by Successive Over Relaxation:
[iter_sor,x_sor,eig_sor,rho_sor] = sor_iter(D,L,U,B,m);

%% Solution of the system:
% Saving a copy of the solution Matrix before changing the diagonal terms
x_jac_old = x_jac; %Solution by Jacobi Method
x_gs_old = x_gs; %Solution by Gauss Siedel method
x_sor_old = x_sor; %Solution by SOR method

% Saving a copy of the spectral radius before changing the diagonal terms
iter_jac_old = iter_jac;
iter_gs_old = iter_gs;
iter_sor_old = iter_sor;

% Saving the values of eigen values before changing the diagonal terms:
eig_jac_old = eig_jacobi;
eig_gs_old = eig_gs;
eig_sor_old = eig_sor;

%% Setting the spectral radius to 1 by modifying diagonal elements:
for k = 1:3
    if k == 1
        rho_jac_old = rho_jacobi;
        m = 0.4602;
        [iter_jac,x_jac,eig_jacobi,rho_jacobi] = jacobi_iter(D,L,U,B,m);
    end
    if k == 2
        rho_gs_old = rho_gs;
        m = 0.5552;
        [iter_gs,x_gs,eig_gs,rho_gs] = siedel_iter(D,L,U,B,m);
    end
    if k == 3
        rho_sor_old = rho_sor;
        m = 0.5337;
        [iter_sor,x_sor,eig_sor,rho_sor] = sor_iter(D,L,U,B,m);
    end
end

2. User Defined Function: Gauss Jacobi Method

function [iter_jac,x_jac,eig_jacobi,rho_jacobi] = jacobi_iter(D,L,U,B,m)
%% Solving the Linear system using Gauss Jacobi Method:
error = 1e9;
tolerance = 1e-5;
iter_count = 0;
D = m*D;
%Iteration Matrix for Jacobi method:
T_jac = inv(D)*(L+U);

x_old = [ 0 ; 0 ; 0 ];
while (error > tolerance)
    x = -T_jac*x_old + inv(D)*B;
    error = max(abs(x-x_old));
    x_old = x;
    iter_count = iter_count + 1;
end
iter_jac = iter_count;
x_jac = x';
%% Calculation of Spectral Radius:
I = eye(3);
syms f(X);
f(X) = det(T_jac-X*I);
Xjac = double(coeffs(f,X,'All'));
eig_jacobi = roots(Xjac);
rho_jacobi = max(real(abs(eig_jacobi))); 
end

3. User Defined Function: Gauss Siedel Method

function [iter_gs,x_gs,eig_gs,rho_gs] = siedel_iter(D,L,U,B,m)
%% Solving the Matrix by Gauss Siedel method
error = 1e9;
tolerance = 1e-5;
iter_count = 0;

D = m*D;
% Iteration Matrix for Gauss Siedel
T_gs = inv(D+L)*U;
x_old = [ 0; 0; 0 ]; 

while error > tolerance
    x = -T_gs*x_old + inv(D+L)*B;
    error = max(abs(x-x_old));
    x_old = x;
    iter_count = iter_count + 1;
end
iter_gs = iter_count;
x_gs = x';
%% Calculation of Spectral Radius:
I = eye(3);
syms f(X);
f(X) = det(T_gs-X*I);
Xgs = double(coeffs(f,X,'All'));
eig_gs = roots(Xgs);
rho_gs = max(abs(real(eig_gs)));
end

4. User Defined Function: Successive Over Relaxation (SOR) method

function [iter_sor,x_sor,eig_sor,rho_sor] = sor_iter(D,L,U,B,m)
%% Solving the system by Succsssive Over Relaxation method:
error = 1e9;
tolerance = 1e-5;
iter_count = 0;
w = 0.9;
D = m*D;
%Iteration Matrix for SOR Method:
T_sor= inv(L+D)*(U-(L+D)*(1-w));

x_old = [0 ; 0; 0 ];
while error > tolerance
    x = x_old - inv(D+w*L)*((D+w*L) - (1-w)*D +w*U)*x_old + w*inv(D+w*L)*B;
    error = max(abs(x - x_old));
    x_old = x;
    iter_count = iter_count + 1;
end
iter_sor = iter_count;
x_sor = x';
%% Calculation of Spectral Radius:
I = eye(3);
syms f(X);
f(X) = det(T_sor-X*I);
Xsor = double(coeffs(f,X,'All'));
eig_sor = roots(Xsor);
rho_sor = max(real(abs(eig_sor)));
end

    

To the solve system iteratively we need to split the coefficient matrix `[A]` into the lower diagonal matrix `[L]`, the upper diagonal matrix `[U]` and the diagonal matrix `[D]` such that

`[A] = [L] +[D] +[U]`

where `L = [[0,0,0],[-3,0,0],[1,2,0]]`

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

          `U = [[0,1,2],[0,0,4],[0,0,0]]`

the following user defined functions solve the system iteratively by computing the iteration matrix for each method. It is updated with each iteration until convergence is achieved.

Iteration Matrices for different iterative schemes: (`[T]` = iteration matrix)

i. Gauss Jacobi method: `[T] = [D]^-1[L+U]`

ii. Gauss Siedel method:  `[T] = [L+D]^-1[U]`

iii. Successive Over Relaxation method: `[T] = [L+D]^-1([U]-(1-omega)[L+D])`

 

(C) RESULTS AND DISCUSSIONS:

1. SOLUTION OF THE SYSTEM:

The exact solution of the system is given by `[X]=[A]^-1[B]`

`[[x1],[x2],[x3]]=[[3.298507462686567],[1.268656716417911],[-3.880597014925373]]`

The solution by iterative methods are given below:

i. Gauss Jacobi Method:

`[[x1],[x2],[x3]]=[[3.298505536984095],[1.268657882614232],[-3.880594095875341]]`

The number of iterations taking place: 21

ii. Gauss Siedel Method:

`[[x1],[x2],[x3]]=[[3.298508561190570],[1.268657661889412],[-3.880596587861515]]`

The number of iterations taking place: 14

iii. Successive Over Relaxation (SOR) method:

`[[x1],[x2],[x3]]=[[3.298507986220588],[1.268656460583268],[-3.880596966137974]]`

The number of iterations taking place: 08

 For the SOR method we have used the relaxation parameter `omega`=0.9 i.e. under-relaxed condition. It is because keeping the system underrelaxed gives us the lowest number of  iterations as compared to `omega`>1 (over-relaxed condition). In that case the system becomes divergent and convergence is achieved after a far larger number of iterations as compared to the under-relaxed condition. The system changes to Gauss-Siedel method at `omega`=1.

2. EIGEN VALUES OF THE SYSTEM:

 i. Gauss Jacobi Method:

`[[lambda1],[lambda2],[lambda3]]=[[-0.0487 + 0.5078i],[-0.0487 - 0.5078i],[0.0975]]`

ii. Gauss Siedel Method:

`[[lambda1],[lambda2],[lambda3]]=[[0],[0.3276],[-0.0387]]`

iii. Successive Over Relaxation Method:

`[[lambda1],[lambda2],[lambda3]]=[[0.2276],[-0.1387],[-0.1]]`

3. SPECTRAL RADIUS OF THE SYSTEM:

By exact method we have `rho` = 8.490

Solving the system iteratively we have the values of spectral radii by different methods:

i. Gauss Jacobi Method: `rho` = 0.5102

ii. Gauss Siedel Method: `rho` = 0.3276

iii. Successive Over Relaxation Method: `rho` = 0.2276

The SOR method gives the lowest value of `rho` with `omega` = 0.9

(D) CHANGING THE SPECTRAL RADIUS TO 1.1 BY MODIFYING THE DIAGONAL ELEMENTS:

In the following program we define m as the Diagonal Magnification Parameter. By changing it suitably we can set the spectral radius to 1.1 each time for different iterative schemes.

i. Gauss Jacobi Method:

Changing m to 0.4602 gives us a value of `rho`=1.1087 and the matrix becomes unsolvable by the above method.

ii. Gauss Siedel Method:

Changing m to 0.5552 gives us a value of `rho`=1.1035 and the matrix becomes unsolvable by the above method.

 iii. Successive Over Relaxation Method:

Changing m to 0.5337 gives us a value of `rho`=1.1007 but surprisingly the matrix is solved by the SOR method.

(E) EFFECT OF DIAGONAL MAGNIFICATION ON SPECTRAL RADIUS AND NUMBER OF ITERATIONS:

Diagonal magnification refers to multiplying the diagonal elements with a constant to change the spectral radius. Here we have used a range of numbers to observe the same.

1. MATLAB CODE:

%% MATLAB Program to observe the effect of Diagonal Magnification on various Iterative solvers;
clear all;
close all;
clc;
%Given Matrix:
A = [5 1 2; -3 9 4; 1 2 -7];
% Defining the Coefficient Matrix:
D = [ A(1,1) 0 0; 0 A(2,2) 0; 0 0 A(3,3) ]; % Diagonal elements
L = [ 0 0 0; A(2,1) 0 0; A(3,1) A(3,2) 0 ]; % Lower triangular elements
U = [ 0 A(1,2) A(1,3); 0 0 A(2,3); 0 0 0 ]; % Upper triangular elements

B = [10 ; -14; 33];

m = linspace(5,0.2,10);
% Computation of Eigen Values:
% Solving the equation |A-xI|=0 
% Characteristic equation for the system: -x^3+7x^2+60x-402=0

coeff = [-1 7 60 -402];
Eigen_value = roots(coeff);
% Checking the Eigen Values of the Matrix
Eigen_actual = eig(A);
%Checking the spectral radius of A:
Spect_rad = max(Eigen_value);

%% Solution of the matrix by Iterative methods:
% Analytical Solution as reference:
exact_sol = AB;

for i = 1:length(m)
% Solution by Gauss Jacobi
[iter_jac(i),x_jac,eig_jacobi,rho_jacobi(i)] = jacobi_iter(D,L,U,B,m(i));
% Solution by Gauss Siedel
[iter_gs(i),x_gs,eig_gs,rho_gs(i)] = siedel_iter(D,L,U,B,m(i));
% Solution by Successive Over Relaxation:
[iter_sor(i),x_sor,eig_sor,rho_sor(i)] = sor_iter(D,L,U,B,m(i));
end

%% Plotting the solution:
% Plotting the Iteration Numbers vs SPectral Radii
fig_1 = figure(1);
hold on;
plot(rho_jacobi,iter_jac);
plot(rho_gs,iter_gs);
plot(rho_sor,iter_sor);
legend('Gauss Jacobi','Gauss Siedel','SOR');
grid minor;
xlabel('Spectral Radius: rho');
ylabel('Number of Iterations: ');
title('Plot of Number of Iterations vs Spectral Radius (rho) for Iterative Schemes');

% Plotting the Diagonal Magnification Values vs Spectral Radii
fig_2 = figure(2);
hold on;
plot(m,rho_jacobi,'--m');
plot(m,rho_gs,'*b');
plot(m,rho_sor,'-.g');
legend('Gauss Jacobi','Gauss Siedel','SOR');
grid on;
xlabel('Diagonal Magnification Parameter');
ylabel('Spectral Radius: (rho)');
title('Plot of the Spectral Radius: (rho) value vs a range of diagonal magnification values');


2. RESULTS AND DISCUSSIONS:

(i) Plotting the Iteration Numbers vs Spectral Radii:

iter_nos

As we increase the spectral radius, the number of iterations also increases. The slope of the curve for the Jacobi method is much more compared to Gauss Siedel and SOR, hence requiring a larger number of iterations for convergence.

(ii) Plotting the Spectral Radii vs Diagonal Magnification Values:

daig_mag

The spectral radius decreases more rapidly for SOR and Gauss Siedel Method as compared to  Jacobi method for the same increase in Diagonal Magnification Parameter.


Projects by Priyotosh Bairagya

LAMINAR INCOMPRESSIBLE FLOW SIMULATION THROUGH A PIPE IN OPENFOAM: [PART 2/2]: SYMMETRY BOUNDARY CONDITION: OBJECTIVES:1. Creating the Mesh-Script (blockMeshDict file) with symmetry boundary condition.2. Simulation Results and post processing the velocity profile at Read more

LAMINAR INCOMPRESSIBLE FLOW SIMULATION THROUGH A PIPE IN OPENFOAM: [PART 1/2]: WEDGE BOUNDARY CONDITION: OBJECTIVES:1. Calculation of Pre-eliminary quantinites related to the flow.2. Creating the Mesh script (blockMeshDict file) for specifying the geometry and bounda Read more

I. Interpolation Schemes in Finite Volume Method: The approximation of surface and volume integrals may require values of the variable at locations other than the computational nodes of the CV. Values at these locations are obtained using interpolation formulae. Some o Read more

ANALYSIS OF NUMERICAL STABILITY OF VARIOUS ITERATIVE SOLVERS FOR TRANSIENT 2D HEAT CONDUCTION: [Part: 3/3] INTRODUCTION: The criterion of stability of a numerical scheme is determined by the way the errors propagate while the solution moves from one time-step to th Read more

MESH GENERATION AND ANALYSIS USING BLOCKMESH FOR FLOW OVER A BACKWARD FACING STEP: The purpose of the following project is to generate the geometry for a variation of the incompressible cavity flow problem in OpenFOAM. For this purpose we have modified the lid-driven c Read more

ANALYSIS OF VARIOUS ITERATIVE SCHEMES FOR THE SOLUTION OF A 2D HEAT CONDUCTION PROBLEM: [PART: 2/3] In the previous part we had explained the problem statement and the MATLAB Program in detail. In this Part we are going to explain the outputs from the 2D Heat Conduct Read more

NUMERICAL SOLUTION OF 1D SUPERSONIC NOZZLE FLOW SIMULATION BY MAC-CORMACK METHOD PROJECT OBJECTIVES: i. Numerical solution of the governing equations in both conservative and non-conservative forms. ii. Creating user defined functions for calculating the flow quan Read more

Simulation of a 2D Heat Conduction problem in steady and unsteady/transient forms using iterative methods. Project Objectives: 1. Solving the 2 Dimensional Heat conduction equation in the generalized form using various iterative techniques: i. Explicit Solver (for Read more

WEEK 4: (Effect of Grid-Size on output for the solution of 1D linear wave equation) 1. Problem Setup: Given Partial Differential Equation: `(delu)/(delt) + c(delu)/(delx) = 0` Numerical Discretization:  `u_(i,n+1) = u_(i,n) + (cDeltat)/(Deltax)*(u_(i-1,n+1) Read more


Loading...

The End