<<

. 2
( 16)



>>

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
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.

1.6 Let „¦ ‚ R2 be a bounded domain. For a su¬ciently smooth function
u : „¦ ’ R determine the di¬erence formula with an order as high as
possible to approximate ‚11 u(x1 , x2 ), using the 9 values u(x1 + γ1 h, x2 +
γ2 h), where γ1 , γ2 ∈ {’1, 0, 1}.

1.7 Let „¦ ‚ R2 be a bounded domain. Show that in (1.21) there exists
no choice of cij such that for an arbitrary smooth function u : „¦ ’ R,
|∆u(x) ’ ∆h u(x)| ¤ Ch3
is valid with a constant C independent of h.

1.8 For the example (1.24), (1.25), investigate the order of consistency
both for the discretization (1.26) and (1.27) in the maximum norm. Are
there improvements possible considering the discrete L2 -norm? (See (1.18).)

1.9 Consider example (1.24) with
“1 := {(a, y) | y ∈ (0, b)} ∪ {(x, b) | x ∈ (0, a]},
“3 := “ \ “1 ,
and discuss the applicability of the one-sided and the symmetric di¬er-
ence quotients for the approximation of the Neumann boundary condition,
in particular with respect to properties (1.15). In which way does the
boundary condition at (a, b), where no unique normal exists, have to be
interpreted?
36 1. Finite Di¬erence Method for the Poisson Equation

1.10 Generalize the discussion concerning the ¬ve-point stencil dis-
cretization (including the order of convergence) to the boundary value
problem
’∆u + ru = f in „¦,
u= g on ‚„¦,
for r > 0 and „¦ := (0, a) — (0, b). To approximate the reactive term ru, the
following schemes in the notation of (1.21) are to be used:
(a) c0,0 = 1, cij = 0 otherwise,
(b) c0,0 > 0, c0,1 , c1,0 , c0,’1 , c’1,0 ≥ 0, cij = 0 otherwise, and
1
i,j=’1 cij = 1 .




1.4 Maximum Principles and Stability
In this section the proof of the stability estimate (1.20), which is still miss-
ing, will be given. For this reason we develop a more general framework, in
which we will then also discuss the ¬nite element method (see Section 3.9)
and the time-dependent problems (see Section 7.5). The boundary value
problem (1.1), (1.2) satis¬es a (weak ) maximum principle in the following
sense: If f is continuous and f (x) ¤ 0 for all x ∈ „¦ (for short f ¤ 0), then
max u(x) ¤ max u(x) .
x∈‚„¦
x∈„¦

This maximum principle is also strong in the following sense: The maxi-
mum of u on „¦ can be attained in „¦ only if u is constant (see, for example,
[13], also for the following assertions). By exchanging u, f, g by ’u, ’f, ’g,
respectively, we see that there is an analogous (strong) minimum principle.
The same holds for more general linear di¬erential equations as in (1.28),
which may also contain convective parts (this means ¬rst-order deriva-
tives). But if the equation contains a reactive part (this means without
derivatives), as in the example
’∆u + ru = f in „¦
with a continuous function r : „¦ ’ R such that r(x) ≥ 0 for x ∈ „¦, there
is a weak maximum principle only in the following form: If f ¤ 0, then

max u(x) ¤ max max u(x), 0 .
x∈‚„¦
x∈„¦

The weak maximum principle directly implies assertions about the de-
pendence of the solution u of the boundary value problem on the data f
and g; this means stability properties. One can also follow this method in
investigating the discretization. For the basic example we have
1.4. Maximum Principles and Stability 37

Theorem 1.8 Let uh be a grid function on „¦h de¬ned by (1.7), (1.8) and
suppose fij ¤ 0 for all i = 1, . . . , l ’ 1, j = 1, . . . , m ’ 1. Then if uh attains
its maximum on „¦h ∪ ‚„¦— at a point (i0 h, j0 h) ∈ „¦h , then the following
h
holds:
uh is constant on „¦h ∪ ‚„¦— .
h

Here
‚„¦— := ‚„¦h \ {(0, 0), (a, 0), (0, b), (a, b)} .
h

In particular, we have
max uh (x, y) ¤ max uh (x, y) .
(x,y)∈‚„¦—
(x,y)∈„¦h h



Proof: Let u := uh (i0 h, j0 h). Then because of (1.7) and fij ¤ 0 we have
¯

4¯ ¤ uh (kh, lh) ¤ 4¯ ,
u u
(k,l)∈N(i0 ,j0 )

