<<

. 4
( 16)



>>

Denoting the matrix B by B = b(1) , b(2) , then it follows that
’1
b(1) · b(1) b(1) · b(2) b(2) · b(2) ’b(1) · b(2)
1
C= =
b(1) · b(2) b(2) · b(2) ’b(1) · b(2) b(1) · b(1)
det(B)2

because det(B T B) = det(B)2 . The previous considerations can be eas-
ily extended to the computation of the sti¬ness matrices of more general
di¬erential operators like

K(x)∇•j (x) · ∇•i (x) dx
„¦

(cf. Section 3.5). For the standard reference element, which we use from
now on, we have b(1) = a(2) ’ a(1) , b(2) = a(3) ’ a(1) . Here a(i) , i = 1, 2, 3,
are the locally numbered nodes of K interpreted as vectors of R2 .
From now on we make also use of the special form of the sti¬ness matrix
and obtain
˜(m)
Aij = γ1 ‚x1 Nj ‚x1 Ni dˆ
x
ˆ ˆ
ˆ
K

+ γ2 ‚x1 Nj ‚x2 Ni + ‚x2 Nj ‚x1 Ni dˆ
x (2.52)
ˆ ˆ ˆ ˆ
ˆ
K

+ γ3 ‚x2 Nj ‚x2 Ni dˆ
x
ˆ ˆ
ˆ
K

with
1
c11 | det(B)| = a(3) ’ a(1) · a(3) ’ a(1) ,
γ1 :=
| det(B)|
1
c12 | det(B)| = ’ a(2) ’ a(1) · a(3) ’ a(1) ,
γ2 :=
| det(B)|
1
c22 | det(B)| = a(2) ’ a(1) · a(2) ’ a(1) .
γ3 :=
| det(B)|
In the implementation it is advisable to compute the values γi just once
from the local geometrical information given in the form of the vertices
a(i) = ari and to store them permanently.
Thus we obtain for the local sti¬ness matrix
˜
A(m) = γ1 S1 + γ2 S2 + γ3 S3 (2.53)
with

S1 := ‚x1 Nj ‚x1 Ni dˆ
x ,
ˆ ˆ
ˆ
K ij
80 2. Finite Element Method for Poisson Equation


S2 := ‚x1 Nj ‚x2 Ni + ‚x2 Nj ‚x1 Ni dˆ
x ,
ˆ ˆ ˆ ˆ
ˆ
K ij

S3 := ‚x2 Nj ‚x2 Ni dˆ
x .
ˆ ˆ
ˆ
K ij

An explicit computation of the matrices Si is possible because the
integrands are constant, and also these matrices can be stored permanently:
«  «  « 
1 ’1 0 2 ’1 ’1 1 0 ’1
1 1 1
S1 = ’1 1 0  , S2 = ’1 0 1  , S3 =  0 0 0  .
2 2 2
’1 1 0 ’1 0 1
000

The right-hand side (q h )i = f (x)•i (x) dx can be treated in a similar
„¦
manner:
N
q (m)
(q h )i = i
m=1

with

(= 0 ’ ai ∈ Km ) .
q (m) = f (x)•i (x) dx
i
Km

(m)
Again, we transform the global numbering qri for the triangle
i=1,2,3
(m)
Km = conv {ar1 , ar2 , ar3 } into the local numbering qi˜ .
Anal- i=1,2,3
ogously to the determination of the entries of the sti¬ness matrix, we
have
(m)
f (F (ˆ)) •ri (F (ˆ)) | det(B)| dˆ
qi
˜ = x x x
ˆ
K
ˆx
f (ˆ) Ni (ˆ) | det(B)| dˆ ,
= x x
ˆ
K

ˆx
where f (ˆ) := f (F (ˆ)) for x ∈ K.
ˆˆ
x
In general, this integral cannot be evaluated exactly. Therefore, it has to
be approximated by a quadrature rule.
A quadrature rule for K g(ˆ) dˆ is of the type
xx
ˆ

R
ωk g ˆ(k)
b
k=1

with certain weights ωk and quadrature points ˆ(k) . As an example, we take
b
the trapezoidal rule (cf. (2.38)), where
ˆ(1) = a1 = (0, 0) , ˆ(2) = a2 = (1, 0) , ˆ(3) = a3 = (0, 1) ,
b ˆ b ˆ b ˆ
1
ωk = , k = 1, 2, 3 .
6
2.4. The Implementation of the Finite Element Method: Part 1 81

Thus for arbitrary but ¬xed quadrature rules, we have
R
(m)
ωk f ˆ(k) Ni ˆ(k) | det(B)| .
ˆb

qi
˜ b (2.54)
k=1

Of course, the application of di¬erent quadrature rules on di¬erent elements
is possible, too. The values Ni ˆ(k) , i = 1, 2, 3, k = 1, . . . , R, should be
b
evaluated just once and should be stored. The discussion on the use of
quadrature rules will be continued in Sections 3.5.2 and 3.6.
In summary, the following algorithm provides the assembling of the
sti¬ness matrix and the right-hand side:

Loop over all elements m = 1, . . . , N :
• Allocating a local numbering to the nodes based on the element-node
table: 1 ’ r1 , 2 ’ r2 , 3 ’ r3 .
• Assembling of the element sti¬ness matrix A(m) according to (2.51)
˜
or (2.53).
Assembling of the right-hand side according to (2.54).
• Loop over i, j = 1, 2, 3:
˜(m)
(Ah )ri rj := (Ah )ri rj + Aij ,
(m)
(q h )ri := (q h )ri + qi
˜ .
For the sake of e¬ciency of this algorithm, it is necessary to adjust the
memory structure to the particular situation; we will see how this can be
done in Section 2.5.


2.4.3 Realization of Dirichlet Boundary Conditions: Part 1
Nodes where a Dirichlet boundary condition is prescribed must be labeled
specially, here, for instance, by the convention M = M1 + M2 , where the
nodes numbered from M1 + 1 to M correspond to the Dirichlet boundary
nodes. In more general cases, other realizations are to be preferred.
In the ¬rst step of assembling of sti¬ness matrix and the load vector, the
Dirichlet nodes are treated like all the other ones. After this, the Dirichlet
nodes are considered separately. If such a node has the number j, the
boundary condition is included by the following procedure:
Replace the jth row and the jth column (for conservation of the sym-
metry) of Ah by the jth unit vector and (q h )j by g(aj ), if u(x) = g(x) is
prescribed for x ∈ ‚„¦. If the jth column is replaced by the unit vector, the
right-hand side (q h )i for i = j must be modi¬ed to (q h )i ’ (Ah )ij g(aj ). In
other words, the contributions caused by the Dirichlet boundary condition
are included into the right-hand side. This is exactly the elimination that
led to the form (1.10), (1.11) in Chapter 1.
82 2. Finite Element Method for Poisson Equation

