<< . .

. 5
( : 59)



. . >>

consisting of corresponding neighbours in the ¬ve-point stencil. The reg-
ularity follows from the irreducible diagonal dominance. From this we
can conclude that (1.10) can be solved by Gaussian elimination without
pivot search. In particular, if the matrix has a band structure, this will be
preserved. This fact will be explained in more detail in Section 2.5.
The matrix Ah has the following further properties:
• Ah is symmetric,
• Ah is positive de¬nite.
It is su¬cient to verify these properties for a ¬xed ordering, for example the
rowwise one, since by a change of the ordering matrix, Ah is transformed
to P Ah P T with some regular matrix P, by which neither symmetry nor
1.2. Derivation and Properties 27

positive de¬niteness is destroyed. Nevertheless, the second assertion is not
obvious. One way to verify it is to compute eigenvalues and eigenvectors
explicitly, but we refer to Chapter 2, where the assertion follows naturally
from Lemma 2.13 and (2.36). The eigenvalues and eigenvectors are speci¬ed
in (5.24) for the special case l = m = n and also in (7.60). Therefore, (1.10)
can be resolved by Cholesky™s method, taking into account the bandedness.

Quality of the Approximation by the Finite Di¬erence Method
We now address the following question: To what accuracy does the grid
function uh corresponding to the solution uh of (1.10) approximate the
solution u of (1.1), (1.2)?
To this end we consider the grid function U : „¦h ’ R, which is de¬ned
by
U (ih, jh) := u(ih, jh). (1.16)
To measure the size of U ’ uh , we need a norm (see Appendix A.4 and also
A.5 for the subsequently used de¬nitions). Examples are the maximum
norm
uh ’ U |(uh ’ U ) (ih, jh)|
:= max (1.17)

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


and the discrete L2 -norm
1/2
l’1 m’1
2
uh ’ U ((uh ’ U )(ih, jh))
:= h . (1.18)
0,h
i=1 j=1

Both norms can be conceived as the application of the continuous norms
· ∞ of the function space L∞ („¦) or · 0 of the function space L2 („¦)
to piecewise constant prolongations of the grid functions (with a special
treatment of the area close to the boundary). Obviously, we have

vh 0,h ¤ ab vh ∞
for a grid function vh , but the reverse estimate does not hold uniformly in
h, so that · ∞ is a stronger norm. In general, we are looking for a norm
· h in the space of grid functions in which the method converges in the
sense
uh ’ U ’0 for h ’ 0
h

or even has an order of convergence p > 0, by which we mean the existence
of a constant C > 0 independent of h such that
uh ’ U ¤ C hp .
h

Due to the construction of the method, for a solution u ∈ C 4 („¦) we have
Ah U = q h + h2 R ,
28 1. Finite Di¬erence Method for the Poisson Equation

where U and R ∈ RM1 are the representations of the grid functions U and
R according to (1.6) in the selected ordering. Therefore, we have:
Ah (uh ’ U ) = ’h2 R
and thus
|Ah (uh ’ U )|∞ = h2 |R|∞ = Ch2
with a constant C(= |R|∞ ) > 0 independent of h.
From Lemma 1.2, 4. we conclude that
1
‚ (4,0) u + ‚ (0,4) u
C= .
∞ ∞
12
That is, for a solution u ∈ C 4 („¦) the method is consistent with the bound-
ary value problem with an order of consistency 2. More generally, the notion
takes the following form:
De¬nition 1.3 Let (1.10) be the system of equations that corresponds to
a (¬nite di¬erence) approximation on the grid points „¦h with a discretiza-
tion parameter h. Let U be the representation of the grid function that
corresponds to the solution u of the boundary value problem according to
(1.16). Furthermore, let · h be a norm in the space of grid functions
on „¦h , and let | · |h be the corresponding vector norm in the space RM1 h ,
where M1 h is the number of grid points in „¦h . The approximation is called
consistent with respect to · h if
|Ah U ’ q h |h ’ 0 for h ’ 0 .
The approximation has the order of consistency p > 0 if
|Ah U ’ q h |h ¤ Chp
with a constant C > 0 independent of h.
Thus the consistency or truncation error Ah U ’ qh measures the quality
of how the exact solution satis¬es the approximating equations. As we have
seen, in general it can be determined easily by Taylor expansion, but at
the expense of unnaturally high smoothness assumptions. But one has to
be careful in expecting the error |uh ’ U |h to behave like the consistency
error. We have
= A’1 Ah (uh ’ U ) ¤ A’1
uh ’ U Ah (uh ’ U ) , (1.19)
h h
h h h h

where the matrix norm · h has to be chosen to be compatible with the
vector norm |·|h . The error behaves like the consistency error asymptotically
in h if A’1 h can be bounded independently of h; that is if the method
h
is stable in the following sense:
De¬nition 1.4 In the situation of De¬nition 1.3, the approximation is
called stable with respect to · h if there exists a constant C > 0
Exercises 29