since in particular uh (kh, lh) ¤ u for (k, l) ∈ N(i0 ,j0 ) . Here we used the
¯
notation
N(i0 ,j0 ) = {((i0 ’ 1), j0 ), ((i0 + 1), j0 ), (i0 , (j0 + 1)), (i0 , (j0 ’ 1))}
for the set of indices of neighbours of (i0 h, j0 h) in the ¬ve-point stencil.
From these inequalities we conclude that
(k, l) ∈ N(i0 ,j0 ) .
uh (kh, lh) = u for
¯
If we apply this argument to the neighbours in „¦h of the grid points (kh, lh)
for (k, l) ∈ N(i0 ,j0 ) and then continue in the same way to the sets of neigh-
bours in „¦h arising in every such step, then ¬nally, for each grid point
(ih, jh) ∈ „¦h ∪ ‚„¦— the claimed identity uh (ih, jh) = u is achieved. 2
¯
h

The exceptional set of vertices ‚„¦h \ ‚„¦— does not participate in any
h
di¬erence stencil, so that the values there are of no relevance for uh .
We want to generalize this result and therefore consider a system of
equations as in (1.10), (1.11):
Ah uh = q h = ’Ah uh + f ,
ˆˆ (1.31)
where Ah ∈ RM1 ,M1 as in (1.10), Ah ∈ RM1 ,M2 as in (1.11), uh , f ∈ RM1 ,
ˆ
and uh ∈ RM2 . This may be interpreted as the discretization of a bound-
ˆ
ary value problem obtained by the ¬nite di¬erence method or any other
approach and without restrictions on the dimensionality of the domain.
At least on one part of the boundary Dirichlet boundary conditions are re-
quired. Then the entries of the vector uh can be interpreted as the unknown
(1) (1)
values at the grid points in „¦h ∪ ‚„¦h , where ‚„¦h correspond to a part
ˆ
of ‚„¦ (with ¬‚ux or mixed boundary condition). Analogously, the vector uh
38 1. Finite Di¬erence Method for the Poisson Equation

(indexed from M1 + 1 to M1 + M2 ) corresponds to the values ¬xed by the
(2)
Dirichlet boundary conditions on ‚„¦h . Again let M = M1 + M2 and

Ah := Ah Ah ∈ RM1 ,M .
˜ ˆ

This means in particular that the dimensions M1 and M2 are not ¬xed,
but are in general unbounded for h ’ 0.
Oriented on (1.15) we require the following general assumptions for the
rest of the section:
(1) (Ah )rr > 0 for all r = 1, . . . , M1 ,
(Ah )rs ¤ 0
(2) for all r, s = 1, . . . , M1 such that r = s ,
M1
(Ah )rs ≥ 0 for all r = 1, . . . , M1 ,
(3) (i)
s=1
(ii) for at least one index the strict inequality holds ,
(4) Ah is irreducible , (1.32)
(Ah )rs ¤ 0 for all r = 1, . . . , M1 , s = M1 + 1, . . . , M ,
ˆ
(5)
M
(Ah )rs ≥ 0 for all r = 1, . . . , M1 ,
˜
(6)
s=1
for every s = M1 + 1, . . . , M there exists r ∈ {1, . . . , M1 }
(7)
ˆ
such that (Ah )rs = 0.
Generalizing the notation above for r ∈ {1, . . . , M1 }, the indices s ∈
{1, . . . , M } \ {r} are called neighbours, for which (Ah )rs = 0, and they are
˜
assembled to form the set Nr . Therefore, the irreducibility of Ah means
that arbitrary r, s ∈ {1, . . . , M1 } can be connected by neighbourhood
relationships.
The condition (7) is not a restriction: It only avoids the inclusion of
known values (ˆ h )s that do not in¬‚uence the solution of (1.31) at all. For
u
the ¬ve-point stencil on the rectangle, these are the values at the corner
points. Because of the condition (7), every index r ∈ {M1 + 1, . . . , M }
is connected to every index s ∈ {1, . . . , M1 } by means of neighbourhood
relationships.
The conditions (2) and (3) imply the weak diagonal dominance of Ah .
Note that the conditions are formulated redundantly: The condition (3)
also follows from (5) through (7).
To simplify the notation we will use the following conventions, where u,
v and A, B are vectors and matrices, respectively, of suitable dimensions:
≥ ≥
u 0 if and only if (u)i 0 for all indices i ,
≥ if u ’ v ≥
u v if and only 0,
(1.33)
≥ ≥
A 0 if and only if (A)ij 0 for all indices (i, j) ,
≥ if A ’ B ≥
A B if and only 0.
1.4. Maximum Principles and Stability 39

Theorem 1.9 We consider (1.31) under the assumptions (1.32). Fur-
thermore, let f ¤ 0. Then a strong maximum principle holds: If the
u
components of uh = uh attain a nonnegative maximum for some in-
˜ ˆh
dex r ∈ {1, . . . , M1 }, then all the components are equal. In particular, a
weak maximum principle is ful¬lled:

(˜ h )r ¤ max 0,
max u max (ˆ h )r
u . (1.34)
r∈{1,...,M} r∈{M1 +1,...,M}


Proof: Let u = maxs∈{1,...,M} (˜ h )s , and u = (uh )r where r ∈
¯ u ¯
{1, . . . , M1 }. Because of (1.32) (2), (5), (6) the rth row of (1.31) implies

¤’ ˜ ˜
(Ah )rr u
¯ Ah (˜ h )s =
u Ah (˜ h )s
u
rs rs
s∈Nr s∈Nr
(1.35)
¤ u ¤ (Ah )rr u ,
˜
Ah ¯ ¯
rs
s∈Nr

where the assumption u ≥ 0 is used in the last estimate. Therefore, ev-
¯
erywhere equality has to hold. Since the second inequality is valid also
˜
for every single term and (Ah )rs = 0 by the de¬nition of Nr , we ¬nally
conclude that
for all s ∈ Nr .
(˜ h )s = u
u ¯
This allows us to apply this argument to all s ∈ Nr © {1, . . . , M1 }, then
to the corresponding sets of neighbours, and so on, until the assertion is
2
proven.

The requirement of irreducibility can be weakened if instead of (1.32) (6)
we have
M
— ˜
(6) Ah = 0 for all r = 1, . . . , M1 .
rs
s=1
Then condition (4) can be replaced by the requirement
(4)— For every r1 ∈ {1, . . . , M1 } such that
M1
(Ah )r1 s = 0 (1.36)
s=1
there are indices r2 , . . . , rl+1 such that
(Ah )ri ri+1 = 0 for i = 1, . . . , l
and
M1
(Ah )rl+1 s > 0 . (1.37)
s=1

These modi¬ed conditions without (7) will be denoted by (1.32) .
40 1. Finite Di¬erence Method for the Poisson Equation

Motivated by the example above we call a point r ∈ {1, . . . , M1 } far from
the boundary if (1.36) holds, and close to the boundary if (1.37) holds, and
the points r ∈ {M1 + 1, . . . , M } are called boundary points.
Theorem 1.10 We consider (1.31) under the assumption (1.32)— .
If f ¤ 0, then
(˜ h )r ¤
max u max (ˆ h )r .
u (1.38)
r∈{1,...,M} r∈{M1 +1,...,M}


Proof: We use the same notation and the same arguments as in the
proof of Theorem 1.9. In (1.35) in the last estimate equality holds, so that
no sign conditions for u are necessary. Because of (4)— the maximum will
¯
also be attained at a point close to the boundary and therefore also at
its neighbours. Because of (6)— a boundary point also belongs to these
2
neighbours, which proves the assertion.
From the maximum principles we immediately conclude a comparison
principle:
Lemma 1.11 We assume (1.32) or (1.32)— .
Let uh1 , uh2 ∈ RM1 be solutions of
Ah uhi = ’Ah uhi + f i
ˆˆ for i = 1, 2
for given f 1 , f 2 ∈ RM1 , uh1 , uh2 ∈ RM2 , which satisfy f 1 ¤ f 2 , uh1 ¤
ˆ ˆ ˆ
ˆ
uh2 . Then
uh1 ¤ uh2 .

Proof: From Ah (uh1 ’ uh2 ) = ’Ah (ˆ h1 ’ uh2 ) + f 1 ’ f 2 we can conclude
ˆu ˆ
with Theorem 1.9 or 1.10 that
(uh1 ’ uh2 )r ¤ 0 .
max
r∈{1,...,M1 }

2
This implies in particular the uniqueness of a solution of (1.31) for
ˆ
arbitrary uh and f and also the regularity of Ah .
In the following we denote by 0 and 0 the zero vector and the zero
matrix, respectively, where all components are equal to 0. An immediate
consequence of Lemma 1.11 is the following
Theorem 1.12 Let Ah ∈ RM1 ,M1 be a matrix with the properties (1.32)
(1)“(3) (i), (4)— , and uh ∈ RM1 . Then
Ah uh ≥ 0 uh ≥ 0 .
implies (1.39)

Proof: To be able to apply Lemma 1.11, one has to construct a matrix
Ah ∈ RM1 ,M2 such that (1.32)* holds. Obviously, this is possible. Then one
ˆ
1.4. Maximum Principles and Stability 41

can choose
ˆ
uh2 := uh , f 2 := Ah uh2 , uh2 := 0 ,
ˆ
uh1 := 0 , f 1 := 0 , uh1 := 0
ˆ
to conclude the assertion. Because of uhi := 0 for i = 1, 2 the speci¬c
2
ˆ
de¬nition of Ah plays no role.
A matrix with the property (1.39) is called inverse monotone. An
equivalent requirement is
A’1 v h ≥ 0 ,
vh ≥ 0 ’ h

