<< . .

. 4
( : 59)



. . >>

nonsingular Jacobi matrix D¦ := ‚¦ .
‚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

<< . .

. 4
( : 59)



. . >>