’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

$$$

rrrr 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