<<

. 12
( 16)



>>

Assume that additionally to our assumptions the bilinear form is symmetric
and the embedding of V into H is compact (see [26]), which is the case here.
Then there are enough eigenvectors in the sense that a sequence (wi , »i ),
292 7. Discretization of Parabolic Problems

0 < »1 ¤ »2 ¤ . . . , exists such that the wi are orthonormal with respect
to ·, · 0 and every v ∈ V has a unique representation (in H) as

v= ci wi . (7.33)
i=1

As in (7.25) the Fourier coe¬cients ci are given by
ci = v, wi . (7.34)
0

In fact, (7.33) gives a rigorous framework to the speci¬c considerations
in (7.16) and subsequent formulas. From (7.33) and (7.34) we conclude
Parseval™s identity

| v, wi 0 |2 .
2
v = (7.35)
0
i=1
’1/2
Furthermore, the sequence vi := »i wi is orthogonal with respect to
a(·, ·), and a representation corresponding to (7.33), (7.34) holds such that
∞ ∞ ∞
»’1 |a(v, wi )|2
|a(v, vi )| = »i | v, wi 0 |2 .
2
a(v, v) = = (7.36)
i
i=1 i=1 i=1

From (7.35) and (7.36) we see that the ellipticity constant can be inter-
preted as the smallest eigenvalue ». In fact, the solution representation
(7.18) (with the sum being in¬nite in H) also holds true under the as-
sumptions mentioned and also leads to the estimate of Lemma 7.4. But
note that the proof there does not require symmetry of the bilinear form.



Exercises
7.1 Consider the initial-boundary value problem
ut ’ uxx = in (0, ∞) — (0, ∞) ,
0
h(t) , t ∈ (0, ∞) ,
u(0, t) =
x ∈ (0, ∞) ,
u(x, 0) = 0,
where h : (0, ∞) ’ R is a di¬erentiable function, the derivative of which
has at most exponential growth.
(a) Show that the function

x2
2 ’s2 /2
h t’ 2
u(x, t) = e ds

π 2s
x/ 2t

is a solution.
(b) Is ut bounded in the domain of de¬nition? If not, give conditions on
h that guarantee the boundedness of ut .
7.2. Semidiscretization by the Vertical Method of Lines 293



7.2 Consider the initial-boundary value problem in one space dimension
ut ’ uxx = in (0, π) — (0, ∞) ,
0
t ∈ (0, ∞) ,
u(0, t) = u(π, t) = 0,
x ∈ (0, π) .
u(x, 0) = u0 (x) ,
(a) Solve it by means of the method of separation.
(b) Give a representation for ut (t) 0 .
(c) Consider the particular initial condition u0 (x) = π’x and investigate,
using the result from subproblem (b), the asymptotic behaviour of
ut (t) 0 near t = 0.

7.3 Let the domain „¦ ‚ Rd be bounded by a su¬ciently smooth boundary
and set V := H0 („¦), H := L2 („¦). Furthermore, a : V — V ’ R is a
1

continuous, V -elliptic, symmetric bilinear form and u0 ∈ H. Prove by using
the so-called energy method (cf. the proof of Lemma 7.4) the following a
priori estimate for the solution u of the initial boundary value problem
+ a(u(t), v) = 0 for all v ∈ V, t ∈ (0, T ) ,
ut (t), v 0
u(0) = u0 .
t t
ds ¤ M
2 2 2
(a) ±t u(t) +2 s ut (s) u(s) 1 ds .
1 0
0 0
M1
¤
(b) ut (t) u0 0 .
0
2± t
Here M and ± denote the corresponding constants in the continuity and
ellipticity conditions, respectively.



7.2 Semidiscretization
by the Vertical Method of Lines
For solving parabolic equations numerically, a wide variety of methods
exists. The most important classes of these methods are the following:
• Full discretizations:
“ Application of ¬nite di¬erence methods to the classical initial
boundary value problem (as of the form (7.1)).
“ Application of so-called space-time ¬nite element methods to a
variational formulation that includes the time variable, too.
• Semidiscretizations:
“ The vertical method of lines: Here the discretization starts with
respect to the spatial variable(s) (e.g., by means of the ¬nite dif-
294 7. Discretization of Parabolic Problems

ference method, the ¬nite element method, or the ¬nite volume
method).
“ The horizontal method of lines (Rothe™s method): Here the
discretization starts with respect to the time variable.

As the name indicates, a semidiscretization has to be followed by a further
discretization step to obtain a full discretization, which may be one of
the above-mentioned or not. The idea behind semidiscretization methods
is to have intermediate problems that are of a well-known structure. In
the case of the vertical method of lines, a system of ordinary di¬erential
equations arises for the solution of which appropriate solvers are often
available. Rothe™s method generates a sequence of elliptic boundary value
problems for which e¬cient solution methods are known, too.
The attributes “vertical” and “horizontal” of the semidiscretizations are
motivated by the graphical representation of the domain of de¬nition of the
unknown function u = u(x, t) in one space dimension (i.e., d = 1), namely,
assigning the abscissa (horizontal axis) of the coordinate system to the
variable x and the ordinate (vertical axis) to the variable t, so that the
spatial discretization yields problems that are setted along vertical lines.
In what follows, the vertical method of lines will be considered in more
detail.
In the following, and similarly in the following sections, we will de-
velop the analogous (semi)discretization approaches for the ¬nite di¬erence
method, the ¬nite element method, and the ¬nite volume method. This
will allow us to analyze these methods in a uniform way, as far as only the
emerging (matrix) structure of the discrete problems will play a role. On
the other hand, di¬erent techniques of analysis as in Chapters 1, 3 and 6
will further elucidate advantages and disadvantages of the methods. Read-
ers who are interested only in a speci¬c approach may skip some of the
following subsections.

The Vertical Method of Lines for the Finite Di¬erence Method
As a ¬rst example we start with the heat equation (7.8) with Dirichlet
boundary conditions on a rectangle „¦ = (0, a) — (0, b). As in Sec-
tion 1.2 we apply the ¬ve-point stencil discretizations at the grid points
x ∈ „¦h (according to (1.5)) for every ¬xed t ∈ [0, T ]. This leads to the
approximation

1
’ ui,j’1 (t) ’ ui’1,j (t) + 4uij (t) ’ ui+1,j (t) ’ ui,j+1 (t)
‚t uij (t) +
h2
i = 1, . . . , l ’ 1, j = 1, . . . , m ’ 1, t ∈ (0, T ) ,
= fij (t) ,
(7.37)

i ∈ {0, l}, j = 0, . . . , m ,
uij (t) = gij (t),
(7.38)
j ∈ {0, m}, i = 0, . . . , l .
7.2. Semidiscretization by Vertical Method of Lines 295

Here we use
fij (t) := f (ih, jh, t) ,
(7.39)
gij (t) := g(ih, jh, t) ,
and the index 3 in the boundary condition is omitted. Additionally, the
initial condition (at the grid points) will be prescribed, that is,
uij (0) = u0 (ih, jh), (ih, jh) ∈ „¦h . (7.40)
The system (7.37), (7.38), (7.40) is a system of (linear) ordinary di¬erential
equations (in the “index” (i, j)). If, as in Section 1.2, we ¬x an ordering of
the grid points, the system takes the form
d
t ∈ (0, T ) ,
uh (t) + Ah uh (t) = q h (t) ,
(7.41)
dt
uh (0) = u0 ,
with Ah , q h as in (1.10), (1.11) (but now q h = q h (t) because of the t-
dependence of f and g).
The unknown is the function
uh : [0, T ] ’ RM1 , (7.42)
which means that the Dirichlet boundary conditions are eliminated as in
Section 1.2.
For a simpli¬cation of the notation we use in the following M instead
of M1 , which also includes the eliminated degrees of freedom. Only in
Sections 7.5 and 7.6 will we return to the original notation.
More generally, if we consider a ¬nite di¬erence approximation, which
applied to the stationary problem (7.6) will lead to the system of equations
Ah uh = q h ,
with uh ∈ RM , then the same method applied to (7.1) for every ¬xed
t ∈ (0, T ) leads to (7.41). In particular, the system (7.41) has a unique
solution due to the theorem of Picard“Lindel¨f (cf. [26]).
o

The Vertical Method of Lines for the Finite Element Method
We proceed as for the ¬nite di¬erence method by now applying the ¬nite
element method to (7.1) in its weak formulation (7.32) for every ¬xed t ∈
(0, T ), using the abbrevation

b(t, v) := f (t), v + g1 (·, t)v dσ + g2 (·, t)v dσ . (7.43)
0
“1 “2

So let Vh ‚ V denote a ¬nite-dimensional subspace with dim Vh = M =
M (h) and let u0h ∈ Vh be some approximation to u0 . Then the semidiscrete
problem reads as follows:
296 7. Discretization of Parabolic Problems

Find uh ∈ L2 ((0, T ), Vh ) with uh ∈ L2 ((0, T ), H) , uh (0) = u0h and
d
+a(uh (t), vh ) = b(t, vh ) for all vh ∈ Vh , t ∈ (0, T ) . (7.44)
uh (t), vh
dt 0

To gain a more speci¬c form of (7.44), again we represent the unknown
uh (t) by its degrees of freedom:
M
Let {•i }M be a basis of Vh , uh (t) = i=1 ξi (t) •i and u0h =
i=1
M
i=1 ξ0i •i . Then for any t ∈ (0, T ), the discrete variational equality (7.44)
is equivalent to
M M
dξj (t)
a(•j , •i ) ξj (t) = b(t, •i ) for all i ∈ {1, . . . , M } .
•j , •i +
0 dt
j=1 j=1

ˆ
Denoting by Ah := (a(•j , •i ))ij the sti¬ness matrix , by Bh := •j , •i 0 ij
the mass matrix, and by
β h (t) := (b(t, •i ))i ,
respectively ξ 0h := (ξ0i )i , the vectors of the right-hand side and of the
initial value, we obtain for ξ h (t) := (ξi (t))i the following system of linear
ordinary di¬erential equations with constant coe¬cients:
d
t ∈ (0, T ) ,
ˆ
Bh ξ (t) + Ah ξ h (t) = β h (t) ,
dt h (7.45)
ξ h (0) = ξ 0h .
Since the matrix Bh is symmetric and positive de¬nite, it can be factored
T
(e.g., by means of Cholesky™s decomposition) as Bh = Eh Eh . Introducing
the new variable uh := Eh ξ h (to maintain the possible de¬niteness of Ah ),
the above system (7.45) can be written as follows:
d
t ∈ (0, T ) ,
uh (t) + Ah uh (t) = q h (t) ,
(7.46)
dt
uh (0) = uh0 ,
’T ˆ ’1 ’T
where Ah := Eh Ah Eh is an RM -elliptic matrix and q h := Eh β h ,
uh0 := Eh ξ 0h .
Thus again the discretization leads us to a system (7.41).
Remark 7.7 By means of the same arguments as in the proof of
Lemma 7.4, an estimate of uh (t) 0 can be derived.

The Vertical Method of Lines for the Finite Volume Method
Based on the ¬nite volume methods introduced in Chapter 6, in this sub-
section a ¬nite volume semidiscretization is given for the problem (7.1)
in its weak formulation (7.32) for every ¬xed t ∈ (0, T ) in the special
case “3 = ‚„¦ and of homogeneous Dirichlet boundary conditions. As in
Chapter 6, the only essential di¬erence to problem (7.1) is that here the
7.2. Semidiscretization by Vertical Method of Lines 297

di¬erential expression L is in divergence form, i.e.,
Lu := ’∇ · (K ∇u ’ c u) + r u = f ,
where the data K, c, r, and f are as in (7.2).
Correspondingly, the bilinear form a in the weak formulation (7.32) is to
be replaced by

[(K ∇u ’ c u) · ∇v + ruv] dx .
a(u, v) = (7.47)
„¦

In order to obtain a ¬nite volume semidiscretization of the problem (7.1)
in divergence form, and of (7.32) with the modi¬cation (7.47), we recall
the way that it was done in the elliptic situation. Namely, comparing the
weak formulation of the elliptic problem (see De¬nition 2.2) with the ¬nite
volume method in the discrete variational formulation (6.11), we see that
the bilinear form a and the linear form b(·) := f, · 0 have been replaced by
certain discrete forms ah and f, · 0,h , respectively. This formal procedure
can be applied to the weak formulation (7.32) of the parabolic problem,
too.
So let Vh ‚ V denote a ¬nite-dimensional subspace as introduced in Sec-
tion 6.2 with dim Vh = M = M (h) and let u0h ∈ Vh be some approximation
to u0 . Then, the semidiscrete ¬nite volume method reads as follows:
Find uh ∈ L2 ((0, T ), Vh ) with uh ∈ L2 ((0, T ), H) , uh (0) = u0h and
d
for all vh ∈ Vh , t ∈ (0, T ) ,
uh (t), vh +ah (uh (t), vh ) = f (t), vh 0,h
dt 0,h
(7.48)
where both the bilinear form ah and the form ·, · 0,h have been formally
de¬ned in Section 6.2. However, to facilitate the comparison of the ¬nite
volume discretization with the previously described methods, here we set
Λ := {1, . . . , M }.
As in Section 6.2 we consider the following discrete L2 -scalar product
·, · 0,h :
M
d d
uh (t), vh = uh (aj , t)vh (aj )mj . (7.49)
dt dt
0,h
j=1

In analogy to the case of the ¬nite element method (cf. Remark 7.7), a
stability estimate for the ¬nite volume method can be obtained. Namely,
under the assumptions of Theorem 6.15, we have that
ah (vh , vh ) ≥ ± vh for all vh ∈ Vh
2
0,h