independent of h such that
A’1 ¤C.
h h

From the above de¬nition we can obviously conclude, with (1.19), the
following result:
Theorem 1.5 A consistent and stable method is convergent, and the order
of convergence is at least equal to the order of consistency.
Therefore, speci¬cally for the ¬ve-point stencil discretization of (1.1),
(1.2) on a rectangle, stability with respect to · ∞ is desirable. In fact, it
follows from the structure of Ah : Namely, we have
12
A’1 ¤ (a + b2 ) . (1.20)

h
16
This follows from more general considerations in Section 1.4 (Theo-
rem 1.14). Putting the results together we have the following theorem:
Theorem 1.6 Let the solution u of (1.1), (1.2) on a rectangle „¦ be
in C 4 („¦). Then the ¬ve-point stencil discretization has an order of
convergence 2 with respect to · ∞ , more precisely,
1
|uh ’ U |∞ ¤ (a2 + b2 ) ‚ (4,0) u + ‚ (0,4) u h2 .
∞ ∞
192


Exercises
1.1 Complete the proof of Lemma 1.2 and also investigate the error of
the respective di¬erence quotients, assuming only u ∈ C 2 [x ’ h, x + h].

1.2 Generalize the discussion concerning the ¬ve-point stencil discretiza-
tion (including the order of convergence) of (1.1), (1.2) on a rectangle for
h1 > 0 in the x1 direction and h2 > 0 in the x2 direction.

1.3 Show that an irreducible weakly row diagonally dominant matrix
cannot have vanishing diagonal elements.



1.3 Generalizations and Limitations of the Finite
Di¬erence Method
We continue to consider the boundary value problem (1.1), (1.2) on a rect-
angle „¦. The ¬ve-point stencil discretization developed may be interpreted
as a mapping ’∆h from functions on „¦h into grid functions on „¦h , which
30 1. Finite Di¬erence Method for the Poisson Equation

is de¬ned by
1
’∆h vh (x1 , x2 ) := cij vh (x1 + ih, x2 + jh) , (1.21)
i,j=’1

where c0,0 = 4/h2 , c0,1 = c1,0 = c0,’1 = c’1,0 = ’1/h2 , and cij = 0 for
all other (i, j). For the description of such a di¬erence stencil as de¬ned
in (1.21) the points of the compass (in two space dimensions) may also
be involved. In the ¬ve-point stencil only the main points of the compass
appear.
The question of whether the weights cij can be chosen di¬erently such
that we gain an approximation of ’∆u with higher order in h has to be
answered negatively (see Exercise 1.7). In this respect the ¬ve-point stencil
is optimal. This does not exclude that other di¬erence stencils with more
entries, but of the same order of convergence, might be worthwhile to con-
sider. An example, which will be derived in Exercise 3.11 by means of the
¬nite element method, has the following form:
8 1
c0,0 = 2 , cij = ’ 2 for all other i, j ∈ {’1, 0, 1} . (1.22)
3h 3h
This nine-point stencil can be interpreted as a linear combination of the
¬ve-point stencil and a ¬ve-point stencil for a coordinate system rotated by

π
(with step size 2 h), using the weights 1 and 2 in this linear combina-
4 3 3
tion. Using a general nine-point stencil a method with order of consistency
greater than 2 can be constructed only if the right-hand side f at the point
(x1 , x2 ) is approximated not by the evaluation f (x1 , x2 ), but by applying
a more general stencil. The mehrstellen method (“Mehrstellenverfahren”)
de¬ned by Collatz is such an example (see, for example, [15, p. 66]).
Methods of higher order can be achieved by larger stencils, meaning
that the summation indices in (1.21) have to be replaced by k and ’k,
respectively, for k ∈ N. But already for k = 2 such di¬erence stencils
cannot be used for grid points close to the boundary, so that there one has
to return to approximations of lower order.
If we consider the ¬ve-point stencil to be a suitable discretization for
the Poisson equation, the high smoothness assumption for the solution in
Theorem 1.6 should be noted. This requirement cannot be ignored, since
in general it does not hold true. On the one hand, for a smoothly bounded
domain (see Appendix A.5 for a de¬nition of a domain with C l -boundary)
the smoothness of the solution is determined only by the smoothness of the
data f and g (see for example [13, Theorem 6.19]), but on the other hand,
corners in the domain reduce this smoothness the more, the more reentrant
the corners are. Let us consider the following examples:
For the boundary value problem (1.1), (1.2) on a rectangle (0, a) — (0, b)
we choose f = 1 and g = 0; this means arbitrarily smooth functions.
Nevertheless, for the solution u, the statement u ∈ C 2 („¦) cannot hold,
because otherwise, ’∆u(0, 0) = 1 would be true, but on the other hand,
1.3. Generalizations and Limitations 31

