method. The de¬ning system of equations has the form

Ch uh = ph , (7.97)

7.5. Maximum Principle for the One-Step-Theta Method 325

where

«

I + „ ˜Ah

¬ ·

..

0

¬ ’I + „ ˜Ah ·

.

¬ ·

¬ ·

.. ..

Ch = ¬ ·,

. .

¬ ·

¬ ·

.. ..

0

. .

’I + „ ˜Ah I + „ ˜Ah

again with ˜ := 1 ’ ˜,

«

„ q h (˜„ ) + (I ’ „ ˜Ah )u0

¬ ·

¬ ·

„ q h ((1 + ˜)„ )

¬ ·

¬ ·

.

ph = ¬ ·.

.

¬ ·

.

¬ ·

.

.

.

„ q h (N ’ 1 + ˜)„ )

Since the spatial discretization is performed as in the stationary case, and

in the nth step the discretization relates to t = tn’1 + ˜„ and also the

approximation

Ah u(tn’1 + ˜„ ) ∼ ˜Ah u(tn’1 ) + ˜Ah u(tn )

enters the formulation (7.66), we can assume to have the following structure

of the right-hand side of (7.66):

q h ((n ’ 1 + ˜)„ ) = ’Ah (˜un’1 + ˜un ) + f ((n ’ 1 + ˜)„ ) for n = 1, . . . , N.

ˆ ˆh ˆh

(7.98)

Here the uh ∈ R

n M2

ˆ are the known spatial boundary values on time level

tn , which have been eliminated from the equation as explained, e.g., in

Chapter 1 for the ¬nite di¬erence method. But as noted, we allow also for

the case where such values do not appear (i.e., M2 = 0) then (7.98) reduces

to

q h ((n ’ 1 + ˜)„ ) = f ((n ’ 1 + ˜)„ ) for n = 1, . . . , N .

For the continuous problem, data are prescribed at the parabolic boundary

„¦ — {0} ∪ ‚„¦ — [0, T ]; correspondingly, the known values ui are collected

ˆh

with the initial data u0 ∈ R M1

to a large vector

«

u0

¬ u0 ·

¬ ˆh ·

¬ ˆ1 ·

uh = ¬ uh · ,

ˆ

¬.·

. .

uN

ˆh

326 7. Discretization of Parabolic Problems

i.e., a vector of dimension M 2 := M1 + (N + 1)M2 , which may reduce to

uh = u0 ∈ RM1 .

ˆ

With this notation we have

ph = ’Ch uh + e

ˆˆ (7.99)

if we de¬ne

«

’I + „ ˜Ah ˆ ˆ

„ ˜Ah „ ˜Ah O

¬ ·

.

.. ..

¬ ·

.

. .

ˆh = ¬ ·

O .

¬ ·,

C . .

¬ ·

.. ..

. .

. .

. .

ˆ ˆ

O „ ˜Ah „ ˜Ah

«

„ f (˜„ )

¬ ·

¬ ·

„ f ((1 + ˜)„ )

¬ ·

¬ ·

.

e=¬ ·.

.

.

¬ ·

¬ ·

.

.

.

„ f ((N ’ 1 + ˜)„ )

In the following the validity of (1.32)— or (1.32) for

˜ ˆ

Ch = (Ch , Ch )

will be investigated on the basis of corresponding properties of

˜ ˆ

Ah = (Ah , Ah ) .

Note that even if Ah is irreducible, the matrix Ch is always reducible,

since un depends only on u1 , . . . , un’1 , but not on the future time levels.

h h h

(Therefore, (7.97) serves only for the theoretical analysis, but not for the

actual computation.)

In the following we assume that

„ ˜(Ah )jj < 1 for j = 1, . . . , M1 , (7.100)

which is always satis¬ed for the implicit Euler method (˜ = 1). Then:

(1) (Ch )rr > 0 f or r = 1, . . . , M 1

holds if (1) is valid for Ah . Actually, also (Ah )jj > ’1/(„ ˜) would

be su¬cient.

(2) (Ch )rs ¤ 0 f or r, s = 1, . . . , M 1 , r = s:

If (2) is valid for Ah , then only the nonpositivity of the diagonal

elements of the o¬-diagonal block of Ch , ’I + „ ˜Ah , is in question.

This is ensured by (7.100) (weakened to “¤”).

M1

≥ 0 f or r = 1, . . . , M 1 :

(3) (i) Cr := Ch

rs

s=1

7.5. Maximum Principle for the One-Step-Theta Method 327

(ii) Cr > 0 f or at least one r ∈ {1, . . . , M 1 }:

We set

M1

Aj := (Ah )jk ,

k=1

so that condition (3) (i) for Ah means that Aj ≥ 0 for j = 1, . . . , M1 .

Therefore, we have

Cr = 1 + „ ˜Aj > 0 (7.101)

