Overview
The investigation of iterative solvers for Ax = b
continues with a look at the Gauss–Seidel method.
Each Gauss–Seidel ite ration requires O(n2) flops. (Jacobi’s method requires O(n)
flops per iteration;
a more detailed comparison of these two methods is the focus of Lab 11.)
New MATLAB commands introduced in this lab include
tril and triu, to extract the lower- and upper-triangular
parts of a matrix, and sparse and full, to work with matrices with lots of zeros
(sparse
matrices).
Part I
As with the Jacobi method, write A = L+D+U where L is the n×n matrix
containing all components
be low the main diagonal of A, D is the n × n matrix containing the diagonal
comp onents of A , and
U is the n × n matrix containing all components above the main diagonal of A.
The Gauss–Seidel method uses S = L +D and T = S −A. Thus, in each step forward
substitution is
needed to solve
A 5 × 5 Example
The following MATLAB commands compute the first ten (10) Gauss–Seidel
iterations with initial guess
x = 0.
| » S = tril( A ) |
% Gauss-Seidel method |
| » T = S - A |
% A = S − T |
| » x = zeros(size(b)) |
% initial guess is x = 0 |
| » for k=1:10, |
|
| » x = S \ (b+T*x) |
% perform one Gauss-Seidel iteration |
| » end |
|
It should be apparent that the vectors x are converging to
the exact solution of Ax = b. (Compare the
iterates and convergence with those obtained with Jacobi’s method for this same
problem in Lab 9.)
Block Construction of Matrices
As matrices become larger it becomes increasingly tedious to enter every
component of the matirx. In
some cases it is possible to decompose a matrix into simpler blocks . For
example,
where I = I5×5 and B is the 5 × 5 diagonal matrix with bkk
= k for k=1, 2, 3, 4, 5.
The following commands can be used to enter A in MATLAB:
| » I = eye(5,5) |
% I5×5 |
| » B = diag(1:5) |
|
| » A = [ 9*I B ; B 9*I ] |
% block definition of A |
Sparse Matrices
A sparse matrix is a matrix with very few non- zero components . Typically, the
number of non -zero
components in a sparse matrix is of order n . Sparse matrices are particularly
well-suited for iterative
methods. The savings typically arises from efficient implementations of matrix
multiplication that
avoids any multiplication by zero. MATLAB can store sparse matrices economically
and can perform
matrix operations involving sparse matrices efficiently.
The sparse commands converts an existing MATLAB matrix into the corresponding
sparse representation.
» Asp = sparse( A ) % create sparse version of A
Note that the matrix Asp is presented in terms of the non-zero components and
their indices. (See
help sparse to learn more about the sparse command and sparse matrices in MATLAB.)
One benefit of using sparse matrices is illustrated by matrix multiplication:
» b = ones(10,1)
» flops(0), Asp*b, flops % sparse matrix- vector product
» flops(0), A*b, flops % full matrix-vector product
» flops(0), Asp*Asp, flops % sparse matrix-matrix product
» flops(0), A*A, flops % full matrix-matrix product
The sparse format of a matrix is preserved by matrix operations such as
transpose, diag, tril, and
lu.
» Asp’ % (sparse) transpose of A
» diag(Asp) % (sparse) diagonal comp of A
» [Lsp,Usp] = lu(Asp) % (sparse) L-U factorization of A
The full command converts from sparse format to the full representation of a
matrix.
» L = full( Lsp ) % full representation of L
