<< . .

. 29
( : 59)



. . >>

e(k) g (k)
¤ κ(A) (0) , (5.16)
e(0) g
where the condition number κ(A) = A A’1 is to be computed with
respect to a matrix norm that is consistent with the chosen vector norm.
Relation (5.16) follows from
A’1 g (k) ¤ A’1
e(k) g (k) ,
=
Ae(0) ¤ A
g (0) e(0) .
=
Therefore, for the selection of δ in (5.15) we have to take into account the
behaviour of the condition number.
For the iteration matrix M , according to (5.8), we have
M = I ’ NA ,
or according to (5.10) with nonsingular W ,
M = I ’ W ’1 A .
5.1. Linear Stationary Iterative Methods 203

To improve the convergence, i.e. to reduce (M ) (or M ), we need
N ≈ A’1 and W ≈ A ,
which is in contradiction to the fast solvability of (5.10).


5.1.2 Classical Methods
The fast solvability of (5.10) (in O(m) operations) is ensured by choosing
W := D , (5.17)
where A = L + D + R is the unique partition of A, with a strictly lower
triangular matrix L, a strictly upper triangular matrix R, and the diagonal
matrix D:
«  « 
0 ··· ··· 0 0 a1,2 · · · a1,m
¬ a2,1 0 ··· 0 · ¬ . .. .. ·
.
¬ · ¬. ·
.
. .
. .
L := ¬ . . · , R := ¬ ·,
.. ..
. .  0 · · · 0 am’1,m 
. .
. .
am,1 · · · am,m’1 0 0 ··· ··· 0
« 
a11
¬ 0·
a22
¬ ·
D := ¬ ·.
..
0 
.
amm
(5.18)
Assume aii = 0 for all i = 1, . . . , m, or equivalently that D is nonsingular,
which can be achieved by row and column permutation.
The choice of (5.17) is called the method of simultaneous displacements
or Jacobi™s method. In the formulation form (5.6) we have
D’1 ,
N =
I ’ N A = I ’ D’1 A = ’D’1 (L + R) .
MJ =
Therefore, the iteration can be written as
D x(k+1) ’ x(k) = ’ Ax(k) ’ b
or
x(k+1) = D’1 ’ Lx(k) ’ Rx(k) + b (5.19)
or
« 
i’1 m
1
aij xj + bi  for all i = 1, . . . , m .
(k+1) (k) (k)
’ aij xj ’
xi =
aii j=1 j=i+1

On the right side in the ¬rst sum it is reasonable to use the new iterate
(k+1)
x where it is already calculated. This leads us to the iteration
x(k+1) = D’1 ’ Lx(k+1) ’ Rx(k) + b (5.20)
204 5. Iterative Methods for Systems of Linear Equations

or
(D + L) x(k+1) = ’Rx(k) + b
or
(D + L) x(k+1) ’ x(k) = ’ Ax(k) ’ b , (5.21)
the so-called method of successive displacements or Gauss“Seidel method.
According to (5.21) we have here a consistent iteration with
W =D+L.
Since D is nonsingular, W is nonsingular. Written in the form (5.6) the
method is de¬ned by
’1
W ’1 = (D + L)
N = ,
’1 ’1
I ’ N A = I ’ (D + L) A = ’ (D + L)
MGS = R.
In contrast to the Jacobi iteration, the Gauss“Seidel iteration depends on
the order of the equations. However, the derivation (5.20) shows that the
number of operations per iteration step is equal,
Jacobi becomes Gauss“Seidel,
if x(k+1) is stored on the same vector as x(k) .
A su¬cient convergence condition is given by the following theorem:
Theorem 5.2 Jacobi™s method and the Gauss“Seidel method converge
· ∞ if the strict row sum
globally and monotonically with respect to
criterion
m
|aij | < |aii | for all i = 1, . . . , m (5.22)
j=1
j=i


is satis¬ed.

Proof : The proof here is given only for the Jacobi iteration. For the other
method see, for example, [16].
The inequality (5.22) is equivalent to MJ ∞ < 1 because of MJ =
’D’1 (L + R) if · ∞ is the matrix norm that is induced by · ∞ , which
2
means the maximum-row-sum norm (see (A3.6)).

It can be shown that the Gauss“Seidel method converges “better” than
Jacobi™s method, as expected: Under the assumption of (5.22) for the
respective iteration matrices,
¤ MJ
MGS <1
∞ ∞

(see, for example, [16]).
5.1. Linear Stationary Iterative Methods 205

Theorem 5.3 If A is symmetric and positive de¬nite, then the Gauss“
Seidel method converges globally. The convergence is monotone in the
1/2
energy norm · A , where x A := xT Ax for x ∈ Rm .

2
Proof: See [16, p. 90].

If the di¬erential operator, and therefore the bilinear form, is symmet-
ric, that is, if (3.12) holds with c = 0, then Theorem 5.3 can be applied.
Concerning the applicability of Theorem 5.2, even for the Poisson equation
with Dirichlet boundary conditions (1.1), (1.2) requirements for the ¬nite
element discretization are necessary in order to satisfy at least a weaker
version of (5.22). This example then satis¬es the weak row sum criterion
only in the following sense:
m
|aij | ¤ |aii | for all i = 1, . . . , m ;
(5.23)
j=1
j=i

