A MATLAB Program to Implement the Jacobi Iteration

A MATLAB Program to Implement Jacobi Iteration to Solve System of Linear Equations:

The following MATLAB codes uses Jacobi iteration formula to solve any system of linear equations where the coefficient matrix is diagonally dominant to achieve desired convergence.

%-----------------------------------------------------------------%
function [x, rel_error] = jacobimethod(A, b, x0, tol, iteration)

% Inputs: A - Coefficient matrix
%         b - Input matrix
%       tol - Defining tolerance for solution
% iteration - Number of iterations

% Outputs: x - Solutions
%  rel_error - Relative error

D = diag(diag(A));   % Making coefficient matrix diagonal
R = A - D;           % Construction of another matrix "R"
N = 1;               % iteration counter
x = x0;
rel_error = tol * 2; % norm(x - x0)/norm(x);

% Implementation of Jacobi method to solve Ax = b
while (rel_error>tol && N <= iteration)
   
    xprev = x;   
    x = inv(D)*(b - R*xprev);
    rel_error = norm(x - xprev)/norm(x);

    Z1=x(1);
    disp(Z1);
    Z2=x(10);
    disp(Z2);
    Z3=x(16);
    disp(Z3);
   
    fprintf('\n Iteration %i: Relative error =%d ',N, rel_error)   
    N = N + 1;       
end
%----------------------------------------------------------------%

Now, define your input matrices in the MATLAB command window and use the above developed function to compute Jacobi iteration.

MATLAB Command Prompt:

>> A = [-2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    1 -2 1 0 0 0 0 0 0 0 0 0 0 0 0 0;
    0 1 -2 1 0 0 0 0 0 0 0 0 0 0 0 0;
    0 0 1 -2 1 0 0 0 0 0 0 0 0 0 0 0;
    0 0 0 1 -2 1 0 0 0 0 0 0 0 0 0 0;
    0 0 0 0 1 -2 1 0 0 0 0 0 0 0 0 0;
    0 0 0 0 0 1 -2 1 0 0 0 0 0 0 0 0;
    0 0 0 0 0 0 1 -2 1 0 0 0 0 0 0 0;
    0 0 0 0 0 0 0 1 -2 1 0 0 0 0 0 0;
    0 0 0 0 0 0 0 0 1 -2 1 0 0 0 0 0;
    0 0 0 0 0 0 0 0 0 1 -2 1 0 0 0 0;
    0 0 0 0 0 0 0 0 0 0 1 -2 1 0 0 0;
    0 0 0 0 0 0 0 0 0 0 0 1 -2 1 0 0;
    0 0 0 0 0 0 0 0 0 0 0 0 1 -2 1 0;
    0 0 0 0 0 0 0 0 0 0 0 0 0 1 -2 1;
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -2];

>> b = [-1; -2; -2; -2; -2; -2; -2; -2; -2; -2; -2; -2; -2; -2; -2; 1];

>> x0 = [0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0];

>> tol=10^-5;

>> iteration=350;

>> jacobimethod(A, b, x0, tol, iteration)


ans =

   14.8499
   28.7009
   40.5541
   50.4104
   58.2708
   64.1360
   68.0066
   69.8830
   69.7653
   67.6537
   63.5477
   57.4473
   49.3516
   39.2600
   27.1715
   13.0852

1 comment:

  1. Don't use inv(), numerically really unstable. Moreover, D is just a diagonal, so inv(D) = ./D.

    ReplyDelete