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ν

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