for the indices r of the ¬rst time level, where the “global” index r

corresponds to the “local” spatial index j. For the following time

levels, the relation

Cr = 1 ’ 1 + „ (˜ + ˜)Aj = „ Aj ≥ 0 (7.102)

holds, i.e., (3) (i) and (ii).

(4)— For every r1 ∈ {1, . . . , M 1 } satisfying

M1

(Ch )rs = 0 (7.103)

r=1

there exist indices r2 , . . . , rl+1 such that

(Ch )ri ri+1 = 0 for i = 1, . . . , l

and

M1

(Ch )rl+1 s > 0 . (7.104)

s=1

To avoid too many technicalities, we adopt the background of a ¬nite

di¬erence method. Actually, only matrix properties enter the reason-

ing. We call (space-time) grid points satisfying (7.103) far from the

boundary, and those satisfying (7.104) close to the boundary. Due to

(7.101), all points of the ¬rst time level are close to the boundary

(consistent with the fact that the grid points for t0 = 0 belong to the

parabolic boundary). For the subsequent time level n, due to (7.102),

a point (xi , tn ) is close to the boundary if xi is close to the bound-

ary with respect to Ah . Therefore, the requirement of (4)— , that a

˜

point far from the boundary can be connected via a chain of neigh-

bours to a point close to the boundary, can be realized in two ways:

Firstly, within the time level n, i.e., the diagonal block of Ch if Ah

satis¬es condition (4)— . Secondly, without this assumption a chain of

neighbours exist by (x, tn ), (x, tn’1 ) up to (x, t1 ), i.e., a point close

to the boundary, since the diagonal element of ’I + „ ˜Ah does not

vanish due to (7.100). This reasoning additionally has established the

following:

328 7. Discretization of Parabolic Problems

(4)# If Ah is irreducible, then a grid point (x, tn ), x ∈ „¦h can be connected

via a chain of neighbours to every grid point (y, tk ), y ∈ „¦h and

0 ¤ k ¤ n.

(5) (Ch )rs ¤ 0 for r = 1, . . . , M 1 , s = M 1 + 1, . . . , M 2 :

ˆ

ˆ

Analogously to (2), this follows from (5) for Ah and (7.100).

M

(6)— Cr :=

˜ ˜

(Ch )rs = 0 for r = 1, . . . , M :

s=1

Analogously to (7.102), we have

M

˜ ˜ ˜

Cr = „ Aj := „ (Ah )jk ,

k=1

˜

so that the property is equivalent to the corresponding one of Ah .

(6) Cr ≥ 0 for r = 1, . . . , M

˜

˜

is equivalent to (6) for Ah by the above argument.

(7) For every s ∈ M 1 + 1, . . . , M there exists an r ∈ {1, . . . , M 1 } such

ˆ

that (Ch )rs = 0:

Every listed boundary value should in¬‚uence the solution: For the

ˆ

values from u0 , . . . , uN this is the case i¬ Ah satis¬es this property.

ˆh ˆh

Furthermore, the “local” indices of the equation, where the boundary

values appear, are the same for each time level. For the values from

u0 ∈ RM1 the assertion follows from (7.100).

From the considerations we have the following theorem:

Theorem 7.28 Consider the one-step-theta method in the form (7.66).

ˆ

Let (7.100) hold. If the spatial discretization Ah satis¬es (1.32) (1), (2),

(3) (i), and (5), then a comparison principle holds:

(1) If for two sets of data f i , u0i and un , n = 0, . . . , N and i = 1, 2, we

ˆ hi

have

f 1 ((n ’ 1 + ˜)„ ) ¤ f 2 ((n ’ 1 + ˜)„ ) for n = 1, . . . , N ,

and

u01 ¤ u02 ; , un ¤ un

ˆ h1 ˆ h1 for n = 0, . . . , N,

then

un ¤ un

ˆ h1 ˆ h2 for n = 1, . . . , N

for the corresponding solutions.

If un = un for n = 1, . . . , N, then condition (1.32) (5) can be

ˆ h1 ˆ h2

omitted.

7.5. Maximum Principle for the One-Step-Theta Method 329

(2) If Ah additionally satis¬es (1.32) (6)— , then the following weak

˜

maximum principle holds:

(ˆ n )r ,

(˜ n )r ¤ max

max uh max (u0 )r , max uh

r∈{1,...,M1 }

r∈{1,...,M } r∈{M1 +1,...,M }

n=0,...,N n=0,...,N

where

un

un := h

˜h .

ˆ

uh

˜

(3) If Ah satis¬es (1.32) (1), (2), (3) (i), (4), (5), (6), (7), then a strong

maximum principle in the following sense holds:

If the components of un , n = 0, . . . , N, attain a nonnegative maxi-

˜h

