Numerical Methods for Differential Equations

Chapter 3: FDM for 2p-BVPs and Sturm–Liouville

Tony Stillfjord, Gustaf Söderlind

Numerical Analysis, Lund University
Contents V4.20

1. Finite difference approximation of derivatives

2. Finite difference methods for 2p-BVPs Lu = f

3. Newton’s method

4. Boundary conditions

5. Adaptive grids

6. Sturm–Liouville eigenvalue problems Lu = λu

7. Toeplitz matrices

8. Convergence: The Lax Principle

2 / 67
1. Approximation of derivatives y 0 = dy /dx

First order approximations

Forward difference
y (x + ∆x) − y (x)
y 0 (x) = + O(∆x)
∆x

Backward difference
y (x) − y (x − ∆x)
y 0 (x) = + O(∆x)
∆x

3 / 67
Spatial symmetric approximation of derivatives

Second order approximations

Symmetric difference quotients

y (x + ∆x) − y (x − ∆x)
y 0 (x) = + O(∆x 2 )
2∆x

y (x + ∆x) − 2y (x) + y (x − ∆x)

y 00 (x) = + O(∆x 2 )
∆x 2

4 / 67
Derivatives → finite differences → matrices

Matrix representation of forward difference

y (x + ∆x) − y (x)
y 0 (x) = + O(∆x)
∆x

Introduce vectors y = {y (xi )} and y 0 = {y 0 (xi )}

y00
    
−1 1 y0
 y10  1  −1 1  y1 
≈
    
.. .. .. ..
 ∆x 
    
 . . .  . 
yN0 −1 1 yN+1

5 / 67
From derivatives to matrices

Note Forward difference ∼ (N + 1) × (N + 2) matrix

y00
    
−1 1 y0
 y10  1  −1 1  y1 
≈
    
.. .. .. ..
 ∆x 
   
 . . .  . 
yN0 −1 1 yN+1

Nullspace spanned by y = (1 1 1 . . . 1)T

Compare nullspace of d/dx , y = 1 ⇒ y 0 ≡ 0

Analogous result for backward difference

6 / 67
From derivatives to matrices. . .

Central difference
y (x + ∆x) − y (x − ∆x)
y 0 (x) ≈
2∆x

Matrix representation
 0   
y1   y0
 y0  −1 0 1  y1
 2  1  . .

 ..  ≈ .. .. ... 
 ..

 .  2∆x
 
 . 
−1 0 1
yN0 yN+1

7 / 67
From derivatives to matrices. . .

Note N × (N + 2) matrix

y10
   
  y0
−1 0 1
 y20  1  . .. ..
 y1 

.. ≈
 .. . .

 ..

 2∆x
  
 .  . 
−1 0 1
yN0 yN+1

Nullspace is now two-dimensional

ȳ = (1 1 1 . . . 1)T and ỹ = (1 −1 1 −1 . . . 1)T

8 / 67
From derivatives to matrices. . .

“False” nullspace

ỹ = (1 −1 1 −1 . . . 1)T does not converge to a C 1 function!

Compare difference equation yn+1 − yn−1 = 0, with characteristic

equation

z2 − 1 = 0 ⇒ z = ±1

and two solutions ȳn = 1 and ỹn = (−1)n

9 / 67
2nd order derivatives → matrices y 00 = d2 y /dx 2

Central difference
y (x + ∆x) − 2y (x) + y (x − ∆x)
y 00 (x) ≈
∆x 2

y100
   
  y0
1 −2 1
 y200  1  . . ..
 y1 

.. ≈
 . . . . .

 ..

 ∆x 2
  
 .  . 
1 −2 1
yN00 yN+1

Note N × (N + 2) matrix with 2D nullspace spanned by

ȳ = (1 1 . . . 1)T and ŷ = (0 1 2 3 . . . N + 1)T

10 / 67
2nd order derivatives. . .

Nullspace of d2 /dx 2

y = 1 and y = x both have y 00 ≡ 0

Compare difference equation yn+1 − 2yn + yn−1 = 0, with

