6.8

(a) Formulate problem (6.11) in terms of an algebraic system of type

(1.31).

(b) Show that for the resulting matrix Ah ∈ RM1 ,M1 , where M1 is the

number of elements of the index set Λ, the following relation is valid:

AT 1 ≥ 0 . Here, as in Section 1.4, 0, respectively 1, denotes a vector

h

of dimension M1 whose components are all equal to 0, respectively 1.

(This is nothing other than the property (1.32)(3)(i) except for the

transpose of Ah .)

7

Discretization Methods for Parabolic

Initial Boundary Value Problems

7.1 Problem Setting and Solution Concept

In this section initial boundary value problems for the linear case of the

di¬erential equation (0.33) are considered. We choose the form (3.12) to-

gether with the boundary conditions (3.18)“(3.20), which have already been

discussed in Section 0.5. In Section 3.2 conditions have been developed to

ensure a unique weak solution of the stationary boundary value problem. In

contrast to Chapter 3, the heterogeneities are now allowed also to depend

on time t, but for the sake of simplicity we do not do so for the coe¬cients

in the di¬erential equations and the boundary conditions, which covers

most of the applications, for example from Chapter 0. Also for the sake

of simplicity, we take the coe¬cient in front of the time derivative to be

constant and thus 1 by a proper scaling. From time to time we will restrict

attention to homogeneous Dirichlet boundary conditions for further ease of

exposition. Thus the problem reads as follows:

The domain „¦ is assumed to be a bounded Lipschitz domain and we

suppose that “1 , “2 , “3 form a disjoint decomposition of the boundary ‚„¦

(cf. (0.39)):

‚„¦ = “1 ∪ “2 ∪ “3 ,

where “3 is a closed subset of the boundary.

In the space-time cylinder QT = „¦ — (0, T ), T > 0, and its boundary

ST = ‚„¦ — (0, T ) there are given functions f : QT ’ R, g : ST ’ R,

g(x, t) = gi (x, t) for x ∈ “i , i = 1, 2, 3, and u0 : „¦ ’ R. The problem is to

284 7. Discretization of Parabolic Problems

¬nd a function u : QT ’ R such that

‚u

+ Lu = f in QT ,

‚t

(7.1)

Ru = g on ST ,

„¦ — {0} ,

u = u0 on

where Lv denotes the di¬erential expression for some function v : „¦ ’ R,

(Lv) (x) := ’∇ · (K(x) ∇v(x)) + c(x) · ∇v(x) + r(x)v(x) (7.2)

with su¬ciently smooth, time-independent coe¬cients

K : „¦ ’ Rd,d , c : „¦ ’ Rd , r : „¦ ’ R.

The boundary condition is expressed by the shorthand notation Ru = g,

which means, for a function ± : “2 ’ R on ‚„¦,

• Neumann boundary condition (cf. (0.41) or (0.36))

K∇u · ν = ‚νK u = g1 on “1 — (0, T ) , (7.3)

• mixed boundary condition (cf. (0.37))

K∇u · ν + ±u = ‚νK u + ±u = g2 on “2 — (0, T ) , (7.4)

• Dirichlet boundary condition (cf. (0.38))

on “3 — (0, T ) .

u = g3 (7.5)

Thus the stationary boundary problem considered so far reads

for x ∈ „¦ ,

Lu(x) = f (x) (7.6)

for x ∈ ‚„¦ .

Ru(x) = g(x)

It is to be expected that both for the analysis and the discretization there

are strong links between (7.6) and (7.1). The formulation (7.1) in particular

includes the heat equation (cf. (0.20))

‚u

’ ∇ · (K∇u) = f in QT , (7.7)

‚t

or for constant scalar coe¬cients in the form (cf. (0.19))

‚u

’ ∆u = f in QT (7.8)

‚t

with appropriate initial and boundary conditions.

Again as in Chapter 1, one of the simplest cases will be, for two space