with some constant ± > 0 independent of h. Then, taking vh = uh (t) in
(7.48), we get
d
uh (t), uh (t) + ah (uh (t), uh (t)) = f (t), uh (t) ,
0,h
dt 0,h
298 7. Discretization of Parabolic Problems

and, after some calculations,
d
¤ f (t)
uh (t) + ± uh (t) 0,h .
0,h 0,h
dt
The subsequent arguments are as in the proof of Lemma 7.4; i.e., we obtain
t
’±t ’±(t’s)
¤ u0
uh (t) 0,h e + f (s) 0,h e ds .
0,h
0

If the right-hand side of (7.48) is a general bounded linear form, i.e., instead
of f (t), vh 0,h we have the term b(t, vh ), where b : (0, T ) — Vh ’ R is such
that
|b(t, v)| ¤ b(t) for all v ∈ Vh , t ∈ (0, T ),
v
— 0,h

< ∞ for all t ∈ (0, T ), then an analogous estimate holds:
with b(t) —
t
’±t
e’±(t’s) ds .
¤ u0
uh (t) 0,h e + b(s) (7.50)

0,h
0

As in the previous subsection, we now want to give a more speci¬c form of
(7.48).
Given a basis {•i }M of the space Vh , such that •i (aj ) = δij for the
i=1
underlying nodes, we have the unique expansions
M M
uh (t) = ξi (t) •i and u0h = ξ0i •i .
i=1 i=1

Then for any t ∈ (0, T ), the discrete variational equality (7.48) is equivalent
to
M
dξi (t)
for all i ∈ {1, . . . , M } ,
mi + ah (•j , •i ) ξj (t) = f (t), •i 0,h
dt j=1

where mi = |„¦i |. Using the notation Ah := (ah (•j , •i ))ij for the ¬nite
ˆ
volume sti¬ness matrix, Bh := diag(mi ) for the ¬nite volume mass matrix,
βh (t) := f (t), •i for the vector of the right-hand side, and ξ 0h :=
0,h
i
(ξ0i )i for the vector of the initial value, we obtain for the unknown vector
function ξ h (t) := (ξi (t))i the following system of linear ordinary di¬erential
equations with constant coe¬cients:
d
t ∈ (0, T ) ,
ˆ
Bh ξ (t) + Ah ξ h (t) = β h (t) ,
dt h (7.51)
ξ h (0) = ξ 0h .
In contrast to the system (7.45) arising in the ¬nite element semidiscretiza-
tion, here the matrix Bh is diagonal. Therefore, it is very easy to introduce

the new variable uh := Eh ξ with Eh := diag( mi ), and the above system
7.2. Semidiscretization by Vertical Method of Lines 299

(7.51) can be written as follows:
d
t ∈ (0, T ) ,
uh (t) + Ah uh (t) = q h (t) ,
(7.52)
dt
uh (0) = uh0 ,
’1 ˆ ’1 ’1
where Ah := Eh Ah Eh is an RM -elliptic matrix and q h := Eh β h , uh0 :=
Eh ξ 0h .
Thus again we have arrived at a system of the form (7.41).

Representation of Solutions in a Special Case
The solution of system (7.41) can be represented explicitly if there is a
basis of RM composed of eigenvectors of Ah . This will be developed in the
following, but is not meant for numerical use, since only in special cases
can eigenvectors and values be given explicitly. Rather, it will serve as a
tool for comparison with the continuous and the fully discrete cases.
Let (wi , »i ), i = 1, . . . , M, be the eigenvectors and real eigenvalues of
Ah . Then the following representation exists uniquely:
M M
i
u0 = ci w i and q h (t) = qh (t)w i . (7.53)
i=1 i=1

Again by a separation of variables approach (cf. (7.18)) we see that
t
M
(ci e’»i t + qh (t)e’»i (t’s) ds)wi .
i
uh (t) = (7.54)
i=1 0

A more compact notation is given by
t
’Ah t
e’Ah (t’s) q h (s)ds
uh (t) = e u0 + (7.55)
0

if we de¬ne for a matrix B ∈ RM,M ,


B
e := .
ν!
ν=0

This can be seen as follows:
Let
T := (w 1 , . . . , wM ) ∈ RM,M ,

Λ := diag(»i ) ∈ RM,M .
Then

T ’1 Ah T = Λ ,
Ah T = T Λ , i.e.,
300 7. Discretization of Parabolic Problems

and therefore
∞ ∞
tν ’1 tν
’1 ’Ah t ν
(’Λ)ν ,
T e T= T (’Ah ) T =
ν! ν!
ν=0 ν=0

since T ’1 (’Ah )ν T = T ’1 (’Ah )T T ’1 (’Ah )T T ’1 . . . T and thus

(’»i t)ν
’1 ’Ah t
= diag e’»i t .
T e T = diag
ν!
ν=0

Then for c = (c1 , . . . , cM )T ∈ RM , because of T c = u0 we conclude for the
case q h = 0 that
uh (t) = T diag(e’»i t )c = T T ’1e’Ah t T c = e’Ah t u0 ,
and similarly in general.
A basis of eigenvalues exists if Ah is self-adjoint with respect to a scalar
product ·, · h in RM , meaning that
for all u, v ∈ RM .
v, Ah u = Ah v, u
h h

Then the eigenvectors even are orthogonal ; that is,
w i , wj =0 for i = j (7.56)
h

because of
»i wi , wj = Ah wi , wj = w i , Ah w j = »j w i , w j ,
h h

and thus (7.56) if »i = »j . But eigenvectors belonging to one eigenvalue
can always be orthonormalized. For orthogonal w i the coe¬cient ci from
(7.53) has the form
u0 , w i h
ci = , (7.57)
wi, wi h
i
and analogously for qh .

Order of Convergence Estimates for the Finite Di¬erence Method
in a Special Case
As an illustrative example, we consider a case where the eigenvectors and
eigenvalues of Ah are known explicitly: the ¬ve-point stencil discretization
of the Poisson equation with Dirichlet conditions in „¦ = (0, a) — (0, b).
Instead of considering a ¬xed ordering of the grid points, we prefer to use
the “natural” two-dimensional indexing; i.e., we regard the eigenvectors as
grid functions. As seen in Section 1.2, Ah is symmetric and thus self-adjoint
with respect to the Euclidean scalar product scaled with hd if „¦ ‚ Rd , i.e.,
d = 2 here:
M
d
u, v =h ui vi . (7.58)
h
i=1
7.2. Semidiscretization by Vertical Method of Lines 301

The norm induced by this scalar product is exactly the discrete L2 -norm
de¬ned in (1.18) (for d = 2 and for the vectors representing the grid func-
tions):
1/2
M
1/2
|u|0,h = u, u |ui |2
= hd/2 . (7.59)
h
i=1