characteristic equation

z 2 − 2z + 1 = 0 ⇒ z = 1, 1

and two solutions ȳn = 1 and ŷn = n , respectively

This corresponds directly to y = 1 and y = x

11 / 67
2. Finite difference methods for 2p-BVP

Consider simplest problem

y 00 = f (x, y )
y (0) = α; y (1) = β

Introduce equidistant grid with ∆x = 1/(N + 1)

FDM discretization

yi+1 − 2yi + yi−1

= f (xi , yi ) i =1:N
∆x 2

y0 = α; yN+1 = β

12 / 67
Discrete 2pBVP Equation system F (y ) = 0

α − 2y1 + y2
F1 (y ) = − f (x1 , y1 )
∆x 2

yi−1 − 2yi + yi+1

Fi (y ) = − f (xi , yi )
∆x 2

yN−1 − 2yN + β
FN (y ) = − f (xN , yN )
∆x 2

A (nonlinear) system F (y ) = 0 for N unknowns y1 , y2 , . . . , yN

Note how boundary values enter

13 / 67
3. Newton’s method F (y ) = 0

Let y (k) approximate the solution y and expand in Taylor series

0 = F (y ) = F (y (k) + y − y (k) ) ≈ F (y (k) ) + F 0 (y (k) ) · (y − y (k) )

Define y (k+1) by 0 =: F (y (k) ) + F 0 (y (k) ) · (y (k+1) − y (k) )

Newton’s method (mathematical formulation)

y (k+1) := y (k) − [F 0 (y (k) )]−1 F (y (k) )

14 / 67
Jacobian matrix F 0 (y ) = {∂Fi /∂yj }

For the FDM the 2p-BVP Jacobian matrix is

∂f
F 0 (y ) = tridiag (1/∆x 2 , −2/∆x 2 − , 1/∆x 2 )
∂yi

Tridiagonal matrix, with

• Super- and subdiagonal elements 1/∆x 2

• Diagonal elements −2/∆x 2 − ∂f /∂yi

• Sparse LU decomposition runs in O(N) time

• Solution effort moderate even when N is large

15 / 67
Newton’s method for F (y ) = 0

Newton iteration

1. Compute Jacobian F 0 (y (k) ) = {∂Fi /∂yj }

2. Factorize Jacobian matrix F 0 (y (k) ) → LU

3. Solve linear system LUδy (k) = −F (y (k) )

4. Update y (k+1) := y (k) + δy (k)

Newton’s method is quadratically convergent

16 / 67
Quadratic convergence

Newton’s method converges if

• kF 0 (y (k) )−1 k ≤ C 0

• kF 00 (y (k) )k ≤ C 00

• ky (0) − y k < ε (close enough starting value)

Then convergence is quadratic

ky (k+1) − y k ≤ C · ky (k) − y k2

17 / 67
4. Boundary conditions come in many types

In many cases the problem is linear, but boundary conditions vary

• Dirichlet conditions
y (0) = α ; y (1) = β straightforward to implement

• Neumann conditions
y 0 (0) = γ ; y (1) = β requires special attention

• Robin conditions
y (0) + c · y 0 (0) = κ ; y (1) = β requires same attention

for the method’s convergence order to be preserved

18 / 67
Neumann problem

Example y 00 = f (x, y )
y (0) = α; y 0 (1) = β

For second-order convergence, we must approximate also y 0 (1) to

second order

Standard symmetric approximation

y (x + ∆x) − y (x − ∆x)
2∆x
requires points to the right of x = 1

What grid to use?

19 / 67
Neumann problem, approach 1

Example y 00 = f (x, y )
y (0) = α; y 0 (1) = β

1
Equidistant grid, xn = n∆x with ∆x = N and xN+1 = 1 + ∆x
yN+1 − yN−1
y 0 (1) = β → =β
2∆x

⇒ yN+1 := 2β ∆x + yN−1 is of second order at x = 1

20 / 67
Neumann problem, approach 2