and therefore by choosing the unit vectors as v h ,
A’1 ≥ 0 .
h

Inverse monotone matrices that also satisfy (1.32) (1), (2) are called M-
matrices.
Finally, we can weaken the assumptions for the validity of the comparison
principle.
Corollary 1.13 Suppose that Ah ∈ RM1 ,M1 is inverse monotone and
(1.32) (5) holds. Let uh1 , uh2 ∈ RM1 be solutions of
Ah uhi = ’Ah uhi + f i
ˆˆ for i = 1, 2
for given f 1 , f 2 ∈ RM1 , uh1 , uh2 ∈ RM2 that satisfy f 1 ¤ f 2 , uh1 ¤ uh2 .
ˆ ˆ ˆ ˆ
Then
uh1 ¤ uh2 .

Proof: Multiplying the equation
Ah (uh1 ’ uh2 ) = ’Ah (ˆ h1 ’ uh2 ) + f 1 ’ f 2
ˆu ˆ
from the left by the matrix A’1 , we get
h

uh1 ’ uh2 = ’ A’1 Ah (ˆ h1 ’ uh2 ) + A’1 (f 1 ’ f 2 ) ¤ 0 .
ˆu ˆ
h h
¤0 ¤0 ¤0
≥0 ≥0
2
The importance of Corollary 1.13 lies in the fact that there exist
˜
discretization methods, for which the matrix Ah does not satisfy, e.g., con-
’1
dition (1.32) (6), or (6)— but Ah ≥ 0. A typical example of such a method
is the ¬nite volume method described in Chapter 6.
In the following we denote by 1 a vector (of suitable dimension) whose
components are all equal to 1.
Theorem 1.14 We assume (1.32) (1)“(3), (4)— , (5). Furthermore, let
(1) (2)
wh , wh ∈ RM1 be given such that
(1) (2)
Ah w h ≥ 1 , Ah w h ≥ ’Ah 1 .
ˆ (1.40)
42 1. Finite Di¬erence Method for the Poisson Equation

Then a solution of Ah uh = ’Ah uh + f satis¬es
ˆˆ
(1) (2) (1) (2)
(1) ’ |f |∞ wh + |ˆ h |∞ wh ¤ uh ¤ |f |∞ wh + |ˆ h |∞ wh ,
u u
(1) (2)
(2) |uh |∞ ¤ wh |f |∞ + wh |ˆ h |∞ .
u
∞ ∞

Under the assumptions (1.32) (1)“(3), (4)— , and (1.40) the matrix norm
· ∞ induced by | · |∞ satis¬es

A’1
(1)
¤ wh .
∞ ∞
h


Proof: Since ’|f |∞ 1 ¤ f ¤ |f |∞ 1 and the analogous statement for uh
ˆ
(1) (2)
is valid, the vector v h := |f |∞ wh + |ˆ h |∞ wh ’ uh satis¬es
u

Ah v h ≥ |f |∞ 1 ’ f ’ Ah (|ˆ h |∞ 1 ’ uh ) ≥ 0 ,
ˆu ˆ
where we have also used ’Ah ≥ 0 in the last estimate. Therefore, the right
ˆ
inequality of (1) implies from Theorem 1.12 that the left inequality can be
proven analogously. The further assertions follow immediately from (1). 2
Because of the inverse monotonicity and from (1.32) (5) the vectors pos-
(i)
tulated in Theorem 1.14 have to satisfy wh ≥ 0 necessarily for i = 1, 2.
Thus stability with respect to · ∞ of the method de¬ned by (1.31) as-
suming (1.32) (1)“(3), (4)* is guaranteed if a vector 0 ¤ wh ∈ RM1 and a
constant C > 0 independent of h can be found such that
Ah wh ≥ 1 and |w h |∞ ¤ C . (1.41)
Finally, this will be proven for the ¬ve-point stencil discretization (1.1),
(1.2) on the rectangle „¦ = (0, a) — (0, b) for C = 16 (a2 + b2 ).
1

