0

u

’0.5

Computed

Exact

’1

’1 0.5

0 1

’0.5 1.5 2

Lax’Friedrichs CFL=0.75 φ=π/3, t=5∆ t

1

0.5

0

u

’0.5

Computed

Exact

’1

1.5 2

’1 1

’0.5 0.5

0

Lax’Friedrichs CFL=0.50 φ=π/3, t=5∆ t

1

0.5

0

u

’0.5

Computed

Exact

’1

’0.5 0 1

0.5 1.5

’1 2

Upwind CFL=0.75 φ=π/3, t=4∆ t

1

0.5

0

u

’0.5

Computed

Exact

’1

2

’0.5 1.5

0 0.5

’1 1

Upwind CFL=0.50 φ=π/3, t=4∆ t

1

0.5

0

u

’0.5

Computed

Exact

’1

2

1.5

1

0 0.5

’1 ’0.5

FIGURE 13.9. Numerical solutions corresponding to the transport of a sinusoidal

wave packet, for di¬erent CFL numbers

so far an equivalent di¬erential equation of the form

vt + avx = µvxx + νvxxx (13.56)

where the terms µvxx and νvxxx represent dissipation and dispersion, re-

spectively. Table 13.2 shows the values of µ and ν for the various methods.

Let us give a proof of this procedure in the case of the upwind scheme.

Let v(x, t) be a smooth function which satis¬es the di¬erence equation

(13.41); then, assuming that a > 0, we have

v(x, t + ∆t) ’ v(x, t) v(x, t) ’ v(x ’ ∆x, t)

+a = 0.

∆t ∆x

Truncating the Taylor expansions of v around (x, t) at the ¬rst and second

order, respectively, we obtain

vt + O(∆t) + avx + O(∆x) = 0 (13.57)

616 13. Parabolic and Hyperbolic Initial Boundary Value Problems

Method µ ν

2

a∆x a ∆t a

’ ’ ∆x2 ’ 3a∆x∆t + 2a2 ∆t2

Upwind

2 2 6

2

a∆x2

∆x

1 ’ (a»)2 1 ’ (a»)2

Lax-Friedrichs

2∆t 3

a∆x2

(a»)2 ’ 1

Lax-Wendro¬ 0

6

TABLE 13.2. Values of dissipation and dispersion coe¬cients for several numer-

ical methods

and

∆t ∆x

vtt + O(∆t2 ) + a vx + vxx + O(∆x2 )

vt + = 0, (13.58)

2 2

where vt = ‚v and vx = ‚x .

‚v

‚t

Di¬erentiating (13.57) with respect to t and then with respect to x, we

get

vtt + avxt = O(∆x + ∆t),

and

vtx + avxx = O(∆x + ∆t).

Thus, it follows that

vtt = a2 vxx + O(∆x + ∆t),

which, after substituting into (13.58), yields the following equation

vt + avx = µvxx (13.59)

where

a∆x a2 ∆t

’

µ= ,

2 2

and having neglected the term O(∆x2 +∆t2 ). Relation (13.59) is the equiv-

alent di¬erential equation up to second order of the upwind scheme.

Following the same procedure and truncating the Taylor expansion at

third order, yields

vt + avx = µvxx + νvxxx (13.60)

where

a22

a ∆t ’ ∆x2 .

ν=

6

We can give a heuristic explanation of the meaning of the dissipative

and dispersive terms in the equivalent equation (13.56) by studying the

following problem

x ∈ R, t > 0

vt + avx = µvxx + νvxxx

(13.61)

(k ∈ Z)

v(x, 0) = eikx ,

13.9 Dissipation and Dispersion 617

Applying the Fourier transform yields, if µ = ν = 0,

v(x, t) = eik(x’at) , (13.62)

while if µ and ν are arbitrary real numbers (with µ > 0) we get

2 2

v(x, t) = e’µk t eik[x’(a+νk )t]

. (13.63)

Comparing (13.62) with (13.63) we can see that the module of the solution

diminishes as µ grows and this becomes more relevant if the frequency k

gets larger. Therefore, the term µvxx in (13.61) has a dissipative e¬ect on

the solution. A further comparison between (13.62) and (13.63) shows that

the presence of the term ν modi¬es the velocity of the propagation of the

solution; the velocity is increased if ν > 0 whereas it is diminuished if ν < 0.

Even in this case the e¬ect is ampli¬ed at high frequencies. Therefore, the

third-order di¬erential term νvxxx introduces a dispersive e¬ect.

Generally speaking, even-order derivatives in the equivalent equation rep-

resent di¬usive terms, while odd-order derivatives mean dispersive e¬ects.

In the case of ¬rst-order schemes (like the upwind method) the dispersive

e¬ect is often only slightly visible since it is hidden by the dissipative one.

Actually, taking ∆t and ∆x of the same order, we have that ν << µ as

∆x ’ 0, since ν = O(∆x2 ) and µ = O(∆x). In particular, for a CFL

number of 1 , the equivalent equation of the upwind method exhibits null

2

dispersion, truncated at second order, according to the results of the pre-

vious section.

On the other hand, the dispersive e¬ect is strikingly visible in the Lax-

Friedrichs and in the Lax-Wendro¬ schemes; the latter, being second-order

