Th g(xj ) = h ° xj (1 ’ xk )»

xk (1 ’ xj ) +

k=1 k=j+1

from which, after straightforward computations, one gets the desired re-

sult.]

8. Prove Young™s inequality (12.40).

¤ vh ∀vh ∈ Vh .

9. Show that vh h h,∞

12.8 Exercises 579

10. Consider the two-point boundary value problem (12.29) with the following

boundary conditions

»0 u(0) + µ0 u(0) = g0 , »1 u(1) + µ1 u(1) = g1 ,

where »j , µj and gj are given data (j = 0, 1). Using the mirror imaging

technique described in Section 12.2.3 write the ¬nite di¬erence discretiza-

tion of the equations corresponding to the nodes x0 and xn .

[Solution:

±1/2 γ0 ± 0 »0 u1 ±0 g0 f0

u0 ’ ±1/2 2 =

node x0 : + + +,

2

h 2 µ0 h h µ0 h 2

±n’1/2 γn ±n »1 un’1 ±n g1 fn

un ’ ±n’1/2 2 =

node xn : + + + .]

2

h 2 µ1 h h µ1 h 2

11. Discretize the fourth-order di¬erential operator Lu(x) = ’u(iv) (x) using

centred ¬nite di¬erences.

[Solution : apply twice the second order centred ¬nite di¬erence operator

Lh de¬ned in (12.9).]

12. Consider problem (12.41) with non homogeneous Neumann boundary con-

ditions ±u (0) = w0 , ±u (1) = w1 . Show that the solution satis¬es prob-

lem (12.43) where V = H1 (0, 1) and the right-hand side is replaced by

(f, v)+w1 v(1)’w0 v(0). Derive the formulation in the case of mixed bound-

ary conditions ±u (0) = w0 , u(1) = u1 .

13. Using Property 1.19 prove that the matrices corresponding to the stabilized

¬nite element method (12.79) using the upwind and SG arti¬cial viscosities

φU P and φSG (see Section 12.5.2) are M-matrices irrespective of h.

[Hint: let us denote respectively by AU P and ASG the two sti¬ness matrices

corrsponding to φU P and φSG . Take v(x) = 1 + x and set vi = 1 + xi ,

i = 0, . . . , n, being xi = ih, h = 1/n. Then, by a direct computation check

that (AU P v)i ≥ β > 0. As for the matrix ASG the same result can be

proved by noting that B(’t) = t + B(t) for any t ∈ R.]

14. Prove that the matrix Afd with entries given by (12.93) is symmetric pos-

itive de¬nite and it is also an M-matrix.

[Solution: to show that Afd is positive de¬nite, proceed as in the corre-

sponding proof in Section 12.2, then proceed as in Exercise 2.]

15. Prove the Green™s formula (12.95).

[Solution: ¬rst, notice that for any u, v su¬ciently smooth, div(v∇u) =

v u + ∇u · ∇v. Then, integrate this relation over „¦ and use the divergence

‚u

theorem div(v∇u) dxdy = v dγ.]

‚n

„¦ ‚„¦

13

Parabolic and Hyperbolic Initial

Boundary Value Problems

The ¬nal chapter of this book is devoted to the approximation of time-

dependent partial di¬erential equations. Parabolic and hyperbolic initial-

boundary value problems will be addressed and either ¬nite di¬erences and

¬nite elements will be considered for their discretization.

13.1 The Heat Equation

The problem we are considering is how to ¬nd a function u = u(x, t) for

x ∈ [0, 1] and t > 0 that satis¬es the partial di¬erential equation

‚u

(13.1)

+ Lu = f 0 < x < 1, t > 0,

‚t

subject to the boundary conditions

u(0, t) = u(1, t) = 0, t>0 (13.2)

and the initial condition

0 ¤ x ¤ 1.

u(x, 0) = u0 (x) (13.3)

The di¬erential operator L is de¬ned as

‚2u

Lu = ’ν . (13.4)

‚x2

582 13. Parabolic and Hyperbolic Initial Boundary Value Problems

Equation (13.1) is called the heat equation. In fact, u(x, t) describes the

temperature at the point x and time t of a metallic bar of unit length

that occupies the interval [0, 1], under the following conditions. Its thermal

conductivity is constant and equal to ν > 0, its extrema are kept at a

constant temperature of zero degrees, at time t = 0 its temperature at

point x is described by u0 (x), and f (x, t) represents the heat production

per unit length supplied at point x at time t. Here we are supposing that

the volumetric density ρ and the speci¬c heat per unit mass cp are both

constant and unitary. Otherwise, the temporal derivative ‚u/‚t should be

multiplied by the product ρcp in (13.1).

A solution of problem (13.1)-(13.3) is provided by a Fourier series. For

instance, if ν = 1 and f ≡ 0, it is given by

∞

2

cn e’(nπ) t sin(nπx)

u(x, t) = (13.5)

n=1

where the coe¬cients cn are the Fourier sine coe¬cients of the initial datum

