¤ κ(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