Understanding linear systems- Stability Analysis

Objective:

Determine:

• Eigenvalues of the iteration matrix.
• The spectral radius of the iteration matrix.
• Solve the system using Jacobi, Gauss-seidel and SOR methods
• Compute the iteration Matrix for the above mentioned iterative methods.
• Increase the Spectral radius to 1.1 by changing the diagonal terms suitably.

System of equations:

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

X=[x1;x2;x3]

B=[10 -14 33]

Introduction:

Iterative methods for AX = B

These techniques can only be applied to square linear systems (n equations in n unknowns), there are mainly 3 methods commonly used:

1. Jacobi Iterative Method
2. Gauss-Seidel Iterative Method
3. Successive Over Relaxation

Jacobi Iterative Method:

The Jacobi method is the simplest iterative method for solving a (square) linear system AX=B.

Gauss-Seidel Iterative Method:

The Gauss-Seidel method is a variant of the Jacobi method that usually improves the rate of convergence.

Successive Over Relaxation:

The Gauss-Seidel method is a variant of the Gauss-Siedel method that usually improves the rate of convergence.

All three iterative schemes can be written analytically in the compact form as:
`X_(n+1) = P*X_n + Q` ---------------(1)
where P is a constant matrix and Q is a constant vector.

To derive expressions for P and Q in terms of the input data: the coefficient matrix A and the right-hand side vector b, let's consider a system of equation AX = B:

`[[A_(11),A_(12),A_(13)],[A_(21),A_(22),A_(23)],[A_(31),A_(32),A_(33)]]*[[X],[Y],[Z]]=[[b_1],[b_2],[b_3]]` ---------------(2)

Each iteration of the Jacobi method, of Eqn(2), takes the form:

`A_(11)*X_(n+1)=b_1 - A_(12)Y_n-A_(13)Z_n`

`A_(22)*Y_(n+1)=b_2- A_(21)X_n-A_(23)Z_n`

`A_(33)*Z_(n+1)=b_3 - A_(31)X_n-A_(32)Y_n`

Or in matrix form:

`[[A_(11),0,0],[0,A_(22),0],[0,0,A_(33)]]*[[X_(n+1)],[Y_(n+1)],[Z_(n+1)]]=-[[0,A_(12),A_(13)],[A_(21),0,A_(23)],[A_(31),A_(32),0]]*[[X_n],[Y_n],[Z_n]]+[[b_1],[b_2],[b_3]]` ---------------(3)

Thus, in general, the Jacobi method can be written as:

`D*X_(n+1)=-(L+U)*X_n +b`

where D is a 3 × 3 matrix containing the diagonal part of A, and L and U are 3 × 3 matrices containing the (strictly) lower and (strictly) upper triangular parts of A

Premultiplying by `D^(-1)` gives:

`X_(n+1)=-D^(-1)*(L+U)*X_n +D^(-1)*b` ---------------(4)

Each iteration of the Gauss-Siedel method, of Eqn(2), takes the form:

`A_(11)*X_(n+1)=b_1 - A_(12)Y_n-A_(13)Z_n`

`A_(22)*Y_(n+1)=b_2- A_(21)X_(n+1)-A_(23)Z_n`

`A_(33)*Z_(n+1)=b_3 - A_(31)X_(n+1)-A_(32)Y_(n+1)`

Or in matrix form:

`[[A_(11),0,0],[0,A_(22),0],[0,0,A_(33)]]*[[X_(n+1)],[Y_(n+1)],[Z_(n+1)]]=-[[0,0,0],[A_(21),0,0],[A_(31),A_(32),0]]*[[X_(n+1)],[Y_(n+1)],[Z_(n+1)]]-[[0,A_(12),A_(13)],[0,0,A_(23)],[0,0,0]]*[[X_n],[Y_n],[Z_n]]+[[b_1],[b_2],[b_3]]` ---------------(5)

Thus, in general, the Gauss-Siedel method can be written as:

`D*X_(n+1)=-L*X_(n+1)-U*X_n +b`

Premultiplying by `D^(-1)` gives:

`X_(n+1)=-D^(-1)*L*X_(n+1)-D^(-1)*U*X_n +D^(-1)*b`

bring all terms related to `X_(n+1)` to the left-hand side and then multiplying whole equation by `(D+L)^(-1)`, we get:

`X_(n+1)=-(D+L)^(-1)*U*X_n +(D+L)^(-1)*b` ---------------(6)

Each iteration of the Successive Over Relaxation method, of Eqn(2), takes the form:

Eqn(4) & Eqn(6) can be said as special cases of the more general formula:

`X_(n+1)=-(D+omega*L)^(-1)*((1-omega)*L+U)*X_n +(D+omega*L)^(-1)*b` ---------------(7)

With `omega=0` we recover the Jacobi method. With `omega=1` we recover the Gauss-Seidel method.

If, `0<omega<1`, then the solution will be called under relaxed and if, `omega>1`, then the solution is over relaxed

Thus, Eqn(4), Eqn(6), Eqn(7) can be summarized as:

To solve a system of equation of type `A*X=B` using the above iterative methods:

We expand square matrix A as `A=D+L+U`, where D is a strictly diagonal matrix of A, L is a strictly lower triangular matrix of A, U is a strictly upper triangular matrix of A. The equation can be solved by using the following formula:

`X_(n+1) = P*X_n + Q` ---------------(8)

where,

`P=-(D+omega*L)^(-1)*((1-omega)*L+U)` ---------------(9)

`Q=(D+omega*L)^(-1)` ---------------(10)

when,

`omega=0text(, we get Jacobi method)`

`0<omega<1text(, we get under relaxed SOR method)`

`omega=1text(, we get Gauss-Siedel method)`

`omega>1text(, we get over relaxed SOR method)`

The eigenvalue of a square matrix:

Eigenvalues are a special set of scalars associated with a linear system of equations (i.e., a matrix equation) that are sometimes also known as characteristic roots, characteristic values, proper values, or latent roots.

For a square matrix, they are calculated as:

Let A be a linear transformation represented by a matrix A. If there is a vector `X in RR^n ne0` such that

`A*X=lambda*X` ---------------(11)

for some scalar `lambda`, then `lambda` is called the eigenvalue of A with corresponding (right) eigenvector X.

Eqn(8) can be written compactly as:

`(A-lambda*I)*X=0` ---------------(12)
where I is the identity matrix. As shown in Cramer's rule, a linear system of equations has nontrivial solutions if the determinant vanishes, so the solutions of Eqn(9) are given by

`det(A-lambda*I)=0` ---------------(13)

This equation is known as the characteristic equation of A, and the left-hand side is known as the characteristic polynomial.

Spectral Radius of a square matrix:

The spectral radius of a square matrix or a bounded linear operator is the largest absolute value of its eigenvalues (i.e. supremum among the absolute values of the elements in its spectrum). It is sometimes denoted by `rho(.)`.

In the program, it is calculated after calculating eigenvalues as:

`bb"SR = max(abs(eigenvalues));"`

A circle drawn of the radius of the spectral radius will include all eigenvalues in a coordinate system. It changes as the values of the diagonal element changes.

Program:

1. Explanation: First we defined matrix A and then divide it into 3 parts such that A = D+L+U. Then for comparison, X is calculated numerically. Error and Tolerance values are initialized, no. of magnification points needed is defined and then mag array is calculated.
2. A combination of 'for' & 'if' loop is used to check if the mag array contains '1' and if it does contain '1' then it will change the condition to 1, a second 'if' loop used to insert '1' in mag array if the condition is still 0. Then the mag array is updated
3. All 3 functions are called one-by-one, function1 calculates the solution using the Jacobi iterative method, function 2 calculates the solution by using the Gauss-Siedel iterative method, function 3 calculates the solution using the SOR method.
4. Figure(1), figure(2), figure(3) is used to plot a subplot of "magnification factor vs spectral radius" & " Spectral radius vs no. of iterations" for the Jacobi, Gauss-Siedel, SOR method respectively.
5. figure(4) shows the plot of "diagonal magnification factor vs convergence" for each method, this helps in visualizing from which mag value the solution starts diverging.
6. figure(5) shows the plot of "diagonal magnification factor vs spectral radius" for each method, this helps in visualizing the change of spectral radius of the iteration matrix due to mag value in each case.
7. figure(6) shows the plot of "spectral radius vs convergence" for each method, this helps in visualizing from which spectral radius value the solution starts diverging.