If we mean the grid function U, we denote the norm by U 0,h .
The eigenvectors, which have already been noted for a special case after
Theorem 5.4, are written as grid functions
π π
for (x, y) ∈ „¦h ,
uνµ (x, y) = sin ν x sin µ y
(7.60)
a b
and ν = 1, . . . , l ’ 1, µ = 1, . . . , m ’ 1
for the eigenvalues
2 π π
»νµ = 2 ’ cos ν h ’ cos µ h .
h
h2 a b
Note that the eigenvectors are the eigenfunctions of the continuous problem
evaluated at the grid points, but the grid points can distinguish only the
maximal frequencies l’1 and m’1 , so that for other indices ν, µ the given
2 2
grid functions would be repeated.
Due to 2 sin2 2 = 1 ’ cos(ξ), an alternative representation is
ξ


4 πh πh
»νµ = sin2 ν + sin2 µ ,
h
h2 a2 b2
so that for h ’ 0,
2
2
νπ πh πh
»νµ = sin ν ν
h
a a2 a2
2
µπ 2 πh πh (7.61)
+ sin µ µ
b b2 b2
νπ 2 µπ 2
’ cos2 (0) + cos2 (0)
a b
holds; i.e., the eigenvalues converge to the eigenvalues (7.27) of the
boundary value problem, with an order of convergence estimate of O(h2 ).
The eigenvectors are orthogonal with respect to ·, · h , since they be-
long to di¬erent eigenvalues (see (7.56)). To specify the Fourier coe¬cients
according to (7.57), we need
ab
uνµ , uνµ = (7.62)
h
4
(see Exercise 7.5).
To investigate the accuracy of the semidiscrete approximation, the so-
lution representations can be compared. To simplify the exposition, we
302 7. Discretization of Parabolic Problems

consider only f = 0, so that because of (7.18), (7.27) we have

π π
cνµ e’»
νµ
t
u(x, y, t) = sin ν x sin µ y ,
a b
ν=1
µ=1


and
b a
4 π π
cνµ = u0 (x, y) sin ν x sin µ y dx dy
ab a b
0 0

because of (7.25) and (7.24) (applied in every space direction), and ¬nally,
νπ 2 µπ 2
»νµ =
+
a b
for the continuous solution. For the semidiscrete approximation at a grid
point (x, y) ∈ „¦h we have, due to (7.54),
l’1 m’1
π π
νµ
ch e’»h t sin ν x sin µ y
uh (x, y, t) = νµ
a b
ν=1 µ=1

and
l’1 m’1
42 π π
ch = h u0 (ih, jh) sin ν ih sin µ jh ,
νµ
ab i=1 a b
j=1
4 πh πh
»νµ = sin2 ν + sin2 µ .
h
h2 a2 b2
Compared at the grid points u has additionally the terms in the in¬nite
series for ν = l, . . ., or µ = m, . . ..
They can be estimated by
∞ ∞ ∞ ∞
π π
cνµ e’»
νµ
t
+ sin ν x sin µ y
a b
ν=1 µ=m
ν=l µ=1
∞ ∞ ∞ ∞
e’»
νµ
¤ t
C1 +
ν=1 µ=m
ν=l µ=1

with C1 := max{|cνµ | , ν, µ ∈ N, ν ∈ {1, . . . , l ’ 1} or µ ∈ {1, . . . , m ’ 1}}
/ /
∞ ∞
2 2
µπ
’( νπ ) t
e ’( )
¤ C1 t
C2 e + C3
a b

µ=m
ν=l
2 π2
∞ ’( µπ ) t
¤ 1’q2 , where q2 := e’( b ) t because of
q2
with C2 := µ=1 e
b


for |q| < 1 , and C3 is de¬ned analogously (µ ←’
q
µ
µ=1 q = 1’q
π2
ν, a ←’ b) with an estimate by 1’q1 , q1 := e’( a ) t for t ≥ t > 0.
q1
7.2. Semidiscretization by Vertical Method of Lines 303

∞ ql
qµ =
Finally, we conclude the estimate because of by
µ=l 1’q
l
qm
q1
¤ C1 C2 + C3 2 .
1 ’ q1 1 ’ q2
Therefore, this error contribution for t ≥ t (for a ¬xed t > 0) approaches
0 for l ’ ∞ and m ’ ∞. The larger t is, the more this error term will
2
decrease. Because of, for example, l = a/h and thus q1 = exp ’ π t h , 1
l
a
the decay in h is exponential and thus much stronger than a term like
O(h2 ). Therefore, we have to compare the terms in the sum only for ν =
1, . . . , l ’ 1, µ = 1, . . . , m ’ 1, i.e., the error in the Fourier coe¬cient and
in the eigenvalue:
νµ νµ
cνµ e’» ’ ch e’»h = cνµ ’ ch e’» + ch e’» ’ e’»h
νµ νµ νµ
t t t t t
.
νµ νµ νµ

Note that ch can be perceived as an approximation of cνµ by the trape-
νµ
zoidal sum with step size h in each spatial direction (see, e.g., [30], p. 129),
since the integrand in the de¬nition of cνµ vanishes for x = 0 or x = a and
y ∈ [0, b], and y = 0 or y = b and x ∈ [0, a]. Thus we have for u0 ∈ C 2 („¦),
¯

|cνµ ’ ch | = O(h2 ).
νµ

Because of
νµ νµ
e’» ’ e’»h = e’» 1 ’ e’(»h ’»νµ )t
νµ νµ
t t t
,

and |»νµ ’ »νµ | = O(h2 ) (see (7.61)), also this term is of order O(h2 ) and
h
will be damped exponentially (depending on t and the size of the smallest
eigenvalue »νµ ).
Summarizing, we expect
O(h2 )
to be the dominating error term in the discrete maximum norm · ∞ at the
grid points (cf. de¬nition (1.17)), which will also be damped exponentially
for increasing t. Note that we have given only heuristic arguments and that
the considerations cannot be transferred to the Neumann case, where the
eigenvalue » = 0 appears.
We now turn to the ¬nite element method.

Order of Convergence Estimates for the Finite Element Method
We will investigate the ¬nite element method on a more abstract level as in
the previous subsection, but we will achieve a result (in di¬erent norms) of
similar character. As worked out at the end of Section 7.1, there is a strong
relation between the V -ellipticity of the bilinear form a with the parameter
± and a positive lower bound of the eigenvalues. Here we rely on the results
already achieved in Section 2.3 and Section 3.4 for the stationary case.
For that, we introduce the so-called elliptic projection of the solution u(t)
of (7.32) as a very important tool in the proof.
304 7. Discretization of Parabolic Problems

De¬nition 7.8 For a V -elliptic, continuous bilinear form a : V — V ’ R,
the elliptic, or Ritz, projection Rh : V ’ Vh is de¬ned by
v ’ Rh v ⇐’ a(Rh v ’ v, vh ) = 0 for all vh ∈ Vh .
Theorem 7.9 Under the assumptions of De¬nition 7.8:
(i) Rh : V ’ Vh is linear and continuous.
(ii) Rh yields quasi-optimal approximations; that is,
M
v ’ Rh v ¤ inf v ’ vh ,
V V
± vh ∈Vh
where M and ± are the Lipschitz and ellipticity constants according to
(2.42) and (2.43).

