Iterative Solvers

Overview: Finite Difference Methods

1. Cover the domain Ω ⊂ Rd with a Cartesian grid.

2. Write down the PDE at the grid point xi . Replace all
derivatives by finite differences.
3. Write down the resulting difference equations for all grid
points, using the boundary conditions.
4. Arrange the difference equations into a system of linear
equations.
5. Solve the system of linear equations.

1
Elliptic Model Problem (d = 2)

Example 3.1 (Poisson Equation)

We consider the unit square domain Ω = (0, 1)2

−∆u(x, y ) = f (x, y ), (x, y ) ∈ Ω,

with the continuous function f : Ω → R and the homogeneous

Dirichlet boundary condition

u(x, y ) = 0 (x, y ) ∈ ∂Ω.

The Laplace operator ∆ is discretized by a five-point stencil

(order of convergence m = 2) on a uniform grid with mesh
size h = 1/(Nh + 1).
Difference equations: Lh u h = f h
2
Content

3. Iterative Solvers
3.1 Properties of the Difference Matrix
3.2 Direct Methods and their Disadvantages
3.3 Linear Iterative Solvers
3.4 Practical Examples
3.5 Eigenvalues and Eigenvectors of the Difference Matrix
3.6 Back to Splitting Methods
3.7 Smoothing Property
3.8 Two-Grid Methods

3
High Dimensionality

The dimension nh ∼ (1/h)d of the difference matrix Lh depends

on the mesh size h and the spatial dimension d!
Convergence analysis (see Section 2.3): Depending on the
desired accuracy, h must be chosen small enough.

Thought Experiment
For quadratic convergence ∥u − u h ∥ ≲ h2 we want to reduce
the error by a factor 10. To do this, h must be reduced by a
√ √ d
factor of 10 and we obtain 10 = 10 times as many
difference equations for d = 2.
With linear convergence ∥u − u h ∥ ≲ h, on the other hand, we
get 10d = 100 times more equations for d = 2.

4
Sparsity, Memory Requirements

Finite difference stencils usually couple only neighboring grid

points. Therefore Lh has at most 5 non-zero entries per line.
The sparse format in Matlab requires ∼ nh memory (zeros are
not stored!)

0 108

10

106
20

30
104

40

102
50

60
100
0 10 20 30 40 50 60 100 101 102 103 104
nz = 288

5
Condition Number

Definition 3.2
Let A ∈ Rn×n be an invertible matrix. Then κ(A) := ∥A∥∥A−1 ∥
is called the condition number of A.

The condition number of a matrix is a measure of the error

amplification when solving the linear system Ax = b.
Example: Consider Ax = b and Ax e=b
e with the perturbed
right-hand side b.
e Then it holds:

∥x − x
e∥ ∥b − b∥
e
≤ κ(A) · (3.1)
∥x∥ ∥b∥

6
Condition Number of Lh

In Section 2.3 it is shown: ∥L−1

h ∥∞ ≤
1
8 independent of h
Row Sum Norm: ∥Lh ∥∞ = (4 · | − 1| + 4)/h2 = 8/h2
In total: κ∞ (Lh ) = ∥Lh ∥∞ ∥L−1 2
h ∥∞ ≤ 8/h · (1/8) = 1/h
2

104

103

102

101

100
10-2 10-1 100
7
Note: The bound in (3.1) is a worst-case bound, i.e., it holds for
all perturbations ∥f h − e
f h ∥.
For Example 3.1 with f (x, y ) = sin(πx) sin(πy ) and
ef (x, y ) = sin(10πx) sin(10πy ) we get:

104

103

102

101

100

10-1
10-2 10-1 100

8
Symmetry, Positive Definiteness

Symmetry Lh = L⊤
h for symmetric stencils, e.g. the five-point
stencil.
Positive Definiteness: x ⊤ Lh x > 0 for all vectors x ̸= 0
Therefore the solution of Lh u h = f h is also the uniquely
determined minimum of the function
1 ⊤
ϕ(x) = x Lh x − x ⊤ f h ,
2
since ∇ϕ(x) = Lh x − f h and the Hessian matrix is H(x) = Lh .
This is important for advanced iterative solvers, e.g., the
conjugate gradient (CG) method.

9
Band Structure

Definition 3.3
A matrix A ∈ Rn×n is called (p, q)-band matrix, if p, q ∈ N exist
such that ai,k = 0 for all k > i + p and i > k + q.
The number w := p + q + 1 is called band width.
A matrix with p = q = 1 is called tridiagonal matrix.

10
Ordering of the Unknowns

Example for Nh = 3, thus 3 · 3 = 9 inner nodes.

(0, 4) (4, 4)

(0, 0) (4, 0)

11
Lexicographic Ordering

Example for Nh = 3, thus 3 · 3 = 9 inner nodes.

(1, 3) (2, 3) (3, 3) 7 8 9

(1, 2) (2, 2) (3, 2) 4 5 6

(1, 1) (2, 1) (3, 1) 1 2 3

Mapping: (i, k ) 7→ ℓ(i, k ) = i + (k − 1)Nh

12
Lexicographic Ordering:
Lh is a (3, 3)-band matrix with a band width of 3 + 3 + 1 = 7:
 
4 −1 −1
 −1 4 −1 −1 
 
−1 4 −1
 
 
 
 −1 4 −1 −1 
1 
 
Lh = 2  −1 −1 4 −1 −1

