Stability analysis in Steady state equation.

Aim : To demonstrate stability analysis of steady state equations.

Procedure : The analysis will be demonstrated by computing spectral radius of the iterative matrix of an equation which is solved by Jacobi , Gauss seidal and SOR methods. And then the data is compared with respect to other parameters to observe the trend.

We first consider a 3x3 matrix and subject it to LDU decomposition as follows :

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

`U = [[0,1,2],[0,0,4],[0,0,0]]`
`D1 = [[5,0,0],[0,9,0],[0,0,-7]]`
A = L+U+D1

Identity matrix ==> `I = [[1,0,0],[0,1,0],[0,0,1]]`

The general iterative function would take the form as : `X^(n ew) = T*X^(old) + C`

The general system of linear equations is solved by taking form : A*X = B. Now expanding this form and replacing it with LDU decomposition terms, we have :

Jacobi method = `X^(n ew) = D^-1*B - D^-1*(L+U)*X^old`

T_jac = `-D^-1*(L+U)*X_(old)`

Gauss Seidal method = `X^(n ew) = (L+D)^-1 *B - (L+D)*U*Xold`

T_gs = `-(L+U)^-1*U`

SOR method = `X^(n ew) = (1-alpha)*X^(old) + alpha*X_(GS)`

T_sor = `I - alpha*((L+U)^-1*U+I))`

where `X^(n ew)` = current iteration values , `X^(old)` = previous iteration values , T = iteration matrix

Note : the eigen values are found out by the 'roots' and 'poly' command and the in-built function is used to check the answer against them.

In the following code, we will be magnifying the values of diagonal elements of the matrix and observe how much iterations it take to converge for the respective methods.

the code is as follows :

close all
clear all
clc

% inputs
alpha = 0.9;

L = [0 0 0;-3 0 0;1 2 0]
U = [0 1 2;0 0 4;0 0 0]
D1 = [5 0 0;0 9 0;0 0 -7]
A = L+U+D1;
B = [1; -2; -1];
I = [1 0 0;0 1 0;0 0 1];

% iteration matrix
%Tjac = -inv(D)*(L+U);
%Tgs = -inv(L+D)*U;
%Tsor = (I - alpha)*(inv(L+U)*U + I);

% spectral radius
%rho_jac = max(abs(eig(Tjac)));
%rho_gs = max(abs(eig(Tgs)));
%rho_sor = max(abs(eig(Tsor)));


% magnification values for diagonal elements
magnify = [0.5,1,2,3,4,5]
%magnify = [0.1:0.1:2];

% iteration input to change the method
iteration = 1

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if iteration == 1   % jacobi method
for i = 1:length(magnify);
    D = D1*magnify(i)

[jacobi_iter(i),Xjac,Tjac] = jacobi(A,B,L,U,D)
rho_jac(i) = max(abs(eig(Tjac)));
eigen_jac = roots(poly(Tjac));     %  eigen value of iteration matrix by jacobi method

end

figure(1)
subplot(2,1,1);
plot(magnify,rho_jac,'-o','linewidth',2);
xlabel('magnification values for diagonal matrix')
ylabel('spectral radius');
title('jacobi method');
grid on
subplot(2,1,2);
plot(rho_jac,jacobi_iter,'-o','linewidth',2);
xlabel('spectral radius');
ylabel('jacobi iterations');
title('jacobi method');
grid on
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


if iteration == 2;        % Gauss seidal method
for i = 1:length(magnify);
    D = D1*magnify(i)

[gs_iter(i),Xgs,Tgs] = gs(A,B,L,U,D)
rho_gs(i) = max(abs(eig(Tgs)));
eigen_gs = roots(poly(Tgs));          %  eigen value of iteration matrix by gauss seidal method

end

figure(1)
subplot(2,1,1);
plot(magnify,rho_gs,'-o','linewidth',2);
xlabel('magnification values for diagonal matrix')
ylabel('spectral radius');
title('gauss seidal method');
grid on
subplot(2,1,2);
plot(rho_gs,gs_iter,'-o','linewidth',2);
xlabel('spectral radius');
ylabel('gs iterations');
grid on
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


if iteration == 3;         % Sucessive over-relaxation method
for i = 1:length(magnify);
    D = D1*magnify(i)

[sor_iter(i),Xsor,Tsor] = sor(A,B,L,U,D,alpha)
rho_sor(i) = max(abs(eig(Tsor)));
eigen_sor = roots(poly(Tsor));       %  eigen value of iteration matrix by sor method


end
figure(1)
subplot(2,1,1);
plot(magnify,rho_sor,'-o','linewidth',2);
xlabel('magnification values for diagonal matrix')
ylabel('spectral radius');
title('sucessive over-relaxation method');
grid on
subplot(2,1,2);
plot(rho_sor,sor_iter,'-o','linewidth',2);
xlabel('spectral radius');
ylabel('sor iterations');
title('sucessive over-relaxation method');

grid on
end

the functions for the code are as follows :

function [jacobi_iter,Xjac,Tjac] = jacobi(A,B,L,U,D)


jacobi_iter = 1;
error = 9e9;
tol = 1e-6;
Xold = zeros(3,1)

while error > tol;

Xjac = inv(D)*B - inv(D)*(L+U)*Xold;
Tjac = -inv(D)*(L+U);

error = max(max(abs(Xold - Xjac)));
Xold = Xjac;
jacobi_iter = jacobi_iter + 1;
end

end

function [gs_iter,Xgs,Tgs] = gs(A,B,L,U,D)


gs_iter = 1;
error = 9e9;
tol = 1e-6;
Xold = zeros(3,1)

while error > tol;

Xgs = (inv(L+D)*B)-(inv(L+D)*U*Xold);
Tgs = -inv(L+D)*U;

error = max(max(abs(Xold - Xgs)));
Xold = Xgs;
gs_iter = gs_iter + 1;
end

end

function [sor_iter,Xsor,Tsor] = sor(A,B,L,U,D,alpha)

I = eye(3)
sor_iter = 1;
error = 9e9;
tol = 1e-6;
Xold = zeros(3,1)

while error > tol;

% Xsor = (alpha*inv(L+D)*B) - (I-alpha)*((inv(L+D)*U+I))*Xold;
Xsor = (1-alpha)*Xold + alpha*((inv(L+D)*B)-(inv(L+D)*U*Xold));
Tsor = (I - alpha)*(inv(L+D)*U + I);

error = max(max(abs(Xold - Xsor)));
Xold = Xsor;
sor_iter = sor_iter + 1;
end

end

the graphs for the acquired data is as follows :

the workspace snipets are as follows :

Conclusion : As we can see from the graph, that as the magnification value for the diagonal matrix increases , the spectral radius decreases.

And as the spectral radius increases the iterations required for that particular method to converge, increases.

the whole point of this data is to show how the spectral radius contributes to the stability of the linear systems.

The variation in diagonal matrix affects the spectral radius. And when the spectral radius goes above 1 , the no. of iterations gets high and the solutions becomes unsolvable.

The diagonal dominance or the spectral radius tells us about the matrix whether it can be solved or not, which gives us an advantage to predetermine the convergence of equation even before undergoing the iteration methods. This can save time and help user understand which method to go for solving the equation.