Matlab Code & Implementation of Jacobi Method or Jacobi Iterative Method For Solving the System Of Linear Equations

In linear algebra, the Jacobi method (or Jacobi iterative method) is an algorithm for determining the solutions of a diagonally dominant system of linear equations.
The process is iterated until it converges to the solution. This algorithm is a stripped-down version of the Jacobi transformation method of matrix diagonalization. The method is named after Carl Gustav Jacob Jacobi. (Source: Wikipedia)

Suppose the system of Linear equation is denoted as 'Ax=b'
In order to solve it we need to find the solution of 'x'
To apply Jacobi Iteration we first need to check for the diagonal dominance of the measurement matrix 'A'. If the matrix 'A' is not diagonally dominant then the solution will not converge.
Example of Diagonal Dominant Matrix is:
A= 5 3 2 
      1 6 2
      2 1 7
Although if the given measurement matrix is not diagonally dominant we need to perform elementary row transform to make it diagonally domimant.
An example of matrix which is not diagonal dominant is:
A= 2 7 3
      4 1 -1
      1 -3 12
In order to make it diagonally dominant we need to exchange 1st & 2nd row, thus the matrix will be:
A= 4 1 -1
      2 7  3
      1 -3 12
Its simple to check for the diagonal dominance of a matrix, please refer this link. Check for the diagonaly Dominance Of a Matrix.
NOTE: We need to apply the same elementary transformation in the vector 'b'

NOTE: We have implemented the Jacobi Iteration method for  a 3x3 measurement matrix, the code will be applicable for higher order system also, but the check for the diagonal dormancy will become tedious.

In normal method we solve for 'x' as x= inv(A)*y (inv(A) is the inverse of matrix A)
But in Jacobi Iterative Method for 'Ax=y' what we do is as Follow:
We define matrix 'A' as the sum of two matrix 'D' & 'E'. Here 'D' is the diagonal matrix & 'E' is the remainder matrix of same order after subtracting 'D' from 'A'.
Then Ax=b will become (D+E)x=y 
=> Dx= -Ex+y
=> x= -inv(D)*E*x+inv(D)*y
*we start with a random value of x (say in MATLAB term it will be rand(3,1) for a 3x3 matrix A)
Then we apply it iteratively for K+1th step the formula, x(K+1)= -inv(D)*E*x(K)+inv(D)*y
The above method is known as Jacobi iteration for solving the system of Linear Equation.

MATLAB Code for Jacobi Method or Jacobi Iterative Method, ( Code in BOLD Text)

function jacobi(A, b, N) 
%function Jacobi with 3 Params A, b from Ax=b & N is the number of iterations
test=all((2*abs(diag(A)))- sum(abs(A),2)>=0);  %1st check of Diagonal Dominance
if test==0
    A([1 2],:) = A([2 1],:);
  
%Elementary Row Transformation of A
    b([1 2]) = b([2 1]);     %Elementary Row Transformation of b
end

test=all((2*abs(diag(A)))- sum(abs(A),2)>=0);
 
%2nd check of Diagonal Dominance after 1st exchange

if test==0
    A([2 1],:) = A([1 2],:); %if test=0 then switch back to original A
    b([2 1]) = b([1 2]);
    A([1 3],:) = A([3 1],:); %Elementary Row Transformation of A   
    b([1 3]) = b([3 1]); %Elementary Row Transformation of b
end
d=diag(A);
%Extracting the diagonal Element of A
D=diag(d); %Making the diagonal matrix D from extracted diagonal elements
D_inv=1/D; %inverse of a diagonal matrix is simply 1/D
E=A-D; %Calculating Remainder Matrix
x=rand(3,1); %Initializing the vector 'x' with a random initial value
T=-D_inv*E;
C=D_inv*b;

for j=1:N
%performing Jacobi Iterations for 'N' times
    x=T*x+C;
end
disp(x)
%Display of result

NOTE: The above code is developed for a 3x3 matrix A & 3x1 vector b, but applicable to any dimension matrix, but in order this method to converge you need to provide A such that it is diagonally dominant.

NOTE: We can also save the content of x to see what happens in each step, just by making the statement x=T*x+C; as x(j+1)=T*x(j)+C;

Test Cases:
Case 1: (When Matrix is diagonally Dominant)
>> A=[5 3 2;1 6 2;2 1 7]

A =

     5     3     2
     1     6     2
     2     1     7

>> b=[14;13;24]

b =

    14
    13
    24

>> jocobi(A,b,40)

x =

    1.0000
    1.0000
    3.0000



Case 2: (When Matrix is not diagonally Dominant, thus not converging)
 A=[5 1 5;6 2 3;7 8 1]

A =

     5     1     5
     6     2     3
     7     8     1

>> b=[11;11;18]

b =

    11
    11
    18

>> jocobi(A,b,40)

%this result is obtained only when the elementary row transformation part was not done to make The matrix A diagonally dominant.

x =

   1.0e+22 *

   -0.3795
   -1.1291
   -0.4246


Although by doing simple mldivide we will get the result
x=b\A

x =   0.4364    0.3127    0.1873 






0 comments: