‚x

Show that the partial di¬erential equation does not change its type if it

is written in the new coordinates ξ = ¦(x).

1

For the Beginning:

The Finite Di¬erence Method for the

Poisson Equation

1.1 The Dirichlet Problem for the Poisson

Equation

In this section we want to introduce the ¬nite di¬erence method using

the Poisson equation on a rectangle as an example. By means of this ex-

ample and generalizations of the problem, advantages and limitations of

the approach will be elucidated. Also, in the following section the Poisson

equation will be the main topic, but then on an arbitrary domain. For the

spatial basic set of the di¬erential equation „¦ ‚ Rd we assume as minimal

requirement that „¦ is a domain, where a domain is a nonempty, open, and

connected set. The boundary of this domain will be denoted by ‚„¦, the

closure „¦ ∪ ‚„¦ by „¦ (see Appendix A.2). The Dirichlet problem for the

Poisson equation is then de¬ned as follows: Given functions g : ‚„¦ ’ R

and f : „¦ ’ R, we are looking for a function u : „¦ ’ R such that

d

‚2

’ u = f in „¦ , (1.1)

‚x2

i

i=1

u = g on ‚„¦ . (1.2)

This di¬erential equation model has already appeared in (0.19) and

(0.38) and beyond this application has an importance in a wide spectrum of

disciplines. The unknown function u can be interpreted as an electromag-

netic potential, a displacement of an elastic membrane, or a temperature.

Similar to the multi-index notation to be introduced in (2.16) (but with

20 1. Finite Di¬erence Method for the Poisson Equation

indices at the top) from now on for partial derivatives we use the following

notation.

Notation: For u : „¦ ‚ Rd ’ R we set

‚

‚i u := u for i = 1, . . . , d ,

‚xi

‚2

for i, j = 1, · · · , d ,

‚ij u := u

‚xi ‚xj

∆u := (‚11 + . . . + ‚dd ) u .

The expression ∆u is called the Laplace operator. By means of this, (1.1)

can be written in abbreviated form as

’∆u = f in „¦ . (1.3)

We could also de¬ne the Laplace operator by

∆u = ∇ · (∇u) ,

T

where ∇u = (‚1 u, . . . , ‚d u) denotes the gradient of a function u, and

∇ · v = ‚1 v1 + · · · + ‚d vd the divergence of a vector ¬eld v. Therefore,

an alternative notation exists, which will not be used in the following:

∆u = ∇2 u. The incorporation of the minus sign in the left-hand side of

(1.3), which looks strange at ¬rst glance, is related to the monotonicity and

de¬niteness properties of ’∆ (see Sections 1.4 and 2.1, respectively).

The notion of a solution for (1.1), (1.2) still has to speci¬ed more pre-

cisely. Considering the equations in a pointwise sense, which will be pursued

in this chapter, the functions in (1.1), (1.2) have to exist, and the equations

have to be satis¬ed pointwise. Since (1.1) is an equation on an open set „¦,

there are no implications for the behaviour of u up to the boundary ‚„¦. To

have a real requirement due to the boundary condition, u has to be at least

continuous up to the boundary, that is, on „¦. These requirements can be

formulated in a compact way by means of corresponding function spaces.

The function spaces are introduced more precisely in Appendix A.5. Some

examples are

u : „¦ ’ R u continuous in „¦ ,

C(„¦) :=

u : „¦ ’ R u ∈ C(„¦) , ‚i u exists in „¦ ,

C 1 („¦) :=

‚i u ∈ C(„¦) for all i = 1, . . . , d .

The spaces C k („¦) for k ∈ N, C(„¦), and C k („¦), as well as C(‚„¦), are

de¬ned analogously. In general, the requirements related to the (contin-

uous) existence of derivatives are called, a little bit vaguely, smoothness

requirements.

In the following, in view of the ¬nite di¬erence method, f and g will also

be assumed continuous in „¦ and ‚„¦, respectively.

1.2. The Finite Di¬erence Method 21

De¬nition 1.1 Assume f ∈ C(„¦) and g ∈ C(‚„¦). A function u is called

a (classical) solution of (1.1), (1.2) if u ∈ C 2 („¦) © C(„¦), (1.1) holds for all

x ∈ „¦, and (1.2) holds for all x ∈ ‚„¦.

1.2 The Finite Di¬erence Method

The ¬nite di¬erence method is based on the following approach: We are

looking for an approximation to the solution of a boundary value problem

at a ¬nite number of points in „¦ (the grid points). For this reason we

substitute the derivatives in (1.1) by di¬erence quotients, which involve

only function values at grid points in „¦ and require (1.2) only at grid

points. By this we obtain algebraic equations for the approximating values

at grid points. In general, such a procedure is called the discretization of the

boundary value problem. Since the boundary value problem is linear, the

system of equations for the approximate values is also linear. In general, for

other (di¬erential equation) problems and other discretization approaches

we also speak of the discrete problem as an approximation of the continuous

problem. The aim of further investigations will be to estimate the resulting

error and thus to judge the quality of the approximative solution.

Generation of Grid Points

In the following, for the beginning, we will restrict our attention to problems

in two space dimensions (d = 2). For simpli¬cation we consider the case

of a constant step size (or mesh width) h > 0 in both space directions.

The quantity h here is the discretization parameter, which in particular

determines the dimension of the discrete problem.

—¦—¦—¦ —¦ —¦ —¦ —¦ —¦ —¦ • : „¦h

—¦ •• • • • • • —¦

3

—¦ : ‚„¦h

—¦ •• • • • • • —¦

l=8 2

—¦ •• • • • • • —¦

2 : far from boundary

m=5

—¦ •• • • • • • —¦

—¦—¦—¦ —¦ —¦ —¦ —¦ —¦ —¦ 3 : close to boundary

Figure 1.1. Grid points in a square domain.

For the time being, let „¦ be a rectangle, which represents the simplest

case for the ¬nite di¬erence method (see Figure 1.1). By translation of the

coordinate system the situation can be reduced to „¦ = (0, a) — (0, b) with

a, b > 0 . We assume that the lengths a, b, and h are such that

for certain l, m ∈ N.

a = lh, b = mh (1.4)

22 1. Finite Di¬erence Method for the Poisson Equation

We de¬ne

(ih, jh) i = 1, . . . , l ’ 1 , j = 1, . . . , m ’ 1

„¦h :=

(1.5)

(x, y) ∈ „¦ x = ih , y = jh with i, j ∈ Z

=

as a set of grid points in „¦ in which an approximation of the di¬erential

equation has to be satis¬ed. In the same way,

(ih, jh) i ∈ {0, l} , j ∈ {0, . . . , m} or i ∈ {0, . . . , l} , j ∈ {0, m}

‚„¦h :=

(x, y) ∈ ‚„¦ x = ih , y = jh with i, j ∈ Z

=

de¬nes the grid points on ‚„¦ in which an approximation of the boundary

condition has to be satis¬ed. The union of grid points will be denoted by

„¦h := „¦h ∪ ‚„¦h .

Setup of the System of Equations

Lemma 1.2 Let „¦ := (x ’ h, x + h) for x ∈ R, h > 0. Then there exists

a quantity R, depending on u and h, the absolute value of which can be

bounded independently of h and such that

(1) for u ∈ C 2 („¦):

u(x + h) ’ u(x) 1

|R| ¤

u (x) = + hR and u ,

∞

h 2

(2) for u ∈ C 2 („¦):

u(x) ’ u(x ’ h) 1

|R| ¤

u (x) = + hR and u ,

∞

h 2

(3) for u ∈ C 3 („¦):

u(x + h) ’ u(x ’ h) 1

|R| ¤

+ h2 R

u (x) = and u ,

∞

2h 6

(4) for u ∈ C 4 („¦):

u(x + h) ’ 2u(x) + u(x ’ h) 1 (4)

+ h2 R and |R| ¤

u (x) = u .

∞

2

h 12

Here the maximum norm · ∞ (see Appendix A.5) has to be taken over

the interval of the involved points (x, x + h), (x ’ h, x), or (x ’ h, x + h).

Proof: The proof follows immediately by Taylor expansion. As an example

we consider statement 3: From

h2 h3

u(x ± h) = u(x) ± hu (x) + u (x) ± u (x ± ξ± ) for certain ξ± ∈ (0, h)

2 6

2

the assertion follows by linear combination.

1.2. Derivation and Properties 23

Notation: The quotient in statement 1 is called the forward di¬erence

quotient, and it is denoted by ‚ + u(x). The quotient in statement 2 is