accurate, does not exhibit a dissipative term of the form µvxx . However,

it ought to be dissipative in order to be stable; actually, the equivalent

equation (truncated at fourth order) for the Lax-Wendro¬ scheme reads

a∆x2 a∆x3

[(a») ’ 1]vxxx ’ a»[1 ’ (a»)2 ]vxxxx

2

vt + avx =

6 6

where the last term is dissipative if |a»| < 1. We thus recover the CFL

condition for the Lax-Wendro¬ method.

618 13. Parabolic and Hyperbolic Initial Boundary Value Problems

13.10 Finite Element Approximation of Hyperbolic

Equations

Let us consider the following ¬rst-order linear, scalar hyperbolic problem

in the interval „¦ = (±, β) ‚ R

±

‚u + a ‚u + a0 u = f in QT = „¦ — (0, T )

‚t

‚x

(13.64)

t ∈ (0, T )

u(±, t) = •(t)

x ∈ „¦,

u(x, 0) = u0 (x)

where a = a(x), a0 = a0 (x, t), f = f (x, t), • = •(t) and u0 = u0 (x) are

given functions.

We assume that a(x) > 0 ∀x ∈ [±, β]. In particular, this implies that

the point x = ± is the in¬‚ow boundary, and the boundary value has to be

speci¬ed there.

13.10.1 Space Discretization with Continuous and

Discontinuous Finite Elements

A semi-discrete approximation of problem (13.64) can be carried out by

means of the Galerkin method (see Section 12.4). De¬ne the spaces

Vh = Xh = vh ∈ C 0 („¦) : vh|Ij ∈ Pr (Ij ), ∀ Ij ∈ Th

r

and

Vh = {vh ∈ Vh : vh (±) = 0} ,

in

where Th is a partition of „¦ (see Section 12.4.5) into n ≥ 2 subintervals

Ij = [xj , xj+1 ], for j = 0, . . . , n ’ 1.

Let u0,h be a suitable ¬nite element approximation of u0 and consider

the problem: for any t ∈ (0, T ) ¬nd uh (t) ∈ Vh such that

± β β

‚uh (t) ‚uh (t)

vh dx + a + a0 (t)uh (t) vh dx

‚t ‚x

± ±

β (13.65)

∀ vh ∈ in

= f (t)vh dx Vh

±

uh (t) = •h (t) at x = ±,

with uh (0) = u0,h ∈ Vh .

13.10 Finite Element Approximation of Hyperbolic Equations 619

If • is equal to zero, uh (t) ∈ Vh , and we are allowed to taking vh = uh (t)

in

and get the following inequality

t t

2

d„ + a(β) u2 („ ) d„

2

uh (t) + µ0 uh („ ) L2 (±,β)

L2 (±,β) h

0 0

t

1

¤ 2

2

u0,h + f („ ) L2 (±,β) d„ ,

L2 (±,β)

µ0

0

for any t ∈ [0, T ], where we have assumed that

1

0 < µ0 ¤ a0 (x, t) ’ a (x). (13.66)

2

Notice that in the special case in which both f and a0 are identically zero,

we obtain

uh (t) L2 (±,β) ¤ u0,h L2 (±,β)

which expresses the conservation of the energy of the system. When (13.66)

does not hold (for example, if a is a constant convective term and a0 = 0),

then an application of Gronwall™s lemma 11.1 yields

t

+ a(β) u2 („ ) d„

2

uh (t) L2 (±,β) h

«

0 (13.67)

t t

¤ u0,h d„ exp [1 + 2µ— („ )] d„,

2

2

+ f („ ) L2 (±,β)

L2 (±,β)

0 0

where µ— (t) = max|µ(x, t)|.

[±,β]

An alternative approach to the semi-discrete approximation of problem

(13.64) is based on the use of discontinuous ¬nite elements. This choice

is motivated by the fact that, as previously pointed out, the solutions of

hyperbolic problems (even in the linear case) may exhibit discontinuities.

The ¬nite element space can be de¬ned as follows

Wh = Yh = vh ∈ L2 (±, β) : vh|Ij ∈ Pr (Ij ), ∀ Ij ∈ Th ,

r

i.e., the space of piecewise polynomials of degree less than or equal to r,

which are not necessarily continuous at the ¬nite element nodes.

Then, the Galerkin discontinuous ¬nite element space discretization reads:

for any t ∈ (0, T ) ¬nd uh (t) ∈ Wh such that

±x

±β

n’1 i+1

‚uh (t) ‚uh (t)

vh dx + a + a0 (x)uh (t) vh dx

‚t ‚x

i=0

± xi

(13.68)

β

+a(u+ ’ U ’ )(x , t)v + (x ) = f (t)v dx ∀v ∈ W ,

i i h h h

h h h

±

620 13. Parabolic and Hyperbolic Initial Boundary Value Problems

where {xi } are the nodes of Th with x0 = ±, and for each node xi , vh (xi ) +

’

denotes the right-value of vh at xi while vh (xi ) is its left-value. Finally,

Uh (xi , t) = u’ (xi , t) if i = 1, . . . , n ’ 1, while Uh (x0 , t) = •(t) ∀t > 0.

’ ’

h

If a is positive, xj is the in¬‚ow boundary of Ij for every j and we set