we have ‚1,1 u(x, 0) = 0 because of the boundary condition and hence also
‚1,1 u(0, 0) = 0 and ‚2,2 u(0, y) = 0 analogously. Therefore, ‚2,2 u(0, 0) = 0.
Consequently, ’∆u(0, 0) = 0, which contradicts the assumption above.
Therefore, Theorem 1.6 is not applicable here.
In the second example we consider the domain with reentrant corner (see
Figure 1.3)
„¦ = (x, y) ∈ R2 x2 + y 2 < 1 , x < 0 or y > 0 .
In general, if we identify R2 and C, this means (x, y) ∈ R2 and z = x + iy ∈
C, we have that if w : C ’ C is analytic (holomorphic), then both the real
and the imaginary parts w, w : C ’ R are harmonic, which means that
they solve ’∆u = 0.

y

„¦

x


Figure 1.3. Domain „¦ with reentrant corner.

We choose w(z) := z 2/3 . Then the function u(x, y) := (x + iy)2/3
solves the equation
’∆u = 0 in „¦ .
In polar coordinates, x = r cos •, y = r sin •, the function u takes the form
2
2/3
rei• = r2/3 sin
u(x, y) = • .
3
Therefore, u satis¬es the boundary conditions
2 3π
for 0 ¤ • ¤
u ei• = sin • , (1.23)
3 2
u(x, y) = 0 otherwise on ‚„¦ .
But note that w (z) = 2 z ’1/3 is unbounded for z ’ 0, so that ‚1 u, ‚2 u
3
are unbounded for (x, y) ’ 0. Therefore, in this case we do not even have
u ∈ C 1 („¦).
The examples do not show that the ¬ve-point stencil discretization is not
suitable for the boundary value problems considered, but they show the ne-
cessity of a theory of convergence, which requires only as much smoothness
as was to be expected.
In the following we discuss some generalizations of the boundary value
problems considered so far.
32 1. Finite Di¬erence Method for the Poisson Equation

General Domains „¦
We continue to consider (1.1), (1.2) but on a general domain in R2 , for
which the parts of the boundary are not necessarily aligned to the coor-
dinate axes. Therefore we can keep the second equation in (1.5) as the
de¬nition of „¦h , but have to rede¬ne the set of boundary grid points ‚„¦h .
For example, if for some point (x, y) ∈ „¦h we have
(x ’ h, y) ∈ „¦ ,
/
then there exists a number s ∈ (0, 1] such that
(x ’ ‘h, y) ∈ „¦ for all ‘ ∈ [0, s) (x ’ sh, y) ∈ „¦ .
and /
Then (x ’ sh, y) ∈ ‚„¦, and therefore we de¬ne
(x ’ sh, y) ∈ ‚„¦h .
The other main points of the compass are treated analogously. In this
way the grid spacing in the vicinity of the boundary becomes variable; in
particular, it can be smaller than h.
For the quality of the approximation we have the following result:
Lemma 1.7 Let „¦ = (x ’ h1 , x + h2 ) for x ∈ R, h1 , h2 > 0.
(1) Then for u ∈ C 3 („¦),
u(x + h2 ) ’ u(x) u(x) ’ u(x ’ h1 )
2

u (x) =
h1 + h2 h2 h1
+ max {h1 , h2 } R ,
where R is bounded independently of h1 , h2 .
(2) There are no ±, β, γ ∈ R such that
u (x) = ± u(x ’ h1 ) + β u(x) + γ u(x + h2 ) + R1 h2 + R2 h2
1 2

for all polynomials u of degree 3 if h1 = h2 .

2
Proof: Exercises 1.4 and 1.5.
This leads to a discretization that is di¬cult to set up and for which the
order of consistency and order of convergence are not easily determined.

Other Boundary Conditions
We want to consider the following example. Let ‚„¦ = “1 ∪ “3 be divided
into two disjoint subsets. We are looking for a function u such that
’∆u =f in „¦ ,
‚ν u := ∇u · ν =g on “1 , (1.24)
u =0 on “3 ,
where ν : ‚„¦ ’ Rd is the outer unit normal, and thus ‚ν u is the normal
derivative of u.
1.3. Generalizations and Limitations 33