h 



 −1 −1 4 −1 


 −1 4 −1 

−1 −1 4 −1 
 

−1 −1 4

Further, Lh is a 3 × 3 block matrix with 3 × 3 blocks.

13
Checkerboard Ordering

Example for Nh = 3, thus 3 · 3 = 9 inner nodes.

(1, 3) (2, 3) (3, 3) 4 9 5

(1, 2) (2, 2) (3, 2) 7 3 8

(1, 1) (2, 1) (3, 1) 1 6 2

14
Checkerboard Ordering:
Lh is a (6, 6)-band matrix with a band width of 6 + 6 + 1 = 13:
 
4 −1 −1

 4 −1 −1 


−1 −1 −1 −1 


 4 
 4 −1 −1 
1  
Lh = 2  4 −1 −1 
 
h  
 −1 −1 −1 4 
 
 −1 −1 −1 4 
 
−1 −1 −1
 
 4 
−1 −1 −1 4

Further, Lh is a 2 × 2 block matrix.

15
 Lexicographic ordering Checkerboard ordering
0 0

100 100

200 200

300 300

400 400

500 500

600 600

700 700

800 800

900 900

0 200 400 600 800 0 200 400 600 800

nz = 4681 nz = 4681

Random ordering
0

100

200

300

400

500

600

700

800

900
16
0 200 400 600 800
nz = 4681
Summary

The matrix Lh in the difference equation in Example 3.1 is

• high dimensional (many rows and columns),

• sparse,
• ill-conditioned (large condition number),
• symmetric and positive-definite,
• a band matrix with small band width (with suitable ordering
of the unknowns).

Question
Which properties can be exploited in the solver? Which
properties are problematic?

17
Content

3. Iterative Solvers
3.1 Properties of the Difference Matrix
3.2 Direct Methods and their Disadvantages
3.3 Linear Iterative Solvers
3.4 Practical Examples
3.5 Eigenvalues and Eigenvectors of the Difference Matrix
3.6 Back to Splitting Methods
3.7 Smoothing Property
3.8 Two-Grid Methods

18
Gaussian Elimination

Task: Solve the linear system Ax = b with the invertible matrix

A ∈ Rn×n .
Idea: Decompose A into the product of two triangular matrices,
where L ∈ Rn×n is a lower triangular matrix with ℓi,i = 1 and
U ∈ Rn×n is an upper triangular matrix such that it holds
PA = LU
with a permutation matrix P ∈ Rn×n , P −1 = P ⊤ .
Solve the linear system as follows:
1st step (Forward elimination): Lz = Pb
2nd step (Backward elimination): Ux = z
Check: Pb = Lz = LUx = PAx, thus it holds b = Ax.
19
Computational Cost

dense matrix (p, q)-band matrix

LU Decomposition ∼ n3 ∼ (2pq + 1)n
Forward Elimination ∼ n2 ∼ pn
Backward Elimination ∼ n2 ∼ qn

Memory Requirements
Without column pivot search: The factor L has the (lower) band
width p and the factor U has the (upper) band width q.
With column pivot search: U has the (upper) band width p + q.
However, the number of non-zero entries in L and U can be
(considerably) larger than the number of non-zero entries in A
(so-called fill-in).

20
 Lower triangular matrix L - Lexicographic ordering Lower triangular matrix L - Checkerboard ordering
0 0

100 100

200 200

300 300

400 400

500 500

600 600

700 700

800 800

900 900

0 200 400 600 800 0 200 400 600 800

nz = 29821 nz = 17009

Lower triangular matrix L - Random ordering

0

100

200

300

400

500

600

700

800

900
21
0 200 400 600 800
nz = 68278
Is There a Better Way to Order the Unknowns?

The permutation matrix P swaps the rows of A (if multiplied

from the left, PA) or the columns of A (if multiplied from the
right, AP).
For all permutation matrices it holds P −1 = P ⊤ .

Ax = b ⇐⇒ PAx = Pb
Ax = b ⇐⇒ AP x
b = b, x = P x
b

Idea: Construct a permutation matrix such that the matrix

PAP ⊤ has the smallest possible bandwidth.
Example: Reverse Cuthill–McKee–Method (e.g. in Matlab)

22
 Lexicographic ordering Reverse Cuthill-McKee ordering
0 0

100 100

200 200

300 300

400 400

500 500

600 600

700 700

800 800

900 900

0 200 400 600 800 0 200 400 600 800

nz = 4681 nz = 4681

Lower triangular matrix L - Lexicographic ordering Lower triangular matrix L - Reverse Cuthill-McKee orderi
0 0

100 100

200 200

300 300

400 400

500 500

600 600

700 700

800 800

900 900
23
0 200 400 600 800 0 200 400 600 800
nz = 29821 nz = 21266
Comparison of Computational Times

0.05

0.04

0.03

0.02

0.01

0
0 0.5 1 1.5 2
4 24
10
Iterative Solvers for Difference Equations

What is the basic idea?

The linear system of equations Lh u h = f h (discretized
boundary value problem) is solved approximately, i.e., we
construct an approximation u ∗ ≈ u h with ∥u ∗ − u h ∥ sufficiently
small.
Only matrix-vector multiplications with Lh are performed (low
memory requirements!)

Which tools are used in the analysis?

• Banach fixed-point theorem
• Eigenvalues and eigenvectors of the matrix Lh

25
Content

