ˆb

ensured by the continuity of v . This implies the same assumption for the

ˆ

coe¬cients, since the shape functions Ni and their derivatives are continu-

ous. In order to ensure the numerical stability of a quadrature formula, it

is usually required that

ωi > 0 for all i = 1, . . . , R ,

ˆ (3.113)

which we will also do. Since all the considered ¬nite elements are such

that their faces with the enclosed degrees of freedom represent again a ¬-

nite element (in Rd’1 ) (see (3.42)), the boundary integrals are included

in a general discussion. In principle, di¬erent quadrature formulas can be

applied for each of the above integrals, but here we will disregard this pos-

sibility (with the exception of distinguishing between volume and boundary

integrals because of their di¬erent dimensions).

ˆ

A quadrature formula on K generates a quadrature formula on a general

element K, recalling

v (ˆ) dˆ | det(B)|

v(x) dx = ˆx x

ˆ

K K

by

R

ωi,K v(bi,K ) ,

i=1

152 3. Finite Element Methods for Linear Elliptic Problems

:= F (ˆi ) are dependent on

where ωi = ωi,K = ωi | det(B)| and bi = bi,K

ˆ b

K. The positivity of the weights is preserved. Here, again F (ˆ) = B x + d

x ˆ

ˆ to K. The errors of the

denotes the a¬ne-linear transformation from K

quadrature formulas

R

ωi v (ˆi ) ,

v (ˆ) dˆ ’

ˆv

E(ˆ) := ˆx x ˆ ˆb

ˆ

K i=1

(3.114)

R

v(x) dx ’

EK (v) := ωi v(bi )

K i=1

are related to each other by

EK (v) = | det(B)|E(ˆ) .

ˆv (3.115)

The accuracy of a quadrature formula will be de¬ned by the requirement

that for l as large as possible,

E(ˆ) = 0 for p ∈ Pl (K)

ˆp ˆ

ˆ

is satis¬ed, which transfers directly to the integration over K. A quadrature

formula should further provide the desired accuracy by using quadrature

nodes as less as possible, since the evaluation of the coe¬cient functions is

often expensive. In contrast, for the shape functions and their derivatives

a single evaluation is su¬cient. In the following we discuss some exam-

ples of quadrature formulas for the elements that have been introduced in

Section 3.3.

The most obvious approach consists in using nodal quadrature formu-

ˆ ˆˆ

las, which have the nodes a1 , . . . , aL of the reference element (K, P , Σ) as

ˆ ˆ

ˆ

quadrature nodes. The requirement of exactness in P is then equivalent to

ωi =

ˆ Ni (ˆ) dˆ ,

xx (3.116)

ˆ

K

so that the question of the validity of (3.113) remains.

ˆ

We start with the unit simplex K de¬ned in (3.47). Here, the weights

of the quadrature formulas can be given directly on a general simplex K: If

the shape functions are expressed by their barycentric coordinates »i , the

integrals can be computed by

±1 !±2 ! · · · ±d+1 ! vol (K)

±

»±1 »±2 · · · »d+1 (x) dx =

d+1

(3.117)

(±1 + ±2 + · · · + ±d+1 + d)! vol (K)

1 2 ˆ

K

(see Exercise 3.28).

If P = P1 (K) and thus the quadrature nodes are the vertices, it follows

that

1

ωi = »i (x) dx = vol (K) for all i = 1, . . . , d + 1 . (3.118)

d+1

K

3.5. The Implementation of the Finite Element Method: Part 2 153

For P = P2 (K) and d = 2 we get, by the shape functions »i (2»i ’ 1), the

weights 0 for the nodes ai and, by the shape functions 4»i »j , the weights

1

ωi = vol (K) for bi = aij , i, j = 1, . . . , 3 , i > j ,

3

so that we have obtained here a quadrature formula that is superior to

(3.118) (for d = 2). However, for d ≥ 3 this ansatz leads to negative weights

and is thus useless. We can also get the exactness in P1 (K) by a single

quadrature node, by the barycentre (see (3.52)):

d+1

1

ω1 = vol (K) and b1 = aS = ai ,

d+1 i=1

which is obvious due to (3.117).

