Effect Of Spectral Radius And Diagonal Magnification In The Convergence Rate Of An Iterative Solver For A Linear System Using MATLAB.

 

(A) Problem defination:- Solve the following Linear equation using Gauss-Jacobi, Gauss-Seidal and SOR iterative method and understand the following parameters:-

  1. Spectral radius and the impact of spectral radius on the convergence of the solution.
  2. Effect of diagonal magnification on spectral radius for each iterative method. 

The coefficient matrix is given as:-

The solution matrix is given as:-

The Residual (RHS) matrix is given as:-

The general matrix form a linear solution is given as:-

The solution matrix is the product of inverse of the coefficient matrix and the residual matrix. The form is given as:-

The characteristic equation is written as:-

The polynomial equation is given as:-

The maximum absolute eigen value of a given matrix is the spectral radius.

The above coefficient matrix [A] can be decomposed as upper matrix [U], lower matrix [L] and diagonal matrix [D]. The matrix [A] can be written as:-

[A] = [U]+[D]+[L]. where

   

The above decomposed matrix form is used to form iteration matrix for each iteration methods.

Now, the iteration matrix for Gauss-jacobi, Seidal and SOR is mentioned below:-

  • Gauss-Jacobi:-

Let,

Ax = b  be a square matrix of n linear equation where:

The A matrix can be decomposed into diagonal matrix and residual matrix.

Then the matrix is solve iteratively via:-

Therefore the iterative matrix [T_jacobi] is given as:-

  • Gauss-Seidal:-

Let,

Ax = b  be a square matrix of n linear equation where:

The A matrix can be decomposed into  A = L + U

Then the matrix can be solve via:-

Therefore the iterative matrix [T_seidal] is given as:-

  • SOR:-

The decomposition matrix is given as:- [A] = [U] +[D] + [L]

The iterative form is given as:-

And therefore the iteration matrix is given as:-

(B) MATLAB CODE

1. Main Syntex to understand the impact of diagonal magnification on the spectral radius .

close all
clear all
clc


% Solving Linear system using iterative matrix
% Defining the iterative matrix A as its decomposed form LDU
% The following linear system is diagonally dominant.

A=  [5 1 2; -3 9 4; 1 2 -7]; % Coefficient Matrix.
L = [0 0 0; -3 0 0; 1 2  0]; % Lower Diagonal Matrix.
D = [5 0 0;  0 9 0; 0 0 -7]; % Diagonal Matrix.
U = [0 1 2;  0 0 4; 0 0  0]; % Upper Diagonal Matrix.

B = [-10; 14; 33];           % Residual Matrix

X = AB;                     % Inverse of Matrix A.

% The Coefficient Matrix is A = L + D + U
% The Eigen value of coefficient matrix is |A-xI| = 0
% The polynimial characteristic equation is given by -x^3+7x^2+60x-402 = 0
Eigen_Value = eig(A);
Spectral_Radius = max(real(Eigen_Value)); % Spectral radius is the maximum absolute eigen value.



Choose_the_iterative_method = input('1. Gauss_Jacobi  2. Gauss_Seidal  3. SOR');
%% Solving the above matrix using iterative solver.
 % Solving using Gauss-Jacobi Method.
 if Choose_the_iterative_method ==1;
      m = 0.4602                  % Diagonal Magnification Factor.
     [spectral_radius,x_jac,iter_jac,jacobi_eigen,error] = Jacobi_iter(L,D,U,B,m);
    
 end
%% Solving using Gauss-Seidal Method.
if Choose_the_iterative_method ==2;
    
       m = 0.5552                   % Diagonal Magnification Factor.
      [spectral_radius,x_seidal,iter_seidal,seidal_eigen,error] = Seidal_iter(L,D,U,B,m);
   
end
%% Solving using Successive-over-Relaxation.
if Choose_the_iterative_method ==3;
        m = 1.0602                   % Diagonal Magnification Factor.
      [spectral_radius,x_sor,sor_iter,sor_eigen,error] = SOR_iter(L,D,U,B,m);
   
end

2. The user defined functions: Gauss-jacobi method.

% Jacobi iteration function

function [spectral_radius,x_jac,iter_jac,jacobi_eigen,error] = Jacobi_iter(L,D,U,B,m)

error     = 1e9;
tolerance = 1e-5;
D = m*D;
x_old = [ 0; 0; 0];
T_jacobi = inv(D)*(L+U); % Iteration matrix for Gauss-Jacobi.
iteration = 1;
  
        while (error > tolerance);
            
            x = (inv(D)*B) - (T_jacobi*x_old);
            error = max(abs(x_old-x));
            x_old = x;
            iteration = iteration + 1;
        
        end
 iter_jac = iteration;
 x_jac = x;
 % Defining spectral Radius
 % Spectral radius is the largest absolute Eigen value of the iteration
 % matrix.
 eye(3);             % Defining 3x3 Identity matrix.
 syms('LAMDA');      % Defining LAMDA as a symbol.
 jacobi_characteristic_eqn = det(T_jacobi - LAMDA*eye(3));              % Characteristic equation.                   
 jacobi_coefficient = double(coeffs(jacobi_characteristic_eqn,'ALL'));  % Coefficients of the the determinant.
 jacobi_eigen = roots(jacobi_coefficient);                              % Eigen values using in-built "roots" command.             
 spectral_radius = max(real(abs(jacobi_eigen)));
 
end

3. User defined function: Gauss-Seidal

% Gauss-Seidal Iterative Method

function [spectral_radius,x_seidal,iter_seidal,seidal_eigen,error] = Seidal_iter(L,D,U,B,m)

error = 1e9;
tolerance = 1e-5;
D = m*D;
x_old = [0;0;0];
T_seidal = (L+D)U;  % Iteration matrix for Gauss-Seidal.
iteration = 1;

        while (error > tolerance);
            
            x = (L+D)(B) - T_seidal*x_old;
            error = max(abs(x_old - x));
            x_old = x;
            iteration = iteration +1;
            
        end
iter_seidal = iteration;
x_seidal = x;
% Defining the Eigen values of the iteration matrix and Spectral Radius.
% Spectral radius is the largest eigen value of the matrix.

eye(3);             % Defining the 3x3 identity matrix.
syms('LAMDA');      % Defining LAMDA as a symbol.
characteristic_equation = det(T_seidal-(LAMDA*eye(3)));
seidal_coefficients = double(coeffs(characteristic_equation,'ALL'));
seidal_eigen = roots(seidal_coefficients);
spectral_radius = max(real(abs(seidal_eigen)));

end

4. User defined function: SOR

% Successive-Over-Relaxation
function [spectral_radius,x_sor,sor_iter,sor_eigen,error] = SOR_iter(L,D,U,B,m)

error = 1e9;
tolerance = 1e-5;
D = m*D;
x_old = [0;0;0];
%prompt = ' what is the Relaxation factor value?';
%omega = input(prompt);
omega = 0.5;
T_SOR = (inv(D+(omega*L))*((omega*U)+((omega-1)*D)));
iteration = 1;

        while (error > tolerance);
            
            x = (inv(D+(omega*L))*(omega*B))-(T_SOR*x_old);
            error = max(abs(x_old - x));
            x_old = x;
            iteration = iteration + 1;
            
        end
sor_iter = iteration;
x_sor = x;
% Defining the Eigen values of the iteration matrix and Spectral Radius.
% Spectral radius is the largest eigen value of the matrix.
eye(3);             % Defining the 3x3 identity matrix.
syms('LAMDA');      % Defining LAMDA as a symbol.
characteristic_equation = det(T_SOR-(LAMDA*eye(3)));
sor_coefficients = double(coeffs(characteristic_equation,'ALL'));
sor_eigen = roots(sor_coefficients);
spectral_radius = max(real(abs(sor_eigen)));
end

    

Results And Discussion:

The exact solution of the matrix is given as 

In matlab the exact soutions result obtained are as follows:-

>> X

X =

    -1
     3
    -4

And the solution by iterative method is given as:-

1. Gauss-Jacobi:-

1. Gauss_Jacobi  2. Gauss_Seidal  3. SOR1

m =

     1

>> x_jac

x_jac =

   -1.0000
    3.0000
   -4.0000

>> Eigen_value
Undefined function or variable 'Eigen_value'.
 
Did you mean:
>> Eigen_Value

Eigen_Value =

    8.4900
    6.1763
   -7.6663

