Solving the set of linear equations using Iterative solvers

Eigen values of Iteration matrix :

From Jacobi method :

-0.0976 + 0.0000i
0.0488 - 0.5079i
0.0488 + 0.5079i

From Gauss-Siedel method :

0
-0.3276
0.0388

From SOR method :

-0.2918 + 0.0000i
-1.9209 - 0.0000i
-0.3854 + 0.0000i

Iteration matrix :

Jacobi iteration matrix:

0	-0.2000	-0.4000
0.3333	0	-0.4444
0.1429	0.2857	0

Gauss-Siedel iteration matrix:

0	-0.2000	-0.4000
0	-0.0667	-0.5778
0	-0.0476	-0.2222

SOR iteration matrix:

-0.6000	-0.3200	-0.6400
-0.3200	-0.7707	-1.0524
-0.2834	-0.4254	-1.2274

Spectral radius :

rho_jac = 0.5102

rho_gs = 0.3276

rho_sor = 1.9209

Function for calculating eigenvalues :

function eigen_values = eigen_value_calc(A)

% Indentity matrix for Eigen value calculations
I = eye(3);

% J represents lamda
syms J

% Solving for eigen values
c = A - J.*I;
determinant = det(c);
coeff_matrix = coeffs(determinant,\'All\');
solution = roots(coeff_matrix);
eigen_values = double(solution);

end

Jacobi method of iteration:

function [rho_jac,iter] = jacobi_method(A)

% Reducing the coefficient matrix
D = diag(diag(A));
L = tril(A-D);
U = triu(A-D);

% Solution array input
b = [10;-14;33];

x = [0;0;0];
iter = 1;

% Error conditions
residual = 100;
tol = 1e-4;

% Eigen value calculator
P_jac = inv(D)*(-L-U);
lamda_jac =  eigen_value_calc(P_jac);
result = max(abs(lamda_jac));
rho_jac = double(result);

% Calculation for final solution
while (residual > tol)
    x = inv(D)*(b - (L+U)*x);
    iter = iter + 1;
    residual = max(abs(b - A*x));
end
end

Gauss Siedal method of iteration

function [rho_gs,iter] = gs_method(A)

% Reducing the coefficient matrix
D = diag(diag(A));
L = tril(A-D);
U = triu(A-D);

% Solution array input
b = [10;-14;33];

x = [0;0;0];
iter = 1;

% Error conditions
residual = 100;
tol = 1e-4;

% Eigen value calculator
P_gs = inv(D+L) * (-U);
lamda_gs = eigen_value_calc(P_gs);
result = max(abs(lamda_gs));
rho_gs = double(result);

% Calculation
while residual > tol
        x = inv(D+L)*(b - (U)*x);
        iter = iter + 1;
        residual = max(abs(b - A*x));
end
end

SOR method :

function [rho_sor,iter] = sor_method(A)

% Reducing the coefficient matrix
D = diag(diag(A));
L = tril(A-D);
U = triu(A-D);

% Solution array input
b = [10;-14;33];

% Initial conditions
x = [0;0;0];
iter = 1;

% Relaxation Factor
w = 1.6;


% Error conditions
residual = 100;
tol = 1e-4;

% Eigen value calculator
P_sor = inv(L+(D/w))*((D/w) - D - U);
lamda_sor = eigen_value_calc(P_sor);
result = max(abs(lamda_sor));
rho_sor = double(result);

% Calculation

while residual > tol
        x = inv(D + w*L)*(w*b - ((w*U + (w-1)*D)*x));
        iter = iter + 1;
        residual = max(abs(b - A*x)); 
end

end

Code for calculation of Iteration and Spectral radius :

close all
clear all
clc

% Coefficient matrix input
A = [5 1 2;-3 9 4;1 2 -7];
[m, ~] = size(A);

% Method of iteration
method =3;


% Solution array input
b = [10;-14;33];

% Jacobi method
if method == 1
    [rho_jac,iter] = jacobi_method(A);
end


% Gauss-Seidel method
if method == 2
    [rho_gs,iter] = gs_method(A);
end

% SOR method
if method == 3
    [rho_sor,iter] = sor_method(A);
end



% Element wise Jacobian
if method == 4
     
    x = zeros(3,1);
    x_old = x;
    D = diag(diag(A));
    L = tril(A-D);
    U = triu(A-D);
    b = [10;-14;33];
    residual = 100;
    tol = 1e-4;
    iter = 1;
    
    while (residual > tol)
        for i = 1:m
            x(i) = 1/A(i,i)*(b(i) - A(i,1:i-1)*x_old(1:i-1) - A(i,i+1:m)*x_old(i+1:m));
        end
        iter = iter + 1;
        residual = max(abs(b - A*x));
        x_old = x;
    end
end


% Element wise Gauss Siedel
if method == 5
    
    x = zeros(3,1);
    x_old = x;
    D = diag(diag(A));
    L = tril(A-D);
    U = triu(A-D);
    b = [10;-14;33];
    residual = 100;
    tol = 1e-4;
    iter = 1;
    
    while (residual > tol)
        for i = 1:m
            x(i) = 1/A(i,i)*(b(i) - A(i,1:i-1)*x(1:i-1) - A(i,i+1:m)*x_old(i+1:m));
        end
        iter = iter + 1;
        residual = max(abs(b - A*x));
        x_old = x;
    end
end

% Element wise SOR
if method == 6
    
    x = zeros(3,1);
    x_old = x;
    D = diag(diag(A));
    L = tril(A-D);
    U = triu(A-D);
    b = [10;-14;33];
    residual = 100;
    tol = 1e-4;
    iter = 1;
    w=1.6;
    
    while (residual > tol)
        for i = 1:m
            x(i) = (1-w)*x_old(i) + (w/A(i,i))*(b(i) - A(i,1:i-1)*x(1:i-1) - A(i,i+1:m)*x_old(i+1:m));
        end
        iter = iter + 1;
        residual = max(abs(b - A*x));
        x_old = x;
    end
end

For increasing the spectral radius :

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

Increased spectral radius :

rho_jac = 1.1163

rho_gs = 1.3333

rho_sor = 2.6483

The condition for convergence is \'spectral radius should be less than 1\'.So, after increasing the spectral radius to more than 1, the solution doesn\'t converge to a particular solution. The residual term becomes NaN (not a number), and iterations increase to quite a large value in the range of 1000\'s.