Proof: The linearity of Rh is obvious. The remaining statements
immediately follow from Lemma 2.16 and Theorem 2.17; see Exercise 7.6. 2

Making use of the elliptic projection, we are able to prove the following
result.
Theorem 7.10 Suppose a is a V -elliptic, continuous bilinear form, f ∈
C([0, T ], H), u0 ∈ V, and u0h ∈ Vh . Then if u(t) is su¬ciently smooth,
e’±t + (I ’ Rh )u(t)
uh (t) ’ u(t) ¤ u0h ’ Rh u0
0 0 0
t
e’±(t’s) ds .
(I ’ Rh )u (s)
+ 0
0


Proof: First, the error is decomposed as follows:
uh (t) ’ u(t) = uh (t) ’ Rh u(t) + Rh u(t) ’ u(t) =: θ(t) + (t) .
We take v = vh ∈ Vh in (7.32) and obtain, by the de¬nition of Rh ,
u (t), vh + a(u(t), vh ) = u (t), vh + a(Rh u(t), vh ) = b(t, vh ) .
0 0

Here b(t, ·) is as de¬ned in (7.43).
Subtracting this equation from (7.44), we get
’ u (t), vh
uh (t), vh + a(θ(t), vh ) = 0 ,
0 0

and thus
d
0’ =’
θ (t), vh 0 + a(θ(t), vh ) = u (t), vh Rh u(t), vh (t), vh .
0
dt 0
The application of Lemma 7.4 yields
t
’±t
e’±(t’s) ds .
¤ θ(0)
θ(t) 0e + (s)
0 0
0

Since the elliptic projection is continuous (Theorem 7.9, (i)) and u(t) is
d
su¬ciently smooth, Rh and the time derivative dt commute; that is, (t) =
7.2. Semidiscretization by Vertical Method of Lines 305

(Rh ’ I)u (t). It remains to apply the triangle inequality to get the stated
2
result.

Theorem 7.10 has the following interpretation:
The error norm uh (t) ’ u(t) 0 is estimated by
• the initial error (exponentially decaying in t), which occurs only if
u0h does not coincide with the elliptic projection of u0 ,
• the projection error of the exact solution u(t) measured in the norm
of H,
• the projection error of u (t) measured in the norm of H and integrally
weighted by the factor e’±(t’s) on (0, t).
Remark 7.11 If the bilinear form a de¬nes an elliptic problem such that
for the elliptic projection an error estimate of the type
(I ’ Rh )w ¤ Ch2 w for all w ∈ V © H 2 („¦)
2
0

is valid, if u0h approximates the elliptic projection Rh u0 of the initial value
u0 at least with the same asymptotic quality, and if the solution u of (7.44)
is su¬ciently smooth, then an optimal L2 -error estimate results:
uh (t) ’ u(t) ¤ C(u(t))h2 .
0

We see that in order to obtain semidiscrete error estimates, we need esti-
mates of the projection error measured in the norm of H = L2 („¦). Due to
· 0 ¤ · V , the quasi-optimality of Rh (Theorem 7.9, (ii)) in conjunction
with the corresponding approximation error estimates (Theorem 3.29) al-
ready yield some error estimate. Unfortunately, this result is not optimal.
However, if the adjoint boundary value problem is regular in the sense of
De¬nition 3.36, the duality argument (Theorem 3.37) can be successfully
used to derive an optimal result.
Theorem 7.12 Suppose the bilinear form a is V -elliptic and continuous,
and the solution of the adjoint boundary value problem is regular.
Furthermore, let the space Vh ‚ V be such that for any function w ∈
V © H 2 („¦),
w ’ vh ¤ C h |w|2 ,
inf V
vh ∈Vh

where the constant C > 0 does not depend on h and w. If u0 ∈ V © H 2 („¦),
then for a su¬ciently smooth solution u of (7.44) we have
e’±t
uh (t) ’ u(t) ¤ u0h ’ u0
0 0
t
’±t
e’±(t’s) ds .
2
+ Ch u0 e + u(t) + u (s)
2 2 2
0
306 7. Discretization of Parabolic Problems

Proof: The ¬rst term in the error bound from Theorem 7.10 is estimated
by means of the triangle inequality:
u0h ’ Rh u0 ¤ u0h ’ u0 + (I ’ Rh )u0 0.
0 0

Then the projection error estimate (Theorem 3.37, (1)) yields the given
bounds of the resulting second term as well as of the remaining two terms
2
in the error bound from Theorem 7.10.

Order of Convergence Estimates for the Finite Volume Method
For simplicity we restrict attention to pure homogeneous Dirichlet condi-
tions (“3 = ‚„¦). The idea is similar to the proof given in the ¬nite element
case. However, here we will meet some additional di¬culties, which are
caused by the use of perturbed bilinear and linear forms.
We take v = vh ∈ Vh in (7.32) and subtract the result from (7.48):
’ u (t), vh + ah (uh (t), vh ) ’ a(u(t), vh )
uh (t), vh 0
0,h
’ f (t), vh
= f (t), vh .
0,h 0

In analogy to the ¬nite element method, we introduce the following aux-
iliary problem: Given some v ∈ V, ¬nd an element Rh v ∈ Vh such
that
ah (Rh v, vh ) = a(v, vh ) for all vh ∈ Vh . (7.63)
With this, the above identity can be rewritten as follows:
’ u (t), vh + ah (uh (t) ’ Rh u(t), vh )
uh (t), vh 0
0,h
’ f (t), vh
= f (t), vh .
0,h 0
d
Subtracting from both sides of this relation the term dt Rh u(t), vh 0,h
and assuming that u (t) is a su¬ciently smooth function of x, a slight
rearrangement yields
+ ah (θ(t), vh ) = ’
θ (t), vh (t), vh + u (t), vh (7.64)
0,h 0,h 0
’ u (t), vh ’ f (t), vh
+ f (t), vh ,
0,h 0,h 0

where, as in the ¬nite element case, θ(t) = uh (t) ’ Rh u(t) and (t) =
Rh u(t) ’ u(t). Furthermore, we de¬ne, for v ∈ Vh , b1 (t, v) := u (t), v 0 ’
u (t), v 0,h and b2 (t, v) := f (t), v 0,h ’ f (t), v 0 .
In order to be able to apply the discrete stability estimate (7.50) to this
situation, we need an error estimate for Rh u (t) as in Remark 7.11 and
bounds (consistency error estimates) for |b1 (t, v)|, |b2 (t, v)|.
So we turn to the ¬rst problem. In fact, the estimate is very similar
to the error estimate for the ¬nite volume method given in the proof of
Theorem 6.18.
For an arbitrary function v ∈ V © H 2 („¦) and vh := Rh v ’ Ih (v), we have
by (7.63) that
ah (vh , vh ) = ah (Rh v, vh ) ’ ah (Ih (v), vh ) = a(v, vh ) ’ ah (Ih (v), vh ).
7.2. Semidiscretization by Vertical Method of Lines 307