>> error

error =

   6.9669e-06

>> iter_jac

iter_jac =

    21

>> spectral_radius

spectral_radius =

    0.5102

>> 

The total iteration taken is  = 21.

2. Gauss-Seidal:-

1. Gauss_Jacobi  2. Gauss_Seidal  3. SOR2

m =

     1

>> x_seidal

x_seidal =

   -1.0000
    3.0000
   -4.0000

>> Eigen_Value

Eigen_Value =

    8.4900
    6.1763
   -7.6663

>> error

error =

   4.2115e-06

>> iter_seidal

iter_seidal =

    15

>> spectral_radius

spectral_radius =

    0.3276

>> 

The total iteration taken is = 15.

3. SOR:-

1. Gauss_Jacobi  2. Gauss_Seidal  3. SOR3

m =

     1

>> x_sor

x_sor =

   -1.0000
    3.0000
   -4.0000

>> Eigen_Value

Eigen_Value =

    8.4900
    6.1763
   -7.6663

>> error

error =

   8.0282e-06

>> sor_iter

sor_iter =

    10

>> spectral_radius

spectral_radius =

    0.1707

>> sor_eigen

sor_eigen =

  -0.1707 + 0.0000i
  -0.0373 + 0.1356i
  -0.0373 - 0.1356i

>> 

The total iteration taken is = 10.

For the SOR method we have used relaxation factor "w=0.85" which is under-relaxed condition.Typically, using under-relaxed condition, the solution is obtained at lesser iterations compared to "w>1" (over-relaxed conditions). In over relaxed condition, the linear equation system becomes divergent and the convergence is obtained at more number of iterations.

Spectral Radius Of The System:

1. For gauss- Jacobi iteration method:- ρ(rho) = 0.5102

2. For gauss-Seidal iteration method:- ρ(rho)  = 0.3276

3. For Successive-Over-relaxation method:- ρ(rho) = 0.1707

 

EFFECT OF SPECTRAL RADIUS AND DIAGONAL MATRIX MAGNIFICATION ON THE CONVERGENCE OF THE SOLUTION:

close all
clear all
clc


% Solving Linear system using iterative matrix
% Defining the iterative matrix A as its decomposed form LDU
% The following linear system is diagonally dominant.

A=  [5 1 2; -3 9 4; 1 2 -7]; % Coefficient Matrix.
L = [0 0 0; -3 0 0; 1 2  0]; % Lower Diagonal Matrix.
D = [5 0 0;  0 9 0; 0 0 -7]; % Diagonal Matrix.
U = [0 1 2;  0 0 4; 0 0  0]; % Upper Diagonal Matrix.

B = [-10; 14; 33];           % Residual Matrix

X = AB;                     % Inverse of Matrix A.

% The Coefficient Matrix is A = L + D + U
% The Eigen value of coefficient matrix is |A-xI| = 0
% The polynimial characteristic equation is given by -x^3+7x^2+60x-402 = 0
Eigen_Value = eig(A);
Spectral_Radius = max(real(Eigen_Value)); % Spectral radius is the maximum absolute eigen value.
              
mag1 =linspace(0.5,5,20);                 % Diagonal Magnification Length.
Choose_the_iterative_method = input('1. Gauss_Jacobi  2. Gauss_Seidal  3. SOR');
%% Solving the above matrix using iterative solver.
 % Solving using Gauss-Jacobi Method.
 if Choose_the_iterative_method ==1;
     for i = 1:20
         D = mag1(i)*D;
         [spectral_radius(i),x_jac,iter_jac(i),jacobi_eigen,error] = Jacobi_iter_1(L,D,U,B);
     end
     
     fig_1 = figure(1);
     loglog(mag1,spectral_radius,'r','linewidth',1)
     xlabel('Diagonal Magnification');
     ylabel('spectral radius');
     grid on
     
     fig_2 = figure(2);
     loglog(spectral_radius,iter_jac,'g','linewidth',1)
     xlabel('spectral radius');
     ylabel('Iteration');
     grid on
     
 end