2.5 Solving Sparse Systems of Linear Equations
by Direct Methods
Let A be an M —M matrix. Given a vector q ∈ RM , we consider the system
of linear equations
Aξ = q .
The matrices arising from the ¬nite element discretization are sparse; i.e.,
they have a bounded number of nonzero entries per row independent of
the dimension of the system of equations. For the simple example of Sec-
tion 2.2, this bound is determined by the number of neighbouring nodes
(see (2.37)). Methods for solving systems of equations should take advan-
tage of the sparse structure. For iterative methods, which will be examined
in Chapter 5, this is easier to reach than for direct methods. Therefore, the
importance of direct methods has decreased. Nevertheless, in adapted form
and for small or medium size problems, they are still the method of choice.

Elimination without Pivoting using Band Structure
In the general case, where the matrix A is assumed only to be nonsingular,
there exist M — M matrices P , L, U such that
P A = LU .
Here P is a permutation matrix, L is a scaled lower triangular matrix, and
U is an upper triangular matrix; i.e., they have the form
«  « 
1 0 u11 uij
¬ · ¬ ·
.. ..
L= , U = .
. .
lij 1 0 uMM
This decomposition corresponds to the Gaussian elimination method with
pivoting. The method is very easy and has favourable properties with re-
spect to the sparse structure, if pivoting is not necessary (i.e., P = I,
A = LU ). Then the matrix A is called LU factorizable.
Denote by Ak the leading principal submatrix of A of dimension k — k,
i.e.,
« 
a11 · · · a1k
¬. . ·,
Ak :=  . .. .
.
. .
· · · akk
ak1
and suppose that it already has been factorized as Ak = Lk Uk . This is
obviously possible for k = 1: A1 = (a11 ) = (1)(a11 ). The matrix Ak+1 can
be represented in the form of a block matrix
Ak b
Ak+1 =
cT d
2.5. Direct Methods for Sparse Linear Systems 83

with b, c ∈ Rk , d ∈ R.
Using the ansatz
Lk 0 Uk u
Lk+1 = , Uk+1 =
lT 1 0s
with unknown vectors u, l ∈ Rk and s ∈ R, it follows that
Ak+1 = Lk+1 Uk+1 ⇐’ Lk u = b , Uk l = c , lT u + s = d .
T
(2.55)
From this, we have the following result:
Let A be nonsingular. Then lower and upper triangular matrices
L, U exist with A = LU if and only if Ak is nonsingular for all (2.56)
1 ¤ k ¤ M . For this case, L and U are determined uniquely.
Furthermore, from (2.55) we have the following important consequences:
If the ¬rst l components of the vector b are equal to zero, then this is valid
for the vector u, too:
0 0
If b = , then u also has the structure u = .
β
Similarly,
0 0
c= implies the structure l = .
γ »
For example, if the matrix A has a structure as shown in Figure 2.16,
then the zeros outside of the surrounded entries are preserved after the
LU factorization. Before we introduce appropriate de¬nitions to generalize
these results, we want to consider the special case of symmetric matrices.

« 
|—|0|—|0 0
¬ 0 | — — | 0 | — |·
¬ ·
A=¬| — — — — — |·
¬ ·
 0 0 | — — 0 |
0|— — 0 —|

Figure 2.16. Pro¬le of a matrix.

If A is as before nonsingular and LU factorizable, then U = DLT with a
diagonal matrix D = diag (di ), and therefore
A = LDLT .
˜
This is true because A has the form A = LDU , where the upper triangular
˜
matrix U satis¬es the scaling condition uii = 1 for all i = 1, . . . , M . Such
˜
a factorization is unique, and thus
˜
A = AT implies LT = U , therefore A = LDLT .
84 2. Finite Element Method for Poisson Equation

If in particular A is symmetric and positive de¬nite, then also di > 0 is
˜
valid. Thus exactly one matrix L of the form
« 
l11 0
˜¬ ·
..
L=  with lii > 0 for all i
.
lij lMM
exists such that
˜˜
A = LLT , the so-called Cholesky decomposition.
We have
√ √
˜
LChol = LGauss D , where D := diag di .
This shows that the Cholesky method for the determination of the Cholesky
˜
factor L also preserves certain zeros of A in the same way as the Gaussian
elimination without pivoting.
In what follows, we want to specify the set of zeros that is preserved by
Gaussian elimination without pivoting. We will not consider a symmetric
matrix; but for the sake of simplicity we will consider a matrix with a
symmetric distribution of its entries.
De¬nition 2.19 Let A ∈ RM—M be a matrix such that aii = 0 for i =
1, . . . , M and
aij = 0 if and only if aji = 0 for all i, j = 1, . . . , M . (2.57)
We de¬ne, for i = 1, . . . , M,
fi (A) := min j aij = 0 , 1 ¤ j ¤ i .
Then
mi (A) := i ’ fi (A)
is called the ith (left-hand side) row bandwidth of A.
The bandwidth of a matrix A that satis¬es (2.57) is the number
m(A) := max mi (A) = max i ’ j aij = 0 , 1 ¤ j ¤ i ¤ M .
1¤i¤M

The band of the matrix A is
B(A) := (i, j), (j, i) i ’ m(A) ¤ j ¤ i , 1 ¤ i ¤ M .
The set
Env (A) := (i, j), (j, i) fi (A) ¤ j ¤ i , 1 ¤ i ¤ M
is called the hull or envelope of A. The number
M
p(A) := M + 2 mi (A)
i=1

is called the pro¬le of A.
2.5. Direct Methods for Sparse Linear Systems 85

The pro¬le is the number of elements of Env(A).
For the matrix A in Figure 2.16 we have (m1 (A), . . . , m5 (A)) =
(0, 0, 2, 1, 3), m(A) = 3, and p(A) = 17.
Summarizing the above considerations, we have proved the following
theorem:
Theorem 2.20 Let A be a matrix with the symmetric structure (2.57).
Then the Cholesky method or the Gaussian elimination without pivoting
preserves the hull and in particular the bandwidth.
The hull may contain zeros that will be replaced by (nonzero) entries during
the decomposition process. Therefore, in order to keep this ¬ll-in small, the
pro¬le should be as small as possible.
Furthermore, in order to exploit the matrix structure for an e¬cient
assembling and storage, this structure (or some estimate of it) should be
known in advance, before the computation of the matrix entries is started.
For example, if A is a sti¬ness matrix with the entries

∇•j · ∇•i dx ,
aij = a(•j , •i ) =
„¦

then the property

aij = 0 ai , aj are neighbouring nodes
can be used for the de¬nition of an (eventually too large) symmetric matrix
structure. This is also valid for the case of a nonsymmetric bilinear form
and thus a nonsymmetric sti¬ness matrix. Also in this case, the de¬nition
of fi (A) can be replaced by
fi (A) := min j 1 ¤ j ¤ i , j is a neighbouring node of i .
Since the characterization (2.56) of the possibility of the Gaussian elim-
ination without pivoting cannot be checked directly, we have to specify
su¬cient conditions. Examples for such conditions are the following (see
[34]):
• A is symmetric and positive de¬nite,
• A is an M-matrix.

Su¬cient conditions for this property were given in (1.32) and (1.32) .
In Section 3.9, geometrical conditions for the family of triangula-
tions (Th )h will be derived that guarantee that the ¬nite element
discretization considered here creates an M-matrix.

