= 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

u(xj , t ) = where (13.54)

k=’∞

Letting

•k = k∆x,

we have k∆t = »•k and thus

gk = e’ia»•k . (13.55)

The real number •k , here expressed in radians, is called the phase angle of

the k-th harmonic. Comparing (13.54) with (13.52) we can see that γk is

the counterpart of gk which is generated by the speci¬c numerical method

at hand. Moreover, |gk | = 1, whereas |γk | ¤ 1, in order to ensure stability.

Thus, γk is a dissipation coe¬cient; the smaller |γk |, the higher the reduc-

tion of the amplitude ±k , and, as a consequence, the higher the numerical

dissipation.

The ratio a (k) = |γk | is called the ampli¬cation error of the k-th harmonic

k|

|g

associated with the numerical scheme (in our case it coincides with the

ampli¬cation coe¬cient). On the other hand, writing

ω

»•k

’i

γk = |γk |e’iω∆t = |γk |e k ,

and comparing this relation with (13.55), we can identify the velocity of

propagation of the numerical solution, relative to its k-th harmonic, as

being ω . The ratio between this velocity and the velocity a of the exact

k

solution is called the dispersion error d relative to the k-th harmonic

ω ω∆x

d (k) = = .

ka •k a

The ampli¬cation and dispersion errors for the numerical schemes exam-

ined so far are functions of the phase angle •k and the CFL number a».

This is shown in Figure 13.7 where we have only considered the interval

0 ¤ •k ¤ π and we have used degrees instead of radians to denote the

values of •k .

In Figure 13.8 the numerical solutions of equation (13.26) with a = 1

and the initial datum u0 given by a packet of two sinusoidal waves of equal

wavelength l and centred at the origin x = 0 are shown. The plots on the

left-side of the ¬gure refer to the case l = 10∆x while on the right-side

13.9 Dissipation and Dispersion 613

Amplification error for Lax’Friedrichs

1

0.8

a

0.6

µ

CFL= 0.25

0.4 CFL= 0.50

CFL= 0.75

0.2 CFL= 1.00

60

20 80

40 100 120

0 140 160 180

φ

Dispersion error for Lax’Friedrichs

5 CFL= 0.25

CFL= 0.50

CFL= 0.75

4

CFL= 1.00

φ

µ

3

2

1

0 20 40 120

60 100 140 180

80 160

φ

Amplification error for Lax’Wendroff

1.2

1

0.8

0.6

a

µ

CFL= 0.25

0.4

CFL= 0.50

0.2 CFL= 0.75

CFL= 1.00

0

140

120 180

100

80 160

40

0 60

20

φ

Dispersion error for Lax’Wendroff

1.4

1.2

1

0.8

φ

µ

0.6

CFL= 0.25

0.4

CFL= 0.50

CFL= 0.75

0.2

CFL= 1.00

0

120 140 180

160

0 20 100

40 80

60

φ

Amplification error for Upwind

1.2

1

0.8

0.6

a

µ

0.4

CFL= 0.25

CFL= 0.50

0.2

CFL= 0.75

0

CFL= 1.00

’0.2

180

160

140

120

100

0 80

40

20 60

φ

Dispersion error for Upwind

1.4

1.2

1

0.8

µφ

0.6

CFL= 0.25

0.4

CFL= 0.50

CFL= 0.75

0.2

CFL= 1.00

0

140 160 180

120

60 80 100

20 40

0

φ

FIGURE 13.7. Ampli¬cation and dispersion errors for several numerical schemes

we have l = 4∆x. Since k = (2π)/l we get •k = ((2π)/l)∆x, so that

•k = π/10 in the left-side pictures and •k = π/4 in the right-side ones.

All numerical solutions have been computed for a CFL number equal to

0.75, using the schemes introduced above. Notice that the dissipation e¬ect

is quite relevant at high frequencies (•k = π/4), especially for ¬rst-order

methods (such as the upwind and the Lax-Friedrichs methods).

In order to highlight the e¬ects of the dispersion, the same computa-

tions have been repeated for •k = π/3 and di¬erent values of the CFL

number. The numerical solutions after 5 time steps are shown in Figure

