Page 111 - Applied Numerical Methods Using MATLAB

100 SYSTEM OF LINEAR EQUATIONS
a m1 x 1 + a m2 x 2 +· · · + a mm x m + ··· + a mN x N = b m
N
a mn (k) b m
(k+1)
x m =− x n + for m = 1, 2,...,N
a mm a mm
n =m
x k+1 = A x k + b for each time stage k (2.5.3)

where
   
0 −a 12 /a 11 ·· · −a 1N /a 11 b 1 /a 11
−a 21 /a 22
 0 ·· ·   
, b =
A N×N =  −a 2N /a 22   b 2 /a 22 

· · · · · · ·
   
−a N1 /a NN −a N2 /a NN ·· · 0 b N /a NN
This scheme is implemented by the following MATLAB routine “jacobi()”.
We run it to solve the above equation.
function X = jacobi(A,B,X0,kmax)
%This function finds a soltuion to Ax=Bby Jacobi iteration.
if nargin < 4, tol = 1e-6; kmax = 100; %called by jacobi(A,B,X0)
elseif kmax < 1, tol = max(kmax,1e-16); kmax = 100; %jacobi(A,B,X0,tol)
else tol = 1e-6; %jacobi(A,B,X0,kmax)
end
if nargin < 3, X0 = zeros(size(B)); end
NA = size(A,1);
X = X0; At = zeros(NA,NA);
for m = 1:NA
for n = 1:NA
if n ~= m, At(m,n) = -A(m,n)/A(m,m); end
end
Bt(m,:) = B(m,:)/A(m,m);
end
fork=1: kmax
X = At*X + Bt; %Eq. (2.5.3)
if nargout == 0, X, end %To see the intermediate results
if norm(X - X0)/(norm(X0) + eps) < tol, break; end
X0=X;
end

>>A=[32;1 2];b=[1 -1]’; %the coefficient matrix and RHS vector
>>x0 = [0 0]’; %the initial value
>>x = jacobi(A,b,x0,20) %to repeat 20 iterations starting from x0
x = 1.0000
-1.0000
>>jacobi(A,b,x0,20) %omit output argument to see intermediate results
X = 0.3333 0.6667 0.7778 0.8889 0.9259 ......
-0.5000 -0.6667 -0.8333 -0.8889 -0.9444 ......

2.5.2 Gauss–Seidel Iteration
Let us take a close look at Eq. (2.5.1). Each iteration of Jacobi method updates
the whole set of N variables at a time. However, so long as we do not use a