+ [ (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

»i , and the time step size „ .

For n time steps and q h = 0 we have

n

e’Ah „ wi = e’»i tn wi ∼ R(’»i „ )n w i . (7.79)

Thus, the restriction to eigenvectors wi with eigenvalues »i has diagonal-

ized the system of ordinary di¬erential equations (7.41) for q h = 0 to the

scalar problems

t ∈ (0, T ) ,

ξ + »i ξ = 0, (7.80)

ξ(0) = ξ0

(for ξ0 = 1) with its solution ξ(t) = e’»i tξ0 , for which the one-step-theta

method gives the approximation

ξn+1 = R(’»i „ )ξn = (R(’»i „ ))n+1 ξ0 (7.81)

at t = tn+1 . A basic requirement for a discretion method is the following:

De¬nition 7.17 A one-step method is called nonexpansive if for two nu-

merical approximations un and un , generated under the same conditions

˜h

h

˜

except for two discrete initial values u0 and u0 , respectively, the following

estimate is valid:

|un+1 ’ un+1 | ¤ |un ’ un | , n ∈ {0, . . . , N ’ 1} .

˜h ˜

h

A recursive application of this estimate immediately results in

|un ’ un | ¤ |u0 ’ u0 | , n ∈ {1, . . . , N } .

˜ ˜

Here a general one-step method has the form

n ∈ {0, . . . , N ’ 1} ,

un+1 = un + „ ¦(„, tn , un ) ,

h h

h

with u0 = u0 and a so-called generating function ¦ : R+ — [0, T ) — RM ’

h

RM that characterizes the particular method. The generating function of

7.4. Stability 317

the one-step-theta method applied to the system (7.41) is