13.9. The Lax-Wendro¬ method is the least dissipative for all the consid-

ered CFL numbers. Moreover, a comparison of the positions of the peaks

of the numerical solutions with respect to the corresponding ones in the

614 13. Parabolic and Hyperbolic Initial Boundary Value Problems

Lax’Wendroff CFL= 0.75 φ=π/10

1

0.5

0

u

Computed at t=1

’0.5

Exact

’1

0

’1

’2

’3

’4 2 4

1 3

Lax’Friedrichs CFL= 0.75 φ=π/10

x

1

0.5

0

u

Computed at t=1

’0.5

Exact

’1

4

3

2

1

’1 0

’3 ’2

’4

Upwind CFL= 0.75 φ=π/10

x

1

0.5

0

u

Computed at t=1

’0.5

Exact

’1

4

’4 3

’2 1 2

0

’1

’3

x

Lax’Wendroff CFL= 0.75 φ=π/4

1

0.5

u

0

Computed at t=1

’0.5

Exact

’1

3 4

2

1

0

’1

’3 ’2

’4

Lax’Friedrichs CFL= 0.75 φ=π/4

x

1

0.5

0

u

Computed at t=1

’0.5

Exact

’1

’3

’4 ’2 0 2

’1 1 3 4

Upwind CFL= 0.75 φ=π/4

x

1

0.5

0

u

Computed at t=1

’0.5

Exact

’1

2 3 4

1

0

’1

’3

’4 ’2

x

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

wave packet with di¬erent wavelengths

exact solution shows that the Lax-Friedrichs scheme is a¬ected by a posi-

tive dispersion error, since the ”numerical” wave advances faster than the

exact one. Also, the upwind scheme exhibits a slight dispersion error for a

CFL number of 0.75 which is absent for a CFL number of 0.5. The peaks

are well aligned with those of the numerical solution, although they have

been reduced in amplitude due to numerical dissipation. Finally, the Lax-

Wendro¬ method exhibits a small negative dispersion error; the numerical

solution is indeed slightly late with respect to the exact one.

13.9.1 Equivalent Equations

Using Taylor™s expansion to the third order to represent the truncation er-

ror, it is possible to associate with any of the numerical schemes introduced

13.9 Dissipation and Dispersion 615

Lax’Wendroff CFL=0.75 φ=π/3, t=4∆ t

1

0.5

u

0

’0.5

Computed

Exact

’1

’0.5

’1 1

0 2

0.5 1.5

Lax’Wendroff CFL=0.50 φ=π/3, t=4∆ t

1

0.5

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

[u]j = u+ (xj ) ’ u’ (xj ), u± (xj ) = lim u(xj + sa), j = 1, . . . , n ’ 1.

± s’0

Then, for any t ∈ [0, T ] the stability estimate for problem (13.68) reads

«

t n’1

uh („ ) a(xj )[uh („ )]2 d„

2

2

uh (t) + +

L2 (±,β)

L2 (±,β) j

j=0

®

0

(13.69)

t

¤ C ° u0,h + a•2 („ ) d„ » .

2

2

+ f („ ) L2 (±,β)

L2 (±,β)

0

As for the convergence analysis, the following error estimate can be proved

for continuous ¬nite elements of degree r, r ≥ 1 (see [QV94], Section 14.3.1)

« 1/2

t

+ a|u(±, „ ) ’ uh (±, „ )|2 d„

max u(t) ’ uh (t) L2 (±,β)

t∈[0,T ]

0

= O( u0 ’ u0,h + hr ),

L2 (±,β)

If, instead, discontinuous ¬nite elements of degree r are used, r ≥ 0, the

convergence estimate becomes (see [QV94], Section 14.3.3 and the refer-

ences therein)

« T

+

max u(t) ’ uh (t) u(t) ’ uh (t) 2

dt

L2 (±,β) L2 (±,β)

t∈[0,T ]

0

1/2

T n’1

a(xj ) [u(t) ’ uh (t)]2 dt = O( u0 ’ u0,h + hr+1/2 ).

+ L2 (±,β)

j

0 j=0

13.10.2 Time Discretization

The time discretization of the ¬nite element schemes introduced in the