Data Structures
For sparse matrices, it is appropriate to store only the components within
the band or the hull. A symmetric matrix A ∈ RM—M with bandwidth
m can be stored in M (m + 1) memory positions. By means of the index
86 2. Finite Element Method for Poisson Equation

conversion aik ; bi,k’i+m+1 for k ¤ i, the matrix
« 
a12 · · · a1,m+1
a11
¬ ·
¬ ·
. ..
.
¬ a21 a22 · · · ·
.
. 0
¬ ·
¬ ·
. . .
.. .. ..
¬ ·
. . .
. . .
¬ ·
. . .
¬ ·
¬ ·
A = ¬ am+1,1 am+1,2 · · · am+1,m+1 . . . .. .. · ∈ RM—M
. .
¬ ·
¬ ·
¬ ·
.. .. .. .. .. ..
¬ ·
. . . . . .
¬ ·
¬ ·
¬ ·
.. .. .. .. ..
. . . . .
 
0
aM,M’m · · · aM,M’1 aM,M
is mapped to the matrix
« 
··· ···
0 0 a11
¬ ·
···
0 0 a21 a22
¬ ·
¬ ·
. . .
¬ ·
. . .
. . .
¬ ·
¬ ·
··· ···
0 am,1 am,m
¬ ·
· ∈ RM—(m+1) .
B = ¬ am+1,1 ··· · · · am+1,m am+1,m+1
¬ ·
¬ ·
. . . . .
¬ ·
. . . . .
. . . . .
¬ ·
¬ ·
. . . . .
 
. . . . .
. . . . .
aM,M’m · · · · · · aM,M’1 aM,M
The unused elements of B, i.e., (B)ij for i = 1, . . . , m, j = 1, . . . , m + 1 ’ i,
are here ¬lled with zeros.
For a general band matrix, the matrix B ∈ RM—(2m+1) obtained by the
above conversion has the following form:
« 
··· ···
0 0 a11 a12 a1,m+1
¬ ·
··· ··· ···
0 a21 a22 a2,m+2
¬ ·
¬ ·
. . . . . .
¬ ·
. . . . . .
¬ ·
. . . . . .
¬ ·
··· ··· ··· ···
¬ ·
0 am,1 am,2m
¬ ·
··· ··· ··· ··· ···
¬ ·
am+1,1 am+1,2m+1
B =¬ ·.
¬ ·
. . . . . . .
. . . . . . .
¬ ·
. . . . . . .
¬ ·
¬ aM ’m,M ’2m ·
··· ··· ··· ··· ··· aM ’m,M
¬ ·
¬ aM ’m+1,M ’2m+1 · · · ·
··· ··· · · · aM ’m+1,M 0
¬ ·
¬ ·
. . . . . .
 
. . . . . .
. . . . . .
··· ··· ···
aM,M ’m aM,M 0 0
Here, in the right lower part of the matrix, a further sector of unused
elements arose, which is also ¬lled with zeros.
If the storage is based on the hull, additionally a pointer ¬eld is needed,
which points to the diagonal elements, for example. If the matrix is sym-
2.5. Direct Methods for Sparse Linear Systems 87

metric, again the storage of the lower triangular matrix is su¬cient. For
the matrix A from Figure 2.16 under the assumption that A is symmetric,
the pointer ¬eld could act as shown in Figure 2.17.


a11 a22 a31 a32 a33 a43 a44 a52 a53 a54
a55
$
X
$$$
‰
rrr‰r T 
  $$
rr rr $$
 
$$
rr rr $$
 
1 2 5 7 11

Figure 2.17. Linear storage of the hull.


Coupled Assembling and Decomposition
A formerly popular method, the so-called frontal method, performs
simultaneously assembling and the Cholesky factorization.
We consider this method for the example of the sti¬ness matrix Ah =
(aij ) ∈ RM—M with bandwidth m (with the original numbering).
The method is based on the kth step of the Gaussian or Cholesky method
(cf. Figure 2.18).
k



0
k

Bk
0


Figure 2.18. kth step of the Cholesky method.

Only the entries of Bk are to be changed, i.e., only those elements aij
with k ¤ i, j ¤ k + m. The corresponding formula is
(k)
aik
(k+1) (k) (k)

aij = aij akj , i, j = k + 1, . . . , k + m . (2.58)
(k)
akk
Here, the upper indices indicate the steps of the elimination method, which
we store in aij . The entries aij are generated by summation of entries of
88 2. Finite Element Method for Poisson Equation

the element sti¬ness matrix of those elements K that contain nodes with
the indices i, j.
(k) (k)
Furthermore, to perform the elimination step (2.58), only aik , akj for
(k)
i, j = k, . . . , k + m must be completely assembled; aij , i, j = k + 1, . . . , k +
(k) (k+1) (k+1) (k+1)
m, can be replaced by aij if aij
˜ is later de¬ned by aij := aij
˜ +
(k) (k)
aij ’ aij . That is, for the present, aij needs to consist of only a few
˜
contributions of elements K with nodes i, j in K.
From these observations, the following algorithm is obtained. The kth
step for k = 1, . . . , M reads as follows:
• Assemble all of the missing contributions of elements K that contain
the node with index k.
• Compute A(k+1) by modi¬cation of the entries of Bk according to
(2.58).
• Store the kth row of A(k+1) , also out of the main memory.
• De¬ne Bk+1 (by a south-east shift).
Here the assembling is node-based and not element-based.
The advantage of this method is that Ah need not be completely assem-
bled and stored in the main memory, but only a matrix Bk ∈ R(m+1)—(m+1) .
Of course, if M is not too large, there may be no advantage.

Bandwidth Reduction
The complexity, i.e., the number of operations, is crucial for the application
of a particular method:
The Cholesky method, applied to a symmetric matrix A ∈ RM—M with
bandwidth m, requires O(m2 M ) operations in order to compute L.
However, the bandwidth m of the sti¬ness matrix depends on the num-
bering of the nodes. Therefore, a numbering is to be found where the
number m is as small as possible.
We want to consider this again for the example of the Poisson equation on
the rectangle with the discretization according to Figure 2.9. Let the inte-
rior nodes have the coordinates (ih, jh) with i = 1, . . . k ’1, j = 1, . . . , l’1.
The discretization corresponds to the ¬nite di¬erence method introduced
beginning with (1.10); i.e., the bandwidth is equal to k ’ 1 for a rowwise
numbering or l ’ 1 for a columnwise numbering.
For k l or k l, this fact results in a large di¬erence of the bandwidth
m or of the pro¬le (of the left triangle), which is of size (k ’ 1)(l ’ 1)(m + 1)
except for a term of m2 . Therefore, the columnwise numbering is preferred
for k l; the rowwise numbering is preferred for k l.
For a general domain „¦, a numbering algorithm based on a given tri-
angulation Th and on a basis {•i } of Vh is necessary with the following
properties:
2.5. Direct Methods for Sparse Linear Systems 89

