Iterative Solution of Linear Equations
Preface to the existing class notes
At the risk of mixing notation a little I want to discuss the general form of iterative methods a general level. We are trying to solve a linear system Ax=b, in a situation where cost of direct solution for matrix A is prohibitive. One approach starts by writing A as the sum of two matrices A=E+F. With E chosen so that a linear system in the form Ey=c can be solved easily by direct means. Now we rearrange the original linear equation to the form:
Ex = b - Fx
Introduce a series of approximations to the solution x starting with an initial guess x(0) and proceeding to the latest available approximation x(j). The next approximation is obtained from the equation:
E x(j+1) = b - F x(j)
or
x(j+1) = E-1 (b - F x(j)).
Understand that we never actually generate the inverse of E, just solve the equation for x(j+1).
To understand convergence properties of this iteration introduce the residual r(j) = b - A x(j). Using this definition and the definition of the iteration, we get the useful expression:
r(j+1) = FE-1 r(j),
or
r(j) = (FE-1)j r(j).
If you know something about linear algebra this tells you that the absolute values of eigenvalues of FE-1 should be less than one.
Another approach is to choose a matrix M such that the product of M-1 and A is close to the identity matrix and the solution of an equation like My=c can be obtained without too much effort. The revised problem can be cast either in the form
M-1Ax = M-1b
Or
AM-1 y = b where x =
Beginning of normal ME540 Notes
a) Introduction
i) Procedure – general
1) Assume initial values for the variable field
2) Use the nodal finite difference equations one at a time or in groups to obtain an improved value for the variables. Repeat the procedure until a converged solution is obtained.
ii) Advantages
Iterative methods are used because:
Requires less computer storage.
Is reasonably computational fast.
There is, however, no guarantee that a converged solution can be obtained.
iii) Symbols to be used
exact solution of differential equation
exact solution of finite difference equations
computer solution to finite difference equations
iv) Convergence
The solution is said to be converged when
Limit
where j is an index for the iteration number.
b) Procedure – details
i) Rearrange coefficient matrix
Solution required H
Modify H matrix by taking each equation and dividing it by diagonal coefficient to form a “C” matrix
,
New Expression
The C Matrix will have all its diagonal elements with the value 1.
ii) Iterative set of equations
initial value of variable
for calculation
residual
Introduce identity matrix, I
Add to both sides of
to
Transpose and let on the right hand
side be
Thus or
The iterative equation can be expressed as
C = I - B
Introduce a lower (L) and upper (U) matrix
This is known as the Jacobi method
Example
Iterative equation for node 5
iii) Error propagation and convergence
1) Error propagation
known
Introduce computational error
Note
So
Back substitute to obtain
Indicating that errors are propagated in an identical manner
To the propagation of
2) Convergence
Limit
residual
Convergence checks
L2 Norm
Limit (with machine accuracy)
Other tests
Maximum
(Local)
Maximum
(Global)
or
Will the iterative solution selected converge?
Error propagation
To answer this question we look at the characteristics roots of
the B matrix.
Latent root of “B” matrix
vector such that the equality is satisfied
Re-arranging
Non-trivial solution
Expand as a polynomial
The characteristic root are the ’s
(Eigenvalues)
Maximum root is called the spectral radius of the matrix
The polynomial will converge to zero if:
<1
Original expression (page 74 of class notes)
if <1 B < 1
i.e.
so general expression
so
Limit
exact solution
so solution converges if
Limit i.e. < 1
c) Influence of search pattern
i) Jacobi
You go through the complete grid pattern with
very slow convergence
ii) Gauss Seidel Method
The new values of
Where the equations are listed in the order of search.
The L and the U depend of the search pattern.
Rewrite expression
Note standard equation
So U is the “B” matrix for the Gauss Seidel method.
i.e. the latent roots of (I – L)-1 U must be less than 1.
Also the error propagation is
Example assume Dx = Dy
L U
Equations for Gauss Seidel
Node 1
Node 5
Change the search pattern
Node 1
Node 5
iii) Extrapolated Liebmann method
Subtract
Multiply that change by (over-under relaxation factor)
Error
Note If we have the Gauss Seidel method.
Example Use original search pattern.
Node 3
Node 6
iv) Other important characteristics of the “C” Matrix
1) Property “A”
If
· Differential equation does not have mixed derivative terms.
· The region is subdivided using a rectangular grid pattern.
And the nodes are labelled in an alternating fashion.
Black Nodes
White Nodes
The matrix will have property “A” if
Consistent ordering
If one has consistent ordering the values of the
Variable field at iteration (j+1) is independent of
Search pattern used.
Simple test: Five point approximation of -
Matrix (“C”) has property A.
1. Select a search pattern.
2. Move node by node through the grid following
the selected search pattern. Look at the four surrounding
nodes. If they are not already connected to the node,
draw an arrow from it with the head of the arrow located at
the surrounding node.
3. Move around a rectangle created by the grid lines. If you
find two arrows going with you and two arrows going
against you, the search pattern has consistent ordering.
If three arrows are with you and one against (or 3 against
and 1 with you) you do not have consistent ordering.
If the “C” matrix is consistent and has property “A”
a relationship can be found between eigenvalues
of the “B” matrix for the Jacobi, Gauss-Seidel and
the over-under relaxation method.
Poisson Equation
Square grid
known on all boundaries
Spectral radius for Jacobi method
Example
If 10 x 10 mesh P = 10 q = 10
For the Gauss Seidel method
thus
4 x 4
and
10 x 10
Asympototic rate of convergence
[A]
so
and
For normal convergence
Using equation [A]
Larger slower convergence
Smaller faster convergence
Note It is assumed that j is very large.
Asympototic rate of convergence
CR
Look at examples
CR
4 x 4 10 x 10
Jacobi 0.346 0.051
Gauss-Seidel 0.695 0.105
Gauss-Seidel method converges twice as fast as the Jacobi
method.
Note The smaller the large CR i.e. the faster the
Convergence
v) Optimum over-under relaxation factor
It has been shown that for an interior node
with the smallest root of
where
Example
4 x 4
10 x 10
For large values of P and q
100 x 100
What happens if you do not have only interior nodes or node
with property “A” and consistent ordering?
Use
Implementation run program with (Gauss Seidel)
Assume
i.e.
Estimate
Limit
where
or
or
Figure 1 Influence of relaxation factor on number iterations to converge