u0 (x), i.e.

1

cn = 2 u0 (x) sin(nπx) dx, n = 1, 2 . . .

0

If instead of (13.2) we consider the Neumann conditions

ux (0, t) = ux (1, t) = 0, t > 0, (13.6)

the corresponding solution (still in the case where f ≡ 0 and ν = 1) would

be

∞

d0 2

dn e’(nπ) t cos(nπx),

u(x, t) = +

2 n=1

where the coe¬cients dn are the Fourier cosine coe¬cients of u0 (x), i.e.

1

dn = 2 u0 (x) cos(nπx) dx, n = 1, 2 . . .

0

These expressions show that the solution decays exponentially fast in time.

A more general result can be stated concerning the behavior in time of the

energy

1

u2 (x, t) dx.

E(t) =

0

13.1 The Heat Equation 583

Indeed, if we multiply (13.1) by u and integrate with respect to x over the

interval [0, 1], we obtain

1 1 1

‚2u ‚u2

‚u 1

(x, t)u(x, t) dx ’ ν (x, t)u(x, t) dx = (x, t) dx

‚x2

‚t 2 ‚t

0 0 0

1

2 x=1

‚u ‚u

(x, t) dx ’ ν

+ν (x, t)u(x, t)

‚x ‚x x=0

0

1

2

1 ‚u

= E (t) + ν (x, t) dx,

2 ‚x

0

having used integration by parts, the boundary conditions (13.2) or (13.6),

and interchanged di¬erentiation and integration.

Using the Cauchy-Schwarz inequality (8.29) yields

1

f (x, t)u(x, t) dx ¤ (F (t))1/2 (E(t))1/2

0

1

f 2 (x, t) dx. Then

where F (t) = 0

1

2

‚u

dx ¤ 2(F (t))1/2 (E(t))1/2 .

E (t) + 2ν (x, t)

‚x

0

Owing to the Poincar´ inequality (12.16) with (a, b) = (0, 1) we obtain

e

ν

E(t) ¤ 2(F (t))1/2 (E(t))1/2 .

E (t) + 2 2

(CP )

By Young™s inequality (12.40) we have

1

2(F (t))1/2 (E(t))1/2 ¤ γE(t) + F (t),

γ

having set γ = ν/CP . Therefore, E (t) + γE(t) ¤ γ F (t), or, equivalently,

1

2

(eγt E(t)) ¤ γ eγt F (t). Then, integrating from 0 to t we get

1

t

1

E(t) ¤ e’γt E(0) + eγ(s’t) F (s)ds. (13.7)

γ

0

In particular, when f ≡ 0, (13.7) shows that the energy E(t) decays expo-

nentially fast in time.

584 13. Parabolic and Hyperbolic Initial Boundary Value Problems

13.2 Finite Di¬erence Approximation of the Heat

Equation

To solve the heat equation numerically we have to discretize both the x

and t variables. We can start by dealing with the x-variable, following the

same approach as in Section 12.2. We denote by ui (t) an approximation of

u(xi , t), i = 0, . . . , n and approximate the Dirichlet problem (13.1)-(13.3)

by the scheme

. (t) ’ ν (u i’1 (t) ’ 2ui (t) + ui+1 (t)) = fi (t), i = 1, . . . , n ’ 1, ∀t > 0,

u i

h2

∀t > 0,

u0 (t) = un (t) = 0,

ui (0) = u0 (xi ), i = 0, . . . , n,

where the upper dot indicates derivation with respect to time, and fi (t) =

f (xi , t). This is actually a semi-discretization of problem (13.1)-(13.3), and

is a system of ordinary di¬erential equations of the following form

u(t) = ’νAfd u(t) + f (t), ∀t > 0,

™

(13.8)

u(0) = u0

where u(t) = (u1 (t), . . . , un’1 (t))T is the vector of unknowns, f (t) =

(f1 (t), . . . , fn’1 (t))T , u0 = (u0 (x1 ), . . . , u0 (xn’1 ))T and Afd is the tridi-

agonal matrix introduced in (12.8). Note that for the derivation of (13.8)

we have assumed that u0 (x0 ) = u0 (xn ) = 0, which is coherent with the

boundary condition (13.2).

A popular scheme for the integration of (13.8) with respect to time is the

so-called θ’method. To construct the scheme, we denote by v k the value

of the variable v at time tk = k∆t, for ∆t > 0; then, the θ-method for the

time-integration of (13.8) is

± k+1

’ uk

u

= ’νAfd (θuk+1 + (1 ’ θ)uk ) + θf k+1 + (1 ’ θ)f k ,

∆t

(13.9)

k = 0, 1, . . .

0

u = u0

or, equivalently,

(I + νθ∆tAfd ) uk+1 = (I ’ ν∆t(1 ’ θ)Afd ) uk + gk+1 , (13.10)

where gk+1 = ∆t(θf k+1 + (1 ’ θ)f k ) and I is the identity matrix of order

n ’ 1.

For suitable values of the parameter θ, from (13.10) we can recover some