dimensions (d = 2), the case of a rectangle „¦ = (0, a) — (0, b) or even the

case d = 1 (with „¦ = (0, a)), for which (7.8) further reduces to

‚2

‚u

’ 2 u = 0 in QT = (0, a) — (0, T ). (7.9)

‚t ‚x

For problem (7.1), the following typical analytical questions arise:

7.1. Problem Setting and Solution Concept 285

• existence of (classical) solutions,

• properties of the (classical) solutions,

• weaker concepts of the solution.

As in the case of elliptic boundary value problems, the theory of classical

solutions requires comparatively strong assumptions on the data of the

initial-boundary value problem. In particular, along the edge ‚„¦ — {0}

of the space-time cylinder initial and boundary conditions meet, so that

additional compatibility conditions have to be taken into account.

Representation of Solutions in a Special Case

To enhance the familiarity with the problem and for further comparison we

brie¬‚y sketch a method, named separation of variables, by which closed-

form solutions in the form of in¬nite series can be obtained for special cases.

Also in these cases, the representations are not meant to be a numerical

method (by its evaluation), but only serve as a theoretical tool.

We start with the case of homogeneous data, i.e., f = 0, gi = 0 (i =

1, 2, 3), so that the process is determined only by the initial data u0 .

We assume a solution of (7.1) to have the form u(x, t) = v(t)w(x) with

v(t) = 0, w(x) = 0. This leads to

’Lw(x)

v (t)

x ∈ „¦, t ∈ (0, T ) .

= , (7.10)

v(t) w(x)

Therefore, the expressions in (7.10) must be constant, for example, equal

to ’» for » ∈ R. Therefore,

v (t) = ’»v(t) , t ∈ (0, T ), (7.11)

which for the initial conditions v(0) = 1 has the solution

v(t) = e’»t .

Furthermore, w has to satisfy

Lw(x) = »w(x) , x ∈ „¦ ,

(7.12)

x ∈ ‚„¦ .

Rw(x) = 0 ,

Such a function w : „¦ ’ R, w = 0, is called an eigenfunction for the eigen-

¯

value » of the boundary value problem (7.6). If (wi , »i ), i = 1, . . . , N, are

eigenfunctions/values for (7.6), then because of the superposition principle,

the function

N

ci e’»i t wi (x)

u(x, t) := (7.13)

i=1

286 7. Discretization of Parabolic Problems

is a solution of the homogeneous initial-boundary value problem for the

initial value

N

u0 (x) := ci wi (x) , (7.14)

i=1

where the ci ∈ R are arbitrary. If there are in¬nitely many eigenfunc-

tions/values (wi , »i ) and if the sums in (7.13) and (7.14) converge in such

a way that also the in¬nite series possesses the derivatives appearing in

(7.6), then also

∞

ci e’»i t wi (x)

u(x, t) = (7.15)

i=1

is a solution to

∞

u0 (x) = ci wi (x) . (7.16)

i=1

For an inhomogeneous right-hand side of the form

N

f (x, t) = fi (t)wi (x) (7.17)

i=1

the solution representation can be extended to (variation of constants

formula)

t

N N

ci e’»i t wi (x) + fi (s)e’»i (t’s) ds wi (x) ,

u(x, t) := (7.18)

i=1 i=1 0

and at least formally the sum can be replaced by the in¬nite series. To

verify (7.18) it su¬ces to consider the case u0 = 0, for which we have

t

N N

fi (s)e’»i (t’s) ds »i wi (x)

fi (t)wi (x) ’

(‚t u)(x, t) =

i=1 i=1 0

t

N (7.19)

fi (s)e’»i (t’s) ds wi (x)

f (x, t) ’ L

=

i=1 0

f (x, t) ’ L(u)(x, t) .

=

From these solution representations we can conclude that initial data (and

thus also perturbances contained in it) and also the in¬‚uence of the right-