Example y 00 = f (x, y )
y (0) = α; y 0 (1) = β

1
Equidistant grid, xn = n∆x with ∆x = N+1/2 and

xN + ∆x/2 = 1 = xN+1 − ∆x/2

yN+1 − yN
y 0 (1) = β → =β
∆x

⇒ yN+1 := β ∆x + yN is of second order at x = 1

yN +yN+1
Need to approximate y (1) ≈ 2

21 / 67
Neumann problem, approach 3

Example y 00 = f (x, y )
y (0) = α; y 0 (1) = β

1
Equidistant grid, xn = n∆x with ∆x = N+1 and xN+1 = 1
yN−1 − 4yN + 3yN+1
y 0 (1) = β → =β
2∆x

1
⇒ yN+1 := 3 (2β∆x + 4yN − yN−1 ) is of second order at x = 1

All three approaches work! It’s a matter of taste.

22 / 67
Robin problem

Example y 00 = f (x, y )
y (0) = α; y (1) + cy 0 (1) = κ

Equidistant grid, with x = 1 between grid points

xN + ∆x/2 = 1 = xN+1 − ∆x/2

yN+1 + yN yN+1 − yN
y (1) + cy 0 (1) = κ → +c =κ
2 ∆x

(2c − ∆x)yN + 2κ∆x

⇒ yN+1 :=
2c + ∆x
Other two approaches analogous
23 / 67
5. FDM on adaptive grids

Left and right divided differences

yi − yi−1 yi − yi−1 yi+1 − yi yi+1 − yi

D − yi = = D + yi = =
xi − xi−1 h− xi+1 − xi h+

Then approximate derivatives by finite differences

h− D + yi + h+ D − yi
yi0 ≈
h+ + h−
D + yi − D − yi
yi00 ≈ 2
h+ + h−

This is 2nd order only on smooth grids with h+ /h− = 1 + O(N −1 )

24 / 67
Nonuniform grids Grid deformation

A smooth nonuniform grid

1

0.9

0.8

0.7

0.6

0.5

0.4

0.3

0.2

0.1

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

25 / 67
Adaptive grids u 00 + 100u 0 = 100

Solution to convection-diffusion boundary layer problem

0

-0.2

-0.4

-0.6
u

-0.8

-1

-1.2
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x

26 / 67
6. Sturm–Liouville eigenvalue problems

Diffusion problem
 
∂u ∂ ∂u
= p(x) ; u(t, a) = u(t, b) = 0
∂t ∂x ∂x

Separation of variables (one space dimension)

v̇ (p(x) y 0 )0
u(t, x) := y (x) · v (t) ⇒ = =: λ
v y

Sturm–Liouville eigenvalue problem

 
d dy
p(x) = λy y (a) = 0 , y (b) = 0
dx dx

27 / 67
Sturm–Liouville eigenvalue problems. . .

Wave equation

∂2u 2
2∂ u
= c ; u(t, a) = u(t, b) = 0
∂t 2 ∂x 2

Express solution as u(t, x) = y (x) eiωt ⇒

−ω 2 y = c 2 y 00 y (a) = y (b) = 0

Sturm–Liouville eigenvalue problem

y 00 = λy with λ = −ω 2 /c 2

28 / 67
Why Sturm–Liouville eigenvalue problems?

Öresund bridge
29 / 67
Fluid–structure interaction

Tacoma Narrows Bridge 1940

30 / 67
Tuned mass dampers

Stockbridge damper (1926) – anti-fatigue devices

31 / 67
Taipei 101

Somewhat larger tuned mass damper (660 000 kg)

32 / 67
It rocks

Gustaf Söderlind: Mathematics rocks too!

33 / 67
Music instruments

34 / 67
Compression loads, buckling and sun kinks

35 / 67
Stationary Schrödinger equation Particle in a box

Quantum mechanics

~2 d2 ψ
− + V (x)ψ = E ψ
2m dx 2

This is a Sturm–Liouville eigenvalue problem, with energy levels Ek

defined by the eigenvalues

