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