For this reason we de¬ne polynomials of second degree w1 , w2 by
1 1
x(a ’ x) y(b ’ y) .
w1 (x) := and w2 (y) := (1.42)
4 4
It is clear that w1 (x) ≥ 0 for all x ∈ [0, a] and w2 (y) ≥ 0 for all y ∈ [0, b].
Furthermore, we have w1 (0) = 0 = w1 (a) and w2 (0) = 0 = w2 (b), and
1 1
w1 (x) = ’and w2 (y) = ’ .
2 2
Therefore w1 and w2 are strictly concave and attain their maximum in a 2
b
and 2 , respectively. Thus the function w(x, y) := w1 (x) + w2 (x) satis¬es
’∆w =1 in „¦ ,
(1.43)
≥0
w on ‚„¦ .
Now let w h ∈ RM1 be, for a ¬xed ordering, the representation of the grid
function wh de¬ned by
for i = 1, . . . , l ’ 1 , j = 1, . . . , m ’ 1 .
(wh )(ih, jh) := w(ih, jh)
1.4. Maximum Principles and Stability 43

Analogously, let w h ∈ RM2 be the representation of the function wh de-
ˆ ˆ

¬ned on ‚„¦h . As can be seen from the error representation in Lemma 1.2,
statement 4, the di¬erence quotient ‚ ’ ‚ + u(x) is exact for polynomials of
second degree. Therefore, we conclude from (1.43) that
Ah w h = ’Ah w h + 1 ≥ 1 ,
ˆˆ

which ¬nally implies
a b 12
|w h |∞ = wh ¤w (a + b2 ) .
= w1 + w2 =
∞ ∞
2 2 16
This example motivates the following general procedure to construct wh ∈
RM1 and a constant C such that (1.41) is ful¬lled.
Assume that the boundary value problem under consideration reads in
an abstract form
x ∈ „¦,
(Lu)(x) = f (x) for
(1.44)
x ∈ ‚„¦ .
(Ru)(x) = g(x) for
Similar to (1.43) we can consider ” in case of existence ” a solution w
of (1.44) for some f, g, such that f (x) ≥ 1 for all x ∈ „¦, g(x) ≥ 0 for all
x ∈ „¦. If w is bounded on „¦, then
(wh )i := w(xi ), i = 1, . . . , M1 ,
for the (non-Dirichlet) grid points xi , is a candidate for wh . Obviously,
|w h |∞ ¤ w .


Correspondingly, we set
(w h )i = w(xi ) ≥ 0 ,
ˆ i = M1 + 1, . . . , M2 ,
for the Dirichlet-boundary grid points.
The exact ful¬llment of the discrete equations by wh cannot be expected
anymore, but in case of consistency the residual can be made arbitrarily
small for small h. This leads to
Theorem 1.15 Assume that a solution w ∈ C(„¦) of (1.44) exists for data
f ≥ 1 and g ≥ 0. If the discretization of the form (1.31) is consistent with
˜
(1.44) (for these data), and there exists H > 0 so that for some ± > 0 :
˜
’Ah wh + f ≥ ±1 h¤H ,
ˆˆ ˜
˜ for (1.45)
then for every 0 < ± < ± there exists H > 0, so that
˜
Ah w h ≥ ±1 h¤H.
for

Proof: Set
„ h := Ah w h + Ah w h ’ f
ˆˆ
44 1. Finite Di¬erence Method for the Poisson Equation

for the consistency error, then
|„ h |∞ ’ 0 for h ’ 0 .
Thus
„ h ’ Ah w h + f
ˆˆ
Ah wh =
≥ ’|„ h |∞ 1 + ±1 for h ¤ H
˜
˜
≥ h¤H
±1 for
2
and some appropriate H > 0.
Thus a proper choice in (1.41) is
1 1
wh and C := w ∞. (1.46)
± ±
The condition (1.45) is not critical: In case of Dirichlet boundary conditions
and (1.32) (5) (for corresponding rows i of Ah ) then, due to (f )i ≥ 1, we
ˆ
can even choose ± = 1. The discussion of Neumann boundary conditions
˜
following (1.24) shows that the same can be expected.
Theorem 1.15 shows that for a discretization with an inverse monotone
system matrix consistency already implies stability.
To conclude this section let us discuss the various ingredients of (1.32)
or (1.32)* that are su¬cient for a range of properties from the inverse
monotonicity up to a strong maximum principle: For the ¬ve-point stencil
on a rectangle all the properties are valid for Dirichlet boundary conditions.
If partly Neumann boundary conditions appear, the situation is the same,
but now close and far from the boundary refers to its Dirichlet part. In
the interpretation of the implications one has to take into account that the
heterogeneities of the Neumann boundary condition are now part of the
right-hand side f , as seen, e.g., in (1.26). If mixed boundary conditions are
applied, as
‚ν u + ±u = g on “2 (1.47)
for some “2 ‚ “ and ± = ±(x) > 0, then the situation is the same again
if ±u is approximated just by evaluation, at the cost that (4)* no longer
holds. The situation is similar if reaction terms appear in the di¬erential
equation (see Exercise 1.10).



Exercises
1.11 Give an example of a matrix Ah ∈ RM1 ,M2 that can be used in the
ˆ
proof of Theorem 1.12.