36 / 67
Sturm–Liouville eigenvalue problem

Find eigenvalues λ and eigenfunctions y (x) with

 
d dy
p(x) + q(x)y = λy ; y (a) = y (b) = 0
dx dx

Discretization Matrix eigenvalue problem

T∆x y = λ∆x y

Note Analytic eigenvalue problem converts to algebraic!

37 / 67
Sturm–Liouville problem Discretization

Symmetric discretizations

p(x)y 0 (x) − p(x)y 0 (x)

 


d dy xi + ∆x
2
xi − ∆x
2
p(x) ≈
dx dx x=xi ∆x
 
∆x 0 ∆x y (xi+1 ) − y (xi )
p(xi + ) y (x) xi + ∆x
≈ p xi +
2 2 2 ∆x
∆x
(similar for xi − 2 ) lead to

pi−1/2 yi−1 − (pi−1/2 + pi+1/2 )yi + pi+1/2 yi+1

+ q(xi )yi = λyi
∆x 2
y0 = yN+1 = 0

38 / 67
Sturm–Liouville problem Discretization II

pi−1/2 yi−1 − (pi−1/2 + pi+1/2 )yi + pi+1/2 yi+1

+ q(xi )yi = λ∆x yi
∆x 2

y0 = yN+1 = 0

Symmetric tridiagonal N × N eigenvalue problem

T∆x y = λ∆x y

There are N eigenvalues λ∆x,n = λn + O(∆x 2 )

39 / 67
Sturm–Liouville problem A simple example

Consider y 00 = λy with boundary conditions y (0) = y (1) = 0

Analytic solution
√ √
y (x) = A sin −λ x + B cos −λ x
√
Boundary values ⇒ B = 0 and A sin −λ = 0

Eigenvalues and eigenfunctions for k = 1, 2, . . .

λk = −(kπ)2
yk (x) = sin kπx

Fourier modes (harmonic analysis) associated with d2 /dx 2

40 / 67
Discrete Sturm–Liouville problem Benchmark test

Discretization of y 00 = λy with BVs ⇒

yi−1 − 2yi + yi+1
= λ∆x yi
∆x 2
y0 = yN+1 = 0 ; ∆x = 1/(N + 1)

Tridiagonal N × N matrix formulation

    
−2 1 y1 y1
1   1 −2 1 y
 2 
  
 y2
 
= λ

. .. . ∆x  ..
∆x 2   .. 
    
 . 
1 −2 yN yN

41 / 67
Discrete Sturm–Liouville problem. . .

Algebraic eigenvalue problem

T∆x y = λ∆x y

Smallest eigenvalue λ∆x = −π 2 + O(∆x 2 )

The first few eigenvalues are well approximated, but the

approximation gradually gets worse

Note There are only N discrete eigenvalues

42 / 67
Discrete Sturm–Liouville problem Computation

First three eigenvectors of T∆x at N = 19

0.4

0.3

0.2

0.1

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0.4

0.2

−0.2

−0.4
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0.4

0.2

−0.2

−0.4
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

43 / 67
Discrete Sturm–Liouville problem Eigenvalues N = 19

Discrete eigenvalues λ∆x −9.8493 −39.1548 −87.1948

Exact eigenvalues λ −9.8696 −39.4784 −88.8264
Relative errors 0.21% 0.82% 1.84%

Note
• Lowest eigenvalues are more accurate
√
• Good approximations for N first eigenvalues

(Here approximately first 4 – 5 modes)

44 / 67
Discrete Sturm–Liouville problem High modes

Eigenvectors 7, 13, 19 of T∆x at N = 19

0.4

0.2

−0.2

−0.4
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0.4

0.2

−0.2

−0.4
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0.4

0.2

−0.2

−0.4
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

45 / 67
7. Toeplitz matrices

A Toeplitz matrix is constant along diagonals

Example (symmetric)
 
−2 1 0 ...
 1 −2 1 
1  1 −2 1

T∆x =
 
∆x 2 
 
 . .. 

... 0 1 −2