previous section can be carried out by resorting either to ¬nite di¬erences

or ¬nite elements. If an implicit ¬nite di¬erence scheme is adopted, both

method (13.65) and (13.68) are unconditionally stable.

As an example, let us use the backward Euler method for the time dis-

cretization of problem (13.65). We obtain for each n ≥ 0: ¬nd un+1 ∈ Vh

h

13.10 Finite Element Approximation of Hyperbolic Equations 621

such that

β β

‚un+1

1

(un+1 ’ un )vh dx + h

a vh dx

h

h

∆t ‚x

± ± (13.70)

β β

∀vh ∈ Vh

an+1 un+1 vh dx = f n+1 vh dx in

+ 0 h

± ±

with un+1 (±) = •n+1 and u0 = u0h . If f ≡ 0 and • ≡ 0, taking vh = un+1

h

h h

in (13.70) we can obtain

1

’ un ¤0

un+1 + a(β)(un+1 (β))2 + µ0 un+1

2 2 2

L2 (±,β) L2 (±,β) L2 (±,β)

h

h h h

2∆t

∀n ≥ 0. Summing for n from 0 to m ’ 1 yields for m ≥ 1

«

m m

+ 2∆t a(β)(uj+1 (β))2 ¤ u0

uj 2 2 (±,β)

um 2 2 (±,β) 2

+ L2 (±,β) .

hL h

hL h

j=1 j=1

In particular, we conclude that

¤ u0 ∀m ≥ 0.

um L2 (±,β) L2 (±,β)

h h

On the other hand, explicit schemes for the hyperbolic equations are

subject to a stability condition: for example, in the case of the forward

Euler method the stability condition is ∆t = O(∆x). In practice, this

restriction is not as severe as happens in the case of parabolic equations

and for this reason explicit schemes are often used in the approximation of

hyperbolic equations.

Programs 102 and 103 provide an implementation of the discontinous

Galerkin-¬nite element method of degree 0 (dG(0)) and 1 (dG(1)) in space

coupled with the backward Euler method in time for the solution of (13.26)

on the space-time domain (±, β) — (t0 , T ).

Program 102 - ipeidg0 : dG(0) implicit Euler

function [u,x]=ipeidg0(I,n,a,u0,bc)

nx = n(1); h = (I(2)-I(1))/nx;

x = [I(1)+h/2:h:I(2)]; t = I(3); u = (eval(u0))™;

nt = n(2); k = (I(4)-I(3))/nt;

lambda = k/h; e = ones(nx,1);

A=spdiags([-a*lambda*e, (1+a*lambda)*e],-1:0,nx,nx);

[L,U]=lu(A);

for t = I(3)+k:k:I(4)

622 13. Parabolic and Hyperbolic Initial Boundary Value Problems

f = u;

if a > 0

f(1) = a*bc(1)+f(1);

elseif a <= 0

f(nx) = a*bc(2)+f(nx);

end

y = L \ f; u = U \ y;

end

Program 103 - ipeidg1 : dG(1) implicit Euler

function [u,x]=ipeidg1(I,n,a,u0,bc)

nx = n(1); h = (I(2)-I(1))/nx;

x = [I(1):h:I(2)]; t = I(3); um = (eval(u0))™;

u = []; xx=[];

for i = 1:nx+1

u = [u, um(i), um(i)]; xx = [xx, x(i), x(i)];

end

u = u™; nt = n(2); k = (I(4)-I(3))/nt;

lambda = k/h; e = ones(2*nx+2,1);

B = spdiags([1/6*e,1/3*e,1/6*e],-1:1,2*nx+2,2*nx+2);

dd = 1/3+0.5*a*lambda; du = 1/6+0.5*a*lambda;

dl = 1/6-0.5*a*lambda; A=sparse([]);

A(1,1) = dd; A(1,2) = du; A(2,1) = dl; A(2,2) = dd;

for i=3:2:2*nx+2

A(i,i-1) =-a*lambda;

A(i,i) = dd; A(i,i+1) = du;

A(i+1,i) = dl; A(i+1,i+1) = A(i,i);

end

[L,U]=lu(A);

for t = I(3)+k:k:I(4)

f = B*u;

if a > 0