hand side act only exponentially damped if all eigenvalues are positive.

For d = 1, „¦ = (0, a) and Dirichlet boundary conditions we have the

eigenfunctions

π

ν ∈ N,

wν (x) = sin ν x , (7.20)

a

7.1. Problem Setting and Solution Concept 287

for the eigenvalues

νπ 2

»ν = . (7.21)

a

If the initial data u0 has the representation

∞

π

u0 (x) = ci sin ν x , (7.22)

a

ν=1

then for example for f = 0 the (formal) solution reads

∞

π

ci e’»ν t sin ν

u(x, t) = x . (7.23)

a

ν=1

The eigenfunctions wν are orthogonal with respect to the scalar product

·, · 0 in L2 („¦), since they satisfy

0 for ν = µ,

π π

· ·

sin ν , sin µ = (7.24)

a

a a for ν = µ,

0

2

which can by checked by means of well-known identities for the trigono-

metric functions.

Therefore (see below (7.57)),

u0 , w ν 0

ci = , (7.25)

wν , wν 0

which is called the Fourier coe¬cient in the Fourier expansion of u0 .

Of course, the (wν , »ν ) depend on the boundary conditions. For Neumann

boundary conditions in x = 0 and x = a we have

wν (x) = cos ν π x , ν = 0, 1, . . . ,

a

(7.26)

2

νπ

ν

»= , ν = 0, 1, . . . .

a

The occurrence of w0 = 1, »0 = 0 re¬‚ects the nontrivial solvability of the

pure Neumann problem (which therefore is excluded by the conditions of

Theorem 3.15).

For Lu = ’∆u and „¦ = (0, a) — (0, b), eigenfunctions and eigenvalues

can be derived from the one-dimensional case because of

˜

’∆(v ν (x)˜µ (y)) = ’v ν (x)˜µ (y) ’ v ν (x)˜µ (y) = (»ν + »µ )v ν (x)˜µ (y).

v v v v

Therefore, for „¦ = (0, a)—(0, b) one has to choose the eigenfunctions/values

(v ν , »ν ) (in x, on (0, a)) for the required boundary conditions at x = 0 and

v˜

x = a, and (˜µ , »µ ) (in y, on (0, b)) for the required boundary conditions

at y = 0, y = b.

For Dirichlet boundary conditions everywhere this leads to

π π

wνµ (x, y) = sin ν x sin µ x (7.27)

a b

288 7. Discretization of Parabolic Problems

for the eigenvalues

2 2

νπ µπ

»νµ = +

a b

2 2

and »νµ ’ ∞ for ν ’ ∞ or

(i.e., the smallest eigenvalue is π + π

a b

µ ’ ∞).

As a further concluding example we note the case

x = 0 or x = a : u(x, y) = 0 for y ∈ [0, b] ,

y = 0 : ∇u · ν(x, y) = ’‚2 u(x, y) = 0 for x ∈ (0, a) ,

y = b : ∇u · ν(x, y) = ‚2 u(x, y) = 0 for x ∈ (0, a) .

Eigenfunctions:

π π

wνµ (x, y) = sin ν x cos µ y , (7.28)

a b

ν = 1, 2, . . . , µ = 0, 1, 2, . . . .

Eigenvalues:

2 2

π π

νµ

» =ν +µ .

a b

A Sketch of the Theory of Weak Solutions

As in the study of the elliptic boundary value problems (3.12), (3.18)“

(3.20), for equation (7.1) a weak formulation can be given that reduces the

requirements with respect to the di¬erentiability properties of the solution.

The idea is to treat time and space variables in a di¬erent way:

• For ¬xed t ∈ (0, T ), the function x ’ u(x, t) is interpreted

(1)

as a parameter-dependent element u(t) of some space V whose

elements are functions of x ∈ „¦. An obvious choice is (see

Subsection 3.2.1, (I)) the space

V = {v ∈ H 1 („¦) : v = 0 on “3 }.