1.12 Show that the transposition of an M-matrix is again an M-matrix.
Exercises 45

1.13 In the assumptions of Theorem 1.9 substitute (1.32) (4) by (4)* and
amend (6) to
#
(6) Condition (1.32) (6) is valid and
M1
(Ah )rs > 0 ’ there exists s ∈ {M1 , . . . , M } such that (Ah )rs < 0.
ˆ
s=1

Under these conditions prove a weak maximum principle as in Theorem 1.9.

1.14 Assuming the existence of wh ∈ RM1 such that Ah wh ≥ 1 and
|wh |∞ ¤ C for some constant C independent of h, show directly (without
Theorem 1.14) a re¬ned order of convergence estimate on the basis of an
order of consistency estimate in which also the shape of wh appears.
2
The Finite Element Method
for the Poisson Equation




The ¬nite element method, frequently abbreviated by FEM, was devel-
oped in the ¬fties in the aircraft industry, after the concept had been
independently outlined by mathematicians at an earlier time. Even today
the notions used re¬‚ect that one origin of the development lies structural
mechanics. Shortly after this beginning, the ¬nite element method was ap-
plied to problems of heat conduction and ¬‚uid mechanics, which form the
application background of this book.
An intensive mathematical analysis and further development was started
in the later sixties. The basics of this mathematical description and analy-
sis are to be developed in this and the following chapter. The homogeneous
Dirichlet boundary value problem for the Poisson equation forms the
paradigm of this chapter, but more generally valid considerations will be
emphasized. In this way the abstract foundation for the treatment of more
general problems in Chapter 3 is provided. In spite of the importance of the
¬nite element method for structural mechanics, the treatment of the linear
elasticity equations will be omitted. But we note that only a small expense
is necessary for the application of the considerations to these equations.
We refer to [11], where this is realized with a very similar notation.



2.1 Variational Formulation for the Model Problem
We will develop a new solution concept for the boundary value problem
(1.1), (1.2) as a theoretical foundation for the ¬nite element method. For
2.1. Variational Formulation 47

such a solution, the validity of the di¬erential equation (1.1) is no longer re-
quired pointwise but in the sense of some integral average with “arbitrary”
weighting functions •. In the same way, the boundary condition (1.2) will
be weakened by the renunciation of its pointwise validity.
For the present, we want to con¬ne the considerations to the case of
homogeneous boundary conditions (i.e., g ≡ 0), and so we consider the
following homogeneous Dirichlet problem for the Poisson equation: Given
a function f : „¦ ’ R, ¬nd a function u : „¦ ’ R such that

’∆u = f in „¦ , (2.1)
u= 0 on ‚„¦ . (2.2)

In the following let „¦ be a domain such that the integral theorem of
Gauss is valid, i.e. for any vector ¬eld q : „¦ ’ Rd with components in
C(„¦) © C 1 („¦) it holds

∇ · q(x) dx = ν(x) · q(x) dσ . (2.3)
„¦ ‚„¦

Let the function u : „¦ ’ R be a classical solution of (2.1), (2.2) in the
sense of De¬nition 1.1, which additionally satis¬es u ∈ C 1 („¦) to facili-

tate the reasoning. Next we consider arbitrary v ∈ C0 („¦) as so-called test
functions. The smoothness of these functions allows all operations of di¬er-

entiation, and furthermore, all derivatives of a function v ∈ C0 („¦) vanish
on the boundary ‚„¦. We multiply equation (2.1) by v, integrate the result
over „¦, and obtain

f (x)v(x) dx = ’ ∇ · (∇u)(x) v(x) dx
f, v =
0
„¦ „¦

∇u(x) · ∇v(x) dx ’ ∇u(x) · ν(x) v(x) dσ
= (2.4)
„¦ ‚„¦

∇u(x) · ∇v(x) dx .
=
„¦

The equality sign at the beginning of the second line of (2.4) is obtained
by integration by parts using the integral theorem of Gauss with q = v∇u .
The boundary integral vanishes because v = 0 holds on ‚„¦.

If we de¬ne, for u ∈ C 1 („¦), v ∈ C0 („¦), a real-valued mapping a by

∇u(x) · ∇v(x) dx ,
a(u, v) :=
„¦

then the classical solution of the boundary value problem satis¬es the
identity

for all v ∈ C0 („¦) .
a(u, v) = f, v (2.5)
0
48 2. Finite Element Method for Poisson Equation


The mapping a de¬nes a scalar product on C0 („¦) that induces the norm
1/2
|∇u| dx
2
u := a(u, u) = (2.6)
a
„¦

(see Appendix A.4 for these notions). Most of the properties of a
scalar product are obvious. Only the de¬niteness (A4.7) requires further
considerations. Namely, we have to show that