As a formula that is exact for P2 (K) and d = 3 (see [53]) we present

R = 4, ωi = 1 vol (K), and the bi are obtained by cyclic exchange of the

4

barycentric coordinates:

√ √ √ √

5’ 5 5’ 5 5’ 5 5+3 5

, , , .

20 20 20 20

ˆ

On the unit cuboid K we obtain nodal quadrature formulas, which are

ˆ

exact for Qk (K), from the Newton“Cˆtes formulas in the one-dimensional

o

situation by

i1 id

for ˆi1 ...id =

ωi1 · · · ωid

ωi1 ...id

ˆ = ˆ ˆ b ,..., (3.119)

k k

for ij ∈ {0, . . . , k} and j = 1, . . . , d .

1

Here the ωij are the weights of the Newton“Cˆtes formula for 0 f (x)dx

ˆ o

(see [30, p. 128]). As in (3.118), for k = 1 we have here a generalization

of the trapezoidal rule (cf. (2.38), (8.31)) with the weights 2’d in the 2d

vertices. From k = 8 on, negative weights arise. This can be avoided and

the accuracy for a given number of points increased if the Newton“Cˆtes o

integration is replaced by the Gauss“(Legendre) integration: In (3.119), ij /k

has to be replaced by the jth node of the kth Gauss“Legendre formula

(see [30, p. 156] there on [’1, 1]) and analogously ωij . In this way, by

ˆ

ˆ

ˆ

(k + 1)d quadrature nodes the exactness in Q2k+1 (K), not only in Qk (K),

is obtained.

Now the question as to which quadrature formula should be chosen arises.

For this, di¬erent criteria can be considered (see also (8.29)). Here, we re-

quire that the convergence rate result that was proved in Theorem 3.29

should not be deteriorated. In order to investigate this question we have

to clarify which problem is solved by the approximation uh ∈ Vh based on

¯

quadrature. To simplify the notation, from now on we do not consider

boundary integrals, that is, only Dirichlet and homogeneous Neumann

154 3. Finite Element Methods for Linear Elliptic Problems

boundary conditions are allowed. However, the generalization should be

clear. Replacing the integrals in (3.111) and (3.111) by quadrature formu-

las R ωi v (ˆi ) leads to some approximation Ah of the sti¬ness matrix

¯

i=1 ˆ ˆ b

¯

and q h of the load vector in the form

¯ ¯

Ah = (ah (•j , •i ))i,j , q h = (bh (•i ))i ,

for i, j = 1, . . . , M . Here the •i are the basis functions of Xh (see (3.101))

without taking into account the Dirichlet boundary condition and

R

ωl,K (K∇v · ∇w)(bl,K )

ah (v, w) :=

K∈Th l=1

R

ωl,K (c · ∇vw)(bl,K ) +

+ ωl,K (rvw)(bl,K )

K∈Th l=1 K∈Th l=1

for v, w ∈ Xh , (3.120)

for v ∈ Xh .

bh (v) := ωl,K (f v)(bl,K )

K∈Th l=1

The above-given mappings ah and bh are well-de¬ned on Xh — Xh and Xh ,

respectively, if the coe¬cient functions can be evaluated in the quadrature

nodes. Here we take into account that for some element K, ∇v for v ∈

Xh can have jump discontinuities on ‚K. Thus, for the quadrature nodes

bl,K ∈ ‚K in ∇v(bl,K ) we have to choose the value “belonging to bl,K ” that

corresponds to the limit of sequences in the interior of K. We recall that

in general ah and bh are not de¬ned for functions of V . Obviously, ah is

bilinear and bh is linear. If we take into account the analysis of incorporating

the Dirichlet boundary conditions in (3.99)“(3.106), we get a system of

¯

equations for the degrees of freedom ξ = (ξ1 , . . . , ξM1 )T , which is equivalent

M1 ¯

to the variational equation on Vh for uh = i=1 ξi •i ∈ Vh :

¯

ah (¯h , v) = bh (v) ’ ah (wh , v) for all v ∈ Vh

u (3.121)

with wh according to (3.105). As has been shown in (3.109), (3.121) is

equivalent, in the sense of the total approximation uh + wh of u + w, to the

¯

¯h ∈ Vh ,

variational equation for u

ah (uh , v) = ¯h (v) := bh (v) ’ ah (Ih (w), v) for all v ∈ Vh ,

¯

¯ b (3.122)

if this system of equations is uniquely solvable.

Exercises

ˆ

3.28 Prove equation (3.117) by ¬rst proving the equation for K = K

and then deducing from this the assertion for the general simplex by

Exercise 3.18.

3.6. Convergence Rate Results in Case of Quadrature and Interpolation 155

3.29 Let K be a triangle with vertices a1 , a2 , a3 . Further, let a12 , a13 , a23

denote the corresponding edge midpoints, a123 the barycenter and |K| the

area of K. Check that the quadrature formula

3

|K|

Qh (u) := 3 u(ai ) + 8 u(aij ) + 27u(a123 )

60 i=1 i<j

computes the integral Q(u) := u dx exactly for polynomials of third

K

degree.

3.6 Convergence Rate Results in the Case of

Quadrature and Interpolation

The purpose of this section is to analyze the approximation quality of a

¯

¯

solution uh + Ih (w) according to (3.122) and thus of uh + wh according to

¯

(3.121) of the boundary value problem (3.12), (3.18)“(3.20).

Hence, we have left the ¬eld of Galerkin methods, and we have to

investigate the in¬‚uence of the errors

a ’ ah , b ’ a(w, ·) ’ bh + ah (Ih (w), ·).

¯

To this end, we consider in general the variational equation in a normed

space (V, · )

u ∈ V satis¬es for all v ∈ V ,

a(u, v) = l(v) (3.123)

and the approximation in subspaces Vh ‚ V for h > 0,

uh ∈ Vh satis¬es for all v ∈ Vh .

ah (uh , v) = lh (v) (3.124)

Here a and ah are bilinear forms on V — V and Vh — Vh , respectively, and

l, lh are linear forms on V and Vh , respectively. Then we have the following

theorem

Theorem 3.38 (First Lemma of Strang)

Suppose there exists some ± > 0 such that for all h > 0 and v ∈ Vh ,

¤ ah (v, v) ,

2

±v (3.125)

and let a be continuous in V — V .

Then, there exists some constant C independent of Vh such that

|a(v, w) ’ ah (v, w)|

u ’ uh ¤ u ’ v + sup

C inf

w

v∈Vh w∈Vh

|l(w) ’ lh (w)|

+ sup .

w

w∈Vh

(3.126)

156 3. Finite Element Methods for Linear Elliptic Problems

Proof: Let v ∈ Vh be arbitrary. Then it follows from (3.123)“(3.125) that

± uh ’ v ¤ ah (uh ’ v, uh ’ v)

2

= a(u ’ v, uh ’ v) + a(v, uh ’ v) ’ ah (v, uh ’ v)

+ lh (uh ’ v) ’ l(uh ’ v)

and moreover, by the continuity of a (cf. (3.2)),

|a(v, w) ’ ah (v, w)|

± uh ’ v ¤ M u ’ v + sup

w

w∈Vh

|lh (w) ’ l(w)|

for v ∈ Vh .

+ sup

w

w∈Vh

By means of u ’ uh ¤ u ’ v + uh ’ v and taking the in¬mum over

all v ∈ Vh , the assertion follows. 2

For ah = a and lh = l the assertion reduces to C´a™s lemma (Theo-

e

rem 2.17), which was the initial point for the analysis of the convergence

rate in Section 3.4. Here we can proceed analogously. For that purpose, the

following conditions must be ful¬lled additionally:

• The uniform Vh -ellipticity of ah according to (3.125) must be ensured.

• For the consistency errors

|a(v, w) ’ ah (v, w)|

Ah (v) := sup (3.127)

w

w∈Vh

for an arbitrarily chosen comparison function v ∈ Vh and for

|l(w) ’ lh (w)|

sup

w

w∈Vh

the behavior in h must be analyzed.

The ¬rst requirement is not crucial if only a itself is V -elliptic and Ah

tends suitably to 0 for h ’ 0 :

Lemma 3.39 Suppose the bilinear form a is V -elliptic and there exists

some function C(h) with C(h) ’ 0 for h ’ 0 such that

Ah (v) ¤ C(h) v for v ∈ Vh .

¯

Then there exists some h > 0 such that ah is uniformly Vh -elliptic for

¯

h ¤ h.

Proof: By assumption, there exists some ± > 0 such that for v ∈ Vh ,

¤ ah (v, v) + a(v, v) ’ ah (v, v)

2

±v

and

|a(v, v) ’ ah (v, v)| ¤ Ah (v) v ¤ C(h) v 2

.

3.6. Convergence Rate Results in Case of Quadrature and Interpolation 157

¯ ¯

Therefore, for instance, choose h such that C(h) ¤ ±/2 for h ¤ h. 2

We concretely address the analysis of the in¬‚uence of numerical quadra-

ture, that is, ah is de¬ned as in (3.120) and lh corresponds to ¯h in (3.122)

b

with the approximate linear form bh according to (3.120). Since this is an

extension of the convergence results (in · 1 ) given in Section 3.4, the as-

sumptions about the ¬nite element discretization are as summarized there

at the beginning. In particular, the triangulations Th consist of elements

that are a¬ne equivalent to each other. Furthermore, for a simpli¬cation of

the notation, let again d ¤ 3 and only Lagrange elements are considered. In

particular, let the general assumptions about the boundary value problems

which are speci¬ed at the end of Section 3.2.1 be satis¬ed.

According to Theorem 3.38, the uniform Vh -ellipticity of ah must be

ensured and the consistency errors (for an appropriate comparison element

v ∈ Vh ) must have the correct convergence behavior. If the step size h is

small enough, the ¬rst proposition is implied by the second proposition

by virtue of Lemma 3.39. Now, simple criteria that are independent of this

restriction will be presented. The quadrature formulas satisfy the properties

(3.112), (3.113) introduced in Section 3.5; in particular, the weights are

positive.

Lemma 3.40 Suppose the coe¬cient function K satis¬es (3.16) and let

c = 0 in „¦, let |“3 |d’1 > 0, and let r ≥ 0 in „¦. If P ‚ Pk (K) for the

ansatz space and if the quadrature formula is exact for P2k’2 (K), then ah

is uniformly Vh -elliptic.

Proof: Let ± > 0 be the constant of the uniform positive de¬niteness of

K(x). Then we have for v ∈ Vh :

R

ah (v, v) ≥ ± ωl,K |∇v|2 (bl,K ) = ± |∇v|2 (x) dx = ±|v|2 ,

1

„¦

K∈Th l=1

since |∇v|2 ∈ P2k’2 (K). The assertion follows from Corollary 3.14. 2

K

Further results of this type can be found in [9, pp. 194]. To investigate

the consistency error we can proceed similarly to the estimation of the

interpolation error in Section 3.4: The error is split into the sum of the errors

over the elements K ∈ Th and there transformed by means of (3.115) into

ˆ

the error over the reference element K. The derivatives (in x) arising in the

ˆ

ˆ

error estimation over K are backtransformed by using Theorem 3.26 and

Theorem 3.27, which leads to the desired hK -factors. But note that powers

of B ’1 or similar terms do not arise. If the powers of det(B) arising in

both transformation steps cancel each other (which will happen), in this

way no condition about the geometric quality of the family of triangulations

arises. Of course, these results must be combined with estimates for the

158 3. Finite Element Methods for Linear Elliptic Problems

approximation error of Vh , for which, in particular, both approaches of

Section 3.4 (either regularity or maximum angle condition) are admissible.

For the sake of simplicity, we restrict our attention in the following to the

case of the polynomial ansatz space P = Pk (K). More general results of

similar type, in particular for triangulations with the cuboid element and

ˆ ˆ

P = Qk (K) as reference element, are summarized in [9, p. 207].

We recall the notation and the relations introduced in (3.114), (3.115)

for the local errors. In the following theorems we make use of the Sobolev

spaces W∞ on „¦ and on K with the norms · l,∞ and · l,∞,K , respectively,

l

and the seminorms | · |l,∞ and | · |l,∞,K , respectively. The essential local

assertion is the following: