<<

. 24
( 26)



>>


uk+1 ’ uk
+ a θuk+1 + (1 ’ θ)uk , vh
h
h
, vh h
h
∆t (13.17)
= θF k+1 (vh ) + (1 ’ θ)F k (vh ) ∀vh ∈ Vh
1
for k ≥ 0 and with u0 = u0 h , F k (vh ) = 0 f (tk )vh (x)dx. Since we are
h
interested in the stability analysis, we can consider the special case where
F = 0; moreover, for the time being, we focus on the case θ = 1 (implicit
Euler scheme), i.e.

uk+1 ’ uk
∀vh ∈ Vh .
+ a uk+1 , vh = 0
h
h
, vh h
∆t

Letting vh = uk+1 , we get
h

uk+1 ’ uk k+1
+ a(uk+1 , uk+1 ) = 0.
h
h
, uh h h
∆t

From the de¬nition of a(·, ·), it follows that
2
‚uk+1
uk+1 , uk+1 h
a =ν . (13.18)
h h
‚x
L2 (0,1)

Moreover, we remark that (see Exercise 3 for the proof of this result)
2
‚uk+1
¤ uk
uk+1 2 2 (0,1) 2
h
+ 2ν∆t L2 (0,1) . (13.19)
h
L
h
‚x
L2 (0,1)

It follows that, ∀n ≥ 1
2
n’1 n’1 n’1
‚uk+1
¤
uk+1 2 2 (0,1) uk 2
h
+ 2ν∆t L2 (0,1) .
h
L
h
‚x
L2 (0,1)
k=0 k=0 k=0
13.3 Finite Element Approximation of the Heat Equation 589

Since these are telescopic sums, we get
2
n’1
‚uk+1
¤ u0h
un 2 2 (0,1) 2
h
+ 2ν∆t L2 (0,1) , (13.20)
hL
‚x
L2 (0,1)
k=0


which shows that the scheme is unconditionally stable. Proceeding similarly
if f = 0, it can be shown that
2
n’1
‚uk+1
h
un 2
+ 2ν∆t
L2 (0,1)
h
‚x
L2 (0,1)
k=0 (13.21)
n
¤ C(n) ∆t f k 2
2
u0 h + ,
L2 (0,1)
L2 (0,1)
k=1

where C(n) is a constant independent of both h and ∆t.

Remark 13.1 The same kind of stability inequalities (13.20) and (13.21)
can be obtained if a(·, ·) is a more general bilinear form provided that it is
continuous and coercive (see Exercise 4).


To carry out the stability analysis of the θ-method for every θ ∈ [0, 1] we
need de¬ning the eigenvalues and eigenvectors of a bilinear form.

De¬nition 13.1 We say that » is an eigenvalue and w ∈ V is the associ-
ated eigenvector for the bilinear form a(·, ·) : V — V ’ R if

∀v ∈ V,
a(w, v) = »(w, v)

where (·, ·) denotes the usual scalar product in L2 (0, 1).

If the bilinear form a(·, ·) is symmetric and coercive, it has in¬nitely many
real positive eigenvalues that form an unbounded sequence; moreover, its
eigenvectors (called also eigenfunctions) form a basis for the space V .
At a discrete level the corresponding pair »h ∈ R, wh ∈ Vh satis¬es

∀vh ∈ Vh .
a(wh , vh ) = »h (wh , vh ) (13.22)

From the algebraic standpoint, problem (13.22) can be formulated as

Afe w = »h Mw

(where w is the vector of the gridvalues of wh ) and can be regarded
as a generalized eigenvalue problem (see Section 5.9). All the eigenvalues
»1 , . . . , »Nh are positive. The corresponding eigenvectors wh , . . . , wh h form
N
1
h h
590 13. Parabolic and Hyperbolic Initial Boundary Value Problems

a basis for the subspace Vh and can be chosen in such a way as to be or-
j
thonormal, i.e., such that (wh , wh ) = δij , ∀i, j = 1, . . . , Nh . In particular,
i

any function vh ∈ Vh can be represented as
Nh
j
vh (x) = vj wh (x).
j=1


Let us now assume that θ ∈ [0, 1] and focus on the case where the bilinear
form a(·, ·) is symmetric. Although the ¬nal stability result still holds in
the nonsymmetric case, the proof that follows cannot apply since in that
i
case the eigenvectors would no longer form a basis for Vh . Let wh be the
eigenvectors of a(·, ·) whose span forms an orthonormal basis for Vh . Since
at each time step uk ∈ Vh , we can express uk as
h h

Nh
j
uk (x) uk wh (x).
=
h j
j=1

i
Letting F = 0 in (13.17) and taking vh = wh , we ¬nd
Nh
1 j
uk+1 ’ uk i
wh , w h
j
j
∆t j=1
Nh
j
θuk+1 + (1 ’ θ)uk a(wh , wh ) = 0,
i
+ i = 1, . . . , Nh .
j
j
j=1

j
Since wh are eigenfunctions of a(·, ·) we obtain

a(wh , wh ) = »j (wh , wh ) = »j δij = »i ,
j j
i i
h
h h

so that
uk+1 ’ uk
+ θuk+1 + (1 ’ θ)uk »i = 0.
i
i
i h
i
∆t
Solving this equation with respect to uk+1 gives
i

1 ’ (1 ’ θ)»i ∆t
h
k+1
uk
ui = .
i i ∆t
1 + θ»h

In order for the method to be unconditionally stable we must have (see
Chapter 11)
1 ’ (1 ’ θ)»i ∆t
h
< 1,
i ∆t
1 + θ»h
that is
2
2θ ’ 1 > ’ .
»i ∆t
h
13.3 Finite Element Approximation of the Heat Equation 591

If θ ≥ 1/2, this inequality is satis¬ed for any value of ∆t. Conversely, if
θ < 1/2 we must have
2
∆t < .
(1 ’ 2θ)»i
h

Since this relation must hold for all the eigenvalues »i of the bilinear form,
h
it su¬ces requiring that it is satis¬ed for the largest of them, which we
assume to be »Nh . h
We therefore conclude that if θ ≥ 1/2 the θ-method is unconditionally
stable (i.e., it is stable ∀∆t), whereas if 0 ¤ θ < 1/2 the θ-method is stable
only if
2
∆t ¤ .
(1 ’ 2θ)»Nh
h

It can be shown that there exist two positive constants c1 and c2 , indepen-
dent of h, such that
c1 h’2 ¤ »Nh = c2 h’2
h

(see for the proof, [QV94], Section 6.3.2). Accounting for this, we obtain
that if 0 ¤ θ < 1/2 the method is stable only if

∆t ¤ C1 (θ)h2 , (13.23)

for a suitable constant C1 (θ) independent of both h and ∆t.

With an analogous proof, it can be shown that if a pseudo-spectral Galerkin
approximation is used for problem (13.12), the θ’method is uncondition-
ally stable if θ ≥ 1 , while for 0 ¤ θ < 1 stability holds only if
2 2

∆t ¤ C2 (θ)N ’4 , (13.24)

for a suitable constant C2 (θ) independent of both N and ∆t. The di¬erence
between (13.23) and (13.24) is due to the fact that the largest eigenvalue
of the spectral sti¬ness matrix grows like O(N 4 ) with respect to the degree
of the approximating polynomial.

Comparing the solution of the globally discretized problem (13.17) with
that of the semi-discrete problem (13.13), by a suitable use of the stability
result (13.21) and of the truncation time discretization error, the following
convergence result can be proved

u(tk ) ’ uk ¤ C(u0 , f, u)(∆tp(θ) + hr+1 ), ∀k ≥ 1
L2 (0,1)
h

where r denotes the piecewise polynomial degree of the ¬nite element space
Vh , p(θ) = 1 if θ = 1/2 while p(1/2) = 2 and C is a constant that depends
on its arguments (assuming that they are su¬ciently smooth) but not on
592 13. Parabolic and Hyperbolic Initial Boundary Value Problems

h and ∆t. In particular, if f ≡ 0 on can obtain the following improved
estimates
r+1 p(θ)
h ∆t

u(t ) ’ ¤C
k
uk L2 (0,1) + u0 L2 (0,1) ,
h
tk
tk

for k ≥ 1, θ = 1 or θ = 1/2. (For the proof of these results, see [QV94], pp.
394-395).