3. Iterative Solvers
3.1 Properties of the Difference Matrix
3.2 Direct Methods and their Disadvantages
3.3 Linear Iterative Solvers
3.4 Practical Examples
3.5 Eigenvalues and Eigenvectors of the Difference Matrix
3.6 Back to Splitting Methods
3.7 Smoothing Property
3.8 Two-Grid Methods

26
Definition 3.4
Let A ∈ Rn×n be an invertible matrix and b ∈ Rn . An iterative
solver for Ax = b constructs a sequence of approximations
x 1 , x 2 , . . . , x k , . . . in Rn starting from a given initial vector
x 0 ∈ Rn .
Definition 3.5
An iterative solver is called convergent if

lim x k = x,
k →+∞

where x = A−1 b is the exact solution of Ax = b.

Definition 3.6
The vector e k = x − x k is called error and the vector
r k = b − Ax k is called residual of the k -th approximation x k ,
k = 0, 1, 2, . . . .
27
Definition 3.7
An iterative solver is called linear, if matrices S, T ∈ Rn×n
exist such that it holds

x k +1 = Sx k + T b, k = 0, 1, 2, . . . . (3.2)

Definition 3.8
An iterative solver is called consistent, if from x k = x it follows
that x k +1 = x holds.
Theorem 3.9
An linear iterative solver is consistent ⇐⇒ S = I − TA.

Proof:
x k +1 = Sx k + T b = Sx + TAx = (S + TA)x = x,
if and only if S + TA = I or, equivalently, S = I − TA.
28
Error Propagation

The behavior of the error sequence {e k }∞

k =0 is crucial for the
convergence of an iterative solver.
Question: How does the error propagate from one iterate to
the next?
Assumption (A1): We consider a consistent, linear iterative
solver. It holds:
e k +1 = x − x k +1
= x − Sx k − T b (linear iterative solver)
k
= x − Sx − TAx (Ax = b)
= (I − TA)x − Sx k
= Sx − Sx k (consistency: S = I − TA)
= Se k = S(Se k −1 ) = · · · = S k +1 e 0 .
29
Definition 3.10
The matrix S = I − TA is called iteration matrix or error
propagation matrix.

Definition 3.11
Let S ∈ Rn×n be a matrix. The number

ρ(S) := max{|λℓ | : λℓ is an eigenvalues of S}

is called spectral radius of S.

30
Theorem 3.12 (Convergence)
Under assumption (A1) the linear iterative solver in (3.2) is
convergent for any initial vector x 0 ⇐⇒ ρ(S) < 1, i.e., the
spectral radius of the iteration matrix is smaller than one.

Sketch of Proof:
Assumption: There is a basis of Rn containing only
eigenvectors of S. (The general case is not considered here.)
Let (λℓ , v ℓ ), ℓ = 1, . . . , n denote the eigenpairs of the matrix S.
In particular, Sv ℓ = λℓ v ℓ , ℓ = 1, . . . , n.
Then the initial error e 0 has the representation

e 0 = c1 v 1 + · · · + cn v n (3.3)

with suitable coefficients c1 , . . . , cn ∈ R.

31
Insertion of e 0 from (3.3) in e k = S k e 0 gives the representation
e k = S k e 0 = S k (c1 v 1 + · · · + cn v n )
= c1 S k v 1 + · · · + cn S k v n
= c1 λk1 v 1 + · · · + cn λkn v n .
With this it follows for ρ(S) < 1,
∥e k ∥ ≤ |c1 ||λk1 |∥v 1 ∥ + · · · + |cn ||λkn |∥v n ∥
≤ ( max |λℓ |)k (|c1 |∥v 1 ∥ + · · · + |cn |∥v n ∥)
ℓ=1,...,n

= ρ(S)k (|c1 |∥v 1 ∥ + · · · + |cn |∥v n ∥) → 0 for k → +∞.

| {z }
<1
Conversely, assume that ρ(S) = |λℓ | ≥ 1 for an index ℓ.
Choosing the initial vector x 0 so that e 0 = v ℓ , it then follows
∥e k ∥ = ∥S k e 0 ∥ = |λℓ |k ∥v ℓ ∥ ≥ ∥e 0 ∥
for all k ∈ N. Thus the iteration does not converge which
32
contradicts the assumption.
A Practical Convergence Criterion

In general, it is difficult to calculate the spectral radius ρ(S). (An

estimation is possible with the Gershgorin circles of S.)
However, from e k = S k e 0 it follows immediately that
∥e k ∥ ≤ ∥S k ∥∥e 0 ∥ ≤ ∥S∥k ∥e 0 ∥.

Conclusion
Under assumption (A1), the iterative solver (3.2) is convergent
for any initial vector x 0 if ∥S∥ < 1 holds for any induced matrix
norm ∥ · ∥, e.g., the row sum norm.

33
Speed of Convergence
Question
How many iteration steps k are needed to reduce the norm of
the initial error by a factor R > 1?

Answer: With
1/k  1/k
∥e k ∥

1
σk := =
∥e 0 ∥ R
it follows
k = − ln(R)/ ln(σk ).

σk is a measure for the average error reduction in k iteration

steps.
But: To calculate σk we need the exact error!
34
Alternative:
From ∥e k ∥ = ∥S k e 0 ∥ ≤ ∥S k ∥∥e 0 ∥ it follows σk ≤ ∥S k ∥1/k .
It can be shown that ∥S k ∥1/k → ρ(S) for k → +∞.
Thus, the spectral radius ρ(S) is a meaningful measure for the
asymptotic error reduction of the iterative solver.
The number
Iter (S) := −⌈ln(R)/ ln(ρ(S))⌉