46 / 67
Toeplitz matrices. . .

Much is known about Toeplitz matrices

• Eigenvalues
• Norms
• Inverses
• etc.
They can be generated in Matlab using the built-in function
toeplitz

47 / 67
Eigenvalues of Toeplitz matrices

Example Solve the eigenvalue problem Ty = λy for

 
−2 1 0 ...
 1 −2 1 
 
T =
 1 −2 1 


 . .. 

... 0 1 −2

Note λ[T ] = −2 + λ[S]

48 / 67
Eigenvalues. . .

. . . the problem gets simplified

  
0 1 0 ... y1
 1 0 1  y2 
  
Sy =  1 0 1 ..
 = λy
  
 .
 ..  
 . 1  
... 0 1 0 yN

Find eigenvalues λ[S], noting that the nth equation of Sy = λy is

yn+1 + yn−1 = λyn

49 / 67
Eigenvalues and difference equations

Linear difference equation yn+1 + yn−1 = λyn with boundary values

y0 = 0 = yN+1

Characteristic equation z 2 − λz + 1 = 0

Two roots z and 1/z (product 1) implies general solution

yn = Az n + Bz −n

Boundary condition y0 = 0 = A + B ⇒ yn = A(z n − z −n )

50 / 67
Eigenvalues and difference equations. . .

Boundary condition yN+1 = 0 = A(z N+1 − z −(N+1) ) ⇒

kπi

z 2(N+1) = 1 ⇒ zk = exp k =1:N
N +1

Sum of the roots of z 2 − λz + 1 = 0 are

kπ
λk [S] = zk + 1/zk = 2 cos
N +1

Hence
kπ kπ
λk [T ] = −2 + 2 cos = −4 sin2
N +1 2(N + 1)

51 / 67
Eigenvalue locations λk [T ]

52 / 67
Eigenvalues of Toeplitz matrices

Theorem The N × N Toeplitz matrix

 
−2 1 0 ...
 1 −2 1 
 
T =
 1 −2 1 

 .. 
 . 
... 0 1 −2

has N real eigenvalues (k = 1 : N)

kπ
λk [T ] = −4 sin2 ∈ (−4, 0)
2(N + 1)

53 / 67
Eigenvalues of Toeplitz matrices. . .

Consider T∆x := T /∆x 2 with ∆x = 1/(N + 1) as an operator

approximation
d2
↔ T∆x
dx 2
on x ∈ [0, 1]

Corollary The eigenvalues of T∆x are

kπ
λk [T∆x ] = −4(N + 1)2 sin2 ≈ −k 2 π 2 = λk [d2 /dx 2 ]
2(N + 1)

for k  N

54 / 67
What are the norms of T ?

Lemma For a symmetric matrix A, it holds

kAk2 = max |λk |

k

Lemma For a symmetric matrix A, it holds

µ2 [A] = max λk
k

(Both results actually hold for normal matrices)

55 / 67
Proofs. Norm

Definition
x T AT Ax
kAk22 = max
x T x6=0 x Tx
Find stationary points of the Rayleigh quotient of AT A, given by
ρ(x) = x T AT Ax/x T x

d
ρ(x) = (2x T AT Ax T x − 2x T x T AT Ax)/(x T x)2 := 0
dx

⇒ x T AT A = x T ρ(x) ⇒ A2 x = ρ(x)x

So ρ(x) = λ2 , therefore kAk2 = max |λ[A]| if AT = A

56 / 67
Proofs. Logarithmic norm

Definition (real A)

x T Ax
µ2 [A] = max
x T x6=0 x T x

Find stationary points of the Rayleigh quotient of A, given by

ρ(x) = x T Ax/x T x

d
ρ(x) = [x T (A + AT )x T x − 2x T x T Ax]/(x T x)2 := 0
dx
1 T
x (A + AT )x = x T ρ(x) ⇒ Ax = ρ(x)x
2

So ρ(x) = λ, therefore µ2 [A] = max λ[A]

57 / 67
What are the norms of T∆x ?

