j=’∞

Note that · ∆,p is an approximation of the norm of Lp (R). For instance, the

implicit Bacward Euler/centred scheme (13.43) is unconditionally stable

with respect to the norm · ∆,2 (see Exercise 7).

13.8.3 The CFL Condition

Courant, Friedrichs and Lewy [CFL28] have shown that a necessary and

su¬cient condition for any explicit scheme of the form (13.36) to be stable

is that the time and space discretization steps must obey the following

condition

∆t

|a»| = a ¤1 (13.48)

∆x

which is known as the CFL condition. The number a», which is an adi-

mensional number since a is a velocity, is commonly referred to as the CFL

number. If a is not constant the CFL condition becomes

∆x

∆t ¤ ,

sup |a(x, t)|

x∈R, t>0

while, in the case of the hyperbolic system (13.30), the stability condition

becomes

∆t

¤ 1 k = 1, . . . , p,

»k

∆x

where {»k : k = 1 . . . , p} are the eigenvalues of A.

The CFL stability condition has the following geometric interpretation.

In a ¬nite di¬erence scheme the value un+1 depends, in general, on the

j

values of un at the three points xj+i , i = ’1, 0, 1. Thus, at the time t = 0

the solution un+1 will depend only on the initial data at the points xj+i ,

j

for i = ’(n + 1), . . . , (n + 1) (see Figure 13.5).

Let us de¬ne numerical domain of dependence D∆t (xj , tn ) to be the set

of values at time t = 0 the numerical solution un depends on, that is

j

tn

D∆t (xj , tn ) ‚ x ∈ R : |x ’ xj | ¤ n∆x = .

»

13.8 Analysis of Finite Di¬erence Methods 607

t

tn+1

tn

t1

x

t0

xj+(n+1)

xj’(n+1) xj’1 xj xj+1

FIGURE 13.5. The numerical domain of dependence D∆t (xj , tn+1 )

Consequently, for any ¬xed point (x, t) we have

t

D∆t (x, t) ‚ x ∈ R : |x ’ x| ¤ .

»

In particular, taking the limit as ∆t ’ 0 for a ¬xed », the numerical domain

of dependence becomes

t

x ∈ R : |x ’ x| ¤

D0 (x, t) = .

»

The condition (13.48) is thus equivalent to the inclusion

D(x, t) ‚ D0 (x, t), (13.49)

where D(x, t) is the domain of dependence de¬ned in (13.31).

In the case of an hyperbolic system, thanks to (13.49), we can conclude

that the CFL condition requires that any straight line x = x ’ »k (t ’ t),

k = 1, . . . , p, must intersect the temporal straight line t = t ’ ∆t at some

point x belonging to the domain of dependence (see Figure 13.6).

Let us analyze the stability properties of some of the methods introduced

in the previous section.

Assuming that a > 0, the upwind scheme (13.41) can be reformulated as

un+1 = un ’ »a(un ’ un ). (13.50)

j j j’1

j

Therefore

¤ ∆x |(1 ’ »a)un | + ∆x |»aun |.

un+1 ∆,1 j j’1

j j

608 13. Parabolic and Hyperbolic Initial Boundary Value Problems

x¯

(¯, t)

¯

t

x¯

(¯, t)

¯

t r1 r2

r1 r2

t ’ ∆t

¯

x ’ ∆x x ’ ∆x

¯ x

¯ x + ∆x

¯ ¯ x

¯ x + ∆x

¯

FIGURE 13.6. Geometric interpretation of the CFL condition for a system with

p = 2, where ri = x ’ »i (t ’ t) i = 1, 2. The CFL condition is satis¬ed for the

¯

¯

left-hand case, while it is violated for the right-hand case

Both »a and 1 ’ »a are nonnegative if (13.48) holds. Thus

¤ ∆x(1 ’ »a) |un | + ∆x»a |un | = un

un+1 ∆,1 .

∆,1 j j’1

j j

Inequality (13.46) is therefore satis¬ed by taking CT = 1 and · =· ∆,1 .

∆

The Lax-Friedrichs scheme is also stable, upon assuming (13.48). Indeed,

from (13.39) we get

1 1

(1 ’ »a)un + (1 + »a)un .

un+1 = j+1 j’1

j

2 2

Therefore,

®

1

∆x ° (1 + »a)un »

¤ (1 ’ »a)un +

un+1 ∆,1 j+1 j’1

2 j j

1 1

¤ (1 ’ »a) un n

= un

∆,1 + (1 + »a) u ∆,1 .

∆,1

2 2

Also the Lax-Wendro¬ scheme is stable under the usual assumption (13.48)

on ∆t (for the proof see, e.g., [QV94] Chapter 14).

13.8.4 Von Neumann Stability Analysis

Let us now show that the condition (13.48) is not su¬cient to ensure that

the forward Euler/centred scheme (13.37) is stable. For this purpose, we

make the assumption that the function u0 (x) is 2π-periodic so that it can

be expanded in a Fourier series as

∞

±k eikx

u0 (x) = (13.51)

k=’∞

13.8 Analysis of Finite Di¬erence Methods 609

where

2π

1

u0 (x) e’ikx dx

±k =

2π

0

is the k-th Fourier coe¬cient of u0 (see Section 10.9). Therefore,

∞

j = 0, ±1, ±2, · · ·

u0 ±k eikjh

= u0 (xj ) =

j

k=’∞

where we have set h = ∆x for ease of notation. Applying (13.37) with n = 0

we get

∞

a∆t ikh

(e ’ e’ikh )

±k eikjh 1 ’

u1 =

j

2h

k=’∞

∞

a∆t

±k eikjh 1 ’

= i sin(kh) .

h

k=’∞

Setting

a∆t

γk = 1 ’

i sin(kh),

h

and proceeding recursively on n yields

∞

j = 0, ±1, ±2, . . . , n ≥ 1.

±k eikjh γk

n

un = (13.52)

j

k=’∞

The number γk ∈ C is said to be the ampli¬cation coe¬cient of the k-th

frequency (or harmonic) at each time step. Since

1

2

2

a∆t

|γk | = 1+ sin(kh) ,

h

we deduce that

mπ

|γk | > 1 if a = 0 and k = m = 0, ±1, ±2, . . .

,

h

Correspondingly, the nodal values |un | continue to grow as n ’ ∞ and the

j

numerical solution ”blows-up” whereas the exact solution satis¬es

|u(x, t)| = |u0 (x ’ at)| ¤ max|u0 (s)| ∀x ∈ R, ∀t > 0.

s∈R

The centred discretization scheme (13.37) is thus unconditionally unstable,

i.e., it is unstable for any choice of the parameters ∆t and ∆x.

610 13. Parabolic and Hyperbolic Initial Boundary Value Problems

The previous analysis is based on the Fourier series expansion and is

called von Neumann analysis. It can be applied to studying the stability of

any numerical scheme with respect to the norm · ∆,2 and for establishing

the dissipation and dispersion of the method.

Any explicit ¬nite di¬erence numerical scheme for problem (13.26) satis-

¬es a recursive relation analogous to (13.52), where γk depends a priori on

∆t and h and is called the k-th ampli¬cation coe¬cient of the numerical

scheme at hand.

Theorem 13.1 Assume that for a suitable choice of ∆t and h, |γk | ¤ 1

∀k; then, the numerical scheme is stable with respect to the · ∆,2 norm.

Proof. Take an initial datum with a ¬nite Fourier expansion

N ’1

2

±k eikx

u0 (x) =

k=’ N

2

where N is a positive integer. Without loss of generality we can assume that

problem (13.26) is well-posed on [0, 2π] since u0 is a 2π-periodic function. Take

in this interval N equally spaced nodes

2π

j = 0, . . . , N ’ 1,

xj = jh with h= ,

N

at which the numerical scheme (13.36) is applied. We get

N N

’1 ’1

2 2

u0 ikjh

un ±k γk eikjh .

n

= u0 (xj ) = ±k e , =

j j

k=’ N k=’ N

2 2

Notice that

N ’1

N ’1 2

un 2 ±k ±m (γk γ m )n ei(k’m)jh .

=h

∆,2

j=0 k,m=’ N

2

Recalling Lemma 10.1 we have

N ’1

N N

’ ¤ k, m ¤ ’ 1,

ei(k’m)jh = 2πδkm ,

h

2 2

j=0

which yields

N ’1

2

|±k |2 |γk |2n .

un 2 = 2π

∆,2

k=’ N

2

As a consequence, since |γk | ¤ 1 ∀k, it turns out that

N ’1

2

¤ 2π ∀n ≥ 0,

2

= u0 2

un 2 ±k ∆,2 ,

∆,2

k=’ N

2

13.9 Dissipation and Dispersion 611

· 3

which proves that the scheme is stable with respect to the norm.

∆,2

In the case of the upwind scheme (13.41), proceeding as was done for

the centred scheme, we ¬nd the following ampli¬cation coe¬cients (see

Exercise 6)

±

1 ’ a ∆t (1 ’ e’ikh ) if a > 0,

h

γk =

1 ’ a ∆t (e’ikh ’ 1) if a < 0.

h

Therefore

h

∀k, |γk | ¤ 1 if ∆t ¤ ,

|a|

which is nothing but the CFL condition.

Thanks to Theorem 13.1, if the CFL condition is satis¬ed, the upwind

scheme is stable with respect to the · ∆,2 norm.

We conclude by noting that the upwind scheme (13.50) satis¬es

un+1 = (1 ’ »a)un + »aun .

j j’1

j

Owing to (13.48), either »a or 1 ’ »a are nonnegative, thus

min(un , un ) ¤ un+1 ¤ max(un , un ).

j j’1 j j’1

j

It follows that

inf u0 ¤ un ¤ sup u0 ∀j ∈ Z, ∀n ≥ 0,

l j l

l∈Z l∈Z

that is,

¤ u0 ∀n ≥ 0,

un (13.53)

∆,∞ ∆,∞

which proves that if (13.48) is satis¬ed, the upwind scheme is stable in the

norm · ∆,∞ . The relation (13.53) is called the discrete maximum principle

(see also Section 12.2.2).

Remark 13.4 For the approximation of the wave equation (13.33) the

Leap-Frog method (13.44) is stable under the CFL restriction ∆t ¤ ∆x/|γ|,

while the Newmark method (13.45) is unconditionally stable if 2β ≥ θ ≥ 1

2

(see [Joh90]).

13.9 Dissipation and Dispersion

The von Neumann analysis on the ampli¬cation coe¬cients enlightens the

study of the stability and dissipation of a numerical scheme.

612 13. Parabolic and Hyperbolic Initial Boundary Value Problems

Consider the exact solution to problem (13.26); the following relation

holds

u(x, tn ) = u0 (x ’ an∆t), ∀n ≥ 0, ∀x ∈ R.

In particular, from applying (13.51) it follows that

∞

gk = e’iak∆t .

n

±k eikjh gk ,

n