mum for some spatial index r ∈ {1, . . . , M1 } and at some time level

k ∈ {1, . . . , N }, then all components for the time levels n = 0, . . . , k

are equal.

Proof: Only part (3) needs further consideration. Theorem 1.9 cannot

be applied directly to (7.97), since Ch is reducible. Therefore, the proof of

Theorem 1.9 has to be repeated: We conclude that the solution is constant

at all points that are connected via a chain of neighbours to the point where

the maximum is attained. According to (4)# these include all grid points

(x, tl ) with x ∈ „¦h and l ∈ {0, . . . , k}. From (7.100) and the discussion of

(7) we see that the connection can also be continued to boundary values

2

up to level k.

The additional condition (7.100), which may be weakened to nonstrict

inequality, as seen above, actually is a time step restriction: Consider again

the example of the ¬ve-point stencil discretization of the heat equation on

a rectangle, for which we have (Ah )jj = 4/h2 . Then the condition takes

the form

„ 1

< (7.105)

4(1 ’ ˜)

2

h

for ˜ < 1. This is very similar to the condition (7.87), (7.88) for the explicit

Euler method, but the background is di¬erent.

As already noted, the results above also apply to the more general form

(7.67) under the assumption (7.69). The condition (7.100) then takes the

form

„ ˜(Ah )jj ¤ bj for j = 1, . . . , M1 .

Exercises

7.14 Formulate the results of this section, in particular condition (7.100),

for the problem in the form (7.67) with Bh according to (7.69) (i.e.

330 7. Discretization of Parabolic Problems

appropriate for ¬nite element discretizations with mass lumping, see

(7.74)).

7.15 Show the validity of (6)# from Exercise 1.13 for Ch if it holds here

for Ah and conclude as in Exercise 1.13 a weak maximum principle for the

one-step-theta method.

7.16 Consider the initial-boundary value problem in one space dimension

±

ut ’ µuxx + cux = f in (0, 1) — (0, T ) ,

u(0, t) = g’ (t), u(1, t) = g+ (t) , t ∈ (0, T ) ,

u(x, 0) = u0 (x) , x ∈ (0, 1) ,

where T > 0 and µ > 0 are constants, and c, f : (0, 1) — (0, T ) ’ R,

u0 : (0, 1) ’ R, and g’ , g+ : (0, T ) ’ R are su¬ciently smooth functions

such that the problem has a classical solution.

De¬ne h := 1/m and „ = T /N for some numbers m, N ∈ N. Then the so-

called full-upwind ¬nite di¬erence method for this problem reads as follows:

Find a sequence of vectors u0 , . . . , uN by

h h

un+1 ’ 2un+1 + un+1 ’ ui+1 ’ ui ’ un+1

n+1 n+1 n+1

un+1 ’ un + ui

i+1 i i’1 i’1

’µ ’c

i

i

+c

2

„ h h h

i = 1, . . . , m ’ 1, n = 0, . . . , N ’ 1,

n+1

= fi ,

where c = c+ ’ c’ with c+ = max{c, 0}, fin = f (ih, n„ ), u0 = u0 (ih),

i

n n

u0 = g’ (n„ ) and um = g+ (n„ ).

Prove that a weak maximum principle holds for this method.

7.6 Order of Convergence Estimates

Based on stability results already derived, we will investigate the (order

of) convergence properties of the one-step-theta method for di¬erent dis-

cretization approaches. Although the results will be comparable, they will

be in di¬erent norms, appropriate for the speci¬c discretization method, as

already seen in Chapters 1, 3, and 6.

Order of Convergence Estimates for the Finite Di¬erence Method

From Section 1.4 we know that the investigation of the (order of)

convergence of a ¬nite di¬erence method consists of two ingredients:

• (order of) convergence of the consistency error

• stability estimates.

The last tool has already been provided by Theorem 7.26 and by Theo-

rem 1.14, which together with the considerations of Section 7.5 allow us to

concentrate on the consistency error. Certain smoothness properties will be

7.6. Order of Convergence Estimates 331

required for the classical solution u of the initial boundary value problem

(7.1), which in particular makes its evaluation possible at the grid points

xi ∈ „¦h at each instance of time t ∈ [0, T ] and also of various derivatives.

The vector representing the corresponding grid function (for a ¬xed order-

ing of the grid points) will be denoted by U (t), or for short by U n := U (tn )

for t = tn . The corresponding grid points depend on the boundary condi-

tion. For a pure Dirichlet problem, the grid points will be from „¦h , but if

Neumann or mixed boundary conditions appear, they are from the enlarged

set

„¦h := „¦h © („¦ ∪ “1 ∪ “2 ) .

˜ (7.106)

Then the error at the grid points and each time level is given by

en := U n ’ un for n = 0, . . . , N , (7.107)

h h

where un is the solution of the one-step-theta method according to (7.66).