familiar methods that have been introduced in Chapter 11. For example, if

θ = 0 the method (13.10) coincides with the forward Euler scheme and we

13.2 Finite Di¬erence Approximation of the Heat Equation 585

can get uk+1 explicitly; otherwise, a linear system (with constant matrix

I + νθ∆tAfd ) needs be solved at each time-step.

Regarding stability, assume that f ≡ 0 (henceforth gk = 0 ∀k > 0),

so that from (13.5) the exact solution u(x, t) tends to zero for every x

as t ’ ∞. Then we would expect the discrete solution to have the same

behaviour, in which case we would call our scheme (13.10) asymptotically

stable, this being coherent with what we did in Chapter 11, Section 11.1

for ordinary di¬erential equations.

If θ = 0, from (13.10) it follows that

uk = (I ’ ν∆tAfd )k u0 , k = 1, 2, . . .

From the analysis of convergent matrices (see Section 1.11.2) we deduce

that uk ’ 0 as k ’ ∞ i¬

ρ(I ’ ν∆tAfd ) < 1. (13.11)

On the other hand, the eigenvalues of Afd are given by (see Exercise 3)

4

i = 1, . . . , n ’ 1.

sin2 (iπh/2),

µi = 2

h

Then (13.11) is satis¬ed i¬

12

∆t < h.

2ν

As expected, the forward Euler method is conditionally stable, and the

time-step ∆t should decay as the square of the grid spacing h.

In the case of the backward Euler method (θ = 1), we would have from

(13.10)

k

uk = (I + ν∆tAfd )’1 u0 , k = 1, 2, . . .

Since all the eigenvalues of the matrix (I + ν∆tAfd )’1 are real, positive and

strictly less than 1 for every value of ∆t, this scheme is unconditionally

stable. More generally, the θ-scheme is unconditionally stable for all the

values 1/2 ¤ θ ¤ 1, and conditionally stable if 0 ¤ θ < 1/2 (see Section

13.3.1).

As far as the accuracy of the θ-method is concerned, its local truncation

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

2

θ = 1 . The method corresponding to θ = 1/2 is frequently called the Crank-

2

Nicolson scheme and is therefore unconditionally stable and second-order

accurate with respect to both ∆t and h.

586 13. Parabolic and Hyperbolic Initial Boundary Value Problems

13.3 Finite Element Approximation of the Heat

Equation

The space discretization of (13.1) can also be accomplished using the Galerkin

¬nite element method by proceeding as in Chapter 12 in the elliptic case.

First, for all t > 0 we multiply (13.1) by a test function v = v(x) and

integrate over (0, 1). Then, we let V = H1 (0, 1) and ∀t > 0 we look for a

0

function t ’ u(x, t) ∈ V (brie¬‚y, u(t) ∈ V ) such that

1

‚u(t)

∀v ∈ V, (13.12)

vdx + a(u(t), v) = F (v)

‚t

0

1

with u(0) = u0 . Here, a(u(t), v) = 0 ν(‚u(t)/‚x) (‚v/‚x) dx and F (v) =

1

f (t)vdx are the bilinear form and the linear functional respectively as-

0

sociated with the elliptic operator L and the right hand side f . Notice that

a(·, ·) is a special case of (12.44) and that the dependence of u and f on

the space variable x will be understood henceforth.

Let Vh be a suitable ¬nite dimensional subspace of V . We consider the

following Galerkin formulation: ∀t > 0, ¬nd uh (t) ∈ Vh such that

1

‚uh (t)

vh dx + a(uh (t), vh ) = F (vh ) ∀vh ∈ Vh (13.13)

‚t

0

where uh (0) = u0 h and u0 h ∈ Vh is a convenient approximation of u0 .

Problem (13.13) is referred to as a semi-discretization of (13.12) since it is

only a space discretization of the heat equation.

Proceeding in a manner similar to that used to obtain the energy estimate

(13.7), we get the following a priori estimate for the discrete solution uh (t)

of (13.13)

t

1

Eh (t) ¤ e’γt Eh (0) + eγ(s’t) F (s)ds,

γ

0

1

where Eh (t) = 0 u2 (x, t) dx.

h

As for the ¬nite element discretization of (13.13), we introduce the ¬nite

element space Vh de¬ned in (12.57) and consequently a basis {•j } for Vh

as already done in Section 12.4.5. Then, the solution uh of (13.13) can be

sought under the form

Nh

uh (t) = uj (t)•j ,

j=1

13.3 Finite Element Approximation of the Heat Equation 587

where {uj (t)} are the unknown coe¬cients and Nh is the dimension of Vh .

Then, from (13.13) we obtain

«

. (t)• • dx + a

1 Nh Nh

uj (t)•j , •i = F (•i ),

u i = 1, . . . , Nh

j j i

0 j=1 j=1

that is,

. (t) 1

Nh Nh

u •j •i dx + uj (t)a(•j , •i ) = F (•i ), i = 1, . . . , Nh .