Gauss-Seidel Method MATLAB Program

Gauss-Seidel method is a popular iterative method of solving linear system of algebraic equations. It is applicable to any converging matrix with non-zero elements on diagonal. The method is named after  two German mathematicians: Carl Friedrich Gauss and Philipp Ludwig von Seidel.

Gauss-Seidel is considered an improvement over Gauss Jacobi Method. In this method, just like any other iterative method, an approximate solution of the given equations is assumed, and iteration is done until the desired degree of accuracy is obtained.

In earlier tutorials, we’ve already gone through the C program and algorithm/flowchart for Gauss-Seidel method. Here, we’re going to write a program code for Gauss-Seidel method in MATLAB, discuss its theoretical background, and analyze the MATLAB program’s result with a numerical example.

Gauss-Seidel Method Theory:

Let’s go through a brief theoretical/mathematical background of Gauss-Seidel method. The matrices, iterations, and the procedure explained below cover the basic guidelines to write the program code for Gauss-Seidel method in MATLAB.

Consider the following system of linear equations:

a₁₁x₁ + a₁₂x₂ + a₁₃x₃ + a₁₄x₄ + a₁₅x₅ + a₁₆x_{6 …….} + a_1nx_n = b₁
a₂₁x₁ + a₂₂x₂ + a₂₃x₃ + a₂₄x₄ + a₂₅x₅ + a₂₆x_{6 …….} + a_2nx_n = b₂
a₃₁x₁ + a₃₂x₂ + a₃₃x₃ + a₃₄x₄ + a₃₅x₅ + a₃₆x_{6 …….} + a_3nx_n = b_{3
. . . . . . . . . . . . . .  . . . . . . . . .   . . . . . . .  . ….. . . . . . . . ..  . . . . . . . . . . . . ..  . . .. .  . .
. . . . . . . . . . . . . .  . . . . . . . . .   . . . . . . .  . ….. . . . . . . . ..  . . . . . . . . . . . . ..  . . .. .  . .}
a_n1x₁ + a_n2x₂ + a_n3x₃ + a_n4x₄ + a_n5x₅ + a_n6x_{6 …….} + a_nnx_n = b_n

where, a_ij represents the coefficient of unknown terms x_i.

The above equations can be presented in matrix form as follows:

Or simply, it can be written as: [A][X] = [B]

Now, decomposing the matrix A into its lower triangular component and upper triangular component, we get:

A = L x U

where,

Further, the system of linear equations can be expressed as:

L x X = B – UX     —–(a)

In Gauss-Seidel method, the equation (a) is solved iteratively by solving the left hand value of x and then using previously found x on right hand side. Mathematically, the iteration process in Gauss-Seidel method can be expressed as:

X^(k+1) = L ^-1( B –UX^(k) )

Applying forward substitution, the elements of X^(k+1) can be computed as follows:

The same procedure aforementioned is followed in the MATLAB program for this method. The process of iteration is continued till the values of unknowns are under the limit of desired tolerance.

Gauss-Seidel Method in MATLAB:

% Gauss-Seidel Method in MATLAB
function x = gauss_siedel( A ,B )
disp ( 'Enter the system of linear equations in the form of AX=B')

%Inputting matrix A
A = input ( 'Enter matrix A :   \n')
% check if the entered matrix is a square matrix
[na , ma ] = size (A);
if na ~= ma
    disp('ERROR: Matrix A must be a square matrix')
    return
end

% Inputting matrix B
B = input ( 'Enter matrix B :   ')
% check if B is a column matrix
[nb , mb ] = size (B);
if nb ~= na || mb~=1
   disp( 'ERROR: Matrix B must be a column matrix')
   return
end

% Separation of matrix A into lower triangular and upper triangular matrices
% A = D + L + U
D = diag(diag(A));
L = tril(A)- D;
U = triu(A)- D

% check for convergence condition for Gauss-Seidel method
e= max(eig(-inv(D+L)*(U)));
if abs(e) >= 1
    disp ('Since the modulus of the largest Eigen value of iterative matrix is not less than 1') 
    disp ('this process is not convergent.')
    return
end

% initial guess for X..?
% default guess is [ 1 1 .... 1]
r = input ( 'Any initial guess for X? (y/n):   ','s');
switch r 
    case 'y'
        % asking for initial guess
    X0 = input('Enter initial guess for X :\n')
        % check for initial guess
    [nx, mx] = size(X0);
        if nx ~= na || mx ~= 1
        disp( 'ERROR: Check input')
        return
    end
    otherwise
    X0 = ones(na,1);
end

% allowable error in final answer
t = input ( 'Enter the error allowed in final answer:  ');
tol = t*ones(na,1);
k= 1;

X( : , 1 ) = X0;
err= 1000000000*rand(na,1);% initial error assumption for looping
while sum(abs(err) >= tol) ~= zeros(na,1)
    X ( : ,k+ 1 ) = -inv(D+L)*(U)*X( : ,k) + inv(D+L)*B;% Gauss-Seidel formula
    err = X( :,k+1) - X( :, k);% finding error
    k = k + 1;
    
end

fprintf ('The final answer obtained after %g iterations is  \n', k)
X( : ,k)

Also see,
Gauss-Seidel C Program
Gauss-Seidel Algorithm/Flowchart

In the above MATLAB program, a function, x = gauss_siedel( A ,B ), is initially defined. Here, A and B are the matrices generated with the coefficients used in the linear system of equations. The elements of A and B are input into the program following the basic syntax of MATLAB programming.

A and B are to be checked: A should be a square matrix and B must be a column matrix to satisfy the criteria of Gauss-Seidel method. Then, as explained in the theory, matrix A is split into its upper triangular and lower triangular parts to get the value of first iteration.

The value of variables obtained from the first iteration are used to start the second iteration, and the program keeps on iterating till the solution are in the desired limit of tolerance as provided by the user.

Here’s a sample output screen of the MATLAB program:

Gauss-Seidel Method Example:

The above MATLAB program of Gauss-Seidel method in MATLAB is now solved here mathematically. The equations given are:

4x₁ – x₂ –x₃ = 3
-2x₁ + 6x₂ + x₃ = 9
-x₁ + x₂ – 7x₃ = -6

In order to get the value of first iteration, express the given equations as follows:

4x₁ – 0 –0 = 3
-2x₁ + 6x₂ + 0 = 9
-x₁ + x₂ – 7x₃ = -6

From the first equation: x₁ = 3/4 = 0.750
Substitute the value of x₁ in the second equation : x₂ = [9 + 2(0.750)] / 6 = 1.750
Substitute the values of x₁ and x₂ in the third equation: x₃ = [-6 + 0.750 − 1.750] / 7 = − 1.000
Thus, the result of first iteration is: ( 0.750, 1.750, -1.000 )

The whole iteration procedure that goes on in Gauss-Seidel method (and the above MATLAB program) is presented below:

where, k is the number of iteration.

The final solution obtained is (1.000, 2.000, -1.000).

If you have any questions regarding Gauss-Seidel method, its theory, or MATLAB program, drop them in the comments. You can find more Numerical methods tutorial using MATLAB here.