The structure of A resulting from the numbering must be such that the
band or the pro¬le of A is as small as possible. Furthermore, the numbering
algorithm should yield the numbers m(A) or fi (A), mi (A) such that the
matrix A can also be assembled using the element matrices A(k) .
Given a triangulation Th and a corresponding basis •i 1 ¤ i ¤ M of
Vh , we start with the assignment of some graph G to this triangulation as
follows:
The nodes of G coincide with the nodes {a1 , . . . , aM } of the triangulation.
The de¬nition of its edges is:

⇐’
(ai , aj ) is an edge of G
there exists a K ∈ Th such that •i |K ≡ 0 , •j |K ≡ 0 .

In Figure 2.19 some examples are given, where the example (2) will be
introduced in Section 3.3.

triangulation:
. . . . . .
(1) (2)


. . . . .
linear ansatz on triangle (bi)linear ansatz on quadrilateral


Graph:
. . . . . .
. . . . .
Figure 2.19. Triangulation and assigned graph.


If several degrees of freedom are assigned to some node of the triangu-
lation Th , then also in G several nodes are assigned to it. This is the case,
for example, if so-called Hermite elements are considered, which will be
introduced in Section 3.3. The costs of administration are small if the same
number of degrees of freedom is assigned to all nodes of the triangulation.
An often-used numbering algorithm is the Cuthill“McKee method. This
algorithm operates on the graph G just de¬ned. Two nodes ai , aj of G are
called neighboured if (ai , aj ) is an edge of G. The degree of a node ai of G
is de¬ned as the number of neighbours of ai .
The kth step of the algorithm for k = 1, . . . , M has the following form:
k = 1: Choose a starting node, which gets the number 1. This starting
node forms the level 1.
k > 1: If all nodes are already numbered, the algorithm is terminated.
Otherwise, the level k is formed by taking all the nodes that are not num-
90 2. Finite Element Method for Poisson Equation

bered yet and that are neighbours of a node of level k ’ 1. The nodes of
level k will be consecutively numbered.
Within a level, we can sort, for example, by the degree, where the node
with the smallest degree is numbered ¬rst.
The reverse Cuthill“McKee method consists of the above method and the
inversion of the numbering at the end; i.e.,
new node number = M + 1 ’ old node number .
This corresponds to a re¬‚ection of the matrix at the counterdiagonal. The
bandwidth does not change by the inversion, but the pro¬le may diminish
drastically for many examples (cf. Figure 2.20).

* *** *** *
* *** *
* *** *
* *** *
* *** *
* *** *** *
Figure 2.20. Change of the hull by re¬‚ection at the counterdiagonal.

The following estimate holds for the bandwidth m of the numbering
created by the Cuthill“McKee algorithm:
D+i
¤ m ¤ max (Nk’1 + Nk ’ 1) .
2 2¤k¤ν

Here D is the maximum degree of a node of G, ν is the number of levels, and
Nk is the number of nodes of level k. The number i is equal to 0 if D is even,
and i is equal to 1 if D is odd. The left-hand side of the above inequality is
easy to understand by means of the following argument: To reach a minimal
bandwidth, all nodes that are neighbours of ai in the graph G should also
be neighbours of ai in the numbering. Then the best situation is given if the
neighboured nodes would appear uniformly immediately before and after
ai . If D is odd, then one side has one node more than the other.
To verify the right-hand side, consider a node ai that belongs to level
k ’ 1 as well as a node aj that is a neighbour of ai in the graph G and that
is not yet numbered in level k ’ 1. Therefore, aj will get a number in the
kth step. The largest bandwidth is obtained if ai is the ¬rst node of the
numbering of level k ’ 1 and if aj is the last node of level k. Hence exactly
(Nk’1 ’ 1) + (Nk ’ 1) nodes lie between both of these; i.e., their distance
in the numbering is Nk’1 + Nk ’ 1.
It is favourable if the number ν of levels is as large as possible and if all
the numbers Nk are of the same size, if possible. Therefore, the starting
node should be chosen “at one end” of the graph G if possible; if all the
2.5. Direct Methods for Sparse Linear Systems 91

˜ ˜
starting nodes are to be checked, the expense will be O(M M ), where M
is the number of edges of G. One possibility consists in choosing a node
with minimum degree for the starting node. Another possibility is to let
the algorithm run once and then to choose the last-numbered node as the
starting node.
If a numbering is created by the (reverse) Cuthill“McKee algorithm, we
can try to improve it “locally”, i.e., by exchanging particular nodes.


Exercise
2.8 Show that the number of arithmetic operations for the Cholesky
method for an M — M matrix with bandwidth m has order M m2 /2;
additionally, M square roots have to be calculated.
3
The Finite Element Method for Linear
Elliptic Boundary Value Problems of
Second Order




3.1 Variational Equations and Sobolev Spaces
We now continue the de¬nition and analysis of the “correct” function spaces
that we began in (2.17)“(2.20). An essential assumption ensuring the exis-
tence of a solution of the variational equation (2.13) is the completeness of
the basic space (V, · ). In the concrete case of the Poisson equation the
“preliminary” function space V according to (2.7) can be equipped with
the norm · 1 , de¬ned in (2.19), which has been shown to be equivalent to
the norm · a , given in (2.6) (see (2.46)). If we consider the minimization
problem (2.14), which is equivalent to the variational equation, the func-
tional F is bounded from below such that the in¬mum assumes a ¬nite
value and there exists a minimal sequence (vn )n in V , that is, a sequence
with the property

lim F (vn ) = inf F (v) v ∈ V .
n’∞


The form of F also implies that (vn )n is a Cauchy sequence. If this sequence
converges to an element v ∈ V , then, due to the continuity of F with respect
to · , it follows that v is a solution of the minimization problem. This
completeness of V with respect to · a , and hence with respect to · 1 , is
not satis¬ed in the de¬nition (2.7), as Example 2.8 has shown. Therefore,
an extension of the basic space V , as formulated in (2.20), is necessary.
This space will turn out to be “correct,” since it is complete with respect
to · 1 .
3.1. Variational Equations and Sobolev Spaces 93

In what follows we use the following general assumption:
V is a vector space with scalar product ·, · and the norm ·
1/2
induced by ·, · (for this, v := v, v for v ∈ V issatisf ied) ;
V is complete with respect to · , i.e. a Hilbert space ; (3.1)
a : V — V ’ R is a (not necessarily symmetric) bilinear form ;
b : V ’ R is a linear form .
The following theorem generalizes the above consideration to nonsym-
metric bilinear forms:
Theorem 3.1 (Lax“Milgram) Suppose the following conditions are sat-
is¬ed:
• a is continuous (cf. (2.42)); that is, there exists some constant
M > 0 such that
|a(u, v)| ¤ Mu for all u, v ∈ V ;
v (3.2)
• a is V -elliptic (cf. (2.43)); that is, there exists some constant ± > 0
such that

a(u, u) ≥ ± u for all u ∈ V ;
2
(3.3)
• b is continuous; that is, there exists some constant C > 0 such that
|b(u)| ¤ C u for all u ∈ V . (3.4)
Then the variational equation (2.21), namely,
¬nd u ∈ V such that for all v ∈ V,
¯ a(¯, v) = b(v)
u (3.5)
has one and only one solution.
Here, one cannot avoid the assumptions (3.1) and (3.2)“(3.4) in general.