f(1) = a*bc(1)+f(1);

elseif a <= 0

f(nx) = a*bc(2)+f(nx);

end

y = L \ f;

u = U \ y;

end

x = xx;

13.11 Applications 623

13.11 Applications

13.11.1 Heat Conduction in a Bar

Consider a homogeneous bar of unit length with thermal conductivity ν,

which is connected at the endpoints to an external thermal source at a ¬xed

temperature, say u = 0. Let u0 (x) be the temperature distribution along

the bar at time t = 0 and f = f (x, t) be a given heat production term.

Then, the initial-boundary value problem (13.1)-(13.4) provides a model of

the time evolution of the temperature u = u(x, t) throughout the bar.

In the following, we study the case where f ≡ 0 and the temperature of

the bar is suddenly raised at the points around 1/2. A rough mathematical

model for this situation is provided, for instance, by taking u0 = K in

a certain subinterval [a, b] ⊆ [0, 1] and equal to 0 outside, where K is a

given positive constant. The initial condition is therefore a discontinuous

function.

We have used the θ-method with θ = 0.5 (Crank-Nicolson method, CN)

and θ = 1 (Backward Euler method, BE). Program 100 has been run

with h = 1/20, ∆t = 1/40 and the obtained solutions at time t = 2 are

shown in Figure 13.10. The results show that the CN method su¬ers a clear

instability due to the low smoothness of the initial datum (about this point,

see also [QV94], Chapter 11). On the contrary, the BE method provides a

stable solution which decays correctly to zero as t grows since the source

term f is null.

1

0.8

0.9

0.8

0.6

0.7

0.6

0.4

0.5

0.4

0.2

0.3

0.2

0

0.1

’0.2 0

0

0

0.5

0.5

1

1

0 0

0.2

1.5 1.5 0.4

0.5 0.6

0.8

2 1 2 1

FIGURE 13.10. Solutions for a parabolic problem with discontinuous initial da-

tum: CN method (left) and BE method (right)

13.11.2 A Hyperbolic Model for Blood Flow Interaction with

Arterial Walls

Let us consider again the problem of the ¬‚uid-structure interaction in a

cylindrical artery considered in Section 11.11.2, where the simple indepen-

dent rings model (11.88) was adopted.

624 13. Parabolic and Hyperbolic Initial Boundary Value Problems

If the axial action due to the tension between the di¬erent rings is

no longer neglected, denoting by z the longitudinal coordinate, equation

(11.88) modi¬es into

‚2· ‚ 2 · HE

ρw H 2 ’ σz 2 + 2 · = P ’ P0 , t > 0, 0<z<L (13.71)

‚t ‚z R0

where σz is the radial component of the axial stress and L is the length of

the cylindrical arterial district which is considered. In particular, neglecting

the third term on the left hand side and letting γ 2 = σz /(ρw H), f =

(P ’ P0 )/(ρw H), we recover the wave equation (13.33).

We have performed two sets of numerical experiments using the Leap-

Frog (LF) and Newmark (NW) methods. In the ¬rst example the space-

time domain of integration is the cylinder (0, 1) — (0, 1) and the source

term is f = (1 + π 2 γ 2 )e’t sin(πx) in such a way that the exact solution is

u(x, t) = e’t sin(πx). Table 13.3 shows the estimated orders of convergence

of the two methods, denoted by pLF and pN W respectively.

To compute these quantities we have ¬rst solved the wave equation on

(k)

four grids with sizes ∆x = ∆t = 1/(2k · 10), k = 0, . . . , 3. Let uh denote

the numerical solution corresponding to the space-time grid at the k-th

(0)

re¬ning level. Moreover, for j = 1, . . . , 10, let tj = j/10 be the time

discretization nodes of the grid at the coarsest level k = 0. Then, for each

level k, the maximum nodal errors ek on the k-th spatial grid have been

j

(0) (k)

evaluated at each time tj in such a way that the convergence order pj

can be estimated as

log(e0 /ek )

j j

(k)

pj = , k = 1, 2, 3.

k)

