Understanding Linear systems
TARGET OF THE PROJECT:
To determine the Eigen values and spectral radius of the iterative matrix using iterative methods such as Jacobi,Gauss Seidel and SOR methods and observing the changes being made to the spectral radius by changing the magnification factor of the diagonal matrix in increments of 1.1.
EIGEN VALUES CALCULATION:
The Eigen values (`lambda`) of a matrix can be computed by taking the roots of the equation|A-`lamda`*I|=0 and the spectral radius is the maximum value of the Eigen values denoted by `rho`(rho).The eigen values of the matrix are determined by finding the roots of the characteristic equation given by
`lambda^3-Alambda^2+Blambda-C=0`
where A is the trace of the coefficient matrix
B is the sum of the determinant of 2x2 matrices with diagonal elements and
C is the determinant of the Coefficient matrix.
The given set of linear equations in matrix form Ax=B is given by
`[[5,1,2],[-3,9,4],[1,2,-7]] [[x1],[x2],[x3]]=[[10],[-14],[33]]`
The matrix A is further split into Lower triangular matrix (L), diagonal matrix (D) and upper triangular matrix (U). The expression is as follows
`[[5,1,2],[-3,9,4],[1,2,-7]]=[[0,0,0],[-3,0,0],[1,2,0]]+[[5,0,0],[0,9,0],[0,0,-7]]+[[0,1,2],[0,0,4],[0,0,0]]`
The Iterative matrix for the Jacobi method is given by `D^-1(L+U)`
The Iterative matrix for the Gauss Seidel method is given by `(D+L)^-1*U`
The Iterative matrix for the SOR method is given by `(omega*L+D)^-1[(1-omega)D-omega*U]`
The diagonal matrix is multiplied with various magnification factors to see its effect on the spectral radius.In order to attain convergence, the spectral radius of the iteration matrices should be less than 1. The code for it is given below.
MAIN SCRIPT
clear all
close all
clc
% coefficient matrix
A=[5 1 2;-3 9 4;1 2 -7];
B=[10; -14; 33];
%diagonal matrix
D=diag(diag(A));
% lower diagonal matrix
L=tril(A)-D;
% upper triangular matrix
U=triu(A)-D;
% magnification factor
m=[0.5 1 2 3 4 6];
omega=1.5;
for i=1:length(m)
%diagonal magnification
D=diag(diag(A));
D=m(i)*D;
%iterative matrix components
X_jac=D;
Y_jac=-(L+U);
X_gs=L+D;
Y_gs=-U;
X_sor=(omega*L+D);
Y_sor=(((1-omega)*D)-omega*U);
%iteration matrices
jac_I=inv(X_jac)*(Y_jac);
gs_I=inv(X_gs)*(Y_gs);
sor_I=inv(X_sor)*(Y_sor);
% computing the eigen values without the inbuilt function
jac_eigen=roots(poly(jac_I));
gs_eigen=roots(poly(gs_I));
sor_eigen=roots(poly(sor_I));
%spectral radius calculation
rho_jac(i)=max(abs(jac_eigen));
rho_gs(i)= max(abs(gs_eigen));
rho_sor(i)= max(abs(sor_eigen));
%calculating no of iterations required to calculate eigen values and
%spectral radius
%jacobi
if rho_jac(i)< 1
[iter_jac(i),x_jac]=iter_sln(jac_I,X_jac,B);
sprintf('No of jacobi iterations with %d as magnification factor and %d as rho= %dn',m(i),rho_jac(i),iter_jac(i));
else
iter_jac(i)=nan;
sprintf('solution is unstable for rho_jac =%d with %d as diagonal magnification factorn ',m(i),rho_jac(i),iter_jac(i));
end
%GS
if rho_gs(i)< 1
[iter_gs(i),x_gs]=iter_sln(gs_I,X_gs,B);
sprintf('No of GS iterations with %d as magnification factor and %d as rho= %dn',m(i),rho_gs(i),iter_gs(i));
else
iter_gs(i)=nan;
sprintf('solution is unstable for rho_gs =%d with %d as diagonal magnification factorn ',m(i),rho_gs(i),iter_gs(i));
end
%SOR
if rho_sor(i)< 1
[iter_sor(i),x_sor]=iter_sln(sor_I,X_sor,B*omega);
sprintf('No of sor iterations with %d as magnification factor and %d as rho= %dn',m(i),rho_sor(i),iter_sor(i));
else
iter_sor(i)=nan;
sprintf('solution is unstable for rho_sor =%d with %d as diagonal magnification factorn ',m(i),rho_sor(i),iter_sor(i));
end
end
% plotting
figure(1)
subplot(2,1,1)
plot(m,rho_jac,'-o','linewidth',2,'markerfacecolor','r','markeredgecolor','r');
grid on;
ylabel('rho')
xlabel('Diagonal magnifier')
title('Jacobi method','fontsize',15);
subplot(2,1,2)
plot(rho_jac,iter_jac,'-ro','linewidth',2,'markerfacecolor','k','markeredgecolor','k');
grid on;
xlabel('rho');
ylabel('no of iterations');
saveas(gca,"Jacobi_rho.jpg");
figure(2)
subplot(2,1,1)
plot(m,rho_gs,'-o','markerfacecolor','r','markeredgecolor','r','linewidth',2);
ylabel('rho')
xlabel('Diagonal magnifier')
title('Gauss Seidel method','fontsize',15);
subplot(2,1,2)
plot(rho_gs,iter_gs,'-ro','markerfacecolor','k','markeredgecolor','k','linewidth',2);
xlabel('rho');
ylabel('no of iterations');
saveas(gca,"GS_rho.jpg");
figure(3)
subplot(2,1,1)
plot(m,rho_sor,'-o','markerfacecolor','r','markeredgecolor','r','linewidth',2);
ylabel('rho')
xlabel('Diagonal magnifier')
title('SOR method','fontsize',15);
subplot(2,1,2)
loglog(rho_sor,iter_sor,'-ro','markerfacecolor','k','markeredgecolor','k','linewidth',2);
xlabel('rho');
ylabel('no of iterations');
saveas(gca,"SOR_rho.jpg");
FUNCTION CODE
function[iter,x_new]= iter_sln(Z,X,B)
% Z is the iteration matrix
error=7e7;
tol=1e-4;
x=[0; 0; 0];
iter=0;
x_old=x;
while error> tol
x = Z*x_old+ inv(X)*B;
error=max(abs(x_old-x));
x_old=x;
iter=iter+1;
end
x_new=x;
end
DIAGONAL MAGNIFIER
m =
0.5000 1.0000 2.0000 3.0000 4.0000 6.0000
SPECTRAL RADIUS FOR JACOBI METHOD
rho_jac
rho_jac =
1.0204 0.5102 0.2551 0.1701 0.1276 0.0850
NO OF ITERATIONS FOR JACOBI CONVERGENCE
iter_jac =
NaN 18 9 7 6 5
SPECTRAL RADIUS FOR GAUSS SEIDEL METHOD
rho_gs
rho_gs =
1.3815 0.3276 0.0859 0.0408 0.0244 0.0120
NO OF ITERATIONS FOR GS CONVERGENCE
iter_gs =
NaN 12 6 5 5 4
SPECTRAL RADIUS FOR SOR METHOD
rho_sor
rho_sor =
4.4467 1.6343 0.9165 0.7500 0.6779 0.6125
NO OF ITERATIONS FOR SOR CONVERGENCE
iter_sor =
NaN NaN 117 35 26 20
Below are the plots obtained for the three Iterative methods
INFERENCE:
- As the diagonal elements multiplied by the magnification factors increase, the spectral radius for all iterative methods decrease,resulting in a stable solution.
- With subsequent increase in spectral radius, the number of iterations required to reach convergence increases drastically, thereby causing errors in the simulation.
- In the given problem , for spectral radius greater than 1, the solution does not converge for all three iterative methods .