Proof: See, for example, [26]; for an alternative proof see Exercise 3.1. 2

Now returning to the example above, the assumptions (3.2) and (3.3) are
obviously satis¬ed for · = · a . However, the “preliminary” de¬nition of
the function space V of (2.7) with norm · a de¬ned in (2.19) is insu¬cient,
since (V, · a ) is not complete. Therefore, the space V must be extended.
Indeed, it is not the norm on V that has been chosen incorrectly, since V
is also not complete with respect to another norm · that satis¬es (3.2)
and (3.3). In this case the norms · and · a would be equivalent (cf.
(2.46)), and consequently,
(V, · ⇐’ (V, · ) complete .
a) complete
Now we extend the space V and thereby generalize de¬nition (2.17).
94 3. Finite Element Methods for Linear Elliptic Problems

De¬nition 3.2 Suppose „¦ ‚ Rd is a (bounded) domain.
The Sobolev space H k („¦) is de¬ned by
H k („¦) := v : „¦ ’ R v ∈ L2 („¦) , the weak derivatives ‚ ± v exist
in L2 („¦) and for all multi-indices ± with |±| ¤ k .
A scalar product ·, · and the resulting norm · in H k („¦) are de¬ned
k
k
as follows:

‚ ± v ‚ ± w dx ,
v, w := (3.6)
k
„¦ ± multi’index
|±|¤k
1/2
2
1/2 ±
v := v, v = ‚ v dx (3.7)
k k
„¦ ± multi’index
|±|¤k
1/2 1/2
2 2
± ±
= ‚ v dx = ‚v .
0
„¦
± multi’index ± multi’index
|±|¤k |±|¤k

Greater ¬‚exibility with respect to the smoothness properties of the func-
tions that are contained in the de¬nition is obtained by requiring that v
and its weak derivatives should belong not to L2 („¦) but to Lp („¦). In the
norm denoted by · k,p the L2 („¦) and 2 norms (for the vector of the
derivative norms) have to be replaced by the Lp („¦) and p norms, respec-
tively (see Appendices A.3 and A.5). However, the resulting space, denoted
k
by Wp („¦), can no longer be equipped with a scalar product for p = 2. Al-
though these spaces o¬er greater ¬‚exibility, we will not use them except in
Sections 3.6, 6.2, and 9.3.
Besides the norms · k , there are seminorms | · |l for 0 ¤ l ¤ k in H k („¦),
de¬ned by
1/2
2
|v|l = ‚±v ,
0
± multi’index
|±|=l

such that
1/2
k
|v|2
v = ,
k l
l=0

In particular, these de¬nitions are compatible with those in (2.18),

vw + ∇v · ∇w dx ,
v, w :=
1
„¦

· for the L2 („¦) norm, giving a meaning to this
and with the notation 0
one.
3.1. Variational Equations and Sobolev Spaces 95

The above de¬nition contains some assertions that are formulated in the
following theorem:
Theorem 3.3 The bilinear form ·, · k is a scalar product on H k („¦); that
is, · k is a norm on H k („¦).
H k („¦) is complete with respect to · k , and is thus a Hilbert space.

2
Proof: See, for example, [37].

Obviously,
H k („¦) ‚ H l („¦) for k ≥ l ,
and the embedding is continuous, since
¤v for all v ∈ H k („¦) .
v (3.8)
l k

In the one-dimensional case (d = 1) v ∈ H 1 („¦) is necessarily continuous:
Lemma 3.4
H 1 (a, b) ‚ C[a, b] ,
and the embedding is continuous, where C[a, b] is equipped with the norm
· ∞ ; that is, there exists some constant C > 0 such that
¤C v for all v ∈ H 1 (a, b) .
v (3.9)
∞ 1


2
Proof: See Exercise 3.2.

Since the elements of H k („¦) are ¬rst of all only square integrable func-
tions, they are determined only up to points of a set of (d-dimensional)
measure zero. Therefore, a result as in Lemma 3.4 means that the func-
tion is allowed to have removable discontinuities at points of such a set of
measure zero that vanish by modifying the function values.
However, in general, H 1 („¦) ‚ C(„¦).
¯
As an example for this, we consider a circular domain in dimension d = 2:
„¦ = BR (0) = x ∈ R2 |x| < R , R < 1.
Then the function
1
v(x) := |log |x| |γ for some γ <
2
¯
is in H 1 („¦), but not in C(„¦) (see Exercise 3.3).
The following problem now arises: In general, one cannot speak of a
value v(x) for some x ∈ „¦ because a set of one point {x} has (Lebesgue)
measure zero. How do we then have to interpret the Dirichlet boundary
conditions? A way out is to consider the boundary (pieces of the boundary,
respectively) not as arbitrary points but as (d ’ 1)-dimensional “spaces”
(manifolds).
96 3. Finite Element Methods for Linear Elliptic Problems

The above question can therefore be reformulated as follows: Is it possible
to interpret v on ‚„¦ as a function of L2 (‚„¦) (‚„¦ “‚” Rd’1 ) ?
It is indeed possible if we have some minimal regularity of ‚„¦ in the
following sense: It has to be possible to choose locally, for some boundary
point x ∈ ‚„¦, a coordinate system in such a way that the boundary is
locally a hyperplane in this coordinate system and the domain lies on one
side. Depending on the smoothness of the parametrisation of the hyperplane
we then speak of Lipschitz, C k - (for k ∈ N), and C ∞ - domains (for an exact
de¬nition see Appendix A.5).
Examples:
(1) A circle „¦ = x ∈ Rd |x ’ x0 | < R is a C k -domain for all k ∈ N,
and hence a C ∞ -domain.
(2) A rectangle „¦ = x ∈ Rd 0 < xi < ai , i = 1, . . . , d is a Lipschitz
domain, but not a C 1 -domain.
(3) A circle with a cut „¦ = x ∈ Rd |x ’ x0 | < R, x = x0 + »e1 for 0 ¤
» < R is not a Lipschitz domain, since „¦ does not lie on one side of
‚„¦ (see Figure 3.1).



„¦ „¦ „¦


Circle Rectangle Circle with cut

Figure 3.1. Domains of di¬erent smoothness.

Hence, suppose „¦ is a Lipschitz domain. Since only a ¬nite number of
overlapping coordinate systems are su¬cient for the description of ‚„¦,
using these, it is possible to introduce a (d ’ 1)-dimensional measure on
‚„¦ and de¬ne the space L2 (‚„¦) of square integrable functions with respect
to this measure (see Appendix A.5 or [37] for an extensive description). In
the following, let ‚„¦ be equipped with this (d ’ 1)-dimensional measure
dσ, and integrals over the boundary are to be interpreted accordingly. This
also holds for Lipschitz subdomains of „¦, since they are given by the ¬nite
elements.
Theorem 3.5 (Trace Theorem) Suppose „¦ is a bounded Lipschitz do-
main. We de¬ne
C ∞ (Rd )|„¦ v : „¦ ’ R v can be extended to v : Rd ’ R and
:= ˜