Program 100 provides an implementation of the θ-method for the solution
of the heat equation on the space-time domain (a, b) — (t0 , T ). The dis-
cretization in space is based on piecewise-linear ¬nite elements. The input
parameters are: the column vector I containing the endpoints of the space
interval (a = I(1), b = I(2)) and of the time interval (t0 = I(3), T = I(4));
the column vector n containing the number of steps in space and time; the
macros u0 and f containing the functions u0h and f , the constant viscosity
nu, the Dirichlet boundary conditions bc(1) and bc(2), and the value of
the parameter theta.

Program 100 - thetameth : θ-method for the heat equation
function [u,x] = thetameth(I,n,u0,f,bc,nu,theta)
nx = n(1); h = (I(2)-I(1))/nx;
x = [I(1):h:I(2)]; t = I(3); uold = (eval(u0))™;
nt = n(2); k = (I(4)-I(3))/nt; e = ones(nx+1,1);
K = spdiags([(h/(6*k)-nu*theta/h)*e, (2*h/(3*k)+2*nu*theta/h)*e, ...
(h/(6*k)-nu*theta/h)*e],-1:1,nx+1,nx+1);
B = spdiags([(h/(6*k)+nu*(1-theta)/h)*e, (2*h/(3*k)-nu*2*(1-theta)/h)*e, ...
(h/(6*k)+nu*(1-theta)/h)*e],-1:1,nx+1,nx+1);
K(1,1) = 1; K(1,2) = 0; B(1,1) = 0; B(1,2) = 0;
K(nx+1,nx+1) = 1; K(nx+1,nx) = 0; B(nx+1,nx+1) = 0; B(nx+1,nx) = 0;
[L,U]=lu(K);
t = I(3); x=[I(1)+h:h:I(2)-h]; fold = (eval(f))™;
fold = h*fold;
fold = [bc(1); fold; bc(2)];
for time = I(3)+k:k:I(4)
t = time; fnew = (eval(f))™;
fnew = h*fnew;
fnew = [bc(1); fnew; bc(2)];
b = theta*fnew+(1-theta)*fold + B*uold;
y = L \ b; u = U \ y;
uold = u;
end
x = [I(1):h:I(2)];


Example 13.1 Let us assess the time-accuracy of the θ-method in the solution
of the heat equation (13.1) on the space-time domain (0, 1) — (0, 1) where f is
13.4 Space-Time Finite Element Methods for the Heat Equation 593

chosen in such a way that the exact solution is u = sin(2πx) cos(2πt). A ¬xed
spatial grid size h = 1/500 has been used while the time step ∆t is equal to
(10k)’1 , k = 1, . . . , 4. Finally, piecewise ¬nite elements are used for the space
discretization. Figure 13.1 shows the convergence behavior in the L2 (0, 1) norm
(evaluated at time t = 1), as ∆t tends to zero, of the backward Euler method
(BE) (θ = 1, solid line) and of the Crank-Nicolson scheme (CN) (θ = 1/2, dashed
line). As expected, the CN method is far more accurate than the BE method. •


’1
10


’2
10


’3
10


’4
10


’5
10


’6
10


’7
10 1 2
10 10


FIGURE 13.1. Convergence analysis of the θ-method as a function of the number
1/∆t of time steps (represented on the x-axis): θ = 1 (solid line) and θ = 0.5
(dashed line)




13.4 Space-Time Finite Element Methods for the
Heat Equation
An alternative approach for time discretization is based on the use of a
Galerkin method to discretize both space and time variables.
Suppose to solve the heat equation for x ∈ [0, 1] and t ∈ [0, T ]. Let us
denote by Ik = [tk’1 , tk ] the k-th time interval for k = 1, . . . , n with ∆tk =
tk ’ tk’1 ; moreover, we let ∆t = maxk ∆tk ; the rectangle Sk = [0, 1] — Ik is
the so called space-time slab. At each time level tk , we consider a partition
Thk of (0, 1) into mk subintervals Kj = [xk , xk ], j = 0, . . . , mk ’ 1. We
k
j j+1
let hk = xk ’ xk and denote by hk = maxj hk and by h = maxk hk .
j j+1 j j
Let us now associate with Sk a space-time partition Sk = ∪mk Rj where k
j=1
Rj = Kj — Ik and Kj ∈ Thk . The space-time slab Sk is thus decomposed
k k k
k
into rectangles Rj (see Figure 13.2).
For each time slab Sk we introduce the space-time ¬nite element space