By partial integration in the ¬rst term of the right-hand side, it follows
that
’ ah (Ih (v), vh ).
ah (vh , vh ) = Lv, vh 0

From [40] an estimate of the right-hand side is known (cf. also (6.22)); thus
1/2
ah (vh , vh ) ¤ Ch v |vh |2 + vh 2
.
2 1 0,h

So Theorem 6.15 yields
1/2
|vh |2 + vh ¤ Ch v
2
2.
1 0,h

By the triangle inequality,
(Rh ’ I)v ¤ Rh v ’ Ih (v) + Ih (v) ’ v 0,h .
0,h 0,h

·
Since the second term vanishes by the de¬nitions of and Ih , we get
0,h
in particular
(Rh ’ I)v ¤ Ch v 2. (7.65)
0,h

Remark 7.13 In contrast to the ¬nite element case (Remark 7.11), this
estimate is not optimal.
To estimate |b1 (t, v)| and |b2 (t, v)|, we prove the following result.
Lemma 7.14 Assume w ∈ C 1 („¦) and v ∈ Vh . Then, if the ¬nite volume
partition of „¦ is a Donald diagram,
| w, v ’ w, v | ¤ Ch|w|1,∞ v 0,h .
0,h 0


Proof: We start with a simple rearrangement of the order of summation:
M
wj vj |„¦j,K | ,
w, v = wj vj mj =
0,h
j=1 K∈Th j:‚K aj

where „¦j,K = „¦j © int K . First, we will consider the inner sum. For any
triangle K ∈ Th with barycentre aS,K , we can write

wj vj |„¦j,K | = [wj ’ w(aS,K )]vj |„¦j,K |
j:‚K aj j:‚K aj


w(aS,K ) vj |„¦j,K | ’
+ v dx
„¦j,K
j:‚K aj

[w(aS,K ) ’ w]v dx +
+ wv dx
„¦j,K „¦j,K
j:‚K aj j:‚K aj

=: I1,K + I2,K + I3,K + wv dx .
K
308 7. Discretization of Parabolic Problems

To estimate I1,K , we apply the Cauchy“Schwarz inequality and get
± 1/2
 
|I1,K | ¤ |wj ’ w(aS,K )| |„¦j,K |
2
v 0,h,K ,
 
j:‚K aj

where
± 1/2
 
vj |„¦j,K |
2
v := .
0,h,K
 
j:‚K aj

Since |wj ’ w(aS,K )| ¤ hK |w|1,∞ , it follows that

|I1,K | ¤ hK |w|1,∞ |K| v 0,h,K .

Similarly, for I3,K we easily get

|I3,K | = [w(aS,K ) ’ w]v dx
„¦K

¤ w(aS,K ) ’ w ¤ hK |w|1,∞ |K| v
v 0,K .
0,K 0,K

So it remains to consider I2,K . Obviously,

[vj ’ v] dx .
I2,K = w(aS,K )
„¦j,K
j:‚K aj

We will show that if „¦j belongs to a Donald diagram, then the sum van-
ishes. To do so, let us suppose that the triangle under consideration has
the vertices ai , aj , and ak . The set „¦j,K can be decomposed into two sub-
triangles by drawing a straight line between aS,K and aj . We will denote
the interior of these triangles by „¦j,K,i and „¦j,K,k ; i.e.,
„¦j,K,i := int conv{aj , aS,K , aij } , „¦j,K,k := int conv{aj , aS,K , akj } .
On each subtriangle, the integral of v can be calculated exactly by means
of the trapezoidal rule. Since |„¦j,K,i | = |„¦j,K,k | = |K|/6 in the case of the
Donald diagram (cf. also (6.4)), we have
|K| vj + vi vj + vi + vk
v dx = vj + +
18 2 3
„¦j,K,i
|K| 11 5 1
= vj + vi + vk ,
18 6 6 3
|K| 11 5 1
v dx = vj + vk + vi .
18 6 6 3
„¦j,K,k

Consequently,
|K| 11 7 7
v dx = vj + vi + vk ,
18 3 6 6
„¦j,K
7.2. Semidiscretization by Vertical Method of Lines 309

and thus
|K|
v dx = vj .
3
„¦j,K
j:‚K aj j:‚K aj

On the other hand, since 3|„¦j,K | = |K| (cf. (6.4)), we have
|K|
vj dx = vj ,
3
„¦j,K
j:‚K aj j:‚K aj

and so I2,K = 0. In summary, we have obtained the following estimate:

|I1,K + I2,K + I3,K | ¤ hK |w|1,∞ |K| v +v .
0,h,K 0,K

So it follows that
| w, v ’ w, v |¤ |I1,K + I2,K + I3,K |
0,h 0
K∈Th

¤ |K| v
h|w|1,∞ +v .
0,h,K 0,K
K∈Th

By the Cauchy“Schwarz inequality,
1/2 1/2

|K| v ¤ |K| |„¦| v
2
v =
0,h,K 0,h
0,h,K
K∈Th K∈Th K∈Th

and, similarly,

|K| v ¤ |„¦| v 0.
0,K
K∈Th

So we ¬nally arrive at
| w, v ’ w, v | ¤ Ch|w|1,∞ v +v .
0,h 0
0,h 0

· ·
Since the norms and are equivalent on Vh (see Remark 6.16),
0,h 0
we get
| w, v ’ w, v | ¤ Ch|w|1,∞ v 0,h .
0,h 0