“<” holds for at least one i ∈ {1, . . . , m} .

In the case of the ¬nite di¬erence method (1.7) for the rectangular do-
main or the ¬nite element method from Section 2.2, which leads to the
same discretization matrix, (5.23) is satis¬ed. For a general triangulation
with linear ansatz functions, conditions for the angles of the elements must
be required (see the angle condition in Section 3.9). The condition (5.23)
is also su¬cient, if A is irreducible (see Appendix A.3).

Theorem 5.4 If A satis¬es the condition (5.23) and is irreducible, then
Jacobi™s method converges globally.

2
Proof: See [28, p. 111].

The qualitative statement of convergence does not say anything about
the usefulness of Jacobi™s and the Gauss“Seidel method for ¬nite element
discretizations. As an example we consider the Dirichlet problem for the
Poisson equation on a rectangular domain as in (1.5), with the ¬ve-point
stencil discretization introduced in Section 1.2. We restrict ourselves to
an equal number of nodes in both space directions for simpli¬city of the
notation. This number is denoted by n + 1, di¬erently than in Chapter 1.
Therefore, A ∈ Rm,m according to (1.14), with m = (n ’ 1)2 being the
number of interior nodes. The factor h’2 can be omitted by multiplying
the equation by h2 .
In the above example the eigenvalues and therefore the spectral radius
can be calculated explicitly. Due to D = 4I we have for Jacobi™s method
1 1
M = ’ (A ’ 4I) = I ’ A ,
4 4
206 5. Iterative Methods for Systems of Linear Equations