%% Solving using Gauss-Seidal Method.
if Choose_the_iterative_method ==2;
    for i = 1:20
        D = mag1(i)*D;
        [spectral_radius(i),x_seidal,iter_seidal(i),seidal_eigen,error] = Seidal_iter_1(L,D,U,B);
    end
    
     fig_2 = figure(1);
     loglog(mag1,spectral_radius,'r','linewidth',2)
     xlabel('Diagonal Magnification');
     ylabel('spectral radius');
     grid on
     
     fig_2 = figure(2);
     loglog(spectral_radius,iter_seidal,'g','linewidth',1)
     xlabel('spectral radius');
     ylabel('Iteration');
     grid on
end
%% Solving using Successive-over-Relaxation.
if Choose_the_iterative_method ==3;
    for i =1:20;
        D = mag1(i)*D;
        [spectral_radius(i),x_sor,sor_iter(i),sor_eigen,error] = SOR_iter_1(L,D,U,B);
    end
    
    fig_1 = figure(1);
    loglog(mag1,spectral_radius,'b','linewidth',1)
    xlabel('Diagonal Magnification');
    ylabel('spectral radius');
    grid on
    
    fig_2 = figure(2);
    loglog(spectral_radius,sor_iter,'g','linewidth',1)
    xlabel('spectral radius');
    ylabel('Iteration');
    grid on
end

PLOTS:

 

  • Spectral radius of the iteration matrix for the SOR method . The plot shows the dependence on the spectral radius of the Jacobi iteration matrix .

  • "*" [in the table] means that the iteration matrix matrix couldn't be solved for both the gauss-jacobi and gauss-seidal method as the ρ(rho)>1 .
  • From the above table it is also clear that as each diagonal matrix element increases, the ρ(rho) [spectral radius] decreases.
  • The decrease in spectral radius increases the convergence rate and thus the total iteration decreases to solve the matrix decreases.
  • For SOR iterative mathod, at ω = 0.9 the lowest value of value of spectral radius is obtained which is the optimum relaxation factor. If considered ω (relaxation factor)>1 i.e., making the condition over-relaxed, then the spectral radius increases and more iterations are required to achieved convergence.

Projects by Ankit Chakrabarty

Title:- Geometry cleanup and meshing of an automotive hood for structural analysis   Objective:  Perform geometry cleanup of automotive hood (both inner and outer panel). Perform mid-surface extraction of both the panels and Hem both the surfaces . Perf Read more

Abstract: This project presents the Modelling of the  behaviour of  fluid and its properties when it passes through a convergent-divergent nozzle using a Finite Difference Method. However, to simplify the understanding of the dynamics of the flow in the nozz Read more

 Turbocharger: Literature Review: A turbocharger is mechanical rotating device which forcebly induces air at high pressure into the internal combustion engine, in turn increases the power and effeciency (thermal and volumetric effeciency) of the engine. A turb Read more

CASE STUDY:-  A Thermodynamic Analysis Of A Tractor Engine For Use In India With Given Engine Performance Targets. A Brief Literature Review: A tractor is an engineered vehicle specifically designed to deliver Power and traction/tractive effort to haul heavy equi Read more

Solving 2-D 2nd Order Linear Transient &Steady  Heat Conduction Equation Using Both Explicit And Implicit Numerical Schemes In MATLAB.   Explicit Scheme : We say that a scheme is explit when the information at time level "n+1" depends on previous time Read more

  CASE STUDY: Effect of Burn Rate on engine performance Of A 6-Cyl-Turbocharged-DI Engine Using A Non-predictive Combustion model (DI WIEBE) Abstract: Analytical functions approximating the burn rate in internal combustion engines are useful and cost-effective to Read more

Solving 1-D Linear Convection Using First-Order Backward Difference And Forward Difference method Using MATLAB LITERATURE REVIEW: WHAT IS CONVECTION? Convection is the sum of bulk transport of the  fluid  and brownian/ osmotic dispersion of fluid constituen Read more

Deriving the 4th order approximation for 2nd order derivative of a function using Taylor table for centered, skewed right-sided stencil and skewed left_sided stencil and compare the error results for the same using MATLAB.   Given function, f(x)    =& Read more

Objective:- (1) Determine the engine parameters running at 1800 rpm and record the same. (2) Further increase the power output obtained at 3600 rpm by 10%.   Model:-  The GT model Of the 1-cylinder engine is shown below:- Observations (@1800rpm):-  Read more


Loading...

The End