Importance of Spectral Radius in linear systems
In this report, the emphasis has been laid about solving the system of linear equations using several iterative methods. The methods implemented are :
1.) Jacobi
2.) Gauss-Seidel
3.) Successive Overrelaxation ( SOR )
But first of all, let's understand why do we need iterative methods to solve linear equations. Consider a Linear system of equations of the form
AX = B
where A is the coefficient matrix, X is the unknown vector, and B is the right-hand side function. One way to solve these types of linear equations is with Direct Methods which include :
1.) Cramers rule
2.) Gaussian Elimination
However, if the matrix is of billions of elements like in the engineering scale problems, it is computationally very expensive to solve these systems using the direct methods. One indirect way to solve these systems using the iterative techniques which uses an initial guess and then iterate it over and over again until the two consecutive iterative differ not more than the tolerance, i.e. we say the solution has converged up to the tolerance level if,
`norm(u^{k+1} - u^k) \leq \epsilon`
where `normu` is the standard defined norm
`u^k` is the solution at the `k^(th)` iterative step
`epsilon` is the specified tolerance limit.
However, these Iterative methods themselves aren't that easy to set up. One very important thing to keep in mind while setting these matrices and the iterative methods is that matrix A must be diagonally dominant or should satisfy the Scarborough criteria i.e.
`|a_{ii}| \geq sum_(j \nei) |a_{ij}|`
meaning that the absolute value of the diagonal element in a particular row must be greater than the absolute value of all the other remaining elements. One more way to think about this is that if A is a positive definite symmetric matrix and if the following condition is satisfied for the iteration matrix R (which can be any of the methods discussed above) then the approximation will converge to the exact solution:
`\rho(R)<1` , where `\rho(R)` is the spectral radius of R
Also, an important fact to note here is that as the spectral radius of A goes near to 0 it will accelerate the convergence. Now, with a basic understanding of the mathematics behind iterative methods we are ready to develop a general form of the iterative methods :
Consider the following basic form of an iteration :
`u^(k+1) = Ru^k +g`
where `u^(k+1)` is the value at the `(k+1)^(th)` iteration and `u^k` is the value at the `k^(th)` iteration and the R is the iteration matrix according to the scheme. According to the different types of schemes the matrix R can have different versions depending upon the method employed.
Consider the system of linear equations in the form:
`\begin{cases} a_1x\: +\:b_1y\:+\:c_1z = d_1 \\a_2x\: +\:b_2y\:+\:c_2z = d_2 \\ a_3x\: +\:b_3y\:+\:c_3z = d_3\end{cases}`
In matrix form these can be written as :
`[[a_1,b_1,c_1],[a_2,b_2,c_2],[a_3,b_3,c_3]]((x_1),(x_2),(x_3)) = ((d_1),(d_2),(d_3))`
Forming an Iteration it results in :
`\begin{cases} \\x_1 = \frac{1}{a_1}[d_1 - b_1x_2 - c_1x_3] \\\\ x_2 = \frac{1}{b_2}[d_2 - a_2x_1 - c_2x_3] \\ \\ x_3 = \frac{1}{c_3}[d_3 - a_3x_1 - b_3x_2] \\\\ \end{cases}`
As per the type of the iteration used the right hand side can be converted to a matrix system as follows:
JACOBI ITERATION
Consider the system of linear equations in the matrix form AX=B.
Let A be decomposed as A = D - L - U where D = Diagonal Matrix, L = Strictly Lower Triangular Matrix, and U = Strictly Upper Triangular Matrix.
We can write the equation in the form :
`(D-L-U)u = F`
`Du = (L+U)u + F`, taking the inverse on both sides we get :
`u = D^(-1)(L+U)u + D^(-1)F`
which can be further simplified to :
`u = R_ju+g`
where `R_j` is the JACOBI ITERATION matrix. In iteration form this can be written as :
`u^(k+1) = R_ju^k + g`
This method will converge only if `rho(R_j) < 1`i.e. the spectral radius of the iteration matrix is less than one. This important fact is very useful in determining the convergence of the iterative scheme.
GAUSS - SEIDEL ITERATION
A similar analysis can be performed for the Gauss-Seidel Iteration. The only difference between the Gauss-Siedel and Jacobi method is the values are used as soon as they are available for the next iteration in Gauss-Siedel as compared to Jacobi where all the values are first calculated and then updated for the next iteration.
We can write :
`(D-L-U)u = F`
`(D-L)u = Uu + F`
`u = (D-L)^(-1)Uu + (D-L)^(-1)F`
which can be simplified to :
`u = R_Gu+g`
where `R_G` is the GAUSS-SEIDEL ITERATION matrix. In Iteration form this can be written as :
`u^(k+1) = R_Gu^k + g`
Similarly this method will converge only if `rho(R_G) < 1`.
SUCCESSIVE OVER-RELAXATION
A more elegant modification can be made to basic iterative techniques. The new iteration can be taken as a weighted average of the new and the old approximation according to the relaxation factor `alpha`.
In equation form :
`u^(k+1) = (alpha)(R_Gu^k + g) + (1-alpha)u^k`
simplifying this a little bit more we get :
`u^(k+1) = [(1-alpha)I + alphaR_G]u^k + alphag`
where the term `[(1-alpha)I + alphaR_G]` can be labelled as a new matrix `R_(alpha)`, so that the iteration can be written as :
`u^(k+1) = R_(alpha)u^k + alphag`
The matrix system given in the question is :
`[[5,1,2],[-3,9,4],[1,2,-7]]((x_1),(x_2),(x_3)) = ((10),(-14),(33))`
To create the Jacobi, Gauss-Seidel, and SOR iteration matrices the following codes are implemented :
To Create the Jacobi Iteration Matrix Following code is Implemented in MATLAB:
function [R_j,D] = create_jacobi_iteration_matrix(A)
%
% This function will take as input a matrix and the weighing parameter and
% will output the weighted Jacobi Iteration matrix
%
% The Matrix A can be written as A = D - L - U
%
% where , D = Diagonal Matrix
% L = lower triangular matrix
% U = Upper Triangular matrix
%
[a,b] = size(A);
if a~=b
R_j = 0;
error('Input a square matrix!!');
else
D = zeros(a);
L = zeros(a);
U = zeros(a);
% Assigning the values to the D matrix
for ii = 1:a
D(ii,ii) = A(ii,ii);
end
% Assigning the value to the L matrix
for ii = 1:a
for jj = 1:ii
L(ii,jj) = A(ii,jj);
end
end
L = -1*(L - D);
% Assigning the value to the L matrix
for ii = 1:a
for jj = ii:a
U(ii,jj) = A(ii,jj);
end
end
U = -1*(U - D);
% creating the weighted Jacobi Matrix
R_j = inv(D) * ( L + U );
end
end
Explanation: The function takes in input as the coefficient matrix A and gives the output two matrices
1.) The Jacobi Iteration Matrix
2.) The Diagonal Matrix Corresponding to A
To Create the Gauss-Seidel Iteration Matrix Following code is Implemented in MATLAB:
function [R_G , D , L] = create_gauss_seidel_matrix(A)
%
% This function will take as input a matrix and the weighing parameter and
% will output the weighted Jacobi Iteration matrix
%
% The Matrix A can be written as A = D - L - U
%
% where , D = Diagonal Matrix
% L = lower triangular matrix
% U = Upper Triangular matrix
%
[a,b] = size(A);
if a~=b
R_G = 0;
error('Input a square matrix');
else
D = zeros(a);
L = zeros(a);
U = zeros(a);
% Assigning the values to the D matrix
for ii = 1:a
D(ii,ii) = A(ii,ii);
end
% Assigning the value to the L matrix
for ii = 1:a
for jj = 1:ii
L(ii,jj) = A(ii,jj);
end
end
L = -1*(L - D);
% Assigning the value to the L matrix
for ii = 1:a
for jj = ii:a
U(ii,jj) = A(ii,jj);
end
end
U = -1*(U - D);
% Creating the Gauss-Seidel Matrix
R_G = inv(D-L)*U;
end
Explanation: The function takes in input as the coefficient matrix A and gives the output two matrices
1.) The Gauss-Seidel Iteration Matrix Rg
2.) The Diagonal Matrix Corresponding to A, D
3.) The lower Triangular matrix corresponding to A, L which will be useful for the iteration.
To Create the Successive-Over_Relaxation Matrix for Jacobi Method using a parameter `omega`, Following code is Implemented in MATLAB:
function [R_w , D] = create_weighted_jacobi(A,w)
%
% This function will take as input a matrix and the weighing parameter and
% will output the weighted Jacobi Iteration matrix
%
% The Matrix A can be written as A = D - L - U
%
% where , D = Diagonal Matrix
% L = lower triangular matrix
% U = Upper Triangular matrix
%
[a,b] = size(A);
if a~=b
R_w = 0;
error('Input a square matrix');
else
D = zeros(a);
L = zeros(a);
U = zeros(a);
% Assigning the values to the D matrix
for ii = 1:a
D(ii,ii) = A(ii,ii);
end
% Assigning the value to the L matrix
for ii = 1:a
for jj = 1:ii
L(ii,jj) = A(ii,jj);
end
end
L = -1*(L - D);
% Assigning the value to the L matrix
for ii = 1:a
for jj = ii:a
U(ii,jj) = A(ii,jj);
end
end
U = -1*(U - D);
I = eye(a);
% creating the weighted Jacobi Matrix
R_w = ( 1 - w ) * I + w * inv(D) * ( L + U );
end
end
Explanation: The function takes in input as the coefficient matrix A and the relaxation parameter `omega` and gives the output two matrices:
1.) The Weighted-Jacobi Iteration Matrix
2.) The Diagonal Matrix Corresponding to A, D
To Create the Successive-Over_Relaxation Matrix for Gauss-Seidel Method using a parameter `omega`, Following code is Implemented in MATLAB:
function [R_G , D , L] = create_weighted_gauss_seidel_matrix(A,w)
%
% This function will take as input a matrix and the weighing parameter and
% will output the weighted Jacobi Iteration matrix
%
% The Matrix A can be written as A = D - L - U
%
% where , D = Diagonal Matrix
% L = lower triangular matrix
% U = Upper Triangular matrix
%
[a,b] = size(A);
if a~=b
R_G = 0;
error('Input a square matrix');
else
D = zeros(a);
L = zeros(a);
U = zeros(a);
% Assigning the values to the D matrix
for ii = 1:a
D(ii,ii) = A(ii,ii);
end
% Assigning the value to the L matrix
for ii = 1:a
for jj = 1:ii
L(ii,jj) = A(ii,jj);
end
end
L = -1*(L - D);
% Assigning the value to the L matrix
for ii = 1:a
for jj = ii:a
U(ii,jj) = A(ii,jj);
end
end
U = -1*(U - D);
I = eye(a);
% Creating the Gauss-Seidel Matrix
R_G = (1-w)*I + w*inv(D-L)*U;
end
Explanation: The function takes in input as the coefficient matrix A and the relaxation parameter `omega` and gives the output two matrices
1.) The Weighted-Gauss-Seidel Iteration Matrix Rg
2.) The Diagonal Matrix Corresponding to A, D
3.) The lower Triangular matrix corresponding to A, L which will be useful for the iteration.
For this report, the relaxation parameter is chosen as `omega = 1.3`
Now after the creation of functions of the iteration matrices we now find each type of iteration matrix for A. The following code is implemented in MATLAB.
clear;
close all;
clc;
% create the coefficient matrix A
A = [ 5 1 2 ; -3 9 4 ; 1 2 -7];
% Defining the right hand side matrix B
B = [ 10 , -14 , 33];
% Create the iteration matrices corresponding to A for all of the following
% methods :
%
% 1.) Jacobi Method
% 2.) Gauss-Seidel Method
% 3.) Jacobi using SOR Method
% 4.) Gauss-Seidel using SOR Method
%
% Already the functions have been created for all the methods discussed
% above
%
% Defining the relaxation parameter w
w = 1.3;
% Jacobi Matrix
[ A_jacobi , D ] = create_jacobi_iteration_matrix(A);
% Gauss-Seidel Matrix
[ A_gauss_seidel , D , L] = create_gauss_seidel_matrix(A) ;
% SOR matrix using the Jacobi Method
[ A_sor_jacobi , D ] = create_weighted_jacobi(A,w);
% SOR matrix using the Gauss-Seidel Method
[ A_sor_gauss_seidel , D , L] = create_weighted_gauss_seidel_matrix(A,w);After the creation of the iteration matrices, we turn to MATLAB to find the Eigenvalues and the spectral radius of the Iteration Matrices. For this, the following code is implemented in MATLAB:
function [ eigenvalues , spectral_radius ] = eig_and_spectral(A)
%
% This function takes in input a matrix and computes its eigen values and
% the spectral radius
%
% First of all we need to check for the square matrix
[a,b] = size(A);
if a~=b
% print an error message
error(' Please input a Square Matrix!!');
else
% Import the symbolic toolbox
syms lambda
% create the equation of the eigenvalue
eqn = det(A-lambda*eye(a));
% Finding the coefficients of the equation
p = sym2poly(eqn);
% solve the equation using an inbuilt
r = roots(p);
% Finding the values in terms of numbers
eigenvalues = vpa(r,6);
% Spectral Radius
spectral_radius = vpa(max(abs(eigenvalues)),3);
end
end
Now with the creation of the function to find the eigenvalues and the spectral radius we now compute the eigenvalues and the spectral radius of the iteration matrices as shown below:
For Jacobi Matrix:
% For Jacobi Matrix
[ eigenvalues_jacobi , spectral_radius_jacobi ] = eig_and_spectral(A_jacobi);The following output is obtained:
For Gauss-Seidel Matrix:
% For Gauss-Seidel Matrix
[ eigenvalues_gauss_seidel , spectral_radius_gauss_seidel ] = eig_and_spectral(A_gauss_seidel);The following output is obtained:
For SOR using Gauss-Seidel Matrix:
% For SOR using Gauss-Seidel
[ eigenvalues_gauss_seidel_sor , spectral_radius_gauss_seidel_sor ] = eig_and_spectral(A_sor_gauss_seidel);
The following output is obtained:
For SOR using Jacobi Matrix:
% For SOR using Jacobi
[ eigenvalues_Jacobi_sor , spectral_radius_Jacobi_sor ] = eig_and_spectral(A_sor_jacobi);The following output is obtained:
Now with the creation of the iteration matrices, we move to solve the linear system using iterative methods. We will be using functions to solve the equation so that the implementations are robust.
Using the Jacobi Method :
The following MATLAB function is implemented to solve the system of equations using the Jacobi Method :
function [ number_of_iter , final_solution ] = iteration_using_jacobi(B,D,A_jacobi)
%
% This function will take in as three inputs
%
% B = right hand-side of the equation AX = B
% D = Diagonal Matrix of A
% A_jacobi = Jacobi iteration matrix for this system
%
% the function will return two outputs
%
% 1.) number of iterations it took to converge the solution
% 2.) the converged solution
%
% Defining the initial Guess
initial_guess = zeros(length(B),1);
% Defining the count of the iteration
number_of_iter = 0;
%
% Specifying the initial tolerance and error
error = 9e9;
tol = 1e-4;
% Main convergence loop
while(error>tol)
% Perform relaxation
x_new = A_jacobi * initial_guess(:) + inv(D) * B(:);
% updating the error
error = max(max(abs(x_new-initial_guess)));
% updating the initial guess
initial_guess = x_new;
% updating the loop count
number_of_iter = number_of_iter + 1;
end
% Return the final converged solution within the tolerance limit
final_solution = x_new;
endThe following inputs are provided in the command prompt window of MATLAB to execute the function using our iteration matrices as described:
% Using Jacobi
[ number_of_iter_jacobi , final_solution_jacobi ] = iteration_using_jacobi(B,D,A_jacobi);The following output is obtained regarding the same:
As it is evident it took 18 Jacobi Iterations for the solution to converge. Now we shall implement the same problem but this time using the Gauss-Seidel Iteration.
Using the Gauss-Seidel Method :
The following MATLAB function is implemented to solve the system of equations using the Gauss-Seidel Method :
function [ number_of_iter , final_solution ] = iteration_using_gauss_seidel(B,D,L,A_gauss_seidel)
%
% This function will take in as three inputs
%
% B = right hand-side of the equation AX = B
% D = Diagonal Matrix of A
% L = Lower Triangular Matrix for A
% A_gauss_seidel = Gauss-Seidel iteration matrix for this system
%
% the function will return two outputs
%
% 1.) number of iterations it took to converge the solution
% 2.) the converged solution
%
% Defining the initial Guess
initial_guess = zeros(length(B),1);
% Defining the count of the iteration
number_of_iter = 0;
%
% Specifying the initial tolerance and error
error = 9e9;
tol = 1e-4;
% Main convergence loop
while(error>tol)
% Perform relaxation
x_new = A_gauss_seidel * initial_guess(:) + inv(D-L) * B(:);
% updating the error
error = max(max(abs(x_new-initial_guess)));
% updating the initial guess
initial_guess = x_new;
% updating the loop count
number_of_iter = number_of_iter + 1;
end
% Return the final converged solution within the tolerance limit
final_solution = x_new;
endThe following inputs are provided in the command prompt window of MATLAB to execute the function using our iteration matrices as described:
% using Gauss-Seidel
[ number_of_iter_GS , final_solution_GS ] = iteration_using_gauss_seidel(B,D,L,A_gauss_seidel);The following output is obtained regarding the same:
As it is evident it took 12 Gauss-Seidel Iterations for the solution to converge. Now we shall implement the same problem but this time using the SOR method.
Using the SOR Jacobi Method :
The following MATLAB function is implemented to solve the system of equations using the SOR Jacobi Method :
function [ number_of_iter , final_solution ] = iteration_using_sor_jacobi(B,D,A_sor_jacobi,w)
%
% This function will take in as three inputs
%
% B = right hand-side of the equation AX = B
% D = Diagonal Matrix of A
% A_sor_jacobi = Jacobi iteration matrix for this system
% w = relaxation paramter
%
% the function will return two outputs
%
% 1.) number of iterations it took to converge the solution
% 2.) the converged solution
%
% Defining the initial Guess
initial_guess = zeros(length(B),1);
% Defining the count of the iteration
number_of_iter = 0;
%
% Specifying the initial tolerance and error
error = 9e9;
tol = 1e-4;
% Main convergence loop
while(error>tol)
% Perform relaxation
x_new = A_sor_jacobi * initial_guess(:) + w * inv(D) * B(:);
% updating the error
error = max(max(abs(x_new-initial_guess)));
% updating the initial guess
initial_guess = x_new;
% updating the loop count
number_of_iter = number_of_iter + 1;
end
% Return the final converged solution within the tolerance limit
final_solution = x_new;
endThe following inputs are provided in the command prompt window of MATLAB to execute the function using our iteration matrices as described:
% Using the SOR jacobi Method
[ number_of_iter_sor_jacobi , final_solution_sor_jacobi ] = iteration_using_sor_jacobi(B,D,A_sor_jacobi,w);
The following output is obtained regarding the same:
Using the SOR Gauss-Seidel Method :
The following MATLAB function is implemented to solve the system of equations using the SOR Gauss-Seidel Method :
function [ number_of_iter_sor_GS , final_solution_sor_GS ] = ...
iteration_using_sor_gauss_seidel(B,D,L,A_sor_gauss_seidel,w)
%
% This function will take in as three inputs
%
% B = right hand-side of the equation AX = B
% D = Diagonal Matrix of A
% L = Lower Triangular Matrix for A
% A_sor_gauss_seidel = the SOR Gauss-Seidel iteration matrix for this system
% w = relaxation parameter
%
% the function will return two outputs
%
% 1.) number of iterations it took to converge the solution
% 2.) the converged solution
%
% Defining the initial Guess
initial_guess = zeros(length(B),1);
% Defining the count of the iteration
number_of_iter_sor_GS = 0;
%
% Specifying the initial tolerance and error
error = 9e9;
tol = 1e-4;
% Main convergence loop
while(error>tol)
% Perform relaxation
x_new = A_sor_gauss_seidel * initial_guess(:) + w * inv(D-L) * B(:);
% updating the error
error = max(max(abs(x_new-initial_guess)));
% updating the initial guess
initial_guess = x_new;
% updating the loop count
number_of_iter_sor_GS = number_of_iter_sor_GS + 1;
end
% Return the final converged solution within the tolerance limit
final_solution_sor_GS = x_new;
endThe following inputs are provided in the command prompt window of MATLAB to execute the function using our iteration matrices as described:
% using the SOR Gauss-Seidel Method
[ number_of_iter_sor_GS , final_solution_sor_GS ] = ...
iteration_using_sor_gauss_seidel(B,D,L,A_sor_gauss_seidel,w)The following output is obtained regarding the same:
A very subtle point to note here is that if we use the relaxation parameter as `omega = 1.1` then we get the following result:
That is by reducing the relaxation parameter we get the convergence faster. This can be explained due to the fact :
In the SOR technique in pushing the approximation towards the solution, it may sometimes push it a little forward and hence the approximation will try to come back to the exact solution and it will keep on oscillating around the solution till the point the frequencies in the error gets damped. It can be though of much like a potential mountain where a particle left from above will settle down at equilibrium ( i.e. the valley ) only if it has the right momentum, if it has too excessive momentum it will bounce off and will never return to the valley. If it has too low a momentum it will take a longer time to come to the valley. If it has slightly more momentum than normal than the particle will oscillate about the valley before stopping at the equilibrium. Much of what described above happens in the iterative techniques. Therefore a careful selection of the relaxation parameter `omega` is a must to avoid over-shooting.
Jacobi - Iteration ( Effect of Different Spectral Radius on Convergence )
Now we will be changing the spectral radius of the iteration matrices to different numbers and will see how the number of iterations for convergence is affected. One way to increase the spectral radius is by magnifying the diagonals of the matrix. The following code is implemented in MATLAB to study the effect of the increase of Spectral Radius on the number of iterations required for the Jacobi matrix to converge.
clear;
close all;
clc;
% create the coefficient matrix A
A = [ 5 1 2 ; -3 9 4 ; 1 2 -7];
% Defining the right hand side matrix B
B = [ 10 , -14 , 33];
% Magnification factors
magnify_diagonal = [0.999 0.9499 0.9099 0.88 0.80 0.75 0.70 0.60 0.50];
% Creating a vector storing all the values for the iteration numbers for
% different magnifications
%
% For jacobi Iteration:
for ii = 1:length(magnify_diagonal)
% The coefficient matrix
A = [ 5 1 2 ; -3 9 4 ; 1 2 -7];
% Changing the diagonal terms
for jj = 1:max(size(A))
A(jj,jj) = A(jj,jj)*magnify_diagonal(ii);
end
% Creating the Jacobi Iteration Matrix
[ A_jacobi , D ] = create_jacobi_iteration_matrix(A);
% Finding the spectral radius
[ ~ , spectral_radius_jacobi ] = eig_and_spectral(A_jacobi);
if spectral_radius_jacobi < 1
% Using Jacobi relaxation
[ number_of_iter_jacobi , ~ ] = iteration_using_jacobi(B,D,A_jacobi);
else
% Assigning the value 0 as they won't converge
number_of_iter_jacobi = 0;
end
% Storing the spectral radius in a vector spec_jacobi
spec_jacobi_different_magnifications(ii) = spectral_radius_jacobi;
% Number of iterations to converge
no_jacobi_iterations_to_converge(ii) = number_of_iter_jacobi;
end
% Plotting the figures
figure(1);clf;
subplot(211)
plot(magnify_diagonal,spec_jacobi_different_magnifications,'k-','linew',2)
set(gca,'fontsize',20);
xlabel('Magnification of the diagonal');
ylabel('Spectral radius');
title('comparison of spectral radius for different diagonal magnification');
subplot(212)
semilogy(spec_jacobi_different_magnifications,no_jacobi_iterations_to_converge,'r-','linew',2)
set(gca,'fontsize',20);
xlabel('Spectral Radius of the Jacobi Iteration Matrix');
ylabel('Jacobi Iterations');
title('comparison of number of Jacobi Iterations for different diagonal magnification');
The code modifies the diagonal terms of A and then calculates the spectral radius for the jacobi_iteration matrix and then computes the number of iterations required to converge to a final solution.
The following plot was obtained:
As is clear from the graph above that as the magnification factor is made less, the spectral radius increases, and the number of iterations required to converge to the solution also increases which is in agreement with the theoretical aspect of the iteration matrices.
Gauss-Seidel - Iteration ( Effect of Different Spectral Radius on Convergence )
Now we will be changing the spectral radius of the iteration matrices to different numbers and will see how the number of iterations for convergence is affected. One way to increase the spectral radius is by magnifying the diagonals of the matrix. The following code is implemented in MATLAB to study the effect of the increase of Spectral Radius on the number of iterations required for the Gauss-Seidel matrix to converge.
clear;
close all;
clc;
% create the coefficient matrix A
A = [ 5 1 2 ; -3 9 4 ; 1 2 -7];
% Defining the right hand side matrix B
B = [ 10 , -14 , 33];
% Magnification factors
magnify_diagonal = [0.999 0.9499 0.9099 0.88 0.80 0.75 0.70 0.60 0.50];
% Creating a vector storing all the values for the iteration numbers for
% different magnifications
%
% For jacobi Iteration:
for ii = 1:length(magnify_diagonal)
% The coefficient matrix
A = [ 5 1 2 ; -3 9 4 ; 1 2 -7];
% Changing the diagonal terms
for jj = 1:max(size(A))
A(jj,jj) = A(jj,jj)*magnify_diagonal(ii);
end
% Gauss-Seidel Matrix
[ A_gauss_seidel , D , L] = create_gauss_seidel_matrix(A) ;
% For Gauss-Seidel Matrix
[ eigenvalues_gauss_seidel , spectral_radius_gauss_seidel ] = eig_and_spectral(A_gauss_seidel);
if spectral_radius_gauss_seidel < 1
% Using Jacobi relaxation
[ number_of_iter_GS , ~ ] = iteration_using_gauss_seidel(B,D,L,A_gauss_seidel);
else
% Assigning the value 0 as they won't converge
number_of_iter_GS = 0;
end
% Storing the spectral radius in a vector spec_jacobi
spec_GS_different_magnifications(ii) = spectral_radius_gauss_seidel;
% Number of iterations to converge
no_GS_iterations_to_converge(ii) = number_of_iter_GS;
end
% Plotting the figures
figure(1);clf;
subplot(211)
plot(magnify_diagonal,spec_GS_different_magnifications,'k-','linew',2)
set(gca,'fontsize',20);
xlabel('Magnification of the diagonal');
ylabel('Spectral radius');
title('comparison of spectral radius for different diagonal magnification');
subplot(212)
semilogy(spec_GS_different_magnifications,no_GS_iterations_to_converge,'r-','linew',2)
set(gca,'fontsize',20);
xlabel('Spectral Radius of the Iteration Matrices');
ylabel('Gauss-Seidel Iterations');
title('comparison of number of Gauss-Seidel Iterations for different diagonal magnification');
The code modifies the diagonal terms of A and then calculates the spectral radius for the Gauss-Seidel iteration matrix and then computes the number of iterations required to converge to a final solution.
The following plot was obtained:
As can be seen here that for the diagonal magnification 0.55 and 0.5 the spectral radius of the Gauss-Seidel Iteration Matrix has crossed over 1. Hence the scheme would never converge. In the rest of the magnification factors, it is clear that as the spectral radius decreases the number of the iterations required also decrease
SOR Using Gauss-Seidel - Iteration ( Effect of Different Spectral Radius on Convergence )
Now we will be changing the spectral radius of the iteration matrices to different numbers and will see how the number of iterations for convergence is affected. One way to increase the spectral radius is by magnifying the diagonals of the matrix. The following code is implemented in MATLAB to study the effect of the increase of Spectral Radius on the number of iterations required for the SOR using the Gauss-Seidel matrix to converge.
clear;
close all;
clc;
% create the coefficient matrix A
A = [ 5 1 2 ; -3 9 4 ; 1 2 -7];
% Defining the right hand side matrix B
B = [ 10 , -14 , 33];
% Defining the relaxation parameter w
w = 1.3;
% Magnification factors
magnify_diagonal = [0.999 0.9499 0.9099 0.88 0.80 0.75 0.70 0.60 0.50];
% Creating a vector storing all the values for the iteration numbers for
% different magnifications
%
% For SOR using Gauss-Seidel Iteration:
for ii = 1:length(magnify_diagonal)
% The coefficient matrix
A = [ 5 1 2 ; -3 9 4 ; 1 2 -7];
% Changing the diagonal terms
for jj = 1:max(size(A))
A(jj,jj) = A(jj,jj)*magnify_diagonal(ii);
end
% SOR matrix using the Gauss-Seidel Method
[ A_sor_gauss_seidel , D , L] = create_weighted_gauss_seidel_matrix(A,w);
% For SOR using Gauss-Seidel
[ ~ , spectral_radius_gauss_seidel_sor ] =...
eig_and_spectral(A_sor_gauss_seidel);
if spectral_radius_gauss_seidel_sor < 1
% Using SOR Gauss-Seidel Iteration
[ number_of_iter__sor_GS , ~ ] = ...
iteration_using_sor_gauss_seidel(B,D,L,A_sor_gauss_seidel,w);
else
% Assigning the value 0 as they won't converge!!
number_of_iter__sor_GS = 0;
end
% Storing the spectral radius in a vector spec_jacobi
spec__sor_GS_different_magnifications(ii) = spectral_radius_gauss_seidel_sor;
% Number of iterations to converge
no__sor_GS_iterations_to_converge(ii) = number_of_iter__sor_GS;
end
% Plotting the figures
figure(1);clf;
subplot(211)
plot(magnify_diagonal,spec__sor_GS_different_magnifications,'k-','linew',2)
set(gca,'fontsize',20);
xlabel('Magnification of the diagonal');
ylabel('Spectral radius');
title('comparison of spectral radius for different diagonal magnification');
subplot(212)
semilogy(spec__sor_GS_different_magnifications,no__sor_GS_iterations_to_converge,'r-','linew',2)
set(gca,'fontsize',20);
xlabel('Spectral Radius of the SOR Using Gauss-Seidel Iteration Matrix');
ylabel('SOR using Gauss-Seidel Iterations');
title('comparison of number of SOR using Gauss-Seidel Iterations for different diagonal magnification');
The code modifies the diagonal terms of A and then calculates the spectral radius for the SOR Using the Gauss-Seidel iteration matrix and then computes the number of iterations required to converge to a final solution.
The following plot was obtained:
As can be seen here that for the diagonal magnification 0.55. 0.5, 0.6 the spectral radius of the SOR using the Gauss-Seidel Iteration Matrix has crossed over 1. Hence the scheme would never converge. In the rest of the magnification factors, it is clear that as the spectral radius decreases the number of the iterations required to converge also decreases.
SOR Using Jacobi - Iteration ( Effect of Different Spectral Radius on Convergence )
Now we will be changing the spectral radius of the iteration matrices to different numbers and will see how the number of iterations for convergence is affected. One way to increase the spectral radius is by magnifying the diagonals of the matrix. The following code is implemented in MATLAB to study the effect of the increase of Spectral Radius on the number of iterations required for the SOR using the Jacobi matrix to converge.
clear;
close all;
clc;
% create the coefficient matrix A
A = [ 5 1 2 ; -3 9 4 ; 1 2 -7];
% Defining the right hand side matrix B
B = [ 10 , -14 , 33];
% Defining the relaxation parameter w
w = 1.3;
% Magnification factors
magnify_diagonal = [0.999 0.9499 0.9099 0.88 0.80 0.75 0.70 0.60 0.50];
% Creating a vector storing all the values for the iteration numbers for
% different magnifications
%
% For SOR using Jacobi Iteration:
for ii = 1:length(magnify_diagonal)
% The coefficient matrix
A = [ 5 1 2 ; -3 9 4 ; 1 2 -7];
% Changing the diagonal terms
for jj = 1:max(size(A))
A(jj,jj) = A(jj,jj)*magnify_diagonal(ii);
end
% SOR matrix using the Jacobi Method
[ A_sor_jacobi , D ] = create_weighted_jacobi(A,w);
% For SOR using Jacobi
[ ~ , spectral_radius_Jacobi_sor ] = ...
eig_and_spectral(A_sor_jacobi);
if spectral_radius_Jacobi_sor < 1
% Using the SOR jacobi Method
[ number_of_iter_sor_jacobi , ~ ] = ...
iteration_using_sor_jacobi(B,D,A_sor_jacobi,w);
else
% Assigning the value 0 as they won't converge!!
number_of_iter_sor_jacobi = 0;
end
% Storing the spectral radius in a vector spec_jacobi
spec__sor_Jacobi_different_magnifications(ii) = spectral_radius_Jacobi_sor;
% Number of iterations to converge
no__sor_Jacobi_iterations_to_converge(ii) = number_of_iter_sor_jacobi;
end
% Plotting the figures
figure(1);clf;
subplot(211)
plot(magnify_diagonal,spec__sor_Jacobi_different_magnifications,'k-','linew',2)
set(gca,'fontsize',20);
xlabel('Magnification of the diagonal');
ylabel('Spectral radius');
title('comparison of spectral radius for different diagonal magnification');
subplot(212)
semilogy(spec__sor_Jacobi_different_magnifications,no__sor_Jacobi_iterations_to_converge,'r-','linew',2)
set(gca,'fontsize',20);
xlabel('Spectral Radius of the SOR Using Jacobi Iteration Matrix');
ylabel('SOR using Jacobi Iterations');
title('comparison of number of SOR using Jacobi Iterations for different diagonal magnification');
The code modifies the diagonal terms of A and then calculates the spectral radius for the SOR Using the Jacobi iteration matrix and then computes the number of iterations required to converge to a final solution.
The following plot was obtained:
As can be seen here that for the diagonal magnification 0.55. 0.5, 0.6 the spectral radius of the SOR using the Jacobi Iteration Matrix has crossed over 1. Hence the scheme would never converge. In the rest of the magnification factors, it is clear that as the spectral radius decreases the number of the iterations required to converge also decreases.
Conclusion: The Spectral radius of the iteration matrix can be suitably modified by magnifying the diagonal elements of the coefficient matrix. As the Spectral radius is decreased the number of iterations required for that particular scheme to converge is also decreased. However, if the spectral radius of an iteration matrix is made higher than 1 than the system of equations will never converge with that iteration method. Hence a proper selection of the iteration method to allow faster convergence.