In each function, divergence is checked by the condition that if the difference between error equals or greater than 60, the solution is considered to diverge and for those mag values, convergence array is updated with 0 which signifies that the solution is no more converging.

For each function, the solution is calculated as a 3Xne matrix, where 'ne' is the no. of mag elements, thus 3 solutions for every mag value is calculated. Eqn(8), Eqn(9), Eqn(10) is used to calculate X using suitable value of `omega`.

Main program:

clear all
close all
clc
A =[5 1 2;-3 9 4;1 2 -7] %coefficiet matrix
D=[A(1,1) 0 0;0 A(2,2) 0;0 0 A(3,3)]; %diagonal Matrix of A
L=[0 0 0;A(2,1) 0 0;A(3,1) A(3,2) 0]; %Lower Triangular Matrix of A
U=[0 A(1,2) A(1,3);0 0 A(2,3);0 0 0]; %Upper Triangular Matrix of A
% A1=D+L+U
% X=[X1 X2 X3]
B=[10;-14;-33]; %RHS
num_x=inv(A)*B  %Numerical solution
error=9e9; %initial error value
tol=1e-8; %initial tolerance value
ne=100; %no. of elements
mag=linspace(0.5,1.3,ne)% magnification value, 100 terms
condition=0;
for o=1:ne %loop to check if mag array include 1 or not
    if (mag(o)==1)
        condition=1;
        break
    end    
end
if (condition==0) %if mag array dose not have 1 then append 1 in it
    [difference,imag]=min(abs(mag-1))
    mag(imag-1:imag+1);
    new_mag=[mag(1:imag) 1 mag(imag+1:end)];
    mag=new_mag; %updating mag array
end
sorfactor=1.2; % relaxation factor
[j_x,j_counter,j_eig,j_sr,j_convergence]=jacobi(D,L,U,B,mag,error,tol); %function for jacobi
[gs_x,gs_counter,gs_eig,gs_sr,gs_convergence]=gs(D,L,U,B,mag,error,tol); %function for gauss siedal
[sorgs_x,sorgs_counter,sorgs_eig,sorgs_sr,sorgs_convergence]=sorgs(D,L,U,B,mag,sorfactor,error,tol); %function for SOR GS

figure(1) 
subplot(2,1,1)
plot(mag,j_sr,'-*','color','r') % comparison of magnification factor and spectral radius for jacobi method
xlabel('Diagonal magnification factor')
ylabel('Spectral radius')
xlim([0.4,1.6])
ylim([0,1.5])
subplot(2,1,2)
plot(j_sr,j_counter,'-^','color','g') % comparison of spectral radius & no. of iteration for jacobi method
xlabel('Spectral radius')
ylabel('No.of iterations')
xlim([0.3,1.1])
% ylim([0,270])
sgtitle('Jacobi iterative method')

figure(2)  
subplot(2,1,1)
plot(mag,gs_sr,'-*','color',[0.8500 0.3250 0.0980])% comparison of magnification factor and spectral radius for Gauss-Siedal method
xlabel('Diagonal magnification factor')
ylabel('Spectral radius')
xlim([0.4,1.6])
ylim([0,1.5])
subplot(2,1,2)
plot(gs_sr,gs_counter,'-^','color',[0.4940 0.1840 0.5560])% comparison of spectral radius & no. of iteration for Gauss-Siedal method
xlabel('Spectral radius')
ylabel('No.of iterations')
xlim([0.1,1.5])
% ylim([0,270])
sgtitle('Gauss-Siedal iterative method')

figure(3)  
subplot(2,1,1)
plot(mag,sorgs_sr,'-*','color','k')% comparison of magnification factor and spectral radius for SOR on GS method
xlabel('Diagonal magnification factor')
ylabel('Spectral radius')
xlim([0.4,1.6])
ylim([0,2.2])
subplot(2,1,2)
plot(sorgs_sr,sorgs_counter,'-^','color','m')% comparison of spectral radius & no. of iteration for SOR on GS method
xlabel('Spectral radius')
ylabel('No.of iterations')
xlim([0.2,2.1])
% ylim([0,270])
sgtitle('SOR on Gauss-Siedal iterative method')