2
Now we are prepared to apply the discrete stability estimate (7.50) to
equation (7.64):
’±t
¤
θ(t) θ(0) 0,h e
0,h
t
+ b2 (s) — ] e’±(t’s) ds ,
+ [ (s) + b1 (s) —
0,h
0

where |bj (t, v)| ¤ bj (t) — v 0,h for all v ∈ Vh , t ∈ (0, T ), and j = 1, 2. The
¬rst term in the integral can be estimated by means of (7.65), whereas the
310 7. Discretization of Parabolic Problems

estimates of b1 (s) — , b2 (s) result from Lemma 7.14:

’±t
¤
θ(t) θ(0) 0,h e
0,h
t
e’±(t’s) ds .
+ |u (s)|1,∞ + f (s)
+Ch [ u (s) 1,∞ ]
2
0

If u0 ∈ V © H 2 („¦), we can write, by (7.65),
¤ uh0 ’ u0 + (I ’ Rh )u0 ¤ uh0 ’ u0
θ(0) + Ch u0 2 .
0,h 0,h 0,h 0,h

So we get
’±t
+ Ch u0 2 e’±t
¤ uh0 ’ u0
θ(t) 0,h e
0,h
t
e’±(t’s) ds .
+ |u (s)|1,∞ + f (s)
+ [ u (s) 1,∞ ]
2
0
Since
uh (t) ’ u(t) ¤ θ(t) + (Rh ’ I)u(t) 0,h ,
0,h 0,h

the obtained estimate and (7.65) yield the following result.
Theorem 7.15 In addition to the assumptions of Theorem 6.15, let f ∈
C([0, T ], C 1 („¦)), u0 ∈ V © H 2 („¦), and u0h ∈ Vh . Then if u(t) is su¬ciently
smooth, the solution uh (t) of the semidiscrete ¬nite volume method (7.48)
on Donald diagrams satis¬es the following estimate:

e’±t + Ch u0 2 e’±t + u(t)
uh (t) ’ u(t) ¤ uh0 ’ u0
0,h 0,h 2
t
e’±(t’s) ds .
+ |u (s)|1,∞ + f (s)
+ [ u (s) 1,∞ ]
2
0
Remark 7.16 In comparison with the ¬nite element method, the result is
not optimal in h. The reason is that, in general, the ¬nite volume method
does not yield optimal L2 -error estimates even in the elliptic case, but this
type of result is necessary to obtain optimal estimates.


Exercises
7.4 Let A ∈ RM,M be an RM-elliptic matrix and let the symmetric positive
de¬nite matrix B ∈ RM,M have the Cholesky decomposition B = E T E.
Show that the matrix A := E ’T AE ’1 is RM-elliptic.
ˆ

7.5 Prove identity (7.62) by ¬rst proving the corresponding identity for
one space dimension:
l’1
π a
sin2 ν ih = .
h
a 2
i=1
7.3. Fully Discrete Schemes 311



7.6 Let V be a Banach space and a : V — V ’ R a V -elliptic, continuous
bilinear form. Show that the Ritz projection Rh : V ’ Vh in a subspace
Vh ‚ V (cf. De¬nition 7.8) has the following properties:
(i) Rh : V ’ Vh is continuous because of Rh u V ¤ M
u ,
V
±
(ii) Rh yields quasi-optimal approximations; that is,
M
u ’ Rh u ¤ inf u ’ vh .
V V
± vh ∈Vh
Here M and ± denote the constants in the continuity and ellipticity
conditions, respectively.

7.7 Let u ∈ C 1 ([0, T ], V ). Show that Rh u ∈ C 1 ([0, T ], V ) and d
dt Rh u(t) =
d
Rh dt u(t).

7.8 Transfer the derivation of the ¬nite volume method given in Sec-
tion 6.2.2 for the case of an elliptic boundary value problem to the parabolic
initial-boundary value problem (7.1) in divergence form; i.e., convince your-
self that the formalism of obtaining (7.48) indeed can be interpreted as a
¬nite volume semidiscretization of (7.1).



7.3 Fully Discrete Schemes
As we have seen, the application of the vertical method of lines results in
the following situation:
• There is a linear system of ordinary di¬erential equations of high
order (dimension) to be solved.
• There is an error estimate for the solution u of the initial-boundary
value problem (7.1) by means of the exact solution uh of the system
(7.41).
A di¬culty in the choice and in the analysis of an appropriate discretiza-
tion method for systems of ordinary di¬erential equations is that many
standard estimates involve the Lipschitz constant of the corresponding
right-hand side, here q h ’ Ah uh (cf. (7.41)). But this constant is typically
large for small spatial parameters h, and so we would obtain nonrealistic
error estimates (cf. Theorem 3.45).
There are two alternatives. For comparatively simple time discretiza-
tions, certain estimates can be derived in a direct way (i.e., without using
standard estimates for systems of ordinary di¬erential equations). The sec-
ond way is to apply speci¬c time discretizations in conjunction with re¬ned
methods of proof.
312 7. Discretization of Parabolic Problems

Here we will explain the ¬rst way for the so-called one-step-theta method.


One-Step Discretizations in Time, in Particular for the Finite
Di¬erence Method
We start from the problem (7.41), which resulted from spatial discretization
techniques. Provided that T < ∞, the time interval (0, T ) is subdivided
into N ∈ N subintervals of equal length „ := T /N. Furthermore, we set
tn := n„ for n ∈ {0, . . . , N } and un ∈ RM for an approximation of uh (tn ).
h
If the time interval is unbounded, the time step „ > 0 is given, and the
number n ∈ N is allowed to increase without bounded; that is, we set
formally N = ∞.
The values t = tn , where an approximation is to be determined, are
called time levels. The restriction to equidistant time steps is only for the
d
sake of simplicity. We approximate dt uh by the di¬erence quotient
d 1
uh (t) ∼ (uh (tn+1 ) ’ uh (tn )).
dt „
If we interpret this approximation to be at t = tn , we take the forward
di¬erence quotient; at t = tn+1 we take the backward di¬erence quotient;
at t = tn + 1 „ we take the symmetric di¬erence quotient. Again we obtain
2
a generalization and uni¬cation by introducing a parameter ˜ ∈ [0, 1] and
interpreting the approximation to be taken at t = tn + ˜„. As for ˜ = 0
or 1, we are not at a time level, and so we need the further approximation
Ah uh ((n + ˜) „ ) ∼ ˜Ah uh (tn ) + ˜Ah uh (tn+1 ).
Here we use the abbreviation ˜ := 1 ’ ˜. The (one-step-)theta method
de¬nes a sequence of vectors u0 , . . . , uN by, for n = 0, 1, . . . , N ’ 1 ,
h h
1 n+1
uh ’ un + ˜Ah un + ˜Ah un+1 = q h ((n + ˜)„ ) , (7.66)
h h h

u0 = u0 .
h

If we apply this discretization to the more general form (7.45), we get
correspondingly
1
Bh un+1 ’ Bh un + ˜Ah un + ˜Ah un+1 = q h ((n + ˜)„ ) .
ˆ ˆ (7.67)
h h
h h

Analogously to (7.45), the more general form can be transformed to (7.66),
’1
assuming that Bh is regular: either by multiplying (7.67) by Bh or in the
T
case of a decomposition Bh = Eh Eh (for a symmetric positive de¬nite Bh )
’T
by multiplying by Eh and a change of variables to Eh un . We will apply
h
two techniques in the following:
One is based on the eigenvector decomposition of Ah ; thus for (7.67),
this means to consider the generalized eigenvalue problem
ˆ
Ah v = »Bh v . (7.68)
7.3. Fully Discrete Schemes 313

Note that the Galerkin approach for the eigenvalue problems according
to De¬nition 7.6 leads to such a generalized eigenvalue problem with the
ˆ
sti¬ness matrix Ah and the mass matrix Bh .
The other approach is based on the matrix properties (1.32)* or (1.32).
For the most important case,

Bh = diag(bi ) , bi > 0 for i = 1, . . . , M , (7.69)

which corresponds to the mass lumping procedure, the above-mentioned
transformation reduces to a diagonal scaling, which does not in¬‚uence any
of their properties.
Having this in mind, in the following we will consider explicitly only the
formulation (7.66).
In the case ˜ = 0, the explicit Euler method, un can be determined
h
explicitly by

un+1 = „ (q h (tn ) ’ Ah un ) + un = (I ’ „ Ah )un + „ q(tn ) .
h h h
h

Thus the e¬ort for one time step consists of a SAXPY operation, a vector
addition, and a matrix-vector operation. For dimension M the ¬rst of these
is of complexity O(M ), and also the last one if the matrix is sparse in the
sense de¬ned at the beginning of Chapter 5. On the other hand, for ˜ = 0,
the method is implicit, since for each time step a system of linear equations
has to be solved with the system matrix I + ˜„ Ah . Here the cases ˜ = 1,
the implicit Euler method, and ˜ = 1 , the Crank“Nicolson method, will
2
be of interest. Due to our restriction to time-independent coe¬cients, the
matrix is the same for every time step (for constant „ ). If direct methods
(see Section 2.5) are used, then the LU factorization has to be computed
only once, and only forward and backward substitutions with changing
right-hand sides are necessary, where computation for ˜ = 1 also requires
a matrix-vector operation. For band matrices, for example, operations of
the complexity bandwidth — dimension are necessary, which means for the
basic example of the heat equation on a rectangle O(M 3/2 ) operations in-
stead of O(M ) for the explicit method. Iterative methods for the resolution
of (7.66) cannot make use of the constant matrix, but with un there is a
h
good initial iterate if „ is not too large.
Although the explicit Euler method ˜ = 0 seems to be attractive, we will
see later that with respect to accuracy or stability one may prefer ˜ = 1 2
or ˜ = 1.
To investigate further the theta method, we resolve recursively the
relations (7.66) to gain the representation
n
’1
I ’ ˜„ Ah
un = (I + ˜„ Ah ) u0 (7.70)
h
n n’k
’1 ’1
I ’ ˜„ Ah q h tk ’ ˜„ .
+„ (I + ˜„ Ah ) (I + ˜„ Ah )
k=1
314 7. Discretization of Parabolic Problems

Here we use the abbreviation A’n = (A’1 )n for a matrix A. Compar-
ing this with the solution (7.55) of the semidiscrete problem, we see the
approximations

e’Ah tn ∼ Eh,„ ,
n


where
’1
I ’ ˜„ Ah
Eh,„ := (I + ˜„ Ah )

and
tn tn
(tn ’s)/„
e’Ah (tn ’s) q h (s)ds e’Ah „
= q h (s)ds
0 0
n
(t ’s)/„
(I + ˜„ Ah )’1 q h s ’ ˜„ .
∼„ n
Eh,„
k=1
s=k„

(7.71)

The matrix Eh,„ thus is the solution operator of (7.66) for one time step
and homogeneous boundary conditions and right-hand side. It is to be
expected that it has to capture the qualitative behaviour of e’Ah „ that it
is approximating. This will be investigated in the next section.

One-Step Discretizations for the Finite Element Method
The fully discrete scheme can be achieved in two ways: Besides apply-
ing (7.66) to (7.41) in the transformed variable or in the form (7.67), the
discretization approach can applied directly to (7.44):
With ‚U n+1 := (U n+1 ’ U n )/„, f n+s := sf (tn+1 ) + (1 ’ s)f (tn ),
bn+s (v) := sb(tn+1 , v) + (1 ’ s)b(tn , v), b according to (7.43), s ∈ [0, 1],
and with a ¬xed number ˜ ∈ [0, 1], the fully discrete method for (7.44)
then reads as follows:
Find a sequence U 0 , . . . , U N ∈ Vh such that for n ∈ {0, . . . , N ’ 1},

‚U n+1 , vh + a(˜U n+1 + ˜U n , vh ) = bn+˜ (vh )
0
for all vh ∈ Vh , (7.72)
U0 = u0h .

An alternative choice for the right-hand side, closer to the ¬nite di¬erence
method, is the direct evaluation at tn + ˜„, e.g., f (tn + ˜„ ). The version
here is chosen to simplify the order of convergence estimate in Section 7.6.
By representing the U n by means of a basis of Vh as after (7.44), again
we get the form (7.67) (or (7.66) in the transformed variable). Note that
also for ˜ = 0 the problem here is implicit if Bh is not diagonal. Therefore,
mass lumping is often applied, and the scalar product ·, · 0 in (7.72) is
7.4. Stability 315

replaced by an approximation due to numerical quadrature, i.e.,
‚U n+1 , vh + a(˜U n+1 + ˜U n , vh ) = bn+˜ (vh )
0,h
for all vh ∈ Vh , (7.73)
U0 = u0h .
As explained in Section 3.5.2, uh , vh 0,h is the sum over all contributions
from elements K ∈ Th , which takes the form (3.112) for the reference
element. In the case of Lagrange elements and a nodal quadrature formula
we have for the nodal basis functions •i :
•j , •i = •i , •i δij =: bi δij for i, j = 1, . . . , M, (7.74)
0,h 0,h

since for i = j the integrand •i •j vanishes at all quadrature points. In this
case we arrive at the form (7.67) with a matrix Bh satisfying (7.69).

One-Step Discretizations for the Finite Volume Method
As in the previous subsection on the ¬nite element approach, the
semidiscrete formulation (7.48) can be discretized in time directly:
Find a sequence U 0 , . . . , U N ∈ Vh such that for n ∈ {0, . . . , N ’ 1},
‚U n+1 , vh + ah (˜U n+1 + ˜U n , vh ) = f n+˜ , vh 0,h
0,h
for all vh ∈ Vh , (7.75)
U0 = u0h ,
where ‚U n+1 , ˜, f n+˜ are de¬ned as before (7.72).
Remember that here we consider only homogeneous boundary conditions.
If the elements U n , U n+1 are represented by means of a basis of Vh , we
recover the form (7.67).
Since the mass matrix Bh is diagonal, the problem can be regarded as
being explicit for ˜ = 0.


Exercise
7.9 Consider linear simplicial elements de¬ned on a general conforming
triangulation of a polygonally bounded domain „¦ ‚ R2 .
(a) Determine the entries of the mass matrix Bh .
(b) Using the trapezoidal rule, determine the entries of the lumped mass
matrix diag(bi ).


7.4 Stability
In Section 7.3 we have seen that at least if a basis of eigenvectors of the
discretization matrix Ah allows for the solution representation (7.55) for
316 7. Discretization of Parabolic Problems

the semidiscrete method, the qualitative behaviour of e’Ah „ u0 should be
preserved by Eh,„ u0 , being one time step „ for homogeneous boundary
conditions and right-hand side (q h = 0) in the semi- and fully discrete
cases. It is su¬cient to consider the eigenvector w i instead of a general u0 .
Thus, we have to compare
e’Ah „ w i = (e’»i „ )w i (7.76)
with
1 ’ ˜„ »i
’1
I ’ ˜„ Ah
(I + ˜„ Ah ) wi = wi . (7.77)
1 + ˜„ »i
We see that the exponential function is approximated by
1 + (1 ’ ˜)z
R(z) = , (7.78)
1 ’ ˜z
the stability function, at the points z = ’»i „ ∈ C, given by the eigenvalues

<<

. 12
( 16)



>>