Eigenvalues of T∆x = T /∆x 2 are

kπ
λk [T∆x ] = −4(N + 1)2 sin2
2(N + 1)

So kT∆x k2 = |λN | and µ2 [T∆x ] = λ1

Theorem The Euclidean norms of T∆x are

4
kT∆x k2 ≈ µ2 [T∆x ] ≈ −π 2
∆x 2

58 / 67
 −1
The norm of T∆x

Recall that µ[A] < 0 ⇒ kA−1 k ≤ −1/µ[A]

Approximate y 00 = f (x) with y (0) = y (1) = 0 by

T∆x u = q

Note µ2 [T∆x ] ≈ −π 2 implies the existence of a unique solution,

as
−1 1
kT∆x k2 /
π2

59 / 67
What norms to use Euclidean and RMS norms

The norm of a function is measured in the L2 norm

Z 1
2
kukL2 = u(x)2 dx
0

A corresponding discrete function (vector) is then measured in the

root mean square (RMS) norm
N N
X 1 X 1
kuk2∆x = 2
u(xi ) ∆x = u(xi )2 = kuk22
N +1 N +1
i=1 i=1

−1 −1
Note For the operator norm, kT∆x k∆x ≡ kT∆x k2

60 / 67
8. Convergence of finite difference methods

Simplest model problem (1D Poisson equation)

y 00 = f (x)
y (0) = α; y (1) = β

Equidistant discretization
yi+1 − 2yi + yi−1
= f (xi )
∆x 2
y0 = α; yN+1 = β

61 / 67
Error definitions Local error

Insert exact continuous solution y (x) into discretization

y (xi−1 ) − 2y (xi ) + y (xi+1 )

= f (xi , y (xi )) − l(xi )
∆x 2

Taylor expansion of local error, using f (xi ) = y 00 (xi )

 2
∆x 4 (6)

∆x (4)
−l(xi ) = 2 y (xi ) + y (xi ) + . . .
4! 6!

Only even powers of ∆x due to symmetry. In particular, we get

1
klk∆x / ∆x 2 max |y (4) (ξ)|
12 ξ∈[0,1]

62 / 67
Error definitions Global error

Definition The global error is defined by e(xi ) = yi − y (xi )

Convergence
Will show that e(x) → 0 as ∆x → 0, or more specifically

kek∆x = c1 ∆x 2 + c2 ∆x 4 + . . .

Again only even powers due to symmetry

63 / 67
Convergence

Consider the problem y 00 = f discretized by 2nd order FDM

T∆x u = f (x)
with T∆x tridiagonal. Then

Numerical solution T∆x u = f (x)

Exact solution T∆x y (x) = f (x) − l(x)

Error equation T∆x e(x) = l(x)

where e(x) = u − y (x) is the global error

64 / 67
Convergence. . .

Solve T∆x u = f formally to get

−1
Numerically u = T∆x · f (x)

−1
Exact y (x) = T∆x · (f (x) − l(x))

−1
Global error e(x) = T∆x · l(x)

−1
Error bound kek∆x ≤ kT∆x k2 · klk∆x

65 / 67
Convergence. . .

Recall
−1
• µ2 [T∆x ] ≈ −π 2 ⇒ kT∆x k2 / 1/π 2

−1
• kek∆x ≤ kT∆x k2 · klk∆x

1
• klk∆x ≤ ∆x 2 12 maxξ∈[0,1] |y (4) (ξ)| + O(∆x 4 )

We therefore have
1 ∆x 2
kek∆x ≤ · klk ∆x = max |y (4) (ξ)| + O(∆x 4 )
π2 12π 2 ξ∈[0,1]

and we have convergence as ∆x → 0

66 / 67
The Lax Principle

Conclusion
Consistency local error l → 0 as ∆x → 0
−1
Stability kT∆x k2 ≤ C as ∆x → 0
Convergence global error e → 0 as ∆x → 0

Theorem (Lax Principle)

Consistency + Stability ⇒ Convergence

“Fundamental theorem of numerical analysis”

67 / 67