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

x¯

(¯, t)

¯

t

x¯

(¯, 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

2π

1

u0 (x) e’ikx dx

±k =

2π

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

mπ

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

2π

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,