<<

. 25
( 26)



>>

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

<<

. 25
( 26)



>>