<< . .

. 88
( : 95)



. . >>

j
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


(¯, t)
¯
t

(¯, 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

1
u0 (x) e’ikx dx
±k =

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

|γ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

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

<< . .

. 88
( : 95)



. . >>