• In a next step, that is, for varying t, a function t ’ u(t) results

with values in the (function) space V.

(2) In addition to V, a further space H = L2 („¦) occurs, from which the

initial value u0 is taken and which contains V as a dense subspace.

A subspace V is called dense in H if the closure of V with respect to

the norm on H coincides with H.

(3) The time derivative is understood in a generalized sense; see (7.29).

(4) The generalized solution t ’ u(t) is sought as an element of a function

space, the elements of which are “function-valued” (cf. (1)).

De¬nition 7.1 Let X denote one of the spaces H or V (in particular, this

means that the elements of X are functions on „¦ ‚ Rd ).

7.1. Problem Setting and Solution Concept 289

(i) The space C l ([0, T ], X), l ∈ N0 , consists of all continuous functions

v : [0, T ] ’ X that have continuous derivatives up to the order l on

[0, T ] with the norm

l

v (i) (t)

sup .

X

i=0 t∈(0,T )

For the sake of simplicity, the notation C([0, T ], X) := C 0 ([0, T ], X)

is used.

(ii) The space Lp ((0, T ), X) with 1 ¤ p ¤ ∞ consists of all functions on

(0, T ) — „¦ with the following properties:

v(t, ·) ∈ X for any t ∈ (0, T ), F ∈ Lp (0, T ) with F (t) := v(t, ·) .

X

Furthermore,

v := F .

Lp ((0,T ),X) Lp (0,T )

Remark 7.2 f ∈ L2 (QT ) ’ f ∈ L2 ((0, T ), H) .

Proof:

2

Basically, the proof is a consequence of Fubini™s theorem (see [1]).

Concerning the interpretation of the time derivative and of the weak

formulation, a comprehensive treatment is possible only within the frame-

work of the theory of distributions; thus a detailed explanation is beyond

the scope of this book. A short but mathematically rigorous introduction

can be found in the book [39, Chapter 23].

The basic idea consists in the following de¬nition:

A function u ∈ L2 ((0, T ), V ) is said to have a weak derivative w if the

following holds:

T T

∞

u(t) Ψ (t) dt = ’ for all Ψ ∈ C0 (0, T ) .

w(t) Ψ(t) dt (7.29)

0 0

du

Usually, this derivative w is denoted by or u .

dt

Remark 7.3 The integrals occurring above are to be understood as so-

called Bochner integrals and are extensions of the Lebesgue integral to

function-valued mappings. Therefore, equation (7.29) is an equality of

functions.

Before we give a weak formulation of (7.1), the following notion is worth

recalling:

u v dx for u, v ∈ H ,

u, v := (7.30)

0

„¦

[K∇u · ∇v + (c · ∇u + ru) v] dx+ ±uv dσ, u, v ∈ V. (7.31)

a(u, v) :=

„¦ “2

290 7. Discretization of Parabolic Problems

Let u0 ∈ H, f ∈ L2 ((0, T ), H), and in case of Dirichlet conditions we

restrict ourselves to the homogeneous case.

An element u ∈ L2 ((0, T ), V ) is called a weak solution of (7.1) if it has

a weak derivative du = u ∈ L2 ((0, T ), H) and the following holds

dt

d

u(t), v + a (u(t), v) = f (t), v + g1 (·, t)v dσ

0

dt “1

0

+ g2 (·, t)v dσ (7.32)

“2

for all v ∈ V and t ∈ (0, T ) ,

u(0) = u0 .

Due to u ∈ L2 ((0, T ), V ) and u ∈ L2 ((0, T ), H) , we also have u ∈

C ([0, T ], H) (see [12, p. 287]), so that the initial condition is meaningful in

the classical sense.

In what follows, the bilinear form a is assumed to be continuous on V —V

(see (3.2)) and V -elliptic (see (3.3)). The latter means that there exists a

number ± > 0 such that

a(v, v) ≥ ± v for all v ∈ V .