MATLAB TUTORIAL for the Second Course. Part 2.5: Iterative Methods

An iterative method to solve the linear system A x = b starts with an initial approximation p₀ to the solution x and generates a sequence of vectors \( \left\{ {\bf p}_k \right\}_{k\ge 0} \) that converges to x. Iterative methods involve a process that converts the system A x = b into an equivalent system of the form x = B x + w, for some fixed matrix B and vector b. After the initial vector p₀ is selected, the sequence of approximate solutions is generated by computing \( {\bf x}^{(k+1)} = {\bf B}\,{\bf x}^{(k)} + {\bf b} , \) for each \( k =0,1,2,\ldots . \)

For large systems containing thousands of equations, iterative methods often have decisive advantages over direct methods in terms of speed and demands on computer memory. Sometimes, if the accuracy requirements are not stringent, a modest number of iterations will suffice to produce an acceptable solution. Also, iterative methods are often very efficient for sparse systems problems. In sparse problems, the nonzero elements of A are sometimes stored in a sparse-storage format. In other case, it is not necessary to store A at all; for example, in problems involving the numerical solution of partial differential equations as, in this case, each row of A might be generated as needed but not retained after use. Another important advantage of iterative methods is that they are usually stable, and they will actually dampen errors, due to roundoff or minor blunders, as the process continues.

function x1 = jacobi1(a,b,x0,tol)
n = length(b);
x = zeros(n,1);
for j = 1 : n
     x(j) = ((b(j) - a(j,[1:j-1,j+1:n]) * x0([1:j-1,j+1:n])) / a(j,j)); % the first iteration
end
x1 = x';
k = 1;
x_ny = zeros(n,1);
while norm(x1-x0,1) > tol
    for j = 1 : n
     x_ny(j) = ((b(j) - a(j,[1:j-1,j+1:n]) * x1([1:j-1,j+1:n])) / a(j,j));
    end
    x0 = x1;
    x1 = x_ny';
    k = k + 1;
end
k
x = x1';

function = jacobi(A,b,p,delta,max1)
% Input:    A is an n x n nonsingular matrix
%           b is an n x 1 matrix (= n-vector)
%           p is an n x 1 matrix: the initial guess of column-vector
%           delta is the tolerance for p
%           max1 is the maximum number of iterations
% Output:   x is an n x 1 matrix (= column vector); the Jacobi approximation
n = length(b);
for k=1:max1
    for j=1:n
        x(j)=(b(j)-A(j,[1:j-1,j+1:n])*p([1:j-1,j+1:n]))/A(j,j);
    end
    err=abs(norm(x'-p));
    relerr=err/(norm(x) +eps);
    p=x';
         if(err < delta)|(relerr < delta)
         break
    end
end
x=x';

Example. The system of equations

We now generalize the Gauss-Seidel iteration procedure. Suppose we are given the linear system

We present matlab code to solve the linear system A x = b by starting with the initial guess \( {\bf x} = {\bf p}_0 \) and generating a sequence \( \left\{ {\bf p}_k \right\}_{k\ge 0} \) that converges to the solution. A sufficient condition for the method to be applicble is that A is either strictly diagonal dominant or symmetric and positive definite.