v ∈ C (R ) .
d
˜
Then, C ∞ (Rd )|„¦ is dense in H 1 („¦); that is, with respect to · 1 an arbitrary
w ∈ H 1 („¦) can be approximated arbitrarily well by some v ∈ C ∞ (Rd )|„¦ .
3.1. Variational Equations and Sobolev Spaces 97

The mapping that restricts v to ‚„¦,
γ0 : C ∞ (Rd )|„¦ , · ’ L2 (‚„¦), · ,
1 0

’ v|‚„¦ ,
v
is continuous.
Thus there exists a unique, linear, and continuous extension
γ0 : H 1 („¦), · ’ L2 (‚„¦), · .
1 0


2
Proof: See, for example, [37].

Therefore, in short form, γ0 (v) ∈ L2 (‚„¦), and there exists some constant
C > 0 such that
¤C v for all v ∈ H 1 („¦) .
γ0 (v) 0 1

Here γ0 (v) ∈ L2 (‚„¦) is called the trace of v ∈ H 1 („¦).
The mapping γ0 is not surjective; that is, γ0 (v) v ∈ H 1 („¦) is a real
subset of L2 (‚„¦). For all v ∈ C ∞ (Rd )|„¦ we have
γ0 (v) = v|‚„¦ .
In the following we will use again v|‚„¦ or “v on ‚„¦” for γ0 (v), but in
the sense of Theorem 3.5. According to this theorem, de¬nition (2.20) is
well-de¬ned with the interpretation of u on ‚„¦ as the trace:
De¬nition 3.6 H0 („¦) := v ∈ H 1 („¦) γ0 (v) = 0 (as a function on ‚„¦) .
1



Theorem 3.7 Suppose „¦ ‚ Rd is a bounded Lipschitz domain. Then
∞ 1
C0 („¦) is dense in H0 („¦).

2
Proof: See [37].

The assertion of Theorem 3.5, that C ∞ (Rd )|„¦ is dense in H 1 („¦), has
severe consequences for the treatment of functions in H 1 („¦) which are in
general not very smooth. It is possible to consider them as smooth functions
if at the end only relations involving continuous expressions in · 1 (and not
requiring something like ‚i v ∞ ) arise. Then, by some “density argument”
the result can be transferred to H 1 („¦) or, as for the trace term, new terms
can be de¬ned for functions in H 1 („¦). Thus, for the proof of Lemma 3.4
it is necessary simply to verify estimate (3.9), for example for v ∈ C 1 [a, b].
1
By virtue of Theorem 3.7, analogous results hold for H0 („¦).
Hence, for v ∈ H 1 („¦) integration by parts is possible:
Theorem 3.8 Suppose „¦ ‚ Rd is a bounded Lipschitz domain. The outer
unit normal vector ν = (νi )i=1,...,d : ‚„¦ ’ Rd is de¬ned almost everywhere
and νi ∈ L∞ (‚„¦).
98 3. Finite Element Methods for Linear Elliptic Problems

For v, w ∈ H 1 („¦) and i = 1, . . . , d,

‚i v w dx = ’ v ‚i w dx + v w νi dσ .
„¦ „¦ ‚„¦


2
Proof: See, for example, [14] or [37].

If v ∈ H 2 („¦), then due to the above theorem, v|‚„¦ := γ0 (v) ∈ L2 (‚„¦)
and ‚i v|‚„¦ := γ0 (‚i v) ∈ L2 (‚„¦), since also ‚i v ∈ H 1 („¦). Hence, the normal
derivative
d
‚ν v|‚„¦ := ‚i v|‚„¦ νi
i=1

is well-de¬ned and belongs to L2 (‚„¦).
Thus, the trace mapping
γ : H 2 („¦) ’ L2 (‚„¦) — L2 (‚„¦) ,
’ (v|‚„¦ , ‚ν v|‚„¦ ) ,
v
is well-de¬ned and continuous. The continuity of this mapping follows from
the fact that it is a composition of continuous mappings:

continuous continuous
v ∈ H 2 („¦) ’ ‚i v ∈ H 1 („¦) ’ ‚i v|‚„¦ ∈ L2 (‚„¦)
continuous
’ ‚i v|‚„¦ νi ∈ L2 (‚„¦) .
Corollary 3.9 Suppose „¦ ‚ Rd is a bounded Lipschitz domain.
(1) Let w ∈ H 1 („¦), qi ∈ H 1 („¦), i = 1, . . . , d. Then

q · ∇w dx = ’ ∇ · q w dx + q · ν w dσ . (3.10)
„¦ „¦ ‚„¦

(2) Let v ∈ H 2 („¦), w ∈ H 1 („¦). Then

∇v · ∇w dx = ’ ∆v w dx + ‚ν v w dσ .
„¦ „¦ ‚„¦
The integration by parts formulas also hold more generally if only it is
ensured that the function whose trace has to be formed belongs to H 1 („¦).
For example, if K = (kij )ij , where kij ∈ W∞ („¦) and v ∈ H 2 („¦), w ∈
1

H 1 („¦), it follows that

K∇v · ∇w dx = ’ ∇ · (K∇v) w dx + K∇v · ν w dσ (3.11)
„¦ „¦ ‚„¦

with conormal derivative (see (0.41))
d
‚νK v := K∇v · ν = ∇v · K ν = T
kij ‚j v νi .
i,j=1
Exercises 99

Here it is important that the components of K∇v belong to H 1 („¦), using
the fact that for v ∈ L2 („¦), k ∈ L∞ („¦),
kv ∈ L2 („¦) and ¤k
kv v .

0 0

Theorem 3.10 Suppose „¦ ‚ Rd is a bounded Lipschitz domain.
If k > d/2, then
H k („¦) ‚ C(„¦) ,
¯

and the embedding is continuous.

2
Proof: See, for example, [37].

For dimension d = 2 this requires k > 1, and for dimension d = 3 we
need k > 3 . Therefore, in both cases k = 2 satis¬es the assumption of the
2
above theorem.



Exercises
3.1 Prove the Lax“Milgram Theorem in the following way:
(a) Show, by using the Riesz representation theorem, the equivalence of
(3.5) with the operator equation

A¯ = f
u
for A ∈ L[V, V ] and f ∈ V .

(b) Show, for Tµ ∈ L[V, V ], Tµ v := v ’ µ(Av ’ f ) and µ > 0, that for
some µ > 0, the operator Tµ is a contraction on V . Then conclude
the assertion by Banach™s ¬xed-point theorem (in the Banach space
setting, cf. Remark 8.5).

3.2 Prove estimate (3.9) by showing that even for v ∈ H 1 (a, b),

|v(x) ’ v(y)| ¤ |v|1 |x ’ y|1/2 for x, y ∈ (a, b) .

3.3 Suppose „¦ ‚ R2 is the open disk with radius 1 and centre 0. Prove
2
±
that for the function u(x) := ln |x| , x ∈ „¦ \ {0}, ± ∈ (0, 1 ) we have
2
u ∈ H 1 („¦), but u cannot be extended continuously to x = 0.

3.4 Suppose „¦ ‚ R2 is the open unit disk. Prove that each u ∈ H 1 („¦)