Qq (Sk ) = v ∈ C 0 (Sk ), v|Rj ∈ P1 (Kj ) — Pq (Ik ), j = 0, . . . , mk ’ 1
k
k
594 13. Parabolic and Hyperbolic Initial Boundary Value Problems


t



tk
j
Rk Sk
k’1
t




1x
0
FIGURE 13.2. Space-time ¬nite element discretization

where, usually, q = 0 or q = 1. Then, the space-time ¬nite element space
over [0, 1] — [0, T ] is de¬ned as follows
Vh,∆t = v : [0, 1] — [0, T ] ’ R : v|Sk ∈ Yh,k , k = 1, . . . , n ,
where
Yh,k = {v ∈ Qq (Sk ) : v(0, t) = v(1, t) = 0 ∀t ∈ Ik } .
The number of degrees of freedom of Vh,∆t is equal to (q + 1)(mk ’ 1).
The functions in Vh,∆t are linear and continuous in space while they are
piecewise polynomials of degree q in time. These functions are in general
discontinuous across the time levels tk and the partitions Th do not match
k

at the interface between contiguous time levels (see Figure 13.2). For this
reason, we adopt henceforth the following notation
v± = lim v(tk ± „ ), [v k ] = v+ ’ v’ .
k k k
„ ’0

The discretization of problem (13.12) using continuous ¬nite elements in
space of degree 1 and discontinuous ¬nite elements of degree q in time
(abbreviated by cG(1)dG(q) method) is: ¬nd U ∈ Vh,∆t such that
n n’1
‚U
([U k ], V+ )
k
,V + a(U, V ) dt +
‚t
k=1I k=1
k
T
0
∀V ∈ V h,∆t ,
0 0
+(U+ , V+ ) = (f, V ) dt,
0

where
0
V h,∆t = {v ∈ Vh,∆t : v(0, t) = v(1, t) = 0 ∀t ∈ [0, T ]} ,
13.4 Space-Time Finite Element Methods for the Heat Equation 595

1
U’ = u0h , U k = U (x, tk ) and (u, v) = 0 uv dx denotes the scalar product
0

of L2 (0, 1). The continuity of U at each point tk is therefore imposed only
in a weak sense.
To construct the algebraic equations for the unknown U we need expand-
ing it over a basis in time and in space. The single space-time basis func-
tion •k (x, t) can be written as •k (x, t) = •k (x)ψl (t), j = 1, . . . , mk ’ 1,
j
jl jl
k
l = 0, . . . , q, where •j is the usual piecewise linear basis function and ψl is
the l-th basis function of Pq (Ik ).
When q = 0 the solution U is piecewise constant in time. In that case
k
Nh
Uj •k (x), x ∈ [0, 1], t ∈ Ik ,
U k (x, t) = k
j
j=1


where Uj = U k (xj , t) ∀t ∈ Ik . Let
k


Ak = (aij ) = (a(•k , •k )), Mk = (mij ) = ((•k , •k )),
j i j i
« 

fk = (fi ) =  f (x, t)•k (x)dx dt , Bk,k’1 = (bij ) = ((•k , •k’1 )),
i j i
Sk

denote the sti¬ness matrix, the mass matrix, the data vector and the pro-
k’1
and Vh , respectively, at the time level tk .
k
jection matrix between Vh
Then, letting Uk = (Uj ), at each k-th time level the cG(1)dG(0) method
k

requires solving the following linear system

Mk + ∆tk Ak Uk = Bk,k’1 Uk’1 + fk ,

which is nothing else than the Euler backward discretization scheme with
a modi¬ed right hand side.
When q = 1, the solution is piecewise linear in time. For ease of notation
we let U k (x) = U’ (x, tk ) and U k’1 (x) = U+ (x, tk’1 ). Moreover, we assume
that the spatial partition Thk does not change with the time level and we
let mk = m for every k = 0, . . . , n. Then, we can write