figure(4)
convergence=[j_convergence;gs_convergence;sorgs_convergence]; %convergence array 3x3
plot(mag,convergence(1,:),'-*',mag,2.*convergence(2,:),'-^',mag,3.*convergence(3,:),'-x','LineWidth',1.2) %convergance for magnification factor
ylim([-0.5,3.5])
xlim([0.4,1.5])
xlabel('Diagonal magnification value')
ylabel({['Convergence'],['0 for no convergence if difference between error>=60'],['1,2,3 for convergence when difference between error<60']})
title('Convergence of iterative solver wrt Diagonal magnification value')
legend('1*Jacobi Method','2*Gauss-Siedal Method','3*SOR on GS Method')

figure(5)
sr=[j_sr;gs_sr;sorgs_sr]; %spectral radius array 3x3
plot(mag,sr(1,:),'-*',mag,sr(2,:),'-^',mag,sr(3,:),'-x','LineWidth',1.2) %comparison of magnification factor and spectral radius
xlim([0.4,1.6])
xlabel('Diagonal magnification value')
ylabel('Spectral Radius')
title('Comparison between Spectral Radius & Diagonal magnification value')
legend('SR Jacobi Method','SR Gauss-Siedal Method','SR SOR on GS Method')

figure(6)
plot(sr(1,:),convergence(1,:),'-*',sr(2,:),2.*convergence(2,:),'-^',sr(3,:),3.*convergence(3,:),'-x','LineWidth',1.2)
xlim([0.4,1.6])
ylim([-0.5,3.5])
xlabel('Spectral Radius')
ylabel({['Convergence'],['Solution converged if y>0'],['Solution diverging if y=0']})
title('Comparison between Spectral Radius & Diagonal magnification value')
legend('SR Jacobi Method','SR Gauss-Siedal Method','SR SOR on GS Method')

Jacobi Method Function:

function[j_x,j_counter,j_eig,j_sr,j_convergence]=jacobi(D,L,U,B,mag,error,tol)
    w=0;    %for jacobi
    j_convergence=ones(1,length(mag)); %convergence true and false matrix intialisation with 1
    counter=zeros(1,length(mag)); %intialising counter with 0, 1xlength(mag) matrix to store no. of itertion for every mag value 
    eigen=zeros(3,length(mag));%intialising eigen with 0, 3xlength(mag) matrix to store 3 eigen values for every mag value 
    sr=zeros(1,length(mag));%intialising counter with 0, 1xlength(mag) matrix to spectral radius for every mag value 
    x_new=zeros(3,length(mag));%initialising solution matrix of 3x3 with 0
    x=x_new;
    for i=1:length(mag) %loop to change magnification value
        e=error; %error re-initialising
        e_old=0;
        d=mag(i).*D;   %new diagonal matrix     
        p = -(inv(d+(w.*L)))*(((1-w).*L)+U);   % jacobi iteration matrix
        q = (inv(d+(w.*L)))*B; %coefficient matrix
        eigen(:,i)=abs(eig(p)); %eigen value
        sr(1,i)=max(eigen(:,i));%spectral radius
        while e>tol %convergence checking
            x_new(:,i) = (p*x(:,i))+q; %new solution
            e=max(abs(x_new(:,i)-x(:,i))); %error calculation
            x=x_new;% update of old solution
            counter(1,i)=counter(1,i)+1; %iteration counter             
            ed=e-e_old %calculating difference in error of current iteration and previous iteration
            e_old=e;        %updating old error values  
            if (ed>=60) %divergence checking
                j_convergence(1,i)=0; % 0 for divergence
                break
            end
        end
    end  
    j_eig=eigen;
    j_x=x;
    j_sr=sr;
    j_counter=counter;
end

Gauss-Siedal Method Function:

function[gs_x,gs_counter,gs_eig,gs_sr,gs_convergence]=gs(D,L,U,B,mag,error,tol)
    w=1;    %for gauss siedal
    gs_convergence=ones(1,length(mag));%convergence true and false matrix intialisation with 1
    counter=zeros(1,length(mag));%intialising counter with 0, 1xlength(mag) matrix to store no. of itertion for every mag value 
    eigen=zeros(3,length(mag));%intialising eigen with 0, 3xlength(mag) matrix to store 3 eigen values for every mag value 
    sr=zeros(1,length(mag));%intialising counter with 0, 1xlength(mag) matrix to spectral radius for every mag value
    x_new=zeros(3,length(mag)); %initialising solution matrix of 3x3 with 0
    x=x_new;
    for i=1:length(mag)%loop to change magnification value
        e=error; %error re-initialising
        e_old=0;
        d=mag(i).*D;   %new diagonal matrix     
        p = -(inv(d+(w.*L)))*(((1-w).*L)+U);   % jacobi iteration matrix
        q = (inv(d+(w.*L)))*B; %coefficient matrix
        eigen(:,i)=abs(eig(p)); %eigen value
        sr(1,i)=max(eigen(:,i));%spectral radius
        while e>tol %convergence checking
            x_new(:,i) = (p*x(:,i))+q; %new solution
            e=max(abs(x_new(:,i)-x(:,i))); %error calculation
            x=x_new;% update of old solution
            counter(1,i)=counter(1,i)+1; %iteration counter             
            ed=e-e_old %calculating difference in error of current iteration and previous iteration
            e_old=e;      %updating old error values    
            if (ed>=60) %divergence checking
                gs_convergence(1,i)=0; % 0 for divergence
                break
            end
        end
    end   
    gs_eig=eigen;
    gs_x=x;
    gs_sr=sr;
    gs_counter=counter;
end

SOR on GS Method:

function[sorgs_x,sorgs_counter,sorgs_eig,sorgs_sr,sorgs_convergence]=sorgs(D,L,U,B,mag,w,error,tol)
    sorgs_convergence=ones(1,length(mag));%convergence true and false matrix intialisation with 1
    counter=zeros(1,length(mag));%intialising counter with 0, 1xlength(mag) matrix to store no. of itertion for every mag value 
    eigen=zeros(3,length(mag));%intialising eigen with 0, 3xlength(mag) matrix to store 3 eigen values for every mag value 
    sr=zeros(1,length(mag));%intialising counter with 0, 1xlength(mag) matrix to spectral radius for every mag value
    x_new=zeros(3,length(mag));%initialising solution matrix of 3x3 with 0
    x=x_new;
    for i=1:length(mag)%loop to change magnification value
        e=error; %error re-initialising
        e_old=0;
        d=mag(i).*D;   %new diagonal matrix     
        p = -(inv(d+(w.*L)))*(((1-w).*L)+U);   % jacobi iteration matrix
        q = (inv(d+(w.*L)))*B; %coefficient matrix
        eigen(:,i)=abs(eig(p)); %eigen value
        sr(1,i)=max(eigen(:,i));%spectral radius
        while e>tol %convergence checking
            x_new(:,i) = (p*x(:,i))+q; %new solution
            e=max(abs(x_new(:,i)-x(:,i))); %error calculation
            x=x_new;% update of old solution
            counter(1,i)=counter(1,i)+1; %iteration counter             
            ed=e-e_old; %calculating difference in error of current iteration and previous iteration
            e_old=e;     %updating old error values   
            if (ed>=60) %divergence checking
                sorgs_convergence(1,i)=0; % 0 for divergence
                break
            end
        end
    end
    sorgs_eig=eigen;
    sorgs_x=x;
    sorgs_sr=sr;
    sorgs_counter=counter;
end

Result:

Figure(1):

Figure(2):

Figure(3):

Figure(4):

Figure(5):

Figure(6):

Observations:

1. We can observe from the above graphs that the value of the spectral radius of the iteration matrix decreases as diagonal magnification value increases for every case.
2. No. of iteration to converge increases when spectral radius proceeds from 0 to 1 and after crossing 1 the solution starts diverging as the difference between error starts increasing rather than decreasing, it is also observed that the solution diverges more quickly as the spectral radius starts increasing beyond 1.
3. In figure(4) it is observed that the solution starts diverging when the mag value reaches a point from which the spectral radius of the iteration matrix reaches 1 and starts increasing.
4. Figure(5) proves point 1 as the trend of the graph tells that the spectral radius of the iteration matrix decreases as the diagonal magnification factor increases.
5. Figure(6) proves point 2 as it clearly shows that the solution starts diverging when the spectral radius reaches 1 and starts increasing, for all value of spectral radius less than 1, the solution converges.