called the backward di¬erence quotient (‚ ’ u(x)), and the one in statement

3 the symmetric di¬erence quotient (‚ 0 u(x)). The quotient appearing in

statement 4 can be written as ‚ ’ ‚ + u(x) by means of the above notation.

In order to use statement 4 in every space direction for the approximation

of ‚11 u and ‚22 u in a grid point (ih, jh), in addition to the conditions of

De¬nition 1.1, the further smoothness properties ‚ (3,0) u, ‚ (4,0) u ∈ C(„¦)

and analogously for the second coordinate are necessary. Here we use, e.g.,

the notation ‚ (3,0) u := ‚ 3 u/‚x3 (see (2.16)).

1

Using these approximations for the boundary value problem (1.1), (1.2),

at each grid point (ih, jh) ∈ „¦h we get

u ((i + 1)h, jh) ’ 2u(ih, jh) + u ((i ’ 1)h, jh)

’

h2

u (ih, (j + 1)h) ’ 2u(ih, jh) + u (ih, (j ’ 1)h)

+ = (1.6)

h2

f (ih, jh) + R(ih, jh)h2 .

=

Here R is as described in statement 4 of Lemma 1.2, a function depending

on the solution u and on the step size h, but the absolute value of which can

be bounded independently of h. In cases where we have less smoothness of

the solution u, we can nevertheless formulate the approximation (1.6) for

’∆u, but the size of the error in the equation is unclear at the moment.

For the grid points (ih, jh) ∈ ‚„¦h no approximation of the boundary

condition is necessary:

u(ih, jh) = g(ih, jh) .

If we neglect the term Rh2 in (1.6), we get a system of linear equations

for the approximating values uij for u(x, y) at points (x, y) = (ih, jh) ∈ „¦h .

They have the form

1

’ ui,j’1 ’ ui’1,j + 4uij ’ ui+1,j ’ ui,j+1 = fij (1.7)

h2

for i = 1, . . . , l ’ 1 , j = 1, . . . , m ’ 1 ,

if i ∈ {0, l}, j = 0, . . . , m or j ∈ {0, m}, i = 0, . . . , l .

uij = gij (1.8)

Here we used the abbreviations

fij := f (ih, jh), gij := g(ih, jh) . (1.9)

Therefore, for each unknown grid value uij we get an equation. The grid

points (ih, jh) and the approximating values uij located at these have a

natural two-dimensional indexing.

In equation (1.7) for a grid point (i, j) only the neighbours at the four

cardinal points of the compass appear, as it is displayed in Figure 1.2. This

24 1. Finite Di¬erence Method for the Poisson Equation

interconnection is also called the ¬ve-point stencil of the di¬erence method

and the method the ¬ve-point stencil discretization.

y

T (i,j+1)

•

(i’1,j) (i,j) (i+1,j)

• • •

(i,j’1)

•

Ex

Figure 1.2. Five-point stencil.

At the interior grid points (x, y) = (ih, jh) ∈ „¦h , two cases can be

distinguished:

(1) (i, j) has a position such that its all neighbouring grid points lie in

„¦h (far from the boundary).

(2) (i, j) has a position such that at least one neighbouring grid point

(r, s) lies on ‚„¦h (close to the boundary). Then in equation (1.7) the

value urs is known due to (1.8) (urs = grs ), and (1.7) can be modi¬ed

in the following way:

Remove the values urs with (rh, sh) ∈ ‚„¦h in the equations for (i, j)

close to the boundary and add the value grs /h2 to the right-hand

side of (1.7). The set of equations that arises by this elimination of

boundary unknowns by means of Dirichlet boundary conditions we

call (1.7)— ; it is equivalent to (1.7), (1.8).

Instead of considering the values uij , i = 1, . . . , l ’ 1, j = 1, . . . , m ’ 1,

one also speaks of the grid function uh : „¦h ’ R, where uh (ih, jh) = uij

for i = 1, . . . , l ’ 1, j = 1, . . . , m ’ 1. Grid functions on ‚„¦h or on „¦h are

de¬ned analogously. Thus we can formulate the ¬nite di¬erence method

in the following way: Find a grid function uh on „¦h such that equations

(1.7), (1.8) hold, or, equivalently ¬nd a grid function uh on „¦h such that

equations (1.7)— hold.

Structure of the System of Equations

