<< . .

. 84
( : 95)



. . >>

j n’1
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.


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 .

<< . .

. 84
( : 95)



. . >>