indicates approximately how many iteration steps we need to

reduce the norm of the initial error by a factor R > 1.
Conclusion: If ρ(S2 ) = ρ(S1 )2 < 1, it holds

Iter (S2 ) = Iter (S1 )/2.

35
Content

3. Iterative Solvers
3.1 Properties of the Difference Matrix
3.2 Direct Methods and their Disadvantages
3.3 Linear Iterative Solvers
3.4 Practical Examples
3.5 Eigenvalues and Eigenvectors of the Difference Matrix
3.6 Back to Splitting Methods
3.7 Smoothing Property
3.8 Two-Grid Methods

36
Let A ∈ Rn×n be an invertible matrix.
Question: How can we construct the matrices S, T in a linear
iterative solver for Ax = b,

x k +1 = Sx k + T b, k = 0, 1, . . .?

Desired properties:

• Consistency: S = I − TA
• Matrix-vector multiplications with S, T are cheap
• Fast convergence, i.e. small spectral radius ρ(S) ≪ 1

37
Splitting Method

Let T be invertible. A consistent method satisfies S = I − TA or,

equivalently, A = T −1 (I − S).
Idea: Splitting of A = M − (M − A), where M is invertible.

A = M − (M − A)
A = M(I − M −1 (M − A))
M (I − (I − M −1 A))
A = |{z}
| {z }
=:T −1 =:S

Conclusion
Choosing T = M −1 and S = I − M −1 A gives a consistent
linear iterative solver.

38
Splitting Method

x k +1 = Sx k + T b
x k +1 = (I − M −1 A)x k + M −1 b
x k +1 = x k + M −1 (b − Ax k )
x k +1 = x k + M −1 r k

Implementation
Given x k , b, A, and M, we carry out the following steps:
1. Calculate r k = b − Ax k
2. Solve the linear system Mz k = r k
3. Update x k +1 = x k + z k
This implementation requires only matrix-vector-multiplications with A and
linear system solves with M, no multiplications with T = M −1 or
S = I − M −1 A. M −1 is not computed, because this is often too expensive. 39
How to choose the matrix M?

Criteria:
• M is invertible and M ≈ A.
• The computational cost to solve a linear system with M is
smaller than the cost to solve a linear system with A.
Example: M is a diagonal or triangular matrix.
Idea of the so called Classical Iterative Solvers:
A=L+D+U
 
a1,1 a1,2 ··· a1,n−1 a1,n

 a2,1 a2,2 ··· a2,n−1 
 a2,n
 .. .. ..  ..
A=
 . . . 
 .
an−1,1 an−1,2 · · · an−1,n−1 an−1,n 
 

an,1 an,2 · · · an−1,n an,n

40
Jacobi Method, Total-Step Process:
M = D contains all entries on the main diagonal of A
Gauss–Seidel Method, Single-Step Process:
M = L + D contains all entries on and below the main diagonal
of A
Backward Gauss–Seidel Method:
M = D + U contains all entries on and above the main diagonal
of A

Note
In all three cases M is invertible if A and D are invertible.

41
Convergence Criteria (Without Proof)

• If both A and 2D − A are symmetric positive-definite, then

ρ(SJac ) = ρ(I − D −1 A) < 1.

• If A is symmetric positive definite, then

ρ(SGS ) = ρ(I − (L + D)−1 A) < 1.

• Other convergence criteria in the literature consider

irreducible diagonally dominant matrices.

42
Matlab Demonstration Example 3.1

f (x, y ) = 2π 2 sin(πx) sin(πy ), 0 < x, y < 1, Nh = 40, x 0 = 0

100 100
Relative norm of residual

Relative norm of residual

10-2 10-2

10-4 10-4

10-6 10-6
0 1000 2000 3000 4000 5000 0 500 1000 1500 2000 2500
Number of iterations Number of iterations

Plot of ∥f h − Lh u kh ∥2 /∥f h − Lh u 0h ∥2 against the iteration number

k for a fixed mesh size h = 1/41

43
Matlab Demonstration Example 3.1

Number of iterations to reduce the norm of the initial residual by a

factor of 10−6 :

h Jacobi Gauss–Seidel
2−2 40 21
2−3 175 89
2−4 713 357
2−5 2863 1433
2−6 11463 5733

Observations:
• The rate of convergence slows down as the mesh size h gets
smaller.
• The number of iterations approximately quadruples when the
mesh size is halved.
• Gauss–Seidel requires about half as many iterations as Jacobi.
44
Matlab-Demonstration Example 3.1

0.9

0.8

0.7

0.6
10-2 10-1
45
Jacobi Method for Lh u h = f h

The iteration matrix is Sh = I − diag(Lh )−1 Lh = I − (h2 /4)Lh .

Example for Nh = 3:
 
0 1/4 1/4
 1/4 0 1/4 1/4 
 

 1/4 0 1/4 

 1/4 0 1/4 1/4

 
Sh =  1/4 1/4 0 1/4 1/4
 

 
 1/4 1/4 0 1/4 
 
 1/4 0 1/4 
 
 1/4 1/4 0 1/4 
1/4 1/4 0

The largest Gerschgorin circle is {z ∈ C : |z| ≤ 1}.