tk ’ t t ’ tk’1
k’1 k
U|Sk = U (x) + U (x) .
∆tk ∆tk
Thus the cG(1)dG(1) method leads to the solution of the following 2 —
2 block-system in the unknowns Uk = (Uik ) and Uk’1 = (Uik’1 ), i =
1, . . . , m ’ 1
± k k
 ’ 1 M + ∆t A Uk’1 + 1 M + ∆t A Uk = f
 k’1
 k’1 + Bk,k’1 U’ ,
k k k k
2 3 2 6
1 ∆tk ∆tk
 1
 Ak Uk’1 + Ak Uk = fk
Mk + Mk +
2 6 2 3
596 13. Parabolic and Hyperbolic Initial Boundary Value Problems

where

f (x, t)•k (x)ψ1 (t)dx dt, fk =
k
f (x, t)•k (x)ψ2 (t)dx dt
k
fk’1 = i i
Sk Sk


and ψ1 (t) = (tk ’ t)/∆tk , ψ2 (t)(t ’ tk’1 )/∆tk are the two basis functions
k k

of P1 (Ik ).

Assuming that Vh,k’1 ‚ Vh,k , it is possible to prove that (see for the
proof [EEHJ96])

u(tn ) ’ U n ¤ C(u0h , f, u, n)(∆t2 + h2 ), (13.25)
L2 (0,1)

where C is a constant that depends on its arguments (assuming that they
are su¬ciently smooth) but not on h and ∆t.
An advantage in using space-time ¬nite elements is the possibility to
perform a space-time grid adaptivity on each time-slab based on a posteriori
error estimates (the interested reader is referred to [EEHJ96] where the
analysis of this method is carried out in detail).

Program 101 provides an implementation of the dG(1)cG(1) method for
the solution of the heat equation on the space-time domain (a, b) — (t0 , T ).
The input parameters are the same as in Program 100.

Program 101 - pardg1cg1 : dG(1)cG(1) method for the heat equation
function [u,x]=pardg1cg1(I,n,u0,f,nu,bc)
nx = n(1); h = (I(2)-I(1))/nx;
x = [I(1):h:I(2)]; t = I(3); um = (eval(u0))™;
nt = n(2); k = (I(4)-I(3))/nt;
e = ones(nx+1,1);
Add = spdiags([(h/12-k*nu/(3*h))*e, (h/3+2*k*nu/(3*h))*e, ...
(h/12-k*nu/(3*h))*e],-1:1,nx+1,nx+1);
Aud = spdiags([(h/12-k*nu/(6*h))*e, (h/3+k*nu/(3*h))*e, ...
(h/12-k*nu/(6*h))*e],-1:1,nx+1,nx+1);
Ald = spdiags([(-h/12-k*nu/(6*h))*e, (-h/3+k*nu/(3*h))*e, ...
(-h/12-k*nu/(6*h))*e],-1:1,nx+1,nx+1);
B = spdiags([h*e/6, 2*h*e/3, h*e/6],-1:1,nx+1,nx+1);
Add(1,1) = 1; Add(1,2) = 0; B(1,1) = 0; B(1,2) = 0;
Aud(1,1) = 0; Aud(1,2) = 0; Ald(1,1) = 0; Ald(1,2) = 0;
Add(nx+1,nx+1)=1; Add(nx+1,nx)=0;
B(nx+1,nx+1)=0; B(nx+1,nx) = 0;
Ald(nx+1,nx+1)=0; Ald(nx+1,nx)=0;
Aud(nx+1,nx+1)=0; Aud(nx+1,nx)=0;
[L,U]=lu([Add Aud; Ald Add]);
x = [I(1)+h:h:I(2)-h]; xx=[I(1),x,I(2)];
for time = I(3)+k:k:I(4)
13.5 Hyperbolic Equations: A Scalar Transport Problem 597

t = time; fq1 = 0.5*k*h*eval(f);
t = time-k; fq0 = 0.5*k*h*eval(f);
rhs0 = [bc(1), fq0, bc(2)]; rhs1 = [bc(1), fq1, bc(2)];
b = [rhs0™; rhs1™] + [B*um; zeros(nx+1,1)];
y = L \ b; u = U \ y; um = u(nx+2:2*nx+2,1);
end
x = [I(1):h:I(2)]; u = um;