For a part of the boundary oriented in a coordinate direction, ‚ν u is
just a positive or negative partial derivative. But if only grid points in
„¦h are to be used, only ±‚ + u and ±‚ ’ u respectively (in the coordinates
orthogonal to the direction of the boundary) are available directly from
the above approximations with a corresponding reduction of the order of
consistency. For a boundary point without these restrictions the question
of how to approximate ‚ν u appropriately is open.
As an example we consider (1.24) for a rectangle „¦ = (0, a)—(0, b), where
“1 := {(a, y) | y ∈ (0, b)} , “3 := “ \ “1 . (1.25)
At the boundary grid points (a, jh), j = 1, . . . , m ’ 1, ‚2 u = ∇u · ν
is prescribed, which can be approximated directly only by ‚ ’ u. Due to
Lemma 1.2, 2 this leads to a reduction in the consistency order (see Ex-
ercise 1.8). The resulting system of equations may include the Neumann
boundary grid points in the set of unknowns, for which an equation with
the entries 1/h in the diagonal and ’1/h in an o¬-diagonal corresponding
to the eastern neighbour (a ’ h, jh) has to be added. Alternatively, those
boundary points can be eliminated, leading for the eastern neighbour to a
modi¬ed di¬erence stencil (multiplied by h2 )
’1
’1 3 (1.26)
’1

for the right-hand side h2 f (a ’ h, jh) + hg(a, jh). In both cases the matrix
properties of the system of equations as collected in (1.15) still hold, with
M1
the exception of s=1 (Ah )rs = 0, both for the Neumann boundary points
and their neighbours, if no Dirichlet boundary point is involved in their
stencil. Thus the term “close to the boundary” has to be interpreted as
“close to the Dirichlet boundary.”
If one wants to take advantage of the symmetric di¬erence quotient ‚ 0 u,
then “arti¬cial” values at new external grid points (a + h, jh) appear.
To keep the balance of unknowns and equations, it can be assumed that
the di¬erential equation also holds at (a, jh), and thus it is discretized
with the ¬ve-point stencil there. If one attributes the discrete boundary
condition to the external grid point, then again the properties (1.15) hold
with the abovementioned interpretation. Alternatively, the external grid
points can be eliminated, leading to a modi¬ed di¬erence stencil (multiplied
by h2 ) at (a, jh):
’1
’2 4 (1.27)
’1

for the right-hand side h2 f (a, jh)+2hg(a, jh), with the same interpretation
of properties (1.15).
34 1. Finite Di¬erence Method for the Poisson Equation

More General Di¬erential Equations
As an example we consider the di¬erential equation
’∇ · (k ∇u) = f on „¦ (1.28)
with a continuous coe¬cient function k : „¦ ’ R, which is bounded from
below by a positive constant on „¦. This equation states the conservation
of an extensive quantity u whose ¬‚ux is ’k∇u (see Section 0.5). This
should be respected by the discretization, and therefore the form of (1.28)
obtained by working out the derivatives is not recommended as a basis for
the discretization. The di¬erential expression in (1.28) can be discretized
by a successive application of central di¬erence quotients, but then again
the order of consistency has to be investigated.
In addition, one has to take into account the fact that the smoothness of
u depends on the smoothness of k. If processes in heterogeneous materials
have to be described, then k is often discontinuous. In the simplest example
k is assumed to take two di¬erent values: Let „¦ = „¦1 ∪ „¦2 and
k|„¦1 = k1 > 0 , k|„¦2 = k2 > 0
with constants k1 = k2 .
As worked out in Section 0.5, on the interior boundary S := „¦1 © „¦2 a
transmission condition has to be imposed:
• u is continuous,
• (k∇u) · ν is continuous, where ν is the outer normal on ‚„¦1 , for
example.
This leads to the following conditions on ui , being the restrictions of u on
„¦i for i = 1, 2:
’k1 ∆u1 =f in „¦1 , (1.29)
’k2 ∆u2 =f in „¦2 ,
u1 = u2 on S , (1.30)
k1 ‚ν u1 = k2 ‚ν u2 on S .
In this case the question of an appropriate discretization is also open.
Summarizing, we have the following catalogue of requirements: We are
looking for a notion of solution for (general) boundary value problems with
nonsmooth coe¬cients and right-hand sides such that, for example, the
transmission condition is ful¬lled automatically.
We are looking for a discretization on general domains such that, for
example, the (order of) convergence can also be assured for less smooth
solutions and also Neumann boundary conditions as in (1.24) can be treated
easily.
The ¬nite element method in the subsequent chapters will ful¬l these
requirements to a large extent.
Exercises 35

Exercises
1.4 Prove Lemma 1.7, 1.

1.5 Under the assumption that u : „¦ ‚ R ’ R is a su¬ciently smooth
function, determine in the ansatz
±u(x ’ h1 ) + βu(x) + γu(x + h2 ) , h1 , h2 > 0 ,
the coe¬cients ± = ±(h1 , h2 ), β = β(h1 , h2 ), γ = γ(h1 , h2 ), such that
(a) for x ∈ „¦, u (x) will be approximated with the order as high as
possible,
(b) for x ∈ „¦, u (x) will be approximated with the order as high as
possible,
and in particular, prove 1.7, 2.
Hint: Determine the coe¬cients such that the formula is exact for
polynomials with the degree as high as possible.

<< . .

. 5
( : 59)



. . >>