has a trace u|‚„¦ ∈ L2 (‚„¦) satisfying u 0,‚„¦ ¤ 4 8 u 1,„¦ .
100 3. Finite Element Methods for Linear Elliptic Problems

3.2 Elliptic Boundary Value Problems of Second
Order
In this section we integrate boundary value problems for the linear, sta-
tionary case of the di¬erential equation (0.33) into the general theory of
Section 3.1.
Concerning the domain we will assume that „¦ is a bounded Lipschitz
domain.
We consider the equation
(Lu)(x) := ’∇ · (K(x)∇u(x)) + c(x) · ∇u(x) + r(x)u(x) = f (x) for x ∈ „¦
(3.12)
with the data
K : „¦ ’ Rd,d , c : „¦ ’ Rd , r, f : „¦ ’ R.

Assumptions about the Coe¬cients and the Right-Hand Side
For an interpretation of (3.12) in the classical sense, we need
‚i kij , ci , r , f ∈ C(„¦) , i, j ∈ {1, . . . , d} ,
¯ (3.13)
and for an interpretation in the sense of L2 („¦) with weak derivatives, and
hence for a solution in H 2 („¦),
‚i kij , ci , r ∈ L∞ („¦) , f ∈ L2 („¦) , i, j ∈ {1, . . . , d} . (3.14)
Once we have obtained the variational formulation, weaker assumptions
about the smoothness of the coe¬cients will be su¬cient for the veri¬ca-
tion of the properties (3.2)“(3.4), which are required by the Lax“Milgram,
namely,
kij , ci , ∇ · c , r ∈ L∞ („¦) , f ∈ L2 („¦) , i, j ∈ {1, . . . , d} ,
(3.15)
ν · c ∈ L∞ (“1 ∪ “2 ) .
and if |“1 ∪ “2 |d’1 > 0 ,
Here we refer to a de¬nition of the boundary conditions as in (0.36)“(0.39)
(see also below). Furthermore, the uniform ellipticity of L is assumed: There
exists some constant k0 > 0 such that for (almost) every x ∈ „¦,
d
kij (x)ξi ξj ≥ k0 |ξ|2 for all ξ ∈ Rd (3.16)
i,j=1

(that is, the coe¬cient matrix K is positive de¬nite uniformly in x).
Moreover, K should be symmetric.
If K is a diagonal matrix, that is, kij (x) = ki (x)δij (this is in particular
the case if ki (x) = k(x) with k : „¦ ’ R, i ∈ {1, . . . , d}, where K∇u
becomes k∇u), this means that
” ki (x) ≥ k0 for (almost) every x ∈ „¦ , i ∈ {1, . . . , d} .
(3.16)
3.2. Elliptic Boundary Value Problems 101

Finally, there exists a constant r0 ≥ 0 such that
1
r(x) ’ ∇ · c(x) ≥ r0 for (almost) every x ∈ „¦ . (3.17)
2

Boundary Conditions
As in Section 0.5, suppose “1 , “2 , “3 is a disjoint decomposition of the
boundary ‚„¦ (cf. (0.39)):
‚„¦ = “1 ∪ “2 ∪ “3 ,
where “3 is a closed subset of the boundary. For given functions gj : “j ’
R , j = 1, 2, 3, and ± : “2 ’ R we assume on ‚„¦
• Neumann boundary condition (cf. (0.41) or (0.36))
K∇u · ν = ‚νK u = g1 on “1 , (3.18)
• mixed boundary condition (cf. (0.37))
K∇u · ν + ±u = ‚νK u + ±u = g2 on “2 , (3.19)
• Dirichlet boundary condition (cf. (0.38))
u = g3 on “3 . (3.20)
Concerning the boundary data the following is assumed: For the classical
approach we need
gj ∈ C(“j ) , ± ∈ C(“2 ) ,
j = 1, 2, 3 , (3.21)
whereas for the variational interpretation,
± ∈ L∞ (“2 )
gj ∈ L2 (“j ) , j = 1, 2, 3 , (3.22)
is su¬cient.


3.2.1 Variational Formulation of Special Cases
The basic strategy for the derivation of the variational formulation of
boundary value problems (3.12) has already been demonstrated in Sec-
tion 2.1. Assuming the existence of a classical solution of (3.12) the
following steps are performed in general:
Step 1: Multiplication of the di¬erential equation by test functions that
are chosen compatible with the type of boundary condition and
subsequent integration over the domain „¦.
Step 2: Integration by parts under incorporation of the boundary condi-
tions in order to derive a suitable bilinear form.
Step 3: Veri¬cation of the required properties like ellipticity and continuity.
102 3. Finite Element Methods for Linear Elliptic Problems

In the following the above steps will be described for some important special
cases.

(I) Homogeneous Dirichlet Boundary Condition
‚„¦ = “3 , g3 ≡ 0 , V := H0 („¦)
1

Suppose u is a solution of (3.12), (3.20); that is, in the sense of classical
solutions let u ∈ C 2 („¦) © C(„¦) and the di¬erential equation (3.12) be
¯
satis¬ed pointwise in „¦ under the assumptions (3.13) as well as u = 0
pointwise on ‚„¦. However, the weaker case in which u ∈ H 2 („¦) © V and
the di¬erential equation is satis¬ed in the sense of L2 („¦), now under the
assumptions (3.14), can also be considered.

Multiplying (3.12) by v ∈ C0 („¦) (in the classical case) or by v ∈ V ,
respectively, then integrating by parts according to (3.11) and taking into
∞ 1
account that v = 0 on ‚„¦ by virtue of the de¬nition of C0 („¦) and H0 („¦),
respectively, we obtain

{K∇u · ∇v + c · ∇u v + r uv} dx
a(u, v) := (3.23)
„¦

for all v ∈ C0 („¦) or v ∈ V .
= b(v) := f v dx
„¦

The bilinear form a is symmetric if c vanishes (almost everywhere).
For f ∈ L2 („¦),
b is continuous on (V, · 1) . (3.24)
This follows directly from the Cauchy“Schwarz inequality, since

|b(v)| ¤ |f | |v| dx ¤ f ¤f for v ∈ V .
v v
0 0 0 1
„¦

Further, by (3.15),
a is continuous (V, · 1) . (3.25)

Proof: First, we obtain

|a(u, v)| ¤ {|K∇u| |∇v| + |c| |∇u||v| + |r| |u| |v|} dx .
„¦

Here | · | denotes the absolute value of a real number or the Euclidean
norm of a vector. Using also · 2 for the (associated) spectral norm, and
· ∞ for the L∞ („¦) norm of a function, we further introduce the following
notations:
< ∞, C2 := |c| < ∞.
C1 := max K 2 ∞, r ∞ ∞

By virtue of
|K(x)∇u(x)| ¤ K(x) |∇u(x)|,
2
3.2. Elliptic Boundary Value Problems 103

we continue to estimate as follows:

|a(u, v)| ¤ C1 {|∇u| |∇v| + |u| |v|} dx + C2 |∇u| |v| dx .
„¦ „¦
=:A1 =:A2