Example 13.2 We assess the accuracy of the dG(1)cG(1) method on the same
test problem considered in Example 13.1. In order to neatly identify both spatial
and temporal contributions in the error estimate (13.25) we have performed the
numerical computations using Program 101 by varying either the time step or
the space discretization step only, having chosen in each case the discretization
step in the other variable su¬ciently small that the corresponding error can be
neglected. The convergence behavior in Figure 13.3 shows perfect agreement with

the theoretical results (second-order accuracy in both space and time).


’1
10




’2
10




’3
10




’4
10




’5
10 1 2 3
10 10 10


FIGURE 13.3. Convergence analysis for the dG(1)cG(1) method. The solid line
is the time discretization error while the dashed line is the space discretization
error. In the ¬rst case the x-axis denotes the number of time steps while in second
case it represents the number of space subintervals




13.5 Hyperbolic Equations: A Scalar Transport
Problem
Let us consider the following scalar hyperbolic problem
±
 ‚u + a ‚u = 0 x ∈ R, t > 0,
‚t ‚x (13.26)

u(x, 0) = u0 (x) x ∈ R,
598 13. Parabolic and Hyperbolic Initial Boundary Value Problems

where a is a positive real number. Its solution is given by

u(x, t) = u0 (x ’ at) t ≥ 0,

and represents a travelling wave with velocity a. The curves (x(t), t) in the
plane (x, t), that satisfy the following scalar ordinary di¬erential equation
±
 dx(t) = a t > 0
dt (13.27)

x(0) = x0 ,

are called characteristic curves. They are the straight lines x(t) = x0 + at,
t > 0. The solution of (13.26) remains constant along them since
du ‚u ‚u dx
= + =0 on (x(t), t).
dt ‚t ‚x dt
For the more general problem
±
 ‚u + a ‚u + a u = f x ∈ R, t > 0,
0
‚t ‚x (13.28)

x ∈ R,
u(x, 0) = u0 (x)

where a, a0 and f are given functions of the variables (x, t), the charac-
teristic curves are still de¬ned as in (13.27). In this case, the solutions of
(13.28) satisfy along the characteristics the following di¬erential equation
du
= f ’ a0 u on (x(t), t).
dt


t
Q
t
t=1

¯
t P

P0 x x
0 1
± β

FIGURE 13.4. Left: examples of characteristics which are straight lines issuing
from the points P and Q. Right: characteristic straight lines for the Burgers
equation

Let us now consider problem (13.26) on a bounded interval. For example,
assume that x ∈ [±, β] and a > 0. Since u is constant along the character-
istics, from Figure 13.4 (left) we deduce that the value of the solution at P
13.6 Systems of Linear Hyperbolic Equations 599

attains the value of u0 at P0 , the foot of the characteristic issuing from P .
On the other hand, the characteristic issuing from Q intersects the straight
¯
line x(t) = ± at a certain time t = t > 0. Thus, the point x = ± is an
in¬‚ow point and it is necessary to assign a boundary value for u at x = ±
for every t > 0. Notice that if a < 0 then the in¬‚ow point is x = β.
Referring to problem (13.26) it is worth noting that if u0 is discontinuous
at a point x0 , then such a discontinuity propagates along the characteristics
issuing from x0 . This process can be rigorously formalized by introducing
the concept of weak solutions of hyperbolic problems, see e.g. [GR96]. An-
other reason for introducing weak solutions is that in the case of nonlinear
hyperbolic problems the characteristic lines can intersect: in this case the
solution cannot be continuous and no classical solution does exist.

Example 13.3 (Burgers equation) Let us consider the Burgers equation

‚u ‚u
x∈R
+u =0 (13.29)
‚t ‚x
which is perhaps the simplest nontrivial example of a nonlinear hyperbolic equa-
tion. Taking as initial condition
±
x ¤ 0,
1
1 ’ x 0 ¤ x ¤ 1,
u(x, 0) = u0 (x) =

x ≥ 1,
0

the characteristic line issuing from the point (x0 , 0) is given by
±
x0 ¤ 0,
 x0 + t
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 ∆,∞
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,

<<

. 24
( 26)



>>