(∇u · ∇u) (x) dx = 0 ⇐’ u ≡ 0 .
a(u, u) =
„¦

To prove this assertion, ¬rst we show that a(u, u) = 0 implies ∇u(x) = 0
for all x ∈ „¦. To do this, we suppose that there exists some point x ∈ „¦
¯
such that ∇u(¯) = 0. Then (∇u · ∇u) (¯) = |∇u| (¯) > 0. Because of
2
x x x
the continuity of ∇u, a small neighbourhood G of x exists with a positive
¯
measure |G| and |∇u|(x) ≥ ± > 0 for all x ∈ G. Since |∇u|2 (x) ≥ 0 for all
x ∈ „¦, it follows that
2
|∇u| (x) dx ≥ ±2 |G| > 0 ,
„¦

which is in contradiction to a(u, u) = 0. Consequently, ∇u(x) = 0 holds
for all x ∈ „¦; i.e., u is constant in „¦. Since u(x) = 0 for all x ∈ ‚„¦, the
assertion follows.

Unfortunately, the space C0 („¦) is too small to play the part of the basic

space because the solution u does not belong to C0 („¦) in general. The
identity (2.4) is to be satis¬ed for a larger class of functions, which include,
as an example for v, the solution u and the ¬nite element approximation
to u to be de¬ned later.

For the present we de¬ne as the basic space V ,
V := u : „¦ ’ R u ∈ C(„¦) , ‚i u exists and is piecewise
¯
(2.7)
continuous for all i = 1, . . . , d, u = 0 on ‚„¦ .
To say that ‚i u is piecewise continuous means that the domain „¦ can be
decomposed as follows:
¯ ¯
„¦= „¦j ,
j

with a ¬nite number of open sets „¦j , with „¦j © „¦k = … for j = k , and ‚i u
¯
is continuous on „¦j and it can continuously be extended on „¦j .
Then the following properties hold:
• a is a scalar product also on V ,

• C0 („¦) ‚ V ,

• C0 („¦) is dense in V with respect to · a ; i.e., for any u ∈ V (2.8)

a sequence (un )n∈N in C0 („¦) exists such that un ’u a ’ 0
for n ’ ∞,
2.1. Variational Formulation 49


• ·
C0 („¦) is dense in V with respect to 0. (2.9)
The ¬rst and second statements are obvious. The two others require a
certain technical e¬ort. A more general statement will be formulated in
Theorem 3.7.
With that, we obtain from (2.5) the following result:
Lemma 2.1 Let u be a classical solution of (2.1), (2.2) and let u ∈ C 1 („¦).
¯
Then
for all v ∈ V .
a(u, v) = f, v (2.10)
0

Equation (2.10) is also called a variational equation.


Proof: Let v ∈ V . Then vn ∈ C0 („¦) exist with vn ’ v with respect
to · 0 and also to · a . Therefore, it follows from the continuity of the
bilinear form with respect to · a (see (A4.22)) and the continuity of the
functional de¬ned by the right-hand side v ’ f, v 0 with respect to · 0
(because of the Cauchy“Schwarz inequality in L2 („¦)) that
’ f, v a(u, vn ) ’ a(u, v) for n ’ ∞ .
f, vn and
0 0

2
Since a(u, vn ) = f, vn 0 , we get a(u, v) = f, v 0 .
The space V in the identity (2.10) can be further enlarged as long as (2.8)
and (2.9) will remain valid. This fact will be used later to give a correct
de¬nition.
De¬nition 2.2 A function u ∈ V is called a weak (or variational) solution
of (2.1), (2.2) if the following variational equation holds:
for all v ∈ V .
a(u, v) = f, v 0

If u models e.g. the displacement of a membrane, this relation is called
the principle of virtual work.
Lemma 2.1 guarantees that a classical solution u is a weak solution.
The weak formulation has the following properties:
• It requires less smoothness: ‚i u has to be only piecewise continuous.
• The validity of the boundary condition is guaranteed by the de¬nition
of the function space V .
We now show that the variational equation (2.10) has exactly the same
solution(s) as a minimization problem:
Lemma 2.3 The variational equation (2.10) has the same solutions u ∈ V
as the minimization problem
F (v) ’ min for all v ∈ V , (2.11)
50 2. Finite Element Method for Poisson Equation

where
1 1
a(v, v) ’ f, v ’ f, v
2
F (v) := = v .
a
0 0
2 2