and therefore A and M have the same eigenvectors, namely,
ikπ jlπ
1 ¤ i, j, k, l ¤ n ’ 1 ,
z k,l = sin sin ,
ij n n
with the eigenvalues
kπ lπ
2 2 ’ cos ’ cos (5.24)
n n
for A and
1 kπ 1 lπ
cos + cos (5.25)
2 n 2 n
for M with 1 ¤ k, l ¤ n ’ 1. This can be proven directly with the help of
trigonometric identities (see, for example, [15, p. 53]). Thus we have
(n ’ 1)π π2
π
= cos = 1 ’ 2 + O n’4 .
(M ) = ’ cos (5.26)
n n 2n
With growing n the rate of convergence becomes worse. The e¬ort to gain
an approximative solution, which means to reduce the error level below a
given threshold µ, is proportional to the number of iterations — operations
for an iteration, as we discussed at the beginning of this chapter. Due to
(5.4) and (5.12) the number of necessary operations is calculated as follows:
ln(1/µ) 1 1
· O(m) = ln · O n2 · O(m) = ln O(m2 ) .
’ ln( (M )) µ µ
Here the well-known expansion ln(1 + x) = x + O(x2 ) is employed in the
determination of the leading term of ’1/(ln( (M )). An analogous result
with better constants holds for the Gauss“Seidel method.
In comparison to this, the elimination or the Cholesky method requires
O band-width2 · m = O(m2 )
operations; i.e., they are of the same complexity. Therefore, both methods
are of use for only moderately large m.
An iterative method has a superior complexity to that of the Cholesky
method if
(M ) = 1 ’ O(n’± ) (5.27)
with ± < 2. In the ideal case (5.5) holds; then the method needs O(m)
operations, which is asymptotically optimal.
In the following we will present a sequence of methods with increasingly
better convergence properties for systems of equations that arise from ¬nite
element discretizations.
The simplest iteration is the Richardson method, de¬ned by
M =I ’A, i.e., N = W = I . (5.28)
5.1. Linear Stationary Iterative Methods 207

For this method we have
(M ) = max {|1 ’ »max (A)|, |1 ’ »min (A)|} ,
with »max (A) and »min (A) being the largest and smallest eigenvalues of A,
respectively. Therefore, this method is convergent for special matrices only.
In the case of a nonsingular D, the Richardson method for the transformed
system of equations
D’1 Ax = D’1 b
is equivalent to Jacobi™s method.
More generally, the following can be shown: If a consistent method is
de¬ned by M, N with I = M +N A, and N nonsingular, then it is equivalent
to the Richardson method applied to
N Ax = N b . (5.29)
The Richardson method for (5.29) has the form

x(k+1) ’ x(k) = ’N N Ax(k) ’ N b
˜

˜
with N = I, which means the form (5.9), and vice versa.
Equation (5.29) can also be interpreted as a preconditioning of the system
of equations (5.1), with the aim to reduce the spectral condition number
κ(A) of the system matrix, since this is essential for the convergence be-
haviour. This will be further speci¬ed in the following considerations (5.33),
(5.73). As already seen in the aforementioned examples, the matrix N A will
not be constructed explicitly, since N is in general densely occupied, even
if N ’1 is sparse. The evaluation of y = N Ax therefore means solving the
auxiliary system of equations
N ’1 y = Ax .
Obviously, we have the following:
Lemma 5.5 If the matrix A is symmetric and positive de¬nite, then for
the Richardson method all eigenvalues of M are real and smaller than 1.


5.1.3 Relaxation
We continue to assume that A is symmetric and positive de¬nite. Therefore,
divergence of the procedure can be caused only by negative eigenvalues of
I ’ A less than or equal to ’1. In general, bad or nonconvergent iterative
methods can be improved in their convergence behaviour by relaxation if
they meet certain conditions.
For an iteration method, given in the form (5.6), (5.7), the corresponding
relaxation method with relaxation parameter ω > 0 is de¬ned by
ω M x(k) + N b + (1 ’ ω)x(k) ,
x(k+1) := (5.30)
208 5. Iterative Methods for Systems of Linear Equations

which means
ωM + (1 ’ ω)I ,
Mω := Nω := ωN , (5.31)
or if the condition of consistency M = I ’ N A holds,
= ω x(k) ’ N Ax(k) ’ b + (1 ’ ω)x(k)
x(k+1)
= x(k) ’ ωN Ax(k) ’ b .
Let us assume for the procedure (5.6) that all eigenvalues of M are real.
For the smallest one »min and the largest one »max we assume
»min ¤ »max < 1 ;
this is, for example, the case for the Richardson method. Then also the
eigenvalues of Mω are real, and we conclude that
»i (Mω ) = ω»i (M ) + 1 ’ ω = 1 ’ ω 1 ’ »i (M )
if the »i (B) are the eigenvalues of B in an arbitrary ordering. Hence
(Mω ) = max |1 ’ ω (1 ’ »min (M ))| , |1 ’ ω (1 ’ »max (M ))| ,
since f (») := 1 ’ ω(1 ’ ») is a straight line for a ¬xed ω (with f (1) = 1
and f (0) = 1 ’ ω).

ρ(Mω )
f
- 1 + ω ( 1 ’ »min )
1 1
» max
ω1 1 - ω ( 1 ’ »max )
ω
ω
» min » max »
» min
1
1
1 - ω ( 1 ’ »min )
ω2
..
f (») for ω1 < 1 and ω2 > 1

Figure 5.1. Calculation of ω .
¯

For the optimal ω , i.e., ω with
¯ ¯
(Mω ) = min (Mω ) ,
¯
ω>0

we therefore have, as can be seen from Figure 5.1,
1 ’ ω (1 ’ »max (M )) = ’1 + ω (1 ’ »min (M ))
¯ ¯
2
⇐’ ω=
¯ .
2 ’ »max (M ) ’ »min (M )
Hence ω > 0 and
¯
(Mω ) = 1 ’ ω (1 ’ »max (M )) < 1 ;
¯
¯
5.1. Linear Stationary Iterative Methods 209

consequently, the method converges with optimal ω even in cases where
it would not converge for ω = 1. But keep in mind that one needs the
eigenvalues of M to determine ω .
¯
Moreover, we have

ω<1
¯ »max (M ) + »min (M ) < 0 .
If »min (M ) = ’»max (M ), that is, ω = 1, we will achieve an improvement
¯
by relaxation:
(Mω ) < (M ) .
¯

The case of ω < 1 is called underrelaxation, whereas in the case of ω > 1
we speak of an overerrelaxation.
In particular, for the Richardson method with the iteration matrix M =
I ’ A, due to »min (M ) = 1 ’ »max (A) and »max (M ) = 1 ’ »min (A), the
optimal ω is given by
¯
2
ω=
¯ . (5.32)
»min (A) + »max (A)
Hence
»max (A) ’ »min (A) κ(A) ’ 1
(Mω ) = 1 ’ ω »min (A) =
¯ = < 1, (5.33)
¯
»min (A) + »max (A) κ(A) + 1
with the spectral condition number of A
»max (A)
κ(A) :=
»min (A)
(see Appendix A.3).
For large κ(A) we have
κ(A) ’ 1 2
≈1’
(Mω ) = ,
¯
κ(A) + 1 κ(A)
the variable of the proportionality being κ(A). For the example of the
¬ve-point stencil discretization, due to (5.24),
n’1 π
»min (A) + »max (A) = 4 2 ’ cos π ’ cos = 8,
n n
and thus due to (5.32),
1
ω=
¯ .
4
Hence the iteration matrix Mω = I ’ 1 A is identical to the Jacobi iteration:
¯ 4
We have rediscovered Jacobi™s method.
By means of (5.33) we can estimate the contraction number, since we
know from (5.24) that
4 1 ’ cos n’1 π 1 + cos π 4n2
≈ 2.
n n
κ(A) = = (5.34)
1 ’ cos π
4 1 ’ cos π π
n
n
210 5. Iterative Methods for Systems of Linear Equations

This shows the stringency of Theorem 3.45, and again we can conclude that
π2
π
≈1’ 2 .
(Mω ) = cos (5.35)
¯
n 2n

<< . .

. 29
( : 59)



. . >>