x0 + t(1 ’ x0 ) 0 ¤ x0 ¤ 1,

x(t) = x0 + tu0 (x0 ) =

x0 ≥ 1.

x0

Notice that the characteristic lines do not intersect only if t < 1 (see Figure 13.4,

•

right).

13.6 Systems of Linear Hyperbolic Equations

Consider the linear hyperbolic systems of the form

‚u ‚u

= 0 x ∈ R, t > 0, (13.30)

+A

‚t ‚x

where u : R — [0, ∞) ’ Rp and A ∈ Rp—p is a matrix with constant

coe¬cients.

This system is said hyperbolic if A is diagonalizable and has real eigen-

values, that is, if there exists a nonsingular matrix T ∈ Rp—p such that

A = TΛT’1 ,

600 13. Parabolic and Hyperbolic Initial Boundary Value Problems

where Λ = diag(»1 , ..., »p ) is the diagonal matrix of the real eigenvalues of

A, while T = (ω 1 , ω 2 , . . . , ω p ) is the matrix whose column vectors are the

right eigenvectors of A (see Section 1.7). The system is said to be strictly

hyperbolic if it is hyperbolic with distinct eigenvalues. Thus

Aω k = »k ω k , k = 1, . . . , p.

Introducing the characteristic variables w = T’1 u, system (13.30) becomes

‚w ‚w

+Λ = 0.

‚t ‚x

This is a system of p independent scalar equations of the form

‚wk ‚wk

+ »k = 0, k = 1, . . . , p.

‚t ‚x

Proceeding as in Section 13.5 we obtain wk (x, t) = wk (x ’ »k t, 0), and thus

the solution u = Tw of problem (13.30) can be written as

p

wk (x ’ »k t, 0)ω k .

u(x, t) =

k=1

The curve (xk (t), t) in the plane (x, t) that satis¬es xk (t) = »k is the k-th

characteristic curve and wk is constant along it. A strictly hyperbolic sys-

tem enjoys the property that p distinct characteristic curves pass through

any point of the plane (x, t), for any ¬xed x and t. Then u(x, t) depends

only on the initial datum at the points x ’ »k t. For this reason, the set of p

points that form the feet of the characteristics issuing from the point (x, t)

D(t, x) = x ∈ R : x = x ’ »k t , k = 1, ..., p , (13.31)

is called the domain of dependence of the solution u(x, t).

If (13.30) is set on a bounded interval (±, β) instead of on the whole real

line, the in¬‚ow point for each characteristic variable wk is determined by the

sign of »k . Correspondingly, the number of positive eigenvalues determines

the number of boundary conditions that can be assigned at x = ±, whereas

at x = β it is admissible to assign a number of conditions which equals the

number of negative eigenvalues. An example is discussed in Section 13.6.1.

Remark 13.2 (The nonlinear case) Let us consider the following non-

linear system of ¬rst-order equations

‚u ‚

+ g(u) = 0 (13.32)

‚t ‚x

where g = (g1 , . . . , gp )T is called the ¬‚ux function. The system is hyperbolic

if the jacobian matrix A(u) whose elements are aij = ‚gi (u)/‚uj , i, j =

1, . . . , p, is diagonalizable and has p real eigenvalues.

13.6 Systems of Linear Hyperbolic Equations 601

13.6.1 The Wave Equation

Consider the second-order hyperbolic equation

‚2u ‚2u

’ γ2 2 = f x ∈ (±, β), t > 0, (13.33)

‚t2 ‚x

with initial data

‚u

(x, 0) = v0 (x), x ∈ (±, β),

u(x, 0) = u0 (x) and

‚t

and boundary data

u(±, t) = 0 and u(β, t) = 0, t > 0. (13.34)

In this case, u may represent the transverse displacement of an elastic

vibrating string of length β ’±, ¬xed at the endpoints, and γ is a coe¬cient

depending on the speci¬c mass of the string and on its tension. The string

is subject to a vertical force of density f .

The functions u0 (x) and v0 (x) denote respectively the initial displace-

ment and the initial velocity of the string.

The change of variables

‚u ‚u

ω1 = , ω2 = ,

‚x ‚t

transforms (13.33) into the following ¬rst-order system

‚ω ‚ω

= 0 x ∈ (±, β), t > 0, (13.35)

+A

‚t ‚x

where

’1

ω1 0

ω= , A= ,

’γ 2

ω2 0

and the initial conditions are ω1 (x, 0) = u0 (x) and ω2 (x, 0) = v0 (x).

Since the eigenvalues of A are the two distinct real numbers ±γ (rep-

resenting the propagation velocities of the wave) we conclude that system

(13.35) is hyperbolic. Moreover, one boundary condition needs to be pre-

scribed at every end-point, as in (13.34). Notice that, also in this case,

smooth solutions correspond to smooth initial data, while any discontinuity

that is present in the initial data will propagate along the characteristics.

2

‚2u

Remark 13.3 Notice that replacing ‚ 2 by t2 ,

u

by x2 and f by 1, the

‚x2

‚t

wave equation becomes

t 2 ’ γ 2 x2 = 1

602 13. Parabolic and Hyperbolic Initial Boundary Value Problems

which represents an hyperbola in the (x, t) plane. Proceeding analogously

in the case of the heat equation (13.1), we end up with

t ’ νx2 = 1

which represents a parabola in the (x, t) plane. Finally, for the Poisson

2 2

equation (12.90), replacing ‚ u by x2 , ‚ u by y 2 and f by 1, we get

2 ‚y 2

‚x

x2 + y 2 = 1

which represents an ellipse in the (x, y) plane.

Due to the geometric interpretation above, the corresponding di¬erential

operators are classi¬ed as hyperbolic, parabolic and elliptic.

13.7 The Finite Di¬erence Method for Hyperbolic

Equations

Let us discretize the hyperbolic problem (13.26) by space-time ¬nite dif-

ferences. To this aim, the half-plane {(x, t) : ’∞ < x < ∞, t > 0} is dis-

cretized by choosing a spatial grid size ∆x, a temporal step ∆t and the

grid points (xj , tn ) as follows

xj = j∆x, j ∈ Z, tn = n∆t, n ∈ N.

Let us set

» = ∆t/∆x,

and de¬ne xj+1/2 = xj + ∆x/2. We look for discrete solutions un which

j

n

approximate the values u(xj , t ) of the exact solution for any j, n.

Quite often, explicit methods are employed for advancing in time in

hyperbolic initial-value problems, even though they require restrictions on

the value of », unlike what typically happens with implicit methods.

Let us focus our attention on problem (13.26). Any explicit ¬nite-di¬erence

method can be written in the form

un+1 = un ’ »(hn

j+1/2 ’ hj’1/2 ),

n

(13.36)

j

j

j+1/2 = h(uj , uj+1 ) for every j and h(·, ·) is a particular function

where hn n n

that is called the numerical ¬‚ux.

13.7.1 Discretization of the Scalar Equation

We illustrate several instances of explicit methods, and provide the corre-

sponding numerical ¬‚ux.

13.7 The Finite Di¬erence Method for Hyperbolic Equations 603

1. Forward Euler/centred

»

un+1 = un ’ a(un ’ un ) (13.37)

j j+1 j’1

j

2

which can be cast in the form (13.36) by setting

1

hj+1/2 = a(uj+1 + uj ). (13.38)

2

2. Lax-Friedrichs

1n »

(uj+1 + un ) ’ a(un ’ un )

un+1 = (13.39)

j’1 j+1 j’1

j

2 2

which is of the form (13.36) with

1

[a(uj+1 + uj ) ’ »’1 (uj+1 ’ uj )].

hj+1/2 =

2

3. Lax-Wendro¬

»2 2 n

»

’ a(uj+1 ’ uj’1 ) + a (uj+1 ’ 2un + un ) (13.40)

un+1 un n n

= j j j’1

j

2 2

which can be written in the form (13.36) provided that

1

[a(uj+1 + uj ) ’ »a2 (uj+1 ’ uj )].

hj+1/2 =

2

4. Upwind (or forward Euler/uncentred)

» »

un+1 = un ’ a(un ’ un ) + |a|(un ’ 2un + un ) (13.41)

j j+1 j’1 j+1 j j’1

j

2 2

which ¬ts the form (13.36) when the numerical ¬‚ux is de¬ned to be

1

[a(uj+1 + uj ) ’ |a|(uj+1 ’ uj )].

hj+1/2 =

2

The last three methods can be obtained from the forward Euler/centred

method by adding a term proportional to a numerical di¬usion, so that

they can be written in the equivalent form

1 un ’ 2un + un

» j+1 j j’1

’ a(uj+1 ’ uj’1 ) + k

un+1 un n n

= , (13.42)

j

j 2

2 2 (∆x)

where the arti¬cial viscosity k is given for the three cases in Table 13.1.