Proof: (2.10) ’ (2.11):
Let u be a solution of (2.10) and let v ∈ V be chosen arbitrarily. We de¬ne
w := v ’ u ∈ V (because V is a vector space), i.e., v = u + w. Then, using
the bilinearity and symmetry, we have
1
a(u + w, u + w) ’ f, u + w 0
F (v) =
2
1 1
a(u, u) + a(u, w) + a(w, w) ’ f, u ’ f, w
= (2.12)
0 0
2 2
1
F (u) + a(w, w) ≥ F (u) ,
=
2
where the last inequality follows from the positivity of a; i.e., (2.11) holds.
(2.10) ⇐ (2.11):
Let u be a solution of (2.11) and let v ∈ V , µ ∈ R be chosen arbitrarily. We
de¬ne g(µ) := F (u + µv) for µ ∈ R. Then
g(µ) = F (u + µv) ≥ F (u) = g(0) for all µ ∈ R ,
because u + µv ∈ V ; i.e., g has a global minimum at µ = 0.
It follows analogously to (2.12):
µ2
1
g(µ) = a(u, u) ’ f, u 0 + µ (a(u, v) ’ f, v 0 ) + a(v, v) .
2 2
Hence the function g is a quadratic polynomial in µ, and in particular,
g ∈ C 1 (R) is valid. Therefore we obtain the necessary condition
0 = g (µ) = a(u, v) ’ f, v 0

for the existence of a minimum at µ = 0. Thus u solves (2.10), because
v ∈ V has been chosen arbitrarily. 2
For applications e.g. in structural mechanics as above, the minimization
problem is called the principle of minimal potential energy.
Remark 2.4 Lemma 2.3 holds for general vector spaces V if a is a sym-
metric, positive bilinear form and the right-hand side f, v 0 is replaced by
b(v), where b : V ’ R is a linear mapping, a linear functional. Then the
variational equation reads as
¬nd u ∈ V for all v ∈ V ,
with a(u, v) = b(v) (2.13)
and the minimization problem as
¬nd u ∈ V F (u) = min F (v) v ∈ V
with , (2.14)
2.1. Variational Formulation 51

1
a(v, v) ’ b(v) .
where F (v) :=
2
Lemma 2.5 The weak solution according to (2.10) (or (2.11)) is unique.

Proof: Let u1 , u2 be two weak solutions, i.e.,
a(u1 , v) = f, v ,
0
for all v ∈ V .
a(u2 , v) = f, v ,
0

By subtraction, it follows that
a(u1 ’ u2 , v) = 0 for all v ∈ V.
Choosing v = u1 ’ u2 implies a(u1 ’ u2 , u1 ’ u2 ) = 0 and consequently
2
u1 = u2 , because a is de¬nite.
Remark 2.6 Lemma 2.5 is generally valid if a is a de¬nite bilinear form
and b is a linear form.
So far, we have de¬ned two di¬erent norms on V : · a and · 0 . The
di¬erence between these norms is essential because they are not equivalent
on the vector space V de¬ned by (2.7), and consequently, they generate
di¬erent convergence concepts, as will be shown by the following example:
Example 2.7 Let „¦ = (0, 1), i.e.
1
a(u, v) := u v dx ,
0

and let vn : „¦ ’ R for n ≥ 2 be de¬ned by (cf. Figure 2.1)
±
for 0 ¤ x ¤ n ,
1
 nx ,

for n ¤ x ¤ 1 ’ n ,
1 1
1,
vn (x) =


n ’ nx , for 1 ’ n ¤ x ¤ 1 .
1




vn
1



1 n-1
n1
n
Figure 2.1. The function vn .

Then
1/2
1
¤
vn 1 dx = 1,
0
0
52 2. Finite Element Method for Poisson Equation

1/2
1

1
n
2n ’ ∞ for n ’ ∞ .
n2 dx + n2 dx
vn = =
a
1
0 1’ n

¤C v
Therefore, there exists no constant C > 0 such that v for all
a 0
v ∈V.
However, as we will show in Theorem 2.18, there exists a constant C > 0
such that the estimate
¤C v for all v ∈ V
v 0 a

holds; i.e., · a is the stronger norm.
It is possible to enlarge the basic space V without violating the previous
statements. The enlargement is also necessary because, for instance, the
proof of the existence of a solution of the variational equation (2.13) or
the minimization problem (2.14) requires in general the completeness of V.
However, the actual de¬nition of V does not imply the completeness, as
the following example shows:
Example 2.8 Let „¦ = (0, 1) again and therefore
1
a(u, v) := u v dx .
0

For u(x) := x± (1’x)± with ± ∈ 1 , 1 we consider the sequence of functions
2
±
for x ∈ n , 1 ’ n ,
1 1
 u(x)

for x ∈ 0, n ,
1 1
un (x) := n u( n ) x


n u(1 ’ n ) (1 ’ x) for x ∈ 1 ’ n , 1 .
1 1


<<

. 2
( 16)



>>