log(2

The results clearly show second-order convergence for both the methods,

as theoretically expected.

In the second example we have taken the following expressions for the

coe¬cient and source term: γ 2 = σz /(ρw H), with σz = 1 [Kgs’2 ], f =

(x∆p · sin(ω0 t))/(ρw H). The parameters ρw , H and the length L of the

vessel are as in Section 11.11.2. The space-time computational domain is

(0, L) — (0, T ), with T = 1 [s].

The Newmark method has been ¬rst used with ∆x = L/10 and ∆t =

T /100; the corresponding value of γ» is 3.6515, where » = ∆t/∆x. Since the

Newmark method is unconditionally stable, no spurious oscillations arise as

is con¬rmed by Figure 13.11, left. Notice the correct periodical behaviour

of the solution with a period corresponding to one heart beat; notice also

that with the present values of ∆t and ∆x the Leap-Frog method cannot

be employed since the CFL condition is not satis¬ed. To overcome this

13.12 Exercises 625

(0) (1) (2) (3) (0) (1) (2) (3)

tj pLF pLF pLF tj pN W pN W pN W

0.1 2.0344 2.0215 2.0151 0.1 1.9549 1.9718 1.9803

0.2 2.0223 2.0139 2.0097 0.2 1.9701 1.9813 1.9869

0.3 2.0170 2.0106 2.0074 0.3 1.9754 1.9846 1.9892

0.4 2.0139 2.0087 2.0061 0.4 1.9791 1.9869 1.9909

0.5 2.0117 2.0073 2.0051 0.5 1.9827 1.9892 1.9924

0.6 2.0101 2.0063 2.0044 0.6 1.9865 1.9916 1.9941

0.7 2.0086 2.0054 2.0038 0.7 1.9910 1.9944 1.9961

0.8 2.0073 2.0046 2.0032 0.8 1.9965 1.9979 1.9985

0.9 2.0059 2.0037 2.0026 0.9 2.0034 2.0022 2.0015

1.0 2.0044 2.0028 2.0019 1.0 2.0125 2.0079 2.0055

TABLE 13.3. Estimated orders of convergence for the Leap-Frog (LF) and New-

mark (NW) methods

problem, we have therefore chosen a much smaller time-step ∆t = T /400,

in such a way that γ» 0.9129 and the Leap-Frog scheme can be applied.

The obtained result is shown in Figure 13.11, right; a similar solution has

been computed by using the Newmark method with the same values of the

discretization parameters.

x 10’4 x 10’4

3

3

2

2

1

1

0

0

’1

’1

’2

’2

0.06

’3

’3 0.04

0.04

1

0.8 0.03 1 0.8

0.6 0.02

0.02 0.6

0.4 0.4

0.01

0.2 0.2 0

0 0

0

FIGURE 13.11. Computed solutions using the NM method on a space-time grid

with ∆t = T /100 and ∆x = L/10 (left) and the LF scheme on a space-time grid

with the same value of ∆x but with ∆t = T /400 (right)

13.12 Exercises

1. Apply the θ-method (13.9) to the approximate solution of the scalar Cauchy

problem (11.1) and using the analysis of Section 11.3 prove that the local

truncation error is of the order of ∆t + h2 if θ = 1 while it is of the order

2

of ∆t2 + h2 if θ = 1 .

2

2. Prove that in the case of piecewise linear ¬nite elements, the mass-lumping

process described in Section 13.3 is equivalent to computing the integrals

626 13. Parabolic and Hyperbolic Initial Boundary Value Problems

1

mij = 0 •j •i dx by the trapezoidal quadrature formula (9.11). This, in

particular, shows that the diagonal matrix M is nonsingular.

[Hint: ¬rst, verify that exact integration yields

±

1 i = j,

h 2

mij =

6

1 i = j.

Then, apply the trapezoidal rule to compute mij recalling that •i (xj ) =

δij .]

3. Prove inequality (13.19).

[Hint: using the Cauchy-Schwarz and Young inequalities, prove ¬rst that

1

1

(u ’ v)u dx ≥ ’v ∀ u, v ∈ L2 (0, 1).

2 2

u ,

L2 (0,1) L2 (0,1)

2

0

Then, use (13.18). ]

4. Assume that the bilinear form a(·, ·) in problem (13.12) is continuous and

coercive over the function space V (see (12.54)-(12.55)) with continuity and

coercivity constants M and ±, respectively. Then, prove that the stability

inequalities (13.20) and (13.21) still hold provided that ν is replaced by ±.

5. Show that the methods (13.39), (13.40) and (13.41) can be written in the

form (13.42). Then, show that the corresponding expressions of the arti¬cial

viscosity K and arti¬cial di¬usion ¬‚ux hdif f are as in Table (13.1).

j+1/2

6. Determine the CFL condition for the upwind scheme.

7. Show that for the scheme (13.43) one has un+1 ∆,2 ¤ un for all

∆,2

n ≥ 0.

[Hint: multiply equation (13.43) by un+1 , and notice that

j

1

(un+1 ’ un )un+1 ≥ |un+1 |2 ’ |un |2 .

j j

j j j

2

Then, sum on j the resulting inequalities, and note that

∞

»a

un+1 ’ un+1 un+1 = 0

2 j=’∞ j+1 j’1 j

since this sum is telescopic.]

8. Show how to ¬nd the values µ and ν in Table 13.2 for Lax-Friedrichs and

Lax-Wendro¬ methods.

9. Prove (13.67).

10. Prove (13.69) when f = 0.

[Hint: take ∀t > 0, vh = uh (t) in (13.68).]

References

[Aas71] Aasen J. (1971) On the Reduction of a Symmetric Matrix to

Tridiagonal Form. BIT 11: 233“242.

[ABB+ 92] Anderson E., Bai Z., Bischof C., Demmel J., Dongarra J., Croz

J. D., Greenbaum A., Hammarling S., McKenney A., Ous-

trouchov S., and Sorensen D. (1992) LAPACK User™s Guide,

Release 1.0. SIAM, Philadelphia.

[Ada75] Adams D. (1975) Sobolev Spaces. Academic Press, New York.

[ADR92] Arioli M., Du¬ I., and Ruiz D. (1992) Stopping Criteria for

Iterative Solvers. SIAM J. Matrix Anal. Appl. 1(13).

[AF83] Alonso M. and Finn E. (1983) Fundamental University

Physics, volume 3. Addison-Wesley, Reading, Massachusetts.

[Arm66] Armijo L. (1966) Minimization of Functions Having Continu-

ous Partial Derivatives. Paci¬c Jour. Math. 16: 1“3.

[Arn73] Arnold V. I. (1973) Ordinary Di¬erential Equations. The MIT

Press, Cambridge, Massachusetts.

[Atk89] Atkinson K. E. (1989) An Introduction to Numerical Analysis.

John Wiley, New York.

[Avr76] Avriel M. (1976) Non Linear Programming: Analysis and

Methods. Prentice-Hall, Englewood Cli¬s, New Jersey.

628 References

[Axe94] Axelsson O. (1994) Iterative Solution Methods. Cambridge

University Press, New York.

[Bar89] Barnett S. (1989) Leverrier™s Algorithm: A New Proof and

Extensions. Numer. Math. 7: 338“352.

[Bat90] Batterson S. (1990) Convergence of the Shifted QR Algorithm

on 3 by 3 Normal Matrices. Numer. Math. 58: 341“352.

[BBC+ 94] Barrett R., Berry M., Chan T., Demmel J., Donato J., Don-

garra J., Eijkhout V., Pozo V., Romine C., and van der Vorst

H. (1994) Templates for the Solution of Linear Systems: Build-

ing Blocks for Iterative Methods. SIAM, Philadelphia.

[BD74] Bj¨rck A. and Dahlquist G. (1974) Numerical Methods.

o

Prentice-Hall, Englewood Cli¬s, N.J.

[BDMS79] Bunch J., Dongarra J., Moler C., and Stewart G. (1979) LIN-

PACK User™s Guide. SIAM, Philadelphia.

[Ber82] Bertsekas D. P. (1982) Constrained Optimization and La-

grange Multiplier Methods. Academic Press. Inc., San Diego,

California.

[Bj¨88]

o Bj¨rck A. (1988) Least Squares Methods: Handbook of Numer-

o

ical Analysis Vol. 1 Solution of Equations in RN . Elsevier

North Holland.

[BM92] Bernardi C. and Maday Y. (1992) Approximations Spectrales

des Probl´mes aux Limites Elliptiques. Springer-Verlag, Paris.

e

[BMW67] Barth W., Martin R. S., and Wilkinson J. H. (1967) Calcula-

tion of the Eigenvalues of a Symmetric Tridiagonal Matrix by

the Method of Bisection. Numer. Math. 9: 386“393.

[BO78] Bender C. M. and Orszag S. A. (1978) Advanced Mathemati-

cal Methods for Scientists and Engineers. McGraw-Hill, New

York.

[Boe80] Boehm W. (1980) Inserting New Knots into B-spline Curves.

Computer Aided Design 12: 199“201.

[Bos93] Bossavit A. (1993) Electromagnetisme, en vue de la modelisa-

tion. Springer-Verlag, Paris.

[BR81] Bank R. E. and Rose D. J. (1981) Global Approximate Newton

Methods. Numer. Math. 37: 279“295.

[Bra75] Bradley G. (1975) A Primer of Linear Algebra. Prentice-Hall,

Englewood Cli¬s, New York.

References 629

[Bre73] Brent R. (1973) Algorithms for Minimization Without Deriva-

tives. Prentice-Hall, Englewood Cli¬s, New York.

[Bri74] Brigham E. O. (1974) The Fast Fourier Transform. Prentice-

Hall, Englewood Cli¬s, New York.

[BS90] Brown P. and Saad Y. (1990) Hybrid Krylov Methods for Non-

linear Systems of equations. SIAM J. Sci. and Stat. Comput.

11(3): 450“481.

[BSG96] B. Smith P. B. and Gropp P. (1996) Domain Decomposition,

Parallel Multilevel Methods for Elliptic Partial Di¬erential

Equations. Univ. Cambridge Press, Cambridge.

[But64] Butcher J. C. (1964) Implicit Runge-Kutta Processes. Math.

Comp. 18: 233“244.

[But66] Butcher J. C. (1966) On the Convergence of Numerical So-

lutions to Ordinary Di¬erential Equations. Math. Comp. 20:

1“10.

[But87] Butcher J. (1987) The Numerical Analysis of Ordinary Di¬er-

ential Equations: Runge-Kutta and General Linear Methods.

Wiley, Chichester.

[CCP70] Cannon M., Cullum C., and Polak E. (1970) Theory and Op-

timal Control and Mathematical Programming. McGraw-Hill,

New York.

¨

[CFL28] Courant R., Friedrichs K., and Lewy H. (1928) Uber die

partiellen di¬erenzengleichungen der mathematischen physik.

Math. Ann. 100: 32“74.

[CHQZ88] Canuto C., Hussaini M. Y., Quarteroni A., and Zang T. A.

(1988) Spectral Methods in Fluid Dynamics. Springer, New

York.

[CI95] Chandrasekaren S. and Ipsen I. (1995) On the Sensitivity of

Solution Components in Linear Systems of equations. SIAM

J. Matrix Anal. Appl. 16: 93“112.

[CL91] Ciarlet P. G. and Lions J. L. (1991) Handbook of Numerical

Analysis: Finite Element Methods (Part 1). North-Holland,

Amsterdam.

[CM94] Chan T. and Mathew T. (1994) Domain Decomposition Algo-

rithms. Acta Numerica pages 61“143.

630 References

[CMSW79] Cline A., Moler C., Stewart G., and Wilkinson J. (1979) An

Estimate for the Condition Number of a Matrix. SIAM J. Sci.

and Stat. Comput. 16: 368“375.

[Col66] Collin R. E. (1966) Foundations for Microwave Engineering.

McGraw-Hill Book Co., Singapore.

[Com95] Comincioli V. (1995) Analisi Numerica Metodi Modelli Appli-

cazioni. McGraw-Hill Libri Italia, Milano.

[Cox72] Cox M. (1972) The Numerical Evaluation of B-splines. Jour-

nal of the Inst. of Mathematics and its Applications 10: 134“

149.

[Cry73] Cryer C. W. (1973) On the Instability of High Order

Backward-Di¬erence Multistep Methods. BIT 13: 153“159.