The integrand of the ¬rst addend is estimated by the Cauchy“Schwarz
inequality for R2 , and then the Cauchy“Schwarz inequality for L2 („¦) is
applied:

1/2 1/2
A1 ¤ C1 |∇u|2 + |u|2 |∇v|2 + |v|2 dx
„¦
1/2 1/2
¤ C1 |u| + |∇u| dx |v| + |∇v| dx
2 2 2 2
= C1 u v 1.
1
„¦ „¦

Dealing with A2 , we can employ the Cauchy“Schwarz inequality for L2 („¦)
directly:
1/2 1/2
¤ |∇u| dx |v| dx
2 2
A2 C2
„¦ „¦
¤ ¤ C2 u for all u, v ∈ V .
C2 u v v
1 0 1 1

2
Thus, the assertion follows.

Remark 3.11 In the proof of the propositions (3.24) and (3.25) it has not
been used that the functions u, v satisfy homogeneous Dirichlet boundary
conditions. Therefore, under the assumptions (3.15) these properties hold
for every subspace V ‚ H 1 („¦).


Conditions for the V -Ellipticity of a
(A) a is symmetric; that is c = 0 (a.e.): Condition (3.17) then has the
simple form r(x) ≥ r0 for almost all x ∈ „¦.
(A1) c = 0, r0 > 0:
Because of (3.16) we directly get

a(u, u) ≥ {k0 |∇u|2 + r0 |u|2 } dx ≥ C3 u for all u ∈ V ,
2
1
„¦

where C3 := min{k0 , r0 }. This also holds for every subspace V ‚ H 1 („¦).
(A2) c = 0, r0 ≥ 0:
According to the Poincar´ inequality (Theorem 2.18), there exists some
e
constant CP > 0, independent of u, such that for u ∈ H0 („¦)
1


1/2
¤ CP |∇u| dx 2
u .
0
„¦
104 3. Finite Element Methods for Linear Elliptic Problems

Taking into account (3.16) and using the simple decomposition k0 =
2
k0 CP
2 + 1 + C 2 k0 we can further conclude that
1 + CP P


a(u, u) ≥ k0 |∇u|2 dx (3.26)
„¦
2
k0 CP 1
≥ |∇u| dx + |u|2 dx = C4 u
2 2
k0 2 1,
2 2
1 + CP 1 + CP CP
„¦ „¦

k0
where C4 := 2 > 0.
1 + CP
For this estimate it is essential that u satis¬es the homogeneous Dirichlet
boundary condition.
(B) |c| ∞ > 0 :

First of all, we consider a smooth function u ∈ C0 („¦). From u∇u = 1 ∇u2
2
we get by integrating by parts
1 1
c · ∇u u dx = c · ∇u2 dx = ’ ∇ · c u2 dx .
2 2
„¦ „¦ „¦

Since according to Theorem 3.7 the space C0 („¦) is dense in V , the above
relation also holds for u ∈ V . Consequently, by virtue of (3.16) and (3.17)
we obtain
1
K∇u · ∇u + r ’ ∇ · c u2 dx
a(u, u) =
2
„¦
(3.27)
≥ {k0 |∇u| + r0 |u| } dx for all u ∈ V .
2 2
„¦

Hence, a distinction concerning r0 as in (A) with the same results
(constants) is possible.
Summarizing, we have therefore proven the following application of the
Lax“Milgram Theorem (Theorem 3.1):
Theorem 3.12 Suppose „¦ ‚ Rd is a bounded Lipschitz domain. Under
the assumptions (3.15)“(3.17) the homogeneous Dirichlet problem has one
and only one weak solution u ∈ H0 („¦).
1



(II) Mixed Boundary Conditions
‚„¦ = “2 , V = H 1 („¦)
Suppose u is a solution of (3.12), (3.19); that is, in the classical sense
let u ∈ C 2 („¦) © C 1 („¦) and the di¬erential equation (3.12) be satis¬ed
¯
pointwise in „¦ and (3.19) pointwise on ‚„¦ under the assumptions (3.13),
(3.21). However, the weaker case can again be considered, now under the
assumptions (3.14), (3.22), that u ∈ H 2 („¦) and the di¬erential equation is
satis¬ed in the sense of L2 („¦) as well as the boundary condition (3.19) in
the sense of L2 (‚„¦).
3.2. Elliptic Boundary Value Problems 105

As in (I), according to (3.11),

{K∇u · ∇v + c · ∇u v + r uv} dx +
a(u, v) := ± uv dσ (3.28)
„¦ ‚„¦

for all v ∈ V .
= b(v) := f v dx + g2 v dσ
„¦ ‚„¦

Under the assumptions (3.15), (3.22) the continuity of b and a, respec-
tively, ((3.24) and (3.25)) can easily be shown. The additional new terms
can be estimated, for instance under the assumptions (3.15), (3.22), by
the Cauchy“Schwarz inequality and the Trace Theorem (Theorem 3.4) as
follows:

g2 v dσ ¤ g2 ¤ C g2 for all v ∈ V
v|‚„¦ v
0,‚„¦ 0,‚„¦ 0,‚„¦ 1
‚„¦

and

±uv dσ ¤ ± ¤ C2 ±
u|‚„¦ v|‚„¦ u v ,
∞,‚„¦ ∞,‚„¦
0,‚„¦ 0,‚„¦ 1 1
‚„¦

respectively, for all u, v ∈ V, where C > 0 denotes the constant appearing
in the Trace Theorem.

Conditions for the V -Ellipticity of a
For the proof of the V -ellipticity we proceed similarly to (I)(B), but now
taking into account the mixed boundary conditions. For the convective
term we have
1 1 1
c · ∇u u dx = c · ∇u2 dx = ’ ∇ · c u2 dx + ν · c u2 dσ ,
2„¦ 2„¦ 2 ‚„¦
„¦

and thus
1 1
K∇u · ∇u + r ’ ∇ · c u2 dx+ ± + ν · c u2 dσ.
a(u, u) =
2 2
„¦ ‚„¦

This shows that ± + 1 ν · c ≥ 0 on ‚„¦ should additionally be assumed. If
2
r0 > 0 in (3.17), then the V -ellipticity of a follows directly. However, if only
r0 ≥ 0 is valid, then the so-called Friedrichs™ inequality, a re¬ned version
of the Poincar´ inequality, helps (see [25, Theorem 1.9]).
e
Theorem 3.13 Suppose „¦ ‚ Rd is a bounded Lipschitz domain and let
the set “ ‚ ‚„¦ have a positive (d ’ 1)-dimensional measure. Then there
˜
exists some constant CF > 0 such that for all v ∈ H 1 („¦),
1/2
¤ CF |∇v| dx
2 2
v v dσ + . (3.29)
1
˜ „¦


If ± + 1 ν · c ≥ ±0 > 0 for x ∈ “ ‚ “2 and “ has a positive (d ’ 1)-
˜ ˜
2
dimensional measure, then r0 ≥ 0 is already su¬cient for the V -ellipticity.
106 3. Finite Element Methods for Linear Elliptic Problems

Indeed, using Theorem 3.13, we have

<<

. 4
( 16)



>>