After choosing an ordering of the uij for i = 0, . . . , l, j = 0, . . . , m, the

system of equations (1.7)— takes the following form:

Ah uh = q h (1.10)

with Ah ∈ RM1 ,M1 and uh , q h ∈ RM1 , where M1 = (l ’ 1)(m ’ 1).

This means that nearly identical notations for the grid function and its

representing vector are chosen for a ¬xed numbering of the grid points.

The only di¬erence is that the representing vector is printed in bold. The

ordering of the grid points may be arbitrary, with the restriction that the

1.2. Derivation and Properties 25

points in „¦h are enumerated by the ¬rst M1 indices, and the points in ‚„¦h

are labelled with the subsequent M2 = 2(l + m) indices. The structure of

Ah is not in¬‚uenced by this restriction.

Because of the described elimination process, the right-hand side q h has

the following form:

q h = ’Ah g + f ,

ˆ (1.11)

where g ∈ RM2 and f ∈ RM1 are the vectors representing the grid functions

f h : „¦h ’ R and gh : ‚„¦h ’ R

according to the chosen numbering with the values de¬ned in (1.9). The

matrix Ah ∈ RM1 ,M2 has the following form:

ˆ

±

’1

if the node i is close to the boundary

h2

ˆh )ij = and j is a neighbour in the ¬ve-point stencil,

(A

0 otherwise .

(1.12)

For any ordering, only the diagonal element and at most four further entries

in a row of Ah , de¬ned by (1.7), are di¬erent from 0; that is, the matrix is

sparse in a strict sense, as is assumed in Chapter 5.

An obvious ordering is the rowwise numbering of „¦h according to the

following scheme:

··· ···

(h,b’h) (2h,b’h) (a’h,b’h)

(l’1)(m’2)+1 (l’1)(m’2)+2 (l’1)(m’1)

··· ···

(h,b’2h) (2h,b’2h) (a’h,b’2h)

(l’1)(m’3)+1 (l’1)(m’3)+2 (l’1)(m’2)

. . .

.. ..

. . . . (1.13)

. .

. . .

··· ···

(h,2h) (2h,2h) (a’h,2h)

l l+1 2l’2

··· ···

(h,h) (2h,h) (a’h,h)

1 2 l’1

Another name of the above scheme is lexicographic ordering. (However,

this name is better suited to the columnwise numbering.)

In this case the matrix Ah has the following form of an (m ’ 1) — (m ’ 1)

block tridiagonal matrix:

«

T ’I

¬ ’I T ’I ·

0

¬ ·

¬ ·

.. .. ..

¬ ·

. . .

Ah = h’2 ¬ · (1.14)

¬ ·

.. .. ..

¬ ·

. . .

¬ ·

’I T ’I

0

’I T

26 1. Finite Di¬erence Method for the Poisson Equation

with the unit matrix I ∈ Rl’1,l’1 and

«

4 ’1

¬ ’1 4 ’1 ·

0

¬ ·

¬ ·

.. .. ..

¬ ·

. . .

T =¬ · ∈ Rl’1,l’1 .

¬ ·

.. .. ..

¬ ·

. . .

¬ ·

4 ’1

’1

0

’1 4

We return to the consideration of an arbitrary numbering. In the fol-

lowing we collect several properties of the matrix Ah ∈ RM1 ,M1 and the

extended matrix

Ah := Ah Ah ∈ RM1 ,M ,

˜ ˆ

˜

where M := M1 + M2 . The matrix Ah takes into account all the grid

points in „¦h . It has no relevance with the resolution of (1.10), but with the

stability of the discretization, which will be investigated in Section 1.4.

• (Ah )rr > 0 for all r = 1, . . . , M1 ,

• (Ah )rs ¤ 0 for all r = 1, . . . , M1 , s = 1, . . . , M such that r = s,

˜

±

≥ 0 for all r = 1, . . . , M1 ,

M1

• (Ah )rs (1.15)

> 0 if r belongs to a grid point close to

s=1 the boundary,

M

• ˜

(Ah )rs = 0 for all r = 1, . . . , M1 ,

s=1

• Ah is irreducible ,

• Ah is regular.

Therefore, the matrix Ah is weakly row diagonally dominant (see Ap-

pendix A.3 for de¬nitions from linear algebra). The irreducibility follows

from the fact that two arbitrary grid points may be connected by a path