Because Sh⊤ = Sh , all eigenvalues of Sh are real and thus in the
interval [−1, 1] (independent of h!)
46
Let (λℓ , v ℓ ) be the eigenvalues and eigenvectors, respectively,
of the difference matrix Lh ∈ Rnh ×nh , i.e., Lh v ℓ = λℓ v ℓ ,
ℓ = 1, . . . , nh .
The iteration matrix
Sh = I − diag(Lh )−1 Lh
of the Jacobi method has the same eigenvectors v ℓ with the
corresponding eigenvalues 1 − 41 h2 λℓ .
In particular, the spectral radius is
 
1 2
ρ(Sh ) = max 1 − h λℓ : λℓ is eigenvalue of Lh .
4

Proof:
 

−1
 1 1 2
Sh v ℓ = I − diag(Lh ) Lh v ℓ = v ℓ − h2 Lh v ℓ = 1 − h λℓ v ℓ
4 4
47
Content

3. Iterative Solvers
3.1 Properties of the Difference Matrix
3.2 Direct Methods and their Disadvantages
3.3 Linear Iterative Solvers
3.4 Practical Examples
3.5 Eigenvalues and Eigenvectors of the Difference Matrix
3.6 Back to Splitting Methods
3.7 Smoothing Property
3.8 Two-Grid Methods

48
Sneak Preview: Eigenvectors of Lh on the Grid

Nh = 63, Arrangement: 0 < λ1 < λ2 < · · · < λℓ < · · ·

49
Eigenfunctions of the Laplace Operator

The functions
eνµ (x, y ) = sin(νπx) sin(µπy ), ν, µ ∈ N
are equal to zero on the boundary of the domain Ω = (0, 1)2 .
Furthermore, for all (x, y ) ∈ Ω it holds
−∆eνµ (x, y ) = −∂xx eνµ (x, y ) − ∂yy eνµ (x, y )
= (νπ)2 sin(νπx) sin(µπy ) + (µπ)2 sin(νπx) sin(µπy )
 
= (νπ)2 + (µπ)2 eνµ (x, y )

= λνµ eνµ (x, y ).

Conclusion
The functions eνµ (x, y ), ν, µ ∈ N, are eigenfunctions of the
operator −∆ on the unit square (with zero-boundary
50
conditions) with associated eigenvalue λνµ .
Trigonometric Identities:

sin(νπ(x ± h)) = sin(νπx) cos(νπh) ± cos(νπx) sin(νπh)

sin(µπ(y ± h)) = sin(µπy ) cos(µπh) ± cos(µπy ) sin(µπh)
Summing Up:

sin(νπ(x + h)) sin(µπy ) + sin(νπ(x − h)) sin(µπy )

= sin(νπx) cos(νπh) sin(µπy ) + cos(νπx) sin(νπh) sin(µπy )
+ sin(νπx) cos(νπh) sin(µπy ) − cos(νπx) sin(νπh) sin(µπy )
=2 cos(νπh) sin(νπx) sin(µπy )
Analogously:

sin(νπx) sin(µπ(y + h)) + sin(νπx) sin(µπ(y − h))

=2 cos(µπh) sin(νπx) sin(µπy )

51
Apply the five-point stencil on eνµ (x, y ):
−eνµ (xi − h, yk ) − eνµ (xi + h, yk ) − eνµ (xi , yk − h)
−∆h eνµ (x, y ) =
h2
νµ νµ
−e (xi , yk + h) + 4e (xi , yk )
+
h2
1
= 2 (4 − 2 cos(νπh) − 2 cos(µπh)) sin(νπx) sin(µπy )
h  
4 1 1 1 1
= 2 − cos(νπh) + − cos(µπh) eνµ (x, y )
h 2 2 2 2
4  2 
= 2 sin (νπh/2) + sin2 (µπh/2) eνµ (x, y )
h
= λνµ νµ
h e (x, y )

Conclusion
The function eνµ (x, y ) is an eigenfunction of the discrete
Laplace operator −∆h on the unit square (with zero-boundary
conditions) with associated eigenvalue λνµ h . 52
Plot of the eigenvalues λνµ und λνµ
h for different mesh sizes

500 500

400 400

300 300

200 200

100 100

0 0
0 5 10 15 20 25 30 0 5 10 15 20 25 30

500 500

400 400

300 300

200 200

100 100

0 0
0 5 10 15 20 25 30 0 5 10 15 20 25 30

53
 νµ νµ
Idea: Define the vector vℓ(i,k ) = e (xi , yk ), i, k = 1, . . . , Nh .

Goal: Show that the vectors v νµ with ν, µ = 1, . . . , Nh are an

2
orthogonal basis of RNh .
Exercise:

Nh
X ′ ′
0, (ν, µ) ̸= (ν ′ , µ′ ),
eνµ (xi , yk )eν µ (xi , yk ) =
 1 h−2 , (ν, µ) = (ν ′ , µ′ ).
i,k =1 4

Because Lh = L⊤ h is symmetric and therefore has an ONB of

Nh · Nh eigenvectors with associated real eigenvalues, we have
thus found all the eigenvectors of Lh .
Note
Since λνµ
h > 0, Lh is a positiv-definite matrix for all h > 0.

54
Content

3. Iterative Solvers
3.1 Properties of the Difference Matrix
3.2 Direct Methods and their Disadvantages
3.3 Linear Iterative Solvers
3.4 Practical Examples
3.5 Eigenvalues and Eigenvectors of the Difference Matrix
3.6 Back to Splitting Methods
3.7 Smoothing Property
3.8 Two-Grid Methods

