Question

How can one write a Matlab code for using Jacobi and Gauss-seidel methods to solve the linear systems in exercise 7.3 question 3(a) and 3(d)? (Numerical Analysis 9th Edition by Burden and Faires)

Answer 1

Since you have not mentioned the linear system of equations, I am providing code for solving any equation of the type Ax = b.

Solution of x in Ax=b using Gauss Seidel Method

A=[5 -2 3 0 6; -3 9 1 -2 7.4; 2 -1 -7 1 6.7; 4 3 -5 7 9; 2 3.5 6.1 -4 -8.1]

b=[-1 2 3 0.5 3.1]'

x=rand(5,1)

n=size(x,1);

normVal=Inf;

tol=1e-3; GaussItr=0;

Algorithm: Gauss Seidel Method

plotGauss=[];

while normVal>tol

x_old=x;

for i=1:n

sigma=0;

for j=1:i-1

sigma=sigma+A(i,j)*x(j);

end

for j=i+1:n

sigma=sigma+A(i,j)*x_old(j);

end

x(i)=(1/A(i,i))*(b(i)-sigma);

end

GaussItr=GaussItr+1;

normVal=norm(x_old-x);

plotGauss=[plotGauss;normVal];

end

fprintf('Solution of the system is : \n%f\n%f\n%f\n%f\n%f in %d iterations',x,GaussItr);

OUTPUT

Solution of the system is : 0.551424 0.469217 -0.595135 -0.649056 -0.171480 in 93 iterations

Solution of x in Ax=b using Jacobi Method

x=rand(5,1)

n=size(x,1);

normVal=Inf;

JacobItr=0;

Algorithm: Jacobi Method

plotJacobi=[];

while normVal>tol

xold=x;

for i=1:n

sigma=0;

for j=1:n

if j~=i

sigma=sigma+A(i,j)*x(j);

end

end

x(i)=(1/A(i,i))*(b(i)-sigma);

end

JacobItr=JacobItr+1;

normVal=norm(xold-x);

plotJacobi=[plotJacobi;normVal];

end

fprintf('Solution of the system is : \n%f\n%f\n%f\n%f\n%f in %d iterations',x,JacobItr);

OUTPUT

Solution of the system is :

0.551539

0.469360

-0.595213

-0.649104

-0.171425 in 86 iterations