d

d d

d d d

d

d d

d

d d

d d d

d

d d

d

d d

d d d

d d d

Figure 2.9. Possibilities of Friedrichs“Keller triangulation.

In both cases, a node aZ belongs to six elements, and consequently, it

has at most six neighbours:

2.2. Linear Elements 65

aNW aN

d II d

for (a): for (b):

d Id

III aZ d

d

aW aE

d IV dVI

d Vd

d d

aS aSE

Figure 2.10. Support of the basis function.

Case (a) becomes case (b) by the transformation x ’ a ’ x, y ’ y.

This transformation leaves the di¬erential equation or the weak formula-

tion, respectively, unchanged. Thus the Galerkin method with the ansatz

space Vh according to (2.27) does not change, because P1 is invariant with

respect to the above transformation. Therefore, the discretization matrices

Ah according to (2.34) are seen to be identical by taking into account the

renumbering of the nodes by the transformation.

Thus it is su¬cient to consider only one case, say (b). A node which is

far away from the boundary has 6 neighbouring nodes in {a1 , . . . , aM1 },

a node close to the boundary has less. The entries of the matrix in the

row corresponding to aZ depend on the derivatives of the basis function

•Z as well as on the derivatives of the basis functions corresponding to the

neighbouring nodes. The values of the partial derivatives of •Z in elements

having the common vertex aZ are listed in Table 2.1, where these elements

are numbered according to Figure 2.10.

I II III IV V VI

’h ’h

1 1 1 1

‚1 •Z 0 0

h h

’h ’h

1 1 1 1

‚2 •Z 0 0

h h

Table 2.1. Derivatives of the basis functions.

Thus for the entries of the matrix in the row corresponding to aZ we

have

2 2 2

|∇•Z | dx = 2

(Ah )Z,Z=a(•Z , •Z ) = (‚1 •Z ) + (‚2 •Z ) dx,

I∪...∪VI I∪II∪III

because the integrands are equal on I and IV, on II and V, and on III and

VI. Therefore

(‚2 •Z )2 dx = 2h’2 h2 + 2h’2 h2 = 4 ,

(‚1 •Z )2 dx + 2

(Ah )Z,Z = 2

I∪III I∪II

∇•N · ∇•Z dx

(Ah )Z,N = a (•N , •Z ) =

I∪II

’h’1 h’1 dx = ’1 ,

= ‚2 •N ‚2 •Z dx =

I∪II I∪II

66 2. Finite Element Method for Poisson Equation

because ‚1 •Z = 0 on II and ‚1 •N = 0 on I. The element I for •N corre-

sponds to the element V for •Z ; i.e., ‚1 •N = 0 on I, analogously, it follows

that ‚2 •N = h’1 on I ∪ II. In the same way we get

(Ah )Z,E = (Ah )Z,W = (Ah )Z,S = ’1

as well as

(Ah )Z,NW = a (•NW , •Z ) = ‚1 •NW ‚1 •Z + ‚2 •NW ‚2 •Z dx = 0 .

II∪III

The last identity is due to ‚1 •NW = 0 on III and ‚2 •NW = 0 on III,

because the elements V and VI for •Z agree with the elements III and II

for •NW , respectively.

Analogously, we obtain for the remaining value

(Ah )Z,SE = 0 ,

such that only 5 (instead of the maximum 7) nonzero entries per row exist.

The way of assembling the sti¬ness matrix described above is called node-

based assembling. However, most of the computer programs implementing

the ¬nite element method use an element-based assembling, which will be

considered in Section 2.4.

If the nodes are numbered rowwise analogously to (1.13) and if the equa-

tions are divided by h2 , then h’2 Ah coincides with the discretization matrix

(1.14), which is known from the ¬nite di¬erence method. But here the

right-hand side is given by

h’2 (qh )i = h’2 f •i dx = h’2 f •i dx

„¦ I∪...∪VI

for aZ = ai and thus it is not identical to f (ai ), the right-hand side of the

¬nite di¬erence method.

However, if the trapezoidal rule, which is exact for g ∈ P1 , is applied to

approximate the right-hand side according to

3

1

g(x) dx ≈ vol (K) g(ai ) (2.38)

3

K i=1

for a triangle K with the vertices ai , i = 1, 2, 3 and with the area vol (K),

then

11 2 1

f •i dx ≈ h (f (aZ ) · 1 + f (aO ) · 0 + f (aN ) · 0) = h2 f (aZ ).

32 6

I

Analogous results are obtained for the other triangles, and thus

h’2 f •i dx ≈ f (aZ ) .

I∪...∪VI

In summary, we have the following result.

2.2. Linear Elements 67

Lemma 2.13 The ¬nite element method with linear ¬nite elements on a

triangulation according to Figure 2.9 and with the trapezoidal rule to ap-

proximate the right-hand side yields the same discretization as the ¬nite

di¬erence method from (1.7), (1.8).

We now return to the general formulation (2.21)“(2.24). The approach

of the Ritz method (2.24), instead of the Galerkin method (2.23), yields an

identical approximation because of the following result.

Lemma 2.14 If a is a symmetric and positive bilinear form and b is a

linear form, then the Galerkin method (2.23) and the Ritz method (2.24)

have identical solutions.

2

Proof: Apply Lemma 2.3 with Vh instead of V .

Hence the ¬nite element method is the Galerkin method (and in our

problem the Ritz method, too) for an ansatz space Vh with the following

properties:

• The coe¬cients have a local interpretation (here as nodal values).

The basis functions have a small support such that:

• the discretization matrix is sparse,

• the entries of the matrix can be assembled locally.

Finally, for the boundary value problem (2.1), (2.2) with the correspond-

ing weak formulation, we consider other ansatz spaces, which to some extent

do not have these properties:

(1) In Section 3.2.1, (3.28), we will show that mixed boundary conditions

need not be included in the ansatz space. Then we can choose the ¬-

nite dimensional polynomial space Vh = span 1, x, y, xy, x2 , y 2 , . . .

for it. But in this case, Ah is a dense matrix and ill-conditioned. Such

ansatz spaces yield the classical Ritz“Galerkin methods.

(2) Let Vh = span {•1 , . . . , •N } and let •i ≡ 0 satisfy, for some »i ,

for all v ∈ V ,

a(•i , v) = »i •i , v 0

i.e., the weak formulation of the eigenvalue problem

’∆u = »u in „¦ ,

u= 0 on ‚„¦ ,

for which eigenvalues 0 < »1 ¤ »2 ¤ . . . and corresponding eigen-

functions •i exist such that •i , •j 0 = δij (e.g., see [12, p. 335]). For

special domains „¦, (»i , •i ) can be determined explicitly, and

(Ah )ij = a(•j , •i ) = »j •j , •i = »j δij

0

68 2. Finite Element Method for Poisson Equation

is obtained. Thus Ah is a diagonal matrix, and the system of equations

Ah ξ = q h can be solved without too great expense. But this kind of

assembling is possible with acceptable costs for special cases only.

(3) The (spectral) collocation method consists in the requirement that

the equations (2.1), (2.2) be satis¬ed only at certain distinct points

xi ∈ „¦, called collocation points, for a special polynomial space Vh .

The above examples describe Galerkin methods without having the typical

properties of a ¬nite element method.

2.3 Stability and Convergence of the Finite

Element Method

We consider the general case of a variational equation of the form (2.21)

and the Galerkin method (2.23). Here let a be a bilinear form, which is not

necessarily symmetric, and let b be a linear form.

Then, if

e := u ’ uh (∈ V )

denotes the error, the important error equation

for all v ∈ Vh

a(e, v) = 0 (2.39)

is satis¬ed. To obtain this equation, it is su¬cient to consider equation

(2.21) only for v ∈ Vh ‚ V and then to subtract from the result the

Galerkin equation (2.23).

If, in addition, a is symmetric and positive de¬nite, i.e.,

a(u, u) ≥ 0 , a(u, u) = 0 ” u = 0

a(u, v) = a(v, u) ,

(i.e., a is a scalar product), then the error is orthogonal to the space Vh

with respect to the scalar product a.

Therefore, the relation (2.39) is often called the orthogonality of the error

(to the ansatz space). In general, the element uh ∈ Vh with minimal distance

to u ∈ V with respect to the induced norm · a is characterized by (2.39):

Lemma 2.15 Let Vh ‚ V be a subspace, let a be a scalar product on V ,

and let u a := a(u, u)1/2 be the norm induced by a. Then for uh ∈ Vh , it

follows that

a(u ’ uh , v) for all v ∈ Vh ”

= 0 (2.40)

u ’ uh u’v v ∈ Vh .

= min (2.41)

a a

Proof: For arbitrary but ¬xed u ∈ V , let b(v) := a(u, v) for v ∈ Vh .

Then b is a linear form on Vh , so (2.40) is a variational formulation on Vh .

2.3. Stability and Convergence 69

According to Lemma 2.14 or Lemma 2.3, this variational formulation has

the same solutions as

min F (v) v ∈ Vh

F (uh ) =

1 1

a(v, v) ’ b(v) = a(v, v) ’ a(u, v) .

with F (v) :=

2 2

Furthermore, F has the same minima as the functional

1/2 1/2

a(v, v) ’ 2a(u, v) + a(u, u)

2F (v) + a(u, u) =

1/2

a(u ’ v, u ’ v) = u’v

= ,

a

because the additional term a(u, u) is a constant. Therefore, F has the

2

same minima as (2.41).

If an approximation uh of u is to be sought exclusively in Vh , then the

element uh , determined by the Galerkin method, is the optimal choice with

respect to · a .

A general, not necessarily symmetric, bilinear form a is assumed to satisfy

the following conditions, where · denotes a norm on V :

• a is continuous with respect to · ; i.e., there exists M > 0 such that

|a(u, v)| ¤ M u for all u, v ∈ V ;

v (2.42)

• a is V -elliptic; i.e., there exists ± > 0 such that

a(u, u) ≥ ± u for u ∈ V .

2

(2.43)

If a is a scalar product, then (2.42) with M = 1 and (2.43) (as equality)

with ± = 1 are valid for the induced norm · := · a due to the

Cauchy“Schwarz inequality.

The V -ellipticity is an essential condition for the unique existence of a

solution of the variational equation (2.21) and of the boundary value prob-

lem described by it, which will be presented in more detail in Sections 3.1

and 3.2. It also implies ” without further conditions ” the stability of the

Galerkin approximation.

Lemma 2.16 The Galerkin solution uh according to (2.23) is stable in

the following sense:

1

uh ¤ b independently of h , (2.44)

±

where

|b(v)|

v∈V ,v=0

b := sup .

v

70 2. Finite Element Method for Poisson Equation

Proof: In the case uh = 0, there is nothing to prove. Otherwise, from

a(uh , v) = b(v) for all v ∈ Vh , it follows that

|b(uh )|

¤ a(uh , uh ) = b(uh ) ¤ uh ¤ b

2

± uh uh .

uh

2

Dividing this relation by ± uh , we get the assertion.

Moreover, the approximation property (2.41) holds up to a constant:

Theorem 2.17 (C´a™s lemma)

e

Assume (2.42), (2.43). Then the following error estimate for the Galerkin

solution holds:

M

u ’ uh ¤ min u ’ v v ∈ Vh . (2.45)

±

Proof: If u ’ uh = 0, then there is nothing to prove. Otherwise, let

v ∈ Vh be arbitrary. Because of the error equation (2.39) and uh ’ v ∈ Vh ,

a(u ’ uh , uh ’ v) = 0 .

Therefore, using (2.43) we have

± u ’ uh ¤ a(u ’ uh , u ’ uh ) = a(u ’ uh , u ’ uh ) + a(u ’ uh , uh ’ v)

2

= a(u ’ uh , u ’ v) .

Furthermore, by means of (2.42) we obtain

± u ’ uh ¤ a(u ’ uh , u ’ v) ¤ M u ’ uh u ’ v for arbitrary v ∈ Vh .

2

Thus the assertion follows by division by ± u ’ uh . 2

Therefore also in general, in order to get an asymptotic error estimate

in h, it is su¬cient to estimate the best approximation error of Vh , i.e.,

u’v v ∈ Vh .

min

However, this consideration is meaningful only in those cases where M/±

is not too large. Section 3.2 shows that this condition is no longer satis¬ed

for convection-dominated problems. Therefore, the Galerkin approach has

to be modi¬ed, which will be described in Chapter 9.

We want to apply the theory developed up to now to the weak formula-

tion of the boundary value problem (2.1), (2.2) with V according to (2.7)

or (2.20) and Vh according to (2.27). According to (2.4) the bilinear form

a and the linear form b read as

∇u · ∇v dx ,

a(u, v) = b(v) = f v dx .

„¦ „¦

In order to guarantee that the linear form b is well-de¬ned on V, it is su¬-

cient to assume that the right-hand side f of the boundary value problem

belongs to L2 („¦).

2.3. Stability and Convergence 71

Since a is a scalar product on V ,

1/2

|∇u| dx 2

u=u =

a

„¦

is an appropriate norm. Alternatively, the norm introduced in (2.19) for

1

V = H0 („¦) can be taken as

1/2

|u(x)| dx + |∇u(x)| dx

2 2

u = .

1

„¦ „¦

In the latter case, the question arises whether the conditions (2.42) and

(2.43) are still satis¬ed. Indeed,

|a(u, v)| ¤ u ¤u for all u, v ∈ V .

v v

a a 1 1

The ¬rst inequality follows from the Cauchy“Schwarz inequality for the

scalar product a, and the second inequality follows from the trivial estimate

1/2

|∇u(x)| dx ¤u for all u ∈ V .

2

u =

a 1

„¦

Thus a is continuous with respect to · 1 with M = 1.

The V -ellipticity of a, i.e., the property

≥± u for some ± > 0 and all u ∈ V ,

2 2

a(u, u) = u a 1

is not valid in general for V = H 1 („¦). However, in the present situation

1

of V = H0 („¦) it is valid because of the incorporation of the boundary

condition into the de¬nition of V :

Theorem 2.18 (Poincar´) Let „¦ ‚ Rn be open and bounded. Then a

e

constant C > 0 exists (depending on „¦) such that

1/2

¤C |∇u(x)| dx for all u ∈ H0 („¦) .

2 1

u 0

„¦

2

Proof: Cf. [13]. For a special case, see Exercise 2.5.

Thus (2.43) is satis¬ed, for instance with

1

±= ,

1 + C2

(see also (3.26) below) and thus in particular