Question

Use Matlab or any other Language /Tool to implement the project work Algorithm for Gauss-Seidel Method to solve the linear (n x n) system Ax = b in matrix form is given by x(0) = initial vector x(4+1) = D-1 (b – Lx(k+1) – Ux(k)), k = 0, 1, 2, ...... where A = L+D+U. Here L, D, U are respectively lower, diag- onal and upper matrices constructed from A.

Use Matlab or any other Language /Tool to implement the project work Algorithm for Gauss-Seidel Method to solve the linear (n x n) system Ax = b in matrix form is given by x(0) = initial vector x(4+1) = D-1 (b – Lx(k+1) – Ux(k)), k = 0, 1, 2, ...... where A = L+D+U. Here L, D, U are respectively lower, diag- onal and upper matrices constructed from A.

I am new to Matlab so details would be appreciated!

Use Matlab or any other Language /Tool to implement the project work Algorithm for Gauss-Seidel Method to solve the linear (n x n) system Ax = b in matrix form is given by x(0) = initial vector x(4+1) = D-1 (b – Lx(k+1) – Ux(k)), k = 0, 1, 2, ...... where A = L+D+U. Here L, D, U are respectively lower, diag- onal and upper matrices constructed from A.

Answer 1

Please find posted the Matlab scripts gaussSeidel.m and gaussSeidelSolver.m. Run the script gaussSeidelSolver.m to test a case.

gaussSeidel.m

function x = gaussSeidel(A, b, tol)
% Gauss Siedel method
% x = gaussSiedel(A, b, tol)

% input:
% A = coefficient matrix
% b = right hand side vector
% tol = tolerance, stopping criterion

% output:
% x = solution vector

[nRows, nCols] = size(A);
% checks for Square Matrix
if nRows ~= nCols
    error('Matrix must be square');
end

% checks for Diagonal Dominance
for m = 1: nRows
    mthRow = abs(A(m, :));
    d = sum(mthRow) - mthRow(m);
    if mthRow(m) < d
        error('Matrix is not diagonally dominant');
    end
end

% split Matrix into three matrices L, U, and D

D = diag(diag(A));
L = tril(A, 1);
U = triu(A, 1);

% sets Initial Guess
x0 = zeros(nRows, 1);
x = x0;

% starts the Iterative Method

iter = 0;
err = inf;
errVector = inf;

while err > tol
  
    disp([num2str(iter) ' iteration,' ' x = ' mat2str(x) ', error = ' mat2str(errVector) '.']);
  
    % update [x]
    x_old = x;
    x = inv(L + D) * (-U*x + b);
    dx = x - x_old;
    % error vector
    errVector = abs(dx./x);
  
    % estimate error
    err = max(abs(dx./x));
    iter = iter + 1;
end

% displays Result
disp(['Gauss Seidel method converges after ' num2str(iter) ' iterations.' ' x = ' mat2str(x) ' error = ' mat2str(errVector) '.']);
end

gaussSeidelSolver.m

clear; clc;

A = [ -7 1 5;
      0 -7 4;
      1   2 -4];

b = [ -3 7 10 ]';

tol = 1e-3;

x = gaussSeidel(A, b, tol);