604 13. Parabolic and Hyperbolic Initial Boundary Value Problems

hdif f

methods k „ (∆t, ∆x)

j+1/2

∆x2

1

’ (uj+1 ’ uj ) O

∆x2

Lax-Friedrichs + ∆t + ∆x

2» ∆t

»a2

’ (uj+1 ’ uj ) O ∆t2 + ∆x2

a2 ∆t2

Lax-Wendro¬

2

|a|

|a|∆x∆t ’ (uj+1 ’ uj ) O(∆t + ∆x)

Upwind

2

TABLE 13.1. Arti¬cial viscosity, arti¬cial ¬‚ux and truncation error for

Lax-Friedrichs, Lax-Wendro¬ and Upwind methods

As a consequence, the numerical ¬‚ux for each scheme can be written equiv-

alently as

hj+1/2 = hF E + hdif f

j+1/2 j+1/2

where hF E is the numerical ¬‚ux of the forward Euler/centred scheme

j+1/2

(which is given in (13.38)) and the arti¬cial di¬usion ¬‚ux hdif f is given

j+1/2

for the three cases in Table 13.1.

An example of an implicit method is the backward Euler/centred scheme

»

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

un+1 + (13.43)

j

j j+1 j’1

2

It can still be written in the form (13.36) provided that hn is replaced by

hn+1 . In the example at hand, the numerical ¬‚ux is the same as for the

Forward Euler/centred method, and so is the arti¬cial viscosity.

Finally, we report the following schemes for the approximation of the

second-order wave equation (13.33):

1. Leap-Frog

un+1 ’ 2un + un’1 = (γ»)2 (un ’ 2un + un ) (13.44)

j j+1 j j’1

j j

2. Newmark

un+1 ’ un = ∆tvj + (γ»)2 βwj + ’ β wj ,

n+1 1

n n

j

j 2

(13.45)

(γ»)2

’ θwj + (1 ’ θ)wj

n+1 n+1 n

n

vj vj =

∆t

with wj = uj+1 ’2uj +uj’1 and where the parameters β and θ satisfy

0 ¤ β ¤ 1 , 0 ¤ θ ¤ 1.

2

13.8 Analysis of Finite Di¬erence Methods 605

13.8 Analysis of Finite Di¬erence Methods

Let us analyze the properties of consistency, stability and convergence, as

well as those of dissipation and dispersion, of the ¬nite di¬erence methods

introduced above.

13.8.1 Consistency

As illustrated in Section 11.3, the local truncation error of a numerical

scheme is the residual that is generated by pretending the exact solution

to satisfy the numerical method itself.

Denoting by u the solution of the exact problem (13.26), in the case of

method (13.37) the local truncation error at (xj , tn ) is de¬ned as follows

u(xj , tn+1 ) ’ u(xj , tn ) u(xj+1 , tn ) ’ u(xj’1 , tn )

’a

n

„j = .

∆t 2∆x

The truncation error is

„ (∆t, ∆x) = max|„j |.

n

j,n

When „ (∆t, ∆x) goes to zero as ∆t and ∆x tend to zero independently,

the numerical scheme is said to be consistent.

Moreover, we say that it is of order p in time and of order q in space (for

suitable integers p and q), if, for a su¬ciently smooth solution of the exact

problem, we have

„ (∆t, ∆x) = O(∆tp + ∆xq ).

Using Taylor™s expansion conveniently we can characterize the truncation

error of the methods previously introduced as indicated in Table 13.1. The

Leap-frog and Newmark methods are both second order accurate if ∆t =

∆x, while the forward (or backward) Euler centred method is O(∆t+∆x2 ).

Finally, we say that a numerical scheme is convergent if

max|u(xj , tn ) ’ un | = 0.

lim j

∆t,∆x’0 j,n

13.8.2 Stability

A numerical method for a hyperbolic problem (linear or nonlinear) is said

to be stable if, for any time T , there exist two constants CT > 0 (possibly

depending on T ) and δ0 > 0, such that

¤ CT u0

un ∆, (13.46)

∆

for any n such that n∆t ¤ T and for any ∆t, ∆x such that 0 < ∆t ¤ δ0 ,

0 < ∆x ¤ δ0 . We have denoted by · ∆ a suitable discrete norm, for

606 13. Parabolic and Hyperbolic Initial Boundary Value Problems

instance one of those indicated below

« p

1

∞

= ∆x |vj |p = sup|vj |.

v p = 1, 2, v (13.47)

∆,p ∆,∞