55
Jacobi Method for Lh u h = f h

Recall: The iteration matrix Sh = I − (h2 /4)Lh has the spectral

radius
 
1 2 νµ νµ
ρ(Sh ) = max 1 − h λh : λh is eigenvalue of Lh .
4
1
1 − h2 λνµ 2 2
h = 1 − (sin (νπh/2) + sin (µπh/2))
4
= 1 − ((1 − cos(νπh))/2 + (1 − cos(µπh))/2)
= (cos(νπh) + cos(µπh))/2
Thus it holds (ν = µ = 1)

ρ(Sh ) = cos(πh) = 1 − π 2 h2 /2 + O(h4 )

and ρ(Sh ) → 1 for h → 0.

56
Matlab Demonstration Example 3.1

f (x, y ) = 2π 2 sin(πx) sin(πy ), 0 < x, y < 1, x 0 = 0

Number of iterations k of the Jacobi method to reduce the
norm of the initial error by the factor R = 106 :

h k Iter (Sh )
2−2 39 40
2−3 174 175
2−4 712 713
2−5 2863 2864
2−6 11462 11463

Iter (Sh ) = −⌈ln(R)/ ln(ρ(Sh ))⌉ = −⌈ln(106 )/ ln(cos(πh))⌉

ln(cos(πh)) ≈ ln(1 − π 2 h2 /2) ≈ −π 2 h2 /2
Conclusion: Iter (Sh/2 ) ≈ 4 × Iter (Sh ) 57
Gauss–Seidel Method for Lh u h = f h

The difference matrix Lh = Uh⊤ + Dh + Uh and the iteration

matrices Sh,Jac = I−Dh−1 Lh and Sh,GS = I−(Uh⊤ + Dh )−1 Lh
satisfy the relation
ρ(Sh,GS ) = ρ(Sh,Jac )2 .

Thus, Iter (Sh,GS ) = Iter (Sh,Jac )/2.

ρ(Sh,GS ) = cos2 (πh) ≈ (1 − π 2 h2 /2)2 ≈ 1 − π 2 h2 + O(h4 )
ln(1 − π 2 h2 ) ≈ −π 2 h2
Conclusion: Iter (Sh/2,GS ) ≈ 4 × Iter (Sh,GS )
Observation
The convergence speed of Gauss-Seidel with respect to the
mesh size h is as bad as for the Jacobi method!
58
Goal: Acceleration of the Gauß–Seidel (GS) method
Matrix Decomposition A = L + D + U
Iterates of the GS method in vector notation:

x k +1,GS = x k + (D + L)−1 (b − Ax k )
(D + L)x k +1,GS = (D + L)x k + (b − (D + L)x k − Ux k )
Dx k +1,GS = b − Lx k +1,GS − Ux k
x k +1,GS = D −1 (b − Lx k +1,GS − Ux k )

Component-wise for i = 1, . . . , n:
!
1
xik +1,GS ai,ℓ xℓk +1,GS
X X
= bi − − ai,ℓ xℓk
ai,i
ℓ<i ℓ>i

59
Relaxation parameter ω ≥ 1, Successive Over Relaxation:
xik +1,SOR = xik + ω(xik +1,GS − xik ), i = 1, . . . , n,
!
ω
xik +1,SOR ai,ℓ xℓk +1,SOR −
X X
= xik + bi − ai,ℓ xℓk − ai,i xik .
ai,i
ℓ<i ℓ>i
Vector notation:
x k +1,SOR = x k + ωD −1 (b − Lx k +1,SOR − Ux k − Dx k )
(D + ωL)x k +1,SOR = Dx k + ωb − ωUx k − ωDx k
= Dx k + ωLx k + ωr k
= (D + ωL)x k + ωr k
x k +1,SOR = x k + ω(D + ωL)−1 r k
This is a splitting method (p. 39) with
M = ω −1 (D + ωL) = ω −1 D + L.

60
SOR Method

The iteration matrix of the SOR method is

SSOR = I − M −1 A = I − (ω −1 D + L)−1 A.

Question: How to chose the parameter ω?

Theorem 3.13
Let A ∈ Rn×n be symmetric and positive-definite. Then
ρ(SSOR ) < 1 for all 0 < ω < 2.

Unfortunately, the optimal value ωopt with the smallest spectral

radius ρ(SSOR ) is strongly problem dependent and often
unknown.

61
SOR Method for Lh u h = f h

Theorem 3.14
Let Lh be the difference matrix of the discretized Poisson
equation and let ρJac be the spectral radius of the iteration
matrix of the Jacobi method. Then the spectral radius of the
iteration matrix of the SOR method is minimal for
 2
2 ρ
ωopt = q =1+ qJac 
2
1 + 1 − ρJac 1 + 1 − ρJac2

and the optimal value given by ρSOR,opt = ωopt − 1.

62
With ρJac = cos(πh) we get
q q
1 − ρ2Jac = 1 − cos2 (πh) = sin(πh)
and thus
2
1 − sin2 (πh)

cos(πh) 1 − sin(πh)
ρSOR,opt = ωopt − 1 = = =
1 + sin(πh) (1 + sin(πh))2 1 + sin(πh)
2 sin(πh)
=1− ≈ 1 − 2πh.
1 + sin(πh)

Conclusion: For ω = ωopt it holds ln(1 − 2πh) ≈ −2πh and

thus
Iter (Sh/2,SOR ) = 2 × Iter (Sh,SOR ),
i.e., the number of iterations doubles when the mesh size is
halved. This is a significant improvement compared to the
Jacobi or GS method, but is still not optimal (h-independent).
63
Matlab Demonstration Example 3.1

f (x, y ) = 2π 2 sin(πx) sin(πy ), 0 < x, y < 1, x 0 = 0

Number of iterations of the SOR method to reduce the norm of

the initial residual by a factor of 106 :

h ω=1 ω = 1.5 ω = ωopt ωopt

2−2 21 22 11 1.1716
2−3 89 23 23 1.4465
2−4 357 112 47 1.6735
2−5 1433 471 95 1.8215
2−6 5733 1905 189 1.9065
2−7 22934 7638 379 1.9521
2−8 – 30573 758 1.9758

64
Computational Cost of the Splitting Method for Lh u h = f h

Total cost = cost per step × number of iterations.

1. What is the cost of an iteration step (p.39)?
2. How many iterations are required?
Remember: The dimension of Lh is nh ∼ (1/h)2 for d = 2.
For all methods (Jacobi, Gauss–Seidel, SOR) the cost per step
is ∼ nh .

Number of Steps Total Cost

Jacobi ∼ nh ∼ h−2 ∼ nh2
Gauss–Seidel ∼ nh ∼ h−2 ∼ nh2
SOR with ωopt ∼ nh0.5 ∼ h−1 ∼ nh1.5

An optimal iterative method has a total cost ∼ nh .

65
Content

3. Iterative Solvers
3.1 Properties of the Difference Matrix
3.2 Direct Methods and their Disadvantages
3.3 Linear Iterative Solvers
3.4 Practical Examples
3.5 Eigenvalues and Eigenvectors of the Difference Matrix
3.6 Back to Splitting Methods
3.7 Smoothing Property
3.8 Two-Grid Methods

66
Error Propagation of the Jacobi Method, Example 3.1

The Jacobi method has the iteration matrix

Sh = I − h2 /4Lh
with the eigenvectors and eigenvalues
v νµ
ℓ(i,k ) = sin(νπxi ) sin(µπyk ), λνµ = (cos(νπh) + cos(µπh))/2,
where ν, µ = 1, . . . , Nh .
The initial error
Nh X
X Nh
0
e = cνµ v νµ
ν=1 µ=1
in the k -th iteration is
Nh X
X Nh
k k 0
e =S e = cνµ (λνµ )k v νµ .
ν=1 µ=1

67
Observation
Error components with different frequencies (ν, µ) are
reduced with different rates.

Example:
1. ν = µ = 1
λνµ = cos(πh) ≈ 1 − π 2 h2 /2 + O(h4 )
Very slow error reduction for small mesh sizes h!
2. ν = µ = Nh /2 = (1 − h)/2h
λνµ = cos((1 − h)π/2)
= cos(π/2) cos(hπ/2) + sin(π/2) sin(hπ/2)
= sin(hπ/2) = hπ/2 − O(h3 )
Very fast error reduction for small mesh sizes h!

68
Eigenvalues of the Iteration Matrix, Example 3.1

Nh = 63, ν = µ, λνν = cos(νπh)

0.5

-0.5

-1
0 0.5 1 1.5 2 2.5 3
69
Damped Jacobi Method

Goal: Fast error reduction for high-frequency error

components
Idea: Change the eigenvalues of the iteration matrix in the
high-frequency range of the eigenvectors by appropriate scaling
Jacobi Method:
x kJac
+1
= x k + D −1 r k

Damped Jacobi Method with relaxation parameter 0 < ω ≤ 1

(see idea of the SOR method):

x kJac,ω
+1
= x k + ω(x kJac
+1
− x k ) = x k + ωD −1 r k

with the iteration matrix S = I − ωD −1 A.

70
Damped Jacobi Method, Example 3.1

Application to the difference equations gives the iteration matrix

Sh,ω = I − ωh2 /4Lh

with the eigenvalues

λνµ
h,ω = (1 − ω) + ω(cos(νπh) + cos(µπh))/2.

71
Nh = 63

0.5

-0.5

-1
0 0.5 1 1.5 2 2.5 3

72
Nh = 63

0.5

-0.5

-1
0 0.5 1 1.5 2 2.5 3

73
Nh = 63, ω = 1/2, error after k = 0, 1, 5, 10 iterations

74
Observation
Linear iterative methods (1) smooth high-frequency errors and
√
(2) require (at least) O(h−1 ) = O( nh ) steps on a grid with
mesh size h to “transport” information from A to B (one-step
method)!

A h
75
Idea
Perform the iterations on a coarser grid with mesh size
H = 2h. Low frequency error components with respect to h
become high frequency error components with respect to H.

A H

76
Content

3. Iterative Solvers
3.1 Properties of the Difference Matrix
3.2 Direct Methods and their Disadvantages
3.3 Linear Iterative Solvers
3.4 Practical Examples
3.5 Eigenvalues and Eigenvectors of the Difference Matrix
3.6 Back to Splitting Methods
3.7 Smoothing Property
3.8 Two-Grid Methods

77
Two-Grid Methods for Lh u h = f h

Initial vector u 0h on a grid with mesh size h.

Perform k ∈ N steps of a smoothing splitting method (e.g.,
damped Jacobi method with ω = 1/2) ⇒ iterate u kh , residual r kh ,
error e kh (low frequency, smooth).
How can we use this information?
Idea: The error e kh = u h − u kh fulfills the system of linear
equations
Lh e kh = r kh . (3.4)

Test: Lh e kh = Lh (u h − u kh ) = f h − Lh u kh = r kh
Problem: The solution of (3.4) is too expensive!

78
Idea: Define grid transfer matrices
R : Ωh → ΩH (fine to coarse, restriction)
P : ΩH → Ωh (coarse to fine, prolongation)
with the properties that
PR ≈ I
for low-frequency grid functions on Ωh , e.g. e kh . Then it holds
RLh e kh = Rr kh
RLh PRe kh ≈ Rr kh
(RLh P)(Re kh ) ≈ Rr kh
eH c H = r k
L (3.5)
H

(3.5) is an linear system on the coarse grid with fewer

unknowns and can therefore (hopefully) be solved with a
smaller cost than the error equation in (3.4).
79
Two-Grid Method

Implementation
Starting from u 0h , f h , Lh , R, P and L
eH = RLh P:

1. Carry out k1 ∈ N steps of a smoothing splitting method.

Partial result: u kh1 , r kh1
2. Grid transfer from fine to coarse: r kH1 = Rr kh1
3. Solve the correction equation: LeH c H = r k1
H
4. Grid transfer from coarse to fine: c h = Pc H
5. Update u kh1 +1 = u kh1 + c h
6. Carry out k2 ∈ N steps of a smoothing splitting method,
starting from u kh1 +1 .
Result: u kh1 +k2 +1 , r kh1 +k2 +1

80
Error Propagation
Initial error: e 0h
After k1 smoothing steps: e kh1 = S k1 e 0h
Coarse Grid Correction (Steps 2.–5., p. 80):
u kh1 +1 = u kh1 + c h = u kh1 + Pc H = u kh1 + P L
e−1 r k1
H H

= u kh1 + P L
e−1 Rr k1
H h

= u kh1 + P L
e−1 RLh e k1
H h

e kh1 +1 = u h − u kh1 +1 = u h − u kh1 − P L

e−1 RLh e k1
H h
e−1 RLh )e k1
= (I − P LH h

After k2 smoothing steps: e kh1 +k2 +1 = S k2 e kh1 +1

In summary:
e kh1 +k2 +1 = S k2 (I − P L
e−1 RLh )S k1 e 0 = S k2 SS
H h
e k1 e 0
h

81
Grid Transfer

Coarse to fine: Bilinear interpolation, matrix P ∈ Rnh ×nH

Example: nH = 4, nh = 9
 
1 0 0 0
1/2 1/2 0 0 
 
 
 0 1 0 0 
 
1/2 0 1/2 0 
 
P = 1/4 1/4 1/4 1/4
 
 
 0 1/2 0 1/2
 
 0 0 1 0 
 
 
 0 0 1/2 1/2
0 0 0 1

82
Grid Transfer

Fine to coarse: Matrix R ∈ RnH ×nh , R = (1/4)P ⊤

Example: nH = 4, nh = 9

 
1/4 1/8 0 1/8 1/16 0 0 0 0
 0 1/8 1/4 0 1/16 1/8 0 0 0 
R=
 

 0 0 0 1/8 1/16 0 1/4 1/8 0 
0 0 0 0 1/16 1/8 0 1/8 1/4

Rule of Thumb
Chose R = cP ⊤ with a constant c ∈ R, c > 0.
Why?
eH = RLh P = cP ⊤ Lh P
L

is then symmetric and positive-definite. 83

Coarse Grid Matrix

0 0

100
500
200
1000
300

1500 400

2000 500

600
2500
700
3000
800

3500 900

1000
0 500 1000 1500 2000 2500 3000 3500 0 200 400 600 800 1000
nz = 19593 nz = 8836

Matrix Lh (Nh = 63) Matrix L

eH = RLh P

84
Matlab Demonstration Example 3.1

f (x, y ) = 2π 2 sin(πx) sin(πy ), 0 < x, y < 1, x 0 = 0

Number of iterations of the two-grid method to reduce the

norm of the initial residual by a factor 106 . k1 = k2 = k
smoothing steps each with Gauss–Seidel method before and
after the coarse grid correction.

h k =1 k =2 k =5
2−2 6 4 3
2−3 6 5 4
2−4 6 5 4
2−5 6 4 4
2−6 5 4 4
2−7 5 4 3
85
Conclusion: The number of iterations of the two-grid method
with mesh sizes (h, 2h) is independent of h! This situation is
optimal and significantly better than using the splitting methods
alone.
Problem: Solving the correction equation (step 3., p. 80) with
direct methods is still too expensive, especially for small mesh
sizes 2h.
Idea: In step 3, p. 80 we use again a two-grid method with
coarser mesh sizes.
Multigrid Methods, e-g. V-Cycle:
h → 2h → 4h → 8h → · → H → · → 8h → 4h → 2h → h

Total cost of V-cylce in 2D (geometric series):

cnd + cnd /4 + cnd /16 + cnd /64 + · · · + ≤ 4/3cnd ∼ nd
86
[Link]

87
Look closely!

88
Image Sources

• Wikimedia Commons: p. 88
• Society for Industrial and Applied Mathematics: p. 87

89