Lectures on The Finite Element Method
By
Ph. Ciarlet
Tata Institute of Fundamental Research Bombay 1975
Lectures on The Finite Element Method
By
Ph. Ciarlet
Notes by
S. Kesavan, Akhil Ranjan M. Vanninathan
Tata Institute of Fundamental Research Bombay 1975
c Tata Institute of Fundamental Research, 1975
No part of this book may be reproduced in any form by print, microfilm or any other means without written permission from the Tata Institute of Fundamental Research, Colaba, Bombay 400005
Printed In India By Anil D. Ved At Prabhat Printers, Bombay 400004 And Published By The Tata Institute of Fundamental Research, Bombay.
Acknowledgements THE AUTHOR EXPRESSES his deep gratitude to Professor Ramanathan for his very kind invitation to deliver these lectures at Bangalore, as part of the T.I.F.R.–I.I.Sc. Mathematics Programme. Also, the author tanks the T.I.F.R. for financial support.
Ph. Ciarlet
v
Preface OUR BASIC AIM has been to present some of the mathematical aspects of the finite element method, as well as some applications of the finite element method for solving problems in Elasticity. This is why important topics, such as curved boundaries, mixed and hybrid methods, timedependent problems, etc..., are not covered here. No attempt has been made to give an exhaustive bibliography.
vii
Contents Acknowledgements
v
Preface
vii
1 The Abstract Problem
1
2 Examples
9
3 The Finite Element Method in its Simplest Form
29
4 Examples of Finite Elements
35
5 General Properties of Finite Elements
53
6 Interpolation Theory in Sobolev Spaces
59
7 Applications to SecondOrder Problems...
67
8 Numerical Integration
77
9 The Obstacle Problem
95
10 Conforming Finite Element Method for the Plate Problem 103 11 NonConforming Methods for the Plate Problem
ix
113
Chapter 1
The Abstract Problem SEVERAL PROBLEMS IN the theory of Elasticity boil down to the 1 solution of a problem described, in an abstract manner, as follows: Let V be a normed linear space over R. Let J : V → R be a functional which can be written in the form (1.1)
J(v) =
1 a(v, v) − f (v) 2
for all
v ∈ V,
where a(·, ·) is a continuous, symmetric bilinear form on V and f is an element of V ′ , the dual of V. Then the problem consists in finding an element u ∈ V such that (1.2)
J(u) = Min J(v). v∈V
Usually J represents the energy of some physical system. More often, instead of minimising J over the entire space V, we do so over a nonempty convex subset K of V and find a element u ∈ K such that (1.3)
J(u) = Min J(v). v∈K
Henceforth we shall denote this abstract problem by the symbol (P). One can ask immediately whether this problem admits of a solution and if so, is the solution unique? We present in this section the essential results regarding existence and uniqueness. 1
1. The Abstract Problem
2
Definition 1.1. Let V be a normed linear space. A bilinear form a(·, ·) on V is said to be Velliptic if there exists a constant α > 0 such that for all v ∈ V. a(v, v) ≥ αv2 .
(1.4) 2
Theorem 1.1. Let V be a Banach space and K a closed convex subset of V. Let a(·, ·) be Velliptic. Then there exists a unique solution for the problem (P). Further this solution is characterised by the property: (1.5)
a(u, v − u) ≥ f (v − u)
for all
v ∈ K.
Remark 1.1. The inequalities (1.5) are known as variational inequalities. Proof. The Vellipticity of a(·, ·) clearly implies that if a(v, v) = 0 then v = 0. This together with the symmetry and bilinearity of a(·, ·) shows that a(·, ·) defines an innerproduct on V. Further the continuity and the Vellipticity of a(·, ·) shows that the norm 1
v ∈ V → a(v, v) 2
(1.6)
defined by the innerproduct is equivalent to the existing norm on V. Thus V acquires the structure of a Hilbert space and we apply the Riesz representation theorem to obtain the following: for all f ∈ V ′ , there exists σ f ∈ V such that (1.7)
f (v) = a(σ f, v)
for all v ∈ V.
The map σ : V ′ → V given by f 7→ σ f is linear. Now, 1 a(v, v) − f (v) 2 1 = a(v, v) − a(σ f, v) 2
J(v) =
1. The Abstract Problem
3
1 1 = a(v − σ f, v − σ f ) − a(σ f, σ f ). 2 2 3
The symmetry of a(·, ·) is essential in obtaining the last equality. For a given f , since σ f is fixed, J is minimised if and only if a(v − σ f, v − σ f ) is minimised. But this being the distance between v and σ f , our knowledge of Hilbert space theory tells us that since K is a closed convex subset, there exists a unique element u ∈ K such that this minimum is obtained. This proves the existence and uniqueness of the solution, which is merely the projection of σ f over K. We know that this projection is characterised by the inequalities: (1.8)
a(σ f − u, v − u) ≤ 0 for all
v ∈ K.
Geometrically, this means that the angle between the vectors (σ f −u) and (v − u) is obtuse. See Fig. 1.1.
Figure 1.1: Thus, a(σ f, v − u) ≤ a(u, v − u) which by virtue of (1.7) is precisely the relation (1.5). This completes the proof. We can state the following Corollary 1.1. (a) If K is a nonempty closed convex cone with vertex at origin 0, then the solution of (P) is characterised by: a(u, v) ≥ f (v) for all v ∈ K (1.9) a(u, u) = f (u).
1. The Abstract Problem
4
(b) If K is a subspace of V then the solution is characterised by (1.10)
a(u, v) = f (v)
for all
v ∈ K.
Remark 1.2. The relations (1.5), (1.9) and (1.10) are all called variational formulations of the problem (P). Proof.
(a) If K is a cone with vertex at 0, then for u, v ∈ K, u + v ∈ K. (cf. Fig. 1.2). If u is the solution to (P), then for all v ∈ K applying (1.5) to (u + v) we get a(u, v) ≥ f (v) for all v ∈ K. In particular this applies to u itself. Setting v = 0 in (1.5) we get −a(u, u) ≥ − f (u) which gives the reverse inequality necessary to complete the proof of (1.9). Conversely, if (1.9) holds, we get (1.5) by just subtracting one inequality from the other.
Figure 1.2: (b) Applying (a) to K, since any subspace is a cone with vertex at 0, we get (b) immediately. For if v ∈ K, then −v ∈ K and applying (1.9) both to v and −v we get (1.10). This completes the proof. Remark 1.3. The solution u of (P) corresponding to f ∈ V ′ (for a fixed a(·, ·)) defines a map V ′ → V. Since this solution is the projection of
4
1. The Abstract Problem
5
σ f on K, it follows that the above map is linear if and only if K is a subspace. The problems associated with variational inequalities are therefore nonlinear in general. Exercise 1.1. Let V be as in Theorem 1.1. For f1 , f2 ∈ V ′ , let u1 , u2 5 be the corresponding solutions of (P). If  · ∗ denotes the norm in V ′ , prove that 1 u1 − u2  ≤  f1 − f2 ∗ . α Remark 1.4. The above exercise shows, in particular, the continuous dependence of the solution of f , in the sense described above. This together with the existence and uniqueness establishes that the problem (P) is “wellposed” in the terminology of partial differential equations. Exercise 1.2. If V is a normed linear space, K a given convex subset of V and J : V → R any functional which is once differentiable everywhere, then (i) if u ∈ K is such that J(u) = Min J(v), u satisfies, J ′ (u)(v − u) ≥ 0 v∈K
for all v ∈ K. (ii) Conversely, if u ∈ K such that J ′ (u)(v − u) ≥ 0 for all v ∈ K, and J is everywhere twice differentiable with J ′′ satisfying J ′′ (v)(w, w) ≥ αw2 , for all v, w ∈ K and some α ≥ 0, then J(u) = Min J(v). v∈K
Exercise 1.3 (1 ). Apply the previous exercise to the functional 1 a(v, v) = f (v) 2 with a(·, ·) and f as in Theorem 1.1. If K is a subspace of V, show that J ′ (u)(v) = 0 for all v ∈ K. In particular if K = V, J ′ (u) = 0. J(v) =
It was essentially the symmetry of the bilinear form which provided the Hilbert space structure in Theorem 1.1. We now drop the symmetry assumption on a(·, ·) but we assume V to be a Hilbert space. In addition we assume that K = V. Theorem 1.2 (LAXMILGRAM LEMMA). Let V be a Hilbert space. 6 1
Exercises (1.2)(i) and 1.3 together give relations (1.5)
1. The Abstract Problem
6
a(·, ·) a continuous, bilinear, Velliptic form, f ∈ V ′ . If (P) is the problem: to find u ∈ V such that for all v ∈ V, (1.11)
a(u, v) = f (v),
then (P) has a unique solution in V 1 . Proof. Since a(·, ·) is continuous and Velliptic, there are constants M, α > 0 such that a(u, v) ≤ Mu v,
(1.12)
a(v, v) ≥ αv2 ,
for all u, v ∈ V. Fix any u ∈ V. Then the map v 7→ a(u, v) is continuous and linear. Let us denote it by Au ∈ V ′ . Thus we have a map A : V → V ′ defined by u 7→ Au. (1.13)
 Au ∗ = sup v∈V v,0
 Au(v) a(u, v) = sup ≤ Mu. v v v∈V v,0
Thus A is continuous and A ≤ M. We are required to solve the equation (1.14)
Au = f.
Let τ be the Riesz isometry, τ : V ′ → V so that (1.15)
f (v) = ((τ f, v)),
where ((·, ·)) denotes the inner product in V. Then, Au = f if and only if τ Au = τ f or equivalently. u = u − ρ(τ Au −τ f ),
(1.16) 7
where ρ > 0 is a constant to be specified. We choose ρ such that g : V → V is a contraction map, where g is defined by (1.17) 1
g(v) = v − ρ(τAv − τ f )
cf. Corollary 1.1(b).
for
v ∈ V.
1. The Abstract Problem
7
Then the solution to (P) will be the unique fixed point of this contraction map, which exists by the contraction mapping theorem. Let v1 , v2 ∈ V. Set v = v1 − v2 . Then g(v1 ) − g(v2 ) = (v1 − v2 ) − ρτA(v1 − v2 ) = v − ρτAv.
But, v − ρτAv2 = ((v − ρτAv, v − ρτAv))
= v2 − 2ρ((τAv, v)) + ρ2 τAv2 2
= v2 − 2ρAv(v) + ρ2 Av∗
≤ v2 − 2ραv2 + ρ2 M 2 v2
= (1 − 2ρα + ρ2 M 2 )v2
2α since Av(v) = a(v, v) ≥ αv2 and A ≤ M. Choosing ρ ∈]0, 2 [, we M get that (1.18)
1 − 2ρα + ρ2 M 2 < 1
and hence g is a contraction, thus completing the proof.
Remark 1.5. The problem (P) of Theorem 1.2 is wellposed. The existence and uniqueness were proved in the theorem. For the continuous dependence of u on f , we have (1.19)
αu2 ≤ a(u, u) = f (u) ≤  f ∗ · u. 8
REFERENCE. For Variational Inequalities, see Lions and Stampacchia [18].
Chapter 2
Examples WE GIVE IN this section several examples of the abstract problem for 9 mulated in Sec. 1. We interpret the solutions of these problems as solutions of classical boundary value problems which often occur in the theory of Elasticity. Before we proceed with the examples, we summarize briefly the results (without proofs) on Sobolev spaces which will prove to be very useful in our discussion. Henceforth Ω ⊂ Rn will denote an open set (more often Ω will be a bounded open set with a specific type of boundary which will be described presently). A multiindex α will denote an ntuple (α1 , . . . , αn ) of nonnegative integers, and we denote (2.1)
α = α1 + · · · + αn ,
and call it the length of the multiindex. If v is a realvalued function on Ω for which all derivatives upto order m exist, for a multiindex α with α ≤ m we define (2.2)
∂αv =
∂αv . ∂αx11 . . . ∂αxnn
The space of test functions on Ω is given by (2.3) D(Ω) = v ∈ C ∞ (Ω); supp(v) is a compact subset of Ω . 9
2. Examples
10 where supp(v) = {x ∈ Ω; v(x) , 0}.
(2.4) 10
Definition 2.1. Let m ≥ 0 be an integer. Then the Sobolev space H m (Ω) is given by n o (2.5) H m (Ω) = v ∈ L2 (Ω); ∂α v ∈ L2 (Ω) for all α ≤ m , where all derivatives are understood in the sense of distributions. On H m (Ω) one can define a norm by means of the formula (2.6)
vm,Ω
21 X Z = ∂α v2 dx ,
v ∈ H m (Ω).
21 X Z = ∂α v2 dx ,
v ∈ H m (Ω).
α≤m
Ω
It is easy to check that  · m,Ω defines a norm on H m (Ω), which makes it a Hilbert space. One can also define a seminorm by (2.7)
vm,Ω
α=m
Ω
Note that since for all m ≥ 0, D(Ω) ⊂ H m (Ω), we may define, (2.8)
H0m (Ω) = D(Ω),
the closure being taken with respect to the topology of H m (Ω). Since H0m (Ω) is a closed subspace of H m (Ω), it is also a Hilbert space under the restriction of the norm  · m,Ω . We also have a stronger result: Theorem 2.1. Assume that Ω is a bounded open set. Then over H0m (Ω) the seminorm  · m,Ω is a norm equivalent to the norm  · m,Ω . This result is a consequence of the following: ´ Theorem 2.2 (POINCAREFRIEDRICHS’ INEQUALITY). If Ω is a bounded open set, there exists a constant C = C(Ω) such that, for all v ∈ H01 (Ω), (2.9)
v0,Ω ≤ Cv1,Ω .
2. Examples 11
11
Henceforth, unless specified to the contrary, the following will be our standing assumptions: Ω is a bounded open subset of Rn . If Γ is the boundary of Ω, then Γ is Lipschitz continuous in the sense of Neˇcas [20]. (Essentially, Γ can be covered by a finite number of local coordinate systems, such that in each, the corresponding portion of Γ is described by a Lipschitz continuous function). If L2 (Γ) is defined in the usual fashion using the Lipschitz continuity of Γ, one has the following result: Theorem 2.3. There exists a constant C = C(Ω) such that, for all v ∈ C ∞ (Ω), (2.10)
vL2 (Ω) ≤ Cv1,Ω .
By virtue of Theorem 2.3 we get that if v ∈ C ∞ (Ω), then its restriction to Γ is an element of L2 (Γ). Thus we have a map from the space C ∞ (Ω) equipped with the norm  · 1,Ω into the space L2 (Γ) which is continuous. We also have: Theorem 2.4. The space C ∞ (Ω) is dense in H 1 (Ω), for domains with Lipschitz continuous boundaries. Consequently, the above map may be extended to a continuous map → L2 (Γ) which we denote by trΓ . It is called the trace operator. An important result on the trace is the characterization: n o (2.11) H01 (Ω) = v ∈ H 1 (Ω); trΓ v = 0 . H 1 (Ω)
When no confusion is likely to occur we will merely write v instead of trΓ v. In fact if v is a “smooth” function then trΓ v is the restriction of v to Γ. −ν is 12 Retaining our assumption on Ω and Γ, the unit outer normal → → − defined a.e. on Γ. Let ν = (ν1 , . . . , νn ). If v is smooth then we may ∂v by define the outer normal derivative ∂ν ∂v X ∂v = . νi ∂ν i=1 ∂xi n
(2.12)
2. Examples
12
We extend this definition to v ∈ H 2 (Ω). If v ∈ H 2 (Ω), then H 1 (Ω) and hence trΓ
∂v ∈ L2 (Γ). We now define ∂xi
∂v ∈ ∂xi
∂v X ∂v = . νi trΓ ∂ν i=1 ∂xi n
(2.13)
However when there is no confusion we write it in the form of (2.12). Then one has the following characterization: ( ) ∂v 2 2 = 0 on Γ . (2.14) H0 (Ω) = v ∈ H (Ω); v = ∂ν Theorem 2.5 (GREEN’S FORMULA IN SOBOLEV SPACES). Let u, v ∈ H 1 (Ω). Then we have Z Z Z ∂v ∂u (2.15) u dx = − v dx + u vνi dγ, Ω ∂xi Ω ∂xi Γ for all 1 ≤ i ≤ n. If we assume u ∈ H 2 (Ω), we may replace u in (2.15) by
ming over all 1 ≤ i ≤ n, we get for u ∈ H 2 (Ω), v ∈ H 1 (Ω),
Z Z Z X n ∂u ∂u ∂v dx = − ∆u v dx + v dγ, Ω Ω ∂ν Ω i=1 ∂xi ∂xi
(2.16)
where ∆ =
∂u ; sum∂xi
Pn
∂2 i=1 ∂x2 i
is the Laplacian.
If both u and v are in H 2 (Ω), we may interchange the roles of u and v in (2.16). Subtracting one formula from the other, we get ! Z Z ∂v ∂u (u∆v − ∆uv)dx = u − v dγ, (2.17) ∂ν ∂ν Γ Ω 13
for u, v ∈ H 2 (Ω).
2. Examples
13
Finally replacing u by ∆u in (2.17) we get, for u ∈ H 4 (Ω), v ∈
H 2 (Ω), (2.18)
Z
Z
Z
∂v ∆u∆v dx = ∆ u v dx + ∆u dγ − ∂ν Ω Ω Γ 2
Z
Γ
∂(∆u) vdγ. ∂ν
The formulae (2.15) through (2.18) are all known as Green’s formulae in Sobolev spaces. We derive two results from these formulae. These results will be useful later. Lemma 2.1. For all v ∈ H02 (Ω), (2.19)
∆v0,Ω = v2,Ω .
Consequently over H02 (Ω), the mapping v 7→ ∆v0,Ω is a norm equivalent to the norm  · 2,Ω . Proof. Since D(Ω) is dense in H02 (Ω), it suffices to prove (2.19) for v ∈ (Ω). Let v ∈ D(Ω). Then n 2 2 X ∂2 v ∂2 v X ∂ v 2 ∆v = . 2 + 2 ∂xi ∂x2i ∂x2j 1≤i< j≤n i=1 By Green’s formula (2.15), !2 Z Z Z 2 2 ∂2 v ∂v ∂3 v ∂ v∂ v dx = − dx = dx (2.20) 2 2 2 Ω ∂xi ∂xi ∂x j Ω ∂xi ∂x j Ω ∂xi ∂x j
for, the integrals over Γ vanish for v ∈ D(Ω). (cf. (2.11)). Now (2.19) follows directly from (2.20). This proves the lemma. Lemma 2.2. Let Ω ⊂ R2 . Then for u ∈ H 3 (Ω), v ∈ H 2 (Ω), Z ∂2 v ∂2 u ∂2 v ∂2 u ∂2 v ∂2 u − 2 2 − 2 2 dx 2 ∂x1 ∂x2 ∂x2 ∂x1 Ω ∂x1 ∂x2 ∂x1 ∂x2 (2.21) ! Z 2 2 ∂ u ∂v ∂ u ∂v = − 2 + dγ, ∂τ ∂ν ∂τ∂ν ∂τ Γ where
∂ denotes the tangential derivative. ∂τ
14
2. Examples
14
−ν = (ν , ν ), → − Proof. Let → 1 2 τ = (τ1 , τ2 ) be the unit vectors along the outer normal and the tangent respectively. Without loss in generality we may assume τ1 = −ν2 , τ2 = ν1 . Also note that ν21 + ν22 = 1. The second derivatives occurring in the right hand side are defined by
(2.22)
2 ∂ u ∂2 u 2 ∂2 u ∂2 u 2 = τ + 2 τ τ τ + 1 2 ∂x1 ∂x2 ∂x22 2 ∂τ2 ∂x21 1 ∂2 u ∂2 u ∂2 u ∂2 u ν τ + ν2 τ2 . = (ν τ + ν τ ) + 2 1 ∂τ∂ν ∂x2 1 1 ∂x1 ∂x2 1 2 ∂x22 1
Using all these relations we get ∂2 u ∂v ∂2 u ∂v ∂2 u ∂v ∂2 u ∂v ν1 − 2 + = − ∂x1 ∂x2 ∂x2 ∂x22 ∂x1 ∂τ ∂ν ∂τ∂ν ∂τ 2 ∂ u ∂v ∂2 u ∂v (2.23) ν2 + − ∂x1 ∂x2 ∂x1 ∂x21 ∂x2 → − − = X ·→ ν, → − where, X = (X1 , X2 ) and ∂2 u ∂v ∂2 u ∂v − X = 1 ∂x1 ∂x2 ∂x2 ∂x22 ∂x1 (2.24) ∂2 u ∂v ∂2 u ∂v − X = 2 ∂x1 ∂x2 ∂x1 ∂x21 ∂x2 Also note that,
(2.25) div X =
∂X1 ∂X2 ∂2 v ∂2 u ∂2 u ∂2 v ∂2 u ∂2 v + =2 − 2 2 − 2 2. ∂x1 ∂x2 ∂x1 ∂x2 ∂x1 ∂x2 ∂x1 ∂x2 ∂x2 ∂x1
15
Now by Green’s formula (2.15) applied to functions vi ∈ H 1 (Ω) and to the constant function 1, Z Z ∂vi vi νi dγ. dx = Γ Ω ∂xi
2. Examples
15
−v = (v , . . . , v ), Hence summing over all i, if → 1 n Z Z −v dx = → −v · → −ν dγ. (2.26) div→ Γ
Γ
(This is known as the Gauss’ Divergence Theorem and also as the Ostrogradsky’s formula). From (2.23), (2.25) and (2.26) the result follows. With this background, we proceed to examples of the abstract problem of Sec. 1. Example 2.1. Let K = V = H01 (Ω). Let a ∈ L∞ (Ω) such that a ≥ 0 a.e. in Ω. Let f ∈ L2 (Ω). Define the bilinear form a(·, ·) and the functional f (·), by ! R Pn ∂u ∂v + auv dx, a(u, v) = Ω i=1 ∂xi ∂xi (2.27) R f (v) = f v dx. Ω The continuity of a(·, ·) and f (·) follows from the CauchySchwarz inequality. For instance (2.28)
 f (v) ≤  f 0,Ω v0,Ω ≤  f 0,Ω v1,Ω .
We now show that a(·, ·) is Velliptic. !2 Z X n ∂v 2 + av dx a(v, v) = Ω i=1 ∂xi !2 Z X n ∂v ≥ dx (since a ≥ 0 a.e. in Ω) Ω i=1 ∂xi = v21,Ω . 16
Since ·1,Ω is equivalent to ·1,Ω over V, this proves the Vellipticity. Hence by our results in Sec. 1 there exists a unique function u ∈ V such that a(u, v) = f (v) for all v ∈ V.
2. Examples
16
Interpretation of this problem: Using the above equation satisfied by u, we get Z X Z n ∂u ∂v (2.29) + auv dx = f v dx, for all v ∈ H01 (Ω). ∂x ∂x i i Ω i=1 Ω
From the inclusion D(Ω) ⊂ H01 (Ω), we get that u satisfies the equation −∆u + au = f in the sense of distributions. If we assume that u is “sufficiently smooth”, for example u ∈ H 2 (Ω), then we may apply Green’s formula (2.16), which gives Z Z Z Z ∂u f v dx. v dγ = ∆uv dx + auv dx − (2.30) Ω Γ ∂ν Ω Ω
Since v ∈ H01 (Ω), trΓ v = 0. Hence the integral over Γ vanishes. Thus we get Z (2.31) (−∆u + au − f )v dx = 0 for all v ∈ H01 (Ω). Ω
17
Varying v over H01 (Ω), we get that u satisfies the equation −∆u+au = f in Ω. Further since u ∈ H01 (Ω), we get the boundary condition u = 0 on Γ. Thus we may interpret u as the solution of the “classical” boundary value problem: −∆u + au = f in Ω, (2.32) u = 0 on Γ.
This is known as the homogeneous Dirichlet problem for the operator −∆u + au. A particular case of this equation arises in the theory of Elasticity, for which Ω ⊂ R2 and a = 0. Thus −∆u = f in Ω and u = 0 on Γ. This corresponds to the membrane problem: Consider an elastic membrane stretched over Ω and kept fixed along Γ. Let Fdx be the density of force acting on an element dx of Ω. Let u(x) be the vertical displacement of the point x ∈ Ω ⊂ R2 , measured
2. Examples
17
in the x3 direction from the (x1 , x2 )plane. If t is the ‘tension’ of the membrane, then u is the solution of the problem (2.33) where f = F/t.
−∆u = f in Ω u = 0 on Γ
Figure 2.1: Remark 2.1. To solve the problem (2.32) by the classical approach, one needs hard analysis involving Schauder’s estimates. By the above procedure viz. the variational method, we have got through with it more easily. The above problem is a typical example of a secondorder problem. 18 Exercise 2.1. The obstacle problem. Let Ω, Γ be as in example (2.1). Let X be an “obstacle” in this region. Let X ≤ 0 on Γ. Let F dx be the density of the force acting on a membrane stretched over Ω, fixed along Γ. The displacement u(x) at x in the vertical direction is the solution of the following problem: n o If V = H01 (Ω), K = v ∈ H01 (Ω); v ≥ X a.e. in Ω ,
2. Examples
18
let a(u, v) =
Z Z X n ∂u ∂v dx, f (v) = f v dx, f ∈ L2 (Ω), ∂x ∂x i i Ω Ω i=1 X ∈ H 2 (Ω), X ≤ 0 on Γ.
Show that K is a closed convex set and hence that this problem admits of a unique solution. Assuming the regularity result u ∈ H 2 (Ω) ∩ H01 (Ω), show that this problem solves the classical problem, u ≥ X in Ω, −∆u = f when u = 0 · Γ
u > X,
( f = F/t)
(We will discuss the Obstacle Problem in Sec. 9).
Figure 2.2: 19
Exercise 2.2. Let V = H 1 (Ω). Define a(·, ·), f (·) as in example 2.1. Assume further that there exists a constant a0 such that a ≥ a0 > 0 in Ω. If u0 is a given function in H 1 (Ω), define n o K = v ∈ H 1 (Ω); v − u0 ∈ H01 (Ω) o n = v ∈ H 1 (Ω); trΓ v = trΓ u0 .
2. Examples
19
Check that K is a closed convex subset. Interpret the solution to be that of the Nonhomogeneous Dirichlet problem, −∆u + au = f in Ω. u = u0 on Γ.
Example 2.2. Let K = V = H 1 (Ω). Let a ∈ L∞ (Ω) such that a ≥ a0 > 0, f ∈ L2 (Ω). Define a(·, ·) and f (·) as in example (2.1). The continuity of a(·, ·) and f (·) follow as usual. For the Vellipticity, we can no longer prove it with the seminorm  · 1,Ω as we did earlier. It is here we use the additional assumption on a, since !2 Z X n ∂v 2 + av dx a(v, v) = Ω i=1 ∂xi ≥ min(1, a0 )v21,Ω .
Thus we have a unique solution u to the abstract problem satisfying a(u, v) = f (v). If we assume again that u is “sufficiently smooth” to apply the Green’s formula (2.16), we get Z Z Z ∂u f v dx. v dγ = (−∆u + au)v dx + (2.34) Γ Γ ∂ν Ω If v ∈ D(Ω), then the integral over Γ will vanish. Thus u satisfies the equation −∆u + au = f as in Example 2.11 . However we now get a 20 different boundary condition. In example (2.1) the boundary condition was built in with the assumption u ∈ V = H01 (Ω). Now from (2.34), we may write: Z Z ∂u v dγ (2.35) (−∆u + au − f )v dx = − Ω Γ ∂ν for all v ∈ H 1 (Ω). But the left hand side of (2.35) is zero since u Rsatisfies the differential equation as above so that for all v ∈ H 1 (Ω), Γ ∂u ∂ν v dγ = 0. Thus
1 As in Example 2.1, the equation −∆u + au = f is always satisfied in the sense of distributions since D(Ω) ⊂ H 1 (Ω).
2. Examples
20
∂u = 0 on Γ, and we may interpret this problem as the classical prob∂ν lem: −∆u + au = f in Ω (2.36) ∂u = 0 on Γ. ∂ν This is a homogeneous Neumann problem.
Exercise 2.3. With K, V, a(·, ·) as in example 2.2., define Z Z f ∈ L2 (Ω), f (v) = f v dx + gv dγ, g ∈ L2 (Ω). Ω Γ
21
Show that the abstract problem leads to a solution of the nonhomogeneous Neumann problem −∆u + au = f in Ω ∂u = g on Γ. ∂ν
Remark 2.2. In these examples one may use the more general bilinear from defined by Z X n ∂u ∂v (2.37) a(u, v) = + auv dx, ai j ∂x j ∂xi Ω i, j=1
where the functions ai j ∈ L∞ (Ω) satisfy the condition that for some ν > 0, (2.38)
n X
i, j=1
ai j ξi ξ j ≥ ν
n X
ξi2
i=1
for all ξ ∈ Rn and a.e. in Ω. This is the classical ellipticity condition for second order partial differential operators. One should check (exercise!) in this case that the abstract problem leads to a solution of the boundary value problem ! n X ∂u ∂ ai j + au = f in Ω (2.39) − ∂xi ∂x j i, j=1
2. Examples
21
with the boundary condition u = 0 on Γ if K = V = H01 (Ω), n P (2.40) ∂u ai j ∂x νi = 0 on Γ if K = V = H 1 (Ω). j i, j=1
The latter boundary operator in (2.40) is called the conormal derivative associated with the partial differential operator, ! n X ∂ ∂ ai j . − ∂xi ∂x j i, j=1
Notice that the term (+ au) contributes nothing. Example 2.3. System of Elasticity. Let Ω ⊂ R3 , with Lipschitz continuous boundary Γ. Further assume that Γ can be partitioned into two portions Γ0 and Γ1 such that the dγmeasure of Γ0 is > 0. Let − −v = (v , v , v ); v ∈ H 1 (Ω), 1 ≤ i ≤ 3 and → −v = → K=V= → 0 on Γ0 . 1 2 3 i Define (2.41)
22
! 1 ∂vi ∂v j → − ǫi j ( v ) = 2 ∂x + ∂x , j i −v ) δ + 2µǫ (→ − −v ) = λ P3 ǫ (→ σi j (→ ij i j v ), k=1 kk
for 1 ≤ i, j ≤ 3. The latter relation is usually known as Hooke’s law. The constants λ(≥ 0) and µ(> 0) are known as Lame’s coefficients. We define the bilinear form a(·, ·) by, −u ,→ −v ) = a(→ (2.42) =
Z X 3
Ω i, j=1
Z
Ω
−u ǫ (→ − σi j (→ i j v )dx
−v + 2µ (λ div →
3 X
−u )ǫ (→ − ǫi j (→ i j v )dx.
i, j=1
→ − −g = (g , g , g ), g ∈ L2 (Γ), be Let f = ( f1 , f2 , f3 ), fi ∈ L2 (Ω), and → 1 2 3 i given.
2. Examples
22 Define the linear functional f (·), by, Z Z → − → → − − → −g · → −v dγ. (2.43) f( v ) = f · v dx + Ω
Γ1
The continuity of a(·, ·) and f (·) follow from the CauchySchwarz inequality. For the Vellipticity of a(·, ·) one uses the inequality −v ,→ −v ) ≥ 2µ a(→
23
Z X 3
Ω i, j=1
−v ))2 dx (ǫi j (→
and the fact that the square root of the integral appearing in the right hand side of the above inequality is a norm over the space V, equivalent −v = (v , v , v ) 7→ (P3 v 2 ) 12 . This is a nontrivial fact to the norm → 1 2 3 i=1 i 1,Ω which uses essentially the fact that Γ0 has measure > 0 and an inequality known as K¨orn’s inequality. We omit the proof here. −u ,→ −v ) = f (→ −v ) admits a unique solution. AsAgain the problem a(→ suming sufficient smoothness, we may apply Green’s formula: Z X 3
1 −u )ǫ (→ − σi j (→ i j v )dx = 2 Ω i, j=1 =
! ∂vi ∂v j → − + dx σi j ( u ) ∂x j ∂xi Ω i, j=1
Z X 3
Z X 3
Ω i, j=1
−u ) σi j (→
∂vi dx ∂x j
(since ǫi j (v) is symmetric in i and j) Z X Z X 3 3 ∂ → − (σi j ( u ))vi dx + =− σi j vi ν j dγ. Γ1 i, j=1 Ω i, j=1 ∂x j Thus, the abstract problem leads to a solution of 3 P ∂ → − − ∂x j (σi j ( u )) = fi (1 ≤ i ≤ 3) in Ω j=1 − → −u = → (2.44) 0 on Γ0 and 3 P → − σi j ( u )ν j = gi on Γ1 (1 ≤ i ≤ 3). j=1
2. Examples
23
Note also that 3 3 X X ∂ ∂ −u )) = − (σi j (→ − ∂x ∂x j j j=1 j=1
=−
3 X j=1
3 X → − → − λ ǫkk ( u )δi j + 2µǫi j ( u ) k=1
3 !! 3 X ∂ui ∂u j ∂ ∂ X ∂uk λ δi j − 2µ + ∂x j ∂xk ∂x j ∂x j ∂xi j=1 k=1
−v ) − µ∆u . = −(λ + µ)(grad div→ i i
Thus the first equation of (2.44) is equivalent to (2.45)
− −u − (λ + µ) grad div→ −u = → −µ∆→ f in Ω.
The equations (2.44) constitute the system of linear Elasticity.
Figure 2.3: 24
If we have an elastic threedimensional body fixed along Γ0 , acted on by an exterior force of density f dx and force of density g dγ along Γ1 and if σi j is the stress tensor, the displacement u satisfies (2.44); cf. Fig. 2.3. −u ,→ −v ) = f (→ −v ), viz., The relation a(→ (2.46)
Z X 3
Ω i, j=1
−u )ǫ (→ − σi j (→ i j v )dx =
Z
Ω
→ − → f · −v dx +
Z
Γ1
→ −g · → −v dγ
24
2. Examples
−v ∈ V is known as the principle of virtual work. The tensor ǫ for all → ij is the strain tensor and the tensor σi j the stress tensor. The expression − → − → − 1 → 2 a( v , v ) is the strain energy, and the functional f ( v ) is the potential energy of exterior forces. This example is of fundamental importance in that the finite element method has been essentially developed for solving this particular problem or some of its special cases (membranes, plates, shells, etc.,) and generalizations (nonlinear elasticity, etc. . . ).
2. Examples
25
Remark 2.3. The above problems are all examples of linear problems2 : 25 The map from the righthand side of the equation and of the boundary conditions to the solution u is linear. The nonlinearity may occur in three ways: (i) When K is not a subspace of V. (e.g. Exercises 2.1 and 2.4); (ii) If in Example 2.3 we have, instead of the first equality in (2.41): ! X 3 ∂vk ∂vk 1 ∂vi ∂v j → − . + + ǫi j ( v ) = 2 ∂x j ∂xi ∂x ∂x i j k=1
This is the case for instance when one derives the socalled Von Karmann’s equations of a clamped plate;
(iii) We may replace Hooke’s law (the second relations in (2.41)) by nonlinear equations connecting ǫi j and σi j , which are known as nonlinear constitutive equations. (e.g. Hencky’s law). Exercise 2.4. Let V = H 1 (Ω), and a(·, ·) and f (·) be as in example 2.2, and let n o 1 K = v ∈ H (Ω); v ≥ 0 a.e. on Γ .
Show that K is a closed convex cone with vertex 0. Using the results of Sec. 1 show that the interpretation is −∆u + au = f in Ω, ∂u ∂u u ≥ 0, ≥ 0, u = 0 on Γ. ∂ν ∂ν (This is called the SIGNORINI problem). We now examine fourthorder problems. Example 2.4. Let K = V = H02 (Ω). Define R ∆u∆v dx, a(u, v) = R Ω (2.47) f (v) = f v dx, f ∈ L2 (Ω), Ω 2
Except in Exercise 2.1.
2. Examples
26 26
for all u, v ∈ V. The continuity follows as usual. For the Vellipticity of a(·, ·) we have, Z (2.48) a(v, v) = (∆v)2 dx = ∆v20,Ω = v22,Ω . Ω
(by Lemma 2.1. Since  · 2,Ω and  · 2,Ω are equivalent on H02 (Ω), the Vellipticity follows from (2.48). Hence there exists a unique function u ∈ H02 (Ω) such that Z Z f v dx for all v ∈ H02 (Ω). ∆u∆v dx = Ω
Ω
Assuming u to be sufficiently smooth (say, u ∈ H 4 (Ω)), then by Green’s formula (2.18), Z Z Z ∂(∆u) ∂v 2 (2.50) (∆ u − f )v dx = v dγ − ∆u dγ, ∂ν Ω Γ ∂ν Γ for all v ∈ H02 (Ω). Hence by varying v over H02 (Ω), we get that u satisfies ∆2 u = f in Ω. Since u ∈ H02 (Ω), the boundary conditions are given by (2.14). Thus we interpret this problem as the classical problem 2 ∆ u = f in Ω, (2.51) ∂u u = = 0 on Γ. ∂ν
27
This is the homogeneous Dirichlet problem for the operator ∆2 . When n = 2, this is an important problem in Hydrodynamics. Here u is known as the stream function and −∆u is the vorticity. A slight modification of a(·, ·) leads to an important problem in Elasticity. Again let n = 2. Let f ∈ L2 (Ω) and if K = V = H02 (Ω), define (2.52) R f (v) = Ω f v dx, !# " R ∂2 v ∂2 u ∂2 v ∂2 u ∂2 v ∂2 u dx. − − ∆u ∆v + (1 − σ) 2 a(u, v) = Ω ∂x1 ∂x2 ∂x1 ∂x2 ∂x21 ∂x22 ∂x22 ∂x21
The integrand occurring in the definition of a(u, v) may also be written as 2 2 2 2 2∂2 u ∂2 v ∂ u ∂ v ∂ u ∂ v . (2.53) σ∆u∆v + (1 − σ) 2 2 + 2 2 + ∂x1 ∂x1 ∂x2 ∂x2 ∂x1 ∂x2 ∂x1 ∂x2
2. Examples
27
1 Usually, from physical considerations, 0 < σ < . 2 Note that (2.54)
a(v, v) = σ∆v20,Ω + (1 − σ)v22,Ω
by (2.53) and this leads to the Vellipticity of a(·, ·). By virtue of (2.21) in Lemma 2.2, we get that the relations a(u, v) = f (v) read as Z Z 2 f v dx, ∆ uv dx = (2.55) Ω
Ω
assuming sufficient smoothness of u. Thus again we get the same equation as in (2.51). Notice that the additional term in the definition of a(·, ·) has contributed nothing towards the differential equation. This latter problem is known as the clamped plate problem: Consider a plate of “small” thickness e lying on the x1 x2 plane. Let E be its Young’s modulus and σ its Poisson coefficient. Let there be a load F acting on the plate. The displacement u is the solution of (2.51), where f is given by (cf. Fig. 2.4): (2.56)
F=
Ee3 f . 12(1 − σ2 )
Figure 2.4: We will return to this problem in Sections 10 and 11.
28
28
2. Examples
Exercise 2.5. Let K = V = {v ∈ H 2 (Ω); v = 0 on Γ} = H 2 (Ω)∩ H01 (Ω), and define a(·, ·) and f (·) as in the case of the clamped plate. Assuming the Vellipticity of a(·, ·) show that the solution of the abstract problem satisfies ∆2 u = f in Ω and u = 0 on Γ. What is the other boundary condition? This is known as the problem of the simply supported plate. Exercise 2.6. Let K = V = H 2 (Ω) ∩ H01 (Ω), and Z a(u, v) = ∆u∆v dx; ZΩ Z ∂v f (v) = f v dx − λ dγ, where f ∈ L2 (Ω), λ ∈ L2 (Γ). Ω Γ ∂ν Show that we may apply the result of Sec. 1 and give an interpretation of this problem. REFERENCES. For details on Sobolev spaces, see Neˇcas [20] and Lions and Magenes [17]. For the theory of Elasticity, one may refer to Duvaut and Lions [10] and Landau and Lipschitz [14].
Chapter 3
The Finite Element Method in its Simplest Form MAINTAINING OUR ASSUMPTIONS as in the LaxMilgram Lemma 29 (Theorem 1.2 we concentrate our attention on the following problem (P): (P): To find u ∈ V such that a(u, v) = f (v) for all v ∈ V. Let Vh be a finitedimensional subspace of V. Then we may state the following problem: (Ph ): To find uh ∈ Vh such that a(uh , vh ) = f (vh ) for all vh ∈ Vh . Vh , being a finitedimensional subspace, is a Hilbert space for the norm of V. Hence by Theorem 1.2, uh exists and is unique. We try to approximate the solution u of (P) by means of solutions uh of the problem (Ph ) for various subspaces Vh . This is known as the internal approximation method. As a first step in this direction, we prove a most fundamental result: Theorem 3.1. There exists a constant C, which is independent of Vh , such that (3.1)
u − uh  ≤ C inf u − vh . vh ∈Vh
29
3. The Finite Element Method in its Simplest Form
30
Proof. If wh ∈ Vh , then a(u, wh ) = f (wh ) = a(uh , wh ). Thus for all wh ∈ Vh a(u − uh , wh ) = 0.
(3.2)
Using this and the Vellipticity of a(·, ·), we get, for all vh ∈ Vh , αu − uh 2 ≤ a(u − uh , u − uh )
= a(u − uh , u − vh ) + a(u − uh , vh − uh ) = a(u − uh , u − vh ) by (3.2) ≤ Mu − uh  u − vh .
30
Hence, u − uh  ≤ M/α u − vh  for all vh ∈ Vh . Setting C = M/α and taking infimum on the righthand side the result follows. The above result estimates the ‘error’ in the solution of (P) when instead we solve (Ph ). To get an upper bound for the error, we only need to compute inf u − vh  which is the distance of u from the subspace vh ∈Vh
Vh . This is a problem in approximation theory. Remark 3.1. If a(·, ·) is also symmetric then we observe the following: (i) J(uh ) = inf J(vh ) by Corollary 1.1 (b) vh ∈Vh
(ii) We saw that a(u − uh , wh ) = 0 for all wh ∈ Vh . Since a(·, ·) is now an inner product, we get that uh is the projection of u to the closed subspace Vh in the sense of this innerproduct. Therefore, p p √ √ αu − uh  ≤ a(u − uh , u − uh ) ≤ a(u − vh , u − vh ) ≤ Mu − vh 
for all vh ∈ √Vh . Hence the constant C in theorem 3.1 can be taken to be here M/α ≤ M/α, since the continuity and Vellipticity imply jointly that M ≥ α.
3. The Finite Element Method in its Simplest Form
31
We may now describe the finite element method (f.e.m.) in its simplest terms. The method consists in making special choices for the subspaces Vh such that the solutions uh of the problems (Ph ) converge to u. We will outline the procedure for obtaining the spaces Vh by consid 31 ering, for example, a secondorder problem. Let V = H01 (Ω) or H 1 (Ω). Let us assume Ω to be a polygonal domain in Rn . That is, Ω is a polygon in Rn . We then have the following stepbystep procedure: (i) We first establish a finite triangulation kh of the domain Ω such S that Ω = K. The sets K are called finite elements. If n = 2, K∈kh
they will be, in general, triangles. They will be tetrahedral in n = 3 and ‘nsimplices’ in any Rn . These have the further property that any side of a finite element K is either a portion of the boundary or the side of an adjacent finite element. (See Fig. 3.1).
Figure 3.1: (ii) The space Vh is such that for each vh ∈ Vh , its restriction vh K to each K belongs to some finitedimensional space Pk of real valued functions over K which are preassigned. In practice we choose PK to be a space of polynomials. (iii) We then need inclusions such as Vh ⊂ H01 (Ω) or H 1 (Ω). We establish a simple criterion to realise this.
32 32
3. The Finite Element Method in its Simplest Form
Theorem 3.2. If for every K = kh , PK ⊂ H 1 (K), and Vh ⊂ C 0 (Ω), then Vh ⊂ H 1 (Ω). If in addition v = 0 on Γ for all v ∈ Vh , then Vh ⊂ H01 (Ω). Proof. Let v ∈ Vh . Since vK ∈ L2 (K) for every K ∈ kh , it follows that v ∈ L2 (Ω). Hence to complete the proof it only remains to show that for 1 ≤ i ≤ n, there exist vi ∈ L2 (Ω) such that for each ϕ = D(Ω), we have, Z Z ∂ϕ (3.3) ϕvi dx = − v dx. (1 ≤ i ≤ n). Ω Ω ∂xi ∂v = vi and hence v ∈ H 1 (Ω). ∂xi ∂(vK) ∈ L2 (K) for However, vK ∈ PK ⊂ H 1 (K) implies that ∂xi 1 ≤ i ≤ n. Let ϕ ∈ D(Ω). Since the boundary ∂K of any K of the triangulation is Lipschitz continuous, we apply the Green’s formula (2.15) to get Z Z Z ∂ϕ ∂(vK) ϕdx = − (vK) dx + (vK)ϕνi,K dγK , (3.4) ∂xi K ∂K K ∂xi Then it will follow that
−ν = (ν , . . . , ν ) is the outer where dγK is the measure on ∂K and → K 1,K n,K normal on ∂K. Summing over all the finite elements K, we get Z X Z ∂(vK) ϕvi dx = ϕ dx ∂xi Ω K∈kh K Z (3.5) XZ ∂ϕ v dx + ϕ(vK)νi,K dγK , =− ∂x i ∂K Ω K∈k h
∂(vK) . ∂xi The summation on the righthand side of the above equation is zero for the following reasons: On the boundary Γ, since ϕ ∈ D(Ω), the integral corresponding to ∂K ∩ Γ is zero. So the problem, if any, is only on the other portions of the boundary of each K. However, these always occur as common boundaries of adjacent finite elements. The value of vK on the common boundary of two adjacent finite elements is the same (Vh ⊂ C 0 (Ω)). But where vi is the function whose restriction to each K is
33
3. The Finite Element Method in its Simplest Form
33
the outer normals are equal and opposite from orientation considerations. (See Fig. 3.2).
Figure 3.2: Hence the contributions from each K along the common boundaries cancel one another. Thus the summation yields only zero. Hence vi satisfies (3.3) for 1 ≤ i ≤ n, and clearly vi ∈ L2 (Ω). The last part of the theorem follows the characterization (2.11). Exercise 3.1. If for all K ∈ kh , PK ⊂ H 2 (K) and Vh ⊂ C 1 (Ω), then ∂v show that Vh ⊂ H 2 (Ω). Also if v = = 0 on Γ, for all v ∈ Vh , then ∂ν Vh ⊂ H02 (Ω). We finally describe the system of linear equations associated with the space Vh . Suppose {w j ; 1 ≤ j ≤ M} is a basis for Vh . Let uh be the solution of (Ph ). If uh is given by (3.6)
uh =
M X
u jw j.
j=1
then we have, since a(uh , wi ) = f (wi ) for 1 ≤ i ≤ M, (3.7)
M X j=1
a(w j , wi )u j = f (wi ),
1 ≤ i ≤ M.
To find uh , the above system of linear equations must be solved. The matrix for this system has for its (i, j)coefficient the value a(w j , wi ). 34
34
3. The Finite Element Method in its Simplest Form
Note that the symmetry of a(·, ·) implies the symmetry of the matrix and the Vellipticity says that the matrix is positive definite. In practical computations these observations are important. Since we have to handle the matrix of the system, it would be ideal of course to have a diagonal matrix. We could in principle achieve this through a GramSchmidt orthogonalisation procedure applied to the basis functions. However such a process is not feasible since it is highly “numerically unstable”. So the best we may hope for is a matrix with “a lot of” zeros in it  what is known as a sparse matrix. For example in the problem given by −∆u + au = f in Ω u = 0 on Γ the (i, j)coefficient of the matrix is Z X n ∂w j ∂wi + aw j wi dx. (3.8) a(w j , wi ) = Ω k=1 ∂xk ∂xk
The matrix will be sparse if the supports of the basis functions are as “small” as possible so that their innerproducts will be most often zero. We will study subsequently methods to achieve this. This trivial criterion extends, of course, to all types of problems.
Chapter 4
Examples of Finite Elements WE SUMMARIZE BELOW our requirements regarding the “finite el 35 ement subspaces” Vh (i) Let Ω ⊂ Rn be a polygonal domain. Let kh be a triangulation of Ω as in Sec. 3. Then Vh is a finitedimensional vector space such that for all v ∈ Vh , vK ∈ PK for every finite element K, where PK is a vector space of finite dimension. Usually, PK is a space of polynomials. This is of practical importance in computing the matrix of the system. We shall see later that it is of theoretical importance as well. Observe for the moment that if PK consists of polynomials, then we automatically have that PK ⊂ H 1 (K) or PK ⊂ H 2 (K). (ii) By Theorem 3.2, Vh ⊂ C 0 (Ω) implies that Vh ⊂ H 1 (Ω) and by Exercise 3.1 Vh ⊂ C 1 (Ω) implies that Vh ⊂ H 2 (Ω). Thus we must choose a proper basis for the “local” spaces PK such that these “global” inclusions hold. (iii) There must exist at least one basis {w j } of Vh which consists of functions with “small” support. We bear these points in mind when constructing examples of finite elements. Before we proceed we need a few definitions. 35
4. Examples of Finite Elements
36
36
Definition 4.1. An nsimplex is the convex hull in Rn of (n + 1) points n {a j }n+1 j=1 such that if a j = (ak j )k=1 and A is the matrix a11 a12 . . . a1,n+1 .. .. .. . . (4.1) A = . an1 an2 . . . an,n+1 1 1 ... 1 then det A , 0.
The above definition generalises the notion of a triangle to n dimensions. Geometrically the condition det A , 0 simply means that the points {a j }n+1 j=1 do not lie in the same hyperplane. For det A is equal to, by elementary column operations, the determinant of the matrix (a11 − a1,n+1 ) . . . (a1,n − a1,n+1 ) .. .. . . (an1 − an,n+1 ) . . . (an,n − an,n+1 )
and that this is nonzero means that (a1 − an+1 ), . . . , (an − an+1 ) are linearly independent vectors in Rn , which is the same as saying that a1 , . . . , an+1 do not lie in the same hyperplane.
n Definition 4.2. Let {a j }n+1 j=1 be (n + 1)points in R satisfying the conditions of definition 4.1. The barycentric coordinates of any x ∈ Rn with respect to these points are numbers {λ j }n+1 j=1 such that n+1 P x= λ ja j, j=1 (4.2) n+1 P λ j. 1 = j=1
The barycentric coordinates exist because they are merely the com→ − ponents of the unique solution vector λ of the system of (n + 1) linear equations in (n + 1) unknowns given by x1 ! → − x → − → − Aλ = where x = ... 1 xn
4. Examples of Finite Elements
37
The functions λ j = λ j (x) are all affine functions of x. Also λ j (ai ) = δi j , where δ is the Kronecker symbol. 37 Remark 4.1. Given points {a j }n+1 j=1 as in the definition (4.1), the corresponding nsimplex is given by n+1 n+1 X X K= x = λ a ; 0 ≤ λ ≤ 1; λ = 1 . j j j j j=1 j=1
Definition 4.3. Let k ≥ 0 be an integer. Then, Pk is the space of all polynomials of degree ≤ k in x1 , . . . , xn . We now proceed with examples of finite elements. Example 4.1. The nsimplex of type (1).
P Let K be an nsimplex. Let PK = P1 . We define a set K = {p(ai ); 1 ≤ i ≤ n + 1} of degrees of freedom for p ∈ PK , where {ai }n+1 i=1 P are the vertices of K: The set K determines every polynomial p ∈ PK uniquely. For, note that dim PK = dim P1 = n + 1. Consider λ1 , . . . , λn+1 ∈ P1 , the barycentric coordinate functions. These are linP early independent since αk λk = 0 implies that its value at each vertex is zero. Since λk (a j ) = δk j we get that α j = 0 for all j. Thus these functions form a basis for P1 . Let us write p=
n+1 X
αi λi .
i=1
Then p(a j ) =
n+1 X
αi λi (a j ) =
i=1
n+1 X i=1
Thus, (4.4)
p=
n+1 X i=1
p(ai )λi .
αi δi j = α j .
4. Examples of Finite Elements
38
Example 4.2. The nsimplex of type (2). Let K be an nsimplex with vertices {ai }n+1 i=1 . Let ai j (i < j) be the 38 midpoints of the line joining ai and a j , i.e. ai j = 12 (ai + a j ).
Figure 4.1: P Let PK = P2 . We define for p ∈ P2 , the set K = {p(ai ), 1 ≤ i ≤ P n + 1; p(ai j ), 1 ≤ i < j ≤ n + 1} of degrees of freedom. Again n+2 K determines p ∈ P2 completely. To see this note that dim PK = 2 and there are as many functions in the set {λi (2λi − 1), 1 ≤ i ≤ n+ 1; λi λ j , 1 ≤ i ≤ j ≤ n + 1}. There are all functions in P2 . Further since 1 2 if i = k or j, λi (a j ) = δi j , λi (ak j ) = 0 otherwise, we see again that these are linearly independent in P2 . Let us write p=
n+1 X i=1
αi λi (2λi − 1) +
Then p(ak ) =
n+1 X i=1
X
βi j λi λ j .
1≤i< j≤n+1
αi δik (2δik − 1) = αk .
Further, p(akl ) =
n X i=1
αi (2λ2i (akl ) − λi (akl ))
4. Examples of Finite Elements +
39
X
βi j λi (akl )λ j (akl ).
1≤i< j≤n+1
But since λi (akl ) = 0 or 12 , the first sum
n P
is zero. Further,
39
i=1
1 4 if (i, j) = (k, l) or (l, k), λi (akl )λ j (akl ) = 0 otherwise.
Hence βkl = 4p(akl ). Thus we have (4.5)
p=
n+1 X i=1
λi (2λi − 1)p(ai ) +
X
4λi λ j p(ai j ).
1≤i< j≤n+1
Example 4.3. The nsimplex of type (3). Let K be an nsimplex with vertices {ai }n+1 i=1 . Let aii j = ai + a j + ak i , j. Let ai jk = for i < j < k. 3
2ai + a j , 3
Figure 4.2: Set PK = P3 . Define the set of degrees of freedom X n = p(ai ), 1 ≤ i ≤ n + 1; p(aii j ), 1 ≤ i , j ≤ n + 1; K
o p(ai jk ), 1 ≤ i < j < k ≤ a + 1 .
4. Examples of Finite Elements
40
2 1 Note that λi (aii j ) = ; λ j (aii j ) = ; λ1 (aii j ) = 0 if 1 , i, 1 , j; 3 3 1 if 1 = i, j or k and 0 otherwise, etc. Using these, one λ1 (ai jk ) = 3 checks the linear independence of the functions n λi (3λi − 1)(3λi − 2), 1 ≤ i ≤ n + 1; λi λ j (3λi − 1), 1 ≤ i , j ≤ n + 1; o λi λ j λk 1 ≤ i < j < k ≤ n + 1 . These then form a basis for P3 , for there are as many functions in the above collection as dim PK . Using the values of λi at the special points described above, we get n+1 X λi (3λi − 1)(3λi − 2) p(ai ) p= 2 i=1 X 9 + λi λ j (3λi − 1)p(aii j ) 2 1≤i, j≤n+1 X 27λi λ j λk p(ai jk ).
(4.6)
1≤i< j
40
Thus
P
K
completely determines p ∈ P3 .
P The points of K at which the polynomials are evaluated to get K P are known as the nodes of the finite element. The set K is the set of degrees of freedom of the finite element. Exercise 4.1. Generalize these ideas and describe the nsimplex of type (k) for any integer k ≥ 1. We now show how these finite elements may be used to define the space Vh . First of all we show the inclusion Vh ⊂ C 0 (Ω). Consider for instance a triangulation by nsimplices of type (1). Number all the nodes of the P triangulation by {b j }. Let us define h = {p(b j ); b j is a node.}: This is the set of degrees of freedom of the space Vh : A function v in the space Vh is, by definition, determined over each K ∈ kh by the values v(b j ) for those nodes b j which belong to K.
4. Examples of Finite Elements
41
Let us examine the twodimensional case, for simplicity. If K1 and K2 are two adjacent triangles with common side K ′ (cf. Fig. 4.3), we need to show that vK1 = vK2 along K ′ for any v ∈ Vh . Let t be an abscissa along K ′ .
Figure 4.3: 41
K′
Now vK1 along is a polynomial of degree 1 in t. So is vK2 along K ′ . But these two agree at the nodes b1 and b3 . Therefore, they must be identical and hence the continuity of v follows. This argument can be extended to any simplex of type (k). These simplices, by Theorem 3.2 yield the inclusion Vh ⊂ H 1 (Ω) and hence we may use them for second order problems. Exercise 4.2. The triangle of type (3′ ). P Let K be a triangle in R2 . Define K to be the values of p at the points {ai , 1 ≤ i ≤ 3}, and the points {aii j , 1 ≤ i , j ≤ 3}. If we define P P P′3 = {p ∈ P3 ; 12p(a123 ) + 2 3i=1 p(ai ) − 3 i, j p(aii j ) = 0}, then show P that K uniquely determines p ∈ P′3 = PK . Further show that P2 ⊂ P′3 .
We now relax our terminological rules about “triangulations” and admit rectangles (and in higher dimensions, hyper rectangles or hypercubes) in triangulations. We describe below some finite elements which are rectangles. We need another space of polynomials.
4. Examples of Finite Elements
42
Definition 4.4. Let k ≥ 1 be an integer. Then X i1 in Qk = p; p(x) = a x . . . x . i1 ,...in 1 n 0≤i j ≤k 1≤ j≤n
42
We have the inclusions Pk ⊂ Qk ⊂ Pnk . Example 4.4. The Rectangle of type (1). Let K be the unit square in R2 , i.e., K = [0, 1]n . Let PK = Q1 . P The set of degrees of freedom is given by K = {p(ai ), 1 ≤ i ≤ 4}; cf. Fig. 4.4 in the case n = 2.
Figure 4.4: P To show that K indeed determines p ∈ Q1 uniquely we adopt a different method now. (There are essentially two methods to show that P K completely determines P K ; the first was used in the previous examples where we exhibited a basis for PK such that the corresponding P coefficients in the expansion of p in terms of this basis came from K ; the second is illustrated now). P Observe first that dim PK = card K = 2n . To determine a polynoP mial completely in terms of the elements of K we must solve 2n linear
4. Examples of Finite Elements
43
equations in as many unknowns. That every polynomial is determined this way is deduced from the existence of a solution to this system. But for such a system the existence and uniqueness of the solution are equiv 43 alent, and one establishes the latter. Thus we show that if p ∈ PK such that all its degrees of freedom are zero, the p ≡ 0. Returning to our example, consider a polynomial p ∈ Q1 such that p(ai ) = 0 for all 1 ≤ i ≤ 4. On each side p is a polynomial of degree 1 either in x1 alone or in x2 alone. Since it vanishes at two points, the polynomial p vanishes on the sides of the square. Now consider various lines parallel to one of the axes. Here too p is a polynomial of degree 1 in one variable only. Since it vanishes at the points where the line meets the side, it also vanishes on this line. Varying the line we get p ≡ 01 . Example 4.5. The Rectangle of type (2). Again consider the unit square (or hypercube in Rn ) to be the finite P element K. Set PK = Q2 , and K = {p(ai ), 1 ≤ i ≤ 9} where the ai are as in the figure below.
Figure 4.5: Here again one can prove the unisolvency as above. Now let p ∈ Q2 be given such that p(ai ) = 0 for all 1 ≤ i ≤ 9; then p = 0 on the four 1 It is not necessary to restrict ourselves to a square. Any rectangle with sides parallel to the coordinate axes would do.
4. Examples of Finite Elements
44
44
sides and on the two central (dotted, in Fig. 4.5) lines. Now take lines parallel to one of the axis and p vanishes on each of these. Thus p ≡ 0 P on K and we get that K uniquely determines p ∈ Q2 . Exercise 4.3. Describe the rectangle of type (3) and generalize to hyperrectangles of type (k).
Exercise 4.4. Prove that in all the preceding examples, we get Vh ⊂ C 0 (Ω). Exercise 4.5. The Rectangle of type (2′ ). Let K be as in example 4.5. However omit the node a9 (the centroid P of K). Let K = {p(ai ), 1 ≤ i ≤ 8} and show that this determines uniquely a function in the space 4 8 X X PK = p ∈ Q ; 4p(a ) + p(a ) − 2 p(a ) = 0 . 2 9 i i i=1
i=5
and that P2 ⊂ PK .
We now turn to different types of finite elements. They differ from the preceding ones in the choice of degrees of freedom as will be seen presently. Example 4.6. The Hermite Triangle of Type (3). Let K ⊂ R2 be a triangle with vertices {a1 , a2 , a3 }. Let λi , 1 ≤ i ≤ 3, be the barycentric coordinate functions. Then one can check that any polynomial p ∈ P3 = PK can be expanded as p=
3 X (−2λ3i + 3λ2i − 7λ1 λ2 λ3 )p(ai ) + 27λ1 λ2 λ3 p(a123 ) i=1
+
3 X X i=1 j=1 j,i
λi λ j (2λi + λ j − 1) Dp(ai )(a j − ai ).
P Thus, K = {p(ai ), 1 ≤ i ≤ 3; Dp(ai )(a j − ai ), 1 ≤ i , j ≤ 3; p(a123 )} is the corresponding set of degrees of freedom.
4. Examples of Finite Elements 45
45
Note that Dp(ai ) is the Frechet derivative of p evaluated at ai : If (e1 , . . . , en ) is the standard basis for Rn , then for v : Rn → R, we have, ∂v (x), the usual partial derivative. (Dv)(x)(ei ) = ∂xi P Notice that we may replace K by the set ′ X K
(
) ∂p ∂p = p(ai ), 1 ≤ i ≤ 3; p(a123 ); (ai ), (ai ), 1 ≤ i ≤ 3 . ∂x1 ∂x2
Remark 4.2. The term “Hermite” means that we assume knowledge of derivatives at some of the nodes. If only the values of p at the nodes appear in the set of degrees of freedom, as was the case upto Example 4.5, we refer to the finite elements as of “Lagrange” type. These ideas will be made precise in Sec. 5. We usually indicate degrees of freedom involving derivatives by circling the nodes  one circle for first derivatives, two for first and second derivatives and so on. Thus the finite element of example 4.6 may be pictured as in Fig. 4.6.
Figure 4.6: Exercise 4.6. The Hermite Triangle of Type (3′ ). This is also known ( as the Zienkiewicz triangle in) Engineering litP P ∂p ∂p erature. Set K = p(ai ), (ai ), (ai ), 1 ≤ i ≤ 3 . Show that K ∂x1 ∂x2 uniquely determines a function in the space 3 3 X X . p ∈ P ; 6p(a ) − 2 p(a ) + Dp(a )(a − a ) = 0 PK = 3 123 i i i 123 i=1
i=1
46
4. Examples of Finite Elements
All examples cited upto now yield the inclusion Vh ⊂ C 0 (Ω) and 46 consequently are useful to solve secondorder problems. In order to solve problems of fourth order, we need the inclusion Vh ⊂ C 1 (Ω). Our subsequent examples will be in this direction. Remark 4.3. Consider a 1simplex K ⊂ R1 . A triangulation is merely a subdivision of Ω into subintervals. In any subinterval K not only vK d(vK) but also must be continuous at both end points. Thus we get dx 4 conditions on vK. Consequently PK must contain all polynomials of degree 3 in it. The analogous result (which is nontrivial) is due ˇ sek [24] that is case of R2 , and K a triangle of R2 , at least to A. Zeniˇ polynomials of degree 5 must be contained in PK . Example 4.7. The Argyris triangle. This is also known as the 21degreeoffreedomtriangle. We set PK = P5 and X ∂p ∂2 p p(a ), (a ), . . . , (ai ), 1 ≤ i ≤ 3; = i i ∂x1 ∂x22 K ) ∂p (ai j ), 1 ≤ i ≤ j ≤ 3 . ∂ν
Figure 4.7: 47
∂p is indicated by a line The knowledge of the normal derivative ∂ν perpendicular to the side at the appropriate point; cf. Fig. 4.7. P We now show that any p ∈ PK is uniquely determined by K . Let p ∈ PK = P5 be given such that all its degrees of freedom are zero. If K ′
4. Examples of Finite Elements
47
is any side of K and t is an abscissa along K ′ then pK ′ is a polynomial dp d2 p , p1 (t) of degree 5. The vanishing of p, at the end points, say b, dt dt2 ′ ′ b , of K imply that all the 6 coefficients of p1 are 0 and hence p1 ≡ 0. ∂p dp on K ′ . The polynomial r(t) = (t) is of degree 4 on Thus p = 0 = dt ∂ν ! b + b′ dr dr ′ ′ ′ K and we also have r(b) = r(b ) = (b) = (b ) = r =0 dt dt 2 ∂p ∂p , all vanish on the which imply that r ≡ 0 on K ′ . Hence p, ∂x1 ∂x2 sides of the triangle K. The sides of K are defined by the equations λi (x1 , x2 ) = 0, (i = 1, 2, 3) where λi are the barycentric coordinate functions. We claim that λ2i divides p for i = 1, 2, 3. To see this it is enough ∂p ∂p to prove that if p is a polynomial such that p, , vanish on any ∂x1 ∂x2 straight line L = {(x1 , x2 ); λ(x1 , x2 ) = 0} then λ2 divides p. In the special P j case, when λ(x1 , x2 ) = x1 writing p(x1 , x2 ) = 5j=0 a j (x2 )x1 (with deg. ∂p = 0 on a j ≤ 5 − j) it follows that a0 (x2 ) = a1 (x2 ) = 0 since p = ∂x1 L. Thus x21 divides p. The general case reduces to this case by an affine transformation. In fact, by translating the origin to a point P, fixed arbitrarily on L and by rotation of the coordinate axes we can assume that L = {(X1 , X2 ); X1 = 0} in the new coordinates. If p′ is the image of p under this transformation then p′ is also a polynomial (of degree 5) and ∂p′ ∂p′ , vanish on L by chain rule for differentiation. Hence X12 dip′ , ∂X1 ∂X2 vides p′ . This is the same thing as saying λ2 divides p which proves the claim. Since λi are mutually coprime we may now write p = qλ21 λ22 λ23 . Then we necessarily have q(x1 , x2 ) ≡ 0 for, otherwise deg. p ≥ 6 which P is impossible since p ∈ P5 . Hence p ≡ 0 on K which proves that K determines p ∈ P5 . To define the corresponding space Vh , we number all the vertices of the triangles by {b j } and all midpoints of the sides by {ck }. The the set of degrees of freedom of the space Vh is 48 X ∂2 v ∂v ∂v ∂2 v ∂2 v ∂v (b ), (b ), (b ), (b ), (b ), (c ) = v(b ), j j j j j k j 2 2 ∂x1 ∂x2 ∂x1 ∂x2 ∂νk ∂x ∂x h
1
2
4. Examples of Finite Elements
48
∂ is one of the two possible normal derivatives at the midpoint where ∂νk ck . We now show that Vh ⊂ C 1 (Ω). Consider two adjacent Argyris triangles K1 and K2 with common boundary K ′ along which t is an abscissa (Fig. 4.8). Let v ∈ Vh . Now vK1 and vK2 are polynomials of degree 5 in t along K ′ and they agree together with their first and second derivatives at the end points. Thus vK1 = vK2 on K ′ , proving continuity. ∂(vK2 ) ∂(vK1 ) and along K ′ are polynomials of degree 4 Now ∂ν ∂ν agreeing in their values with first derivatives at end points and agree at ∂(vK1 ) ∂(vK2 ) the midpoint in their values. Thus = on K ′ . Similarly, ∂ν ∂ν ∂(vK1 ) ∂(vK2 ) = on K ′ and hence v ∈ C 1 (Ω). Thus Vh ⊂ C 1 (Ω). ∂t ∂t
Figure 4.8: Exercise 4.7. The 18DegreeofFreedomTriangle Let K be a triangle in R2 . Let PK consist of those polynomials of degree ≤ 5 for which, along each side of K, the normal derivative is a polynomial of degree ≤ 3, in one variable of course. Show that a polynomial in PK is uniquely determined by the following set of degrees of freedom: X ∂2 p ∂p (a ), . . . , (a ), 1 ≤ i ≤ 3 = p(a ), i i i . 2 ∂x ∂x K
49
1
Note that P4 ⊂ PK and dim PK = 18.
2
4. Examples of Finite Elements
49
Exercise 4.8. The HCTTriangle. This element is due to Hsieh, Clough and Tocher. Let a be any interior point of the triangle K with vertices a1 , a2 , a3 . With a as common vertex subdivide the triangle into triangles K1 , K2 , K3 ; cf. Fig. 4.9. Define n o PK = p ∈ C 1 (K); PKi ∈ P3 , 1 ≤ i ≤ 3 . Obviously, P3 ⊂ PK . The degrees of freedom are given by ) X ( ∂p ∂p ∂p (ai ), (ai ), 1 ≤ i ≤ 3; (ai j ), 1 ≤ i < j ≤ 3 . = p(ai ), ∂x1 ∂x2 ∂ν K
Figure 4.9: Show that
P
K
uniquely determines p ∈ PK .
Note: Since we have to determine 3 polynomials pi = pKi each of degree ≤ 3, we need to determine 30 coefficients on the whole. For this we have the following conditions: (i) The values at the vertices together with first derivatives and also the normal derivative at the mid points give 7 conditions for each pi = pKi . Thus we have 21 conditions from these. (ii) p1 (a) = p2 (a) = p3 (a) gives 2 conditions. (iii)
∂p1 ∂p2 ∂p3 (a) = (a) = (a) for i = 1, 2, gives 4 more conditions. ∂xi ∂xi ∂x1
50
4. Examples of Finite Elements
50 (iv)
∂p2 ∂p1 = along a1 a and two more similar conditions give 3 ∂ν ∂ν conditions.
Thus we have 30 conditions to determine the 30 coefficients. But, of course this is no proof, which is left as an exercise! Exercise 4.9. The BognerLoxSchmidt Rectangle; cf. Fig. 4.10.
Figure 4.10: Let PK = Q3 , the degrees of freedom being given by ) X ( ∂p ∂2 p ∂p (ai ), (ai ), (ai ), 1 ≤ i ≤ 4 . = p(ai ), ∂x1 ∂x2 ∂x1 ∂x2 K P Show that K determines uniquely a polynomial p ∈ Q3 (a double dotted arrow indicates that the mixed second derivative is a degree of freedom). Show also that in this case Vh ⊂ C 1 (Ω).
51
So far, we have verified requirements (i) and (ii) mentioned at the beginning of this Section. Let us now examine requirement (iii), which P will be fulfilled by a “canonical” choice for the basis functions. Let h be the set of degrees of freedom of the space Vh derived in an obvious P P way from the sets K , K ∈ kh ; Examples of such sets h have been given for nsimplices of type (k) and for Argyris triangles. Then if X n o = ϕ jh , 1 ≤ j ≤ M , h
4. Examples of Finite Elements
51
we let the basis functions w j , 1 ≤ j ≤ M, be those functions in the space Vh which satisfies ϕi (w j ) = δi j , 1 ≤ i ≤ M. Then it is easily seen that this choice will result in functions with “small” support: in Fig. 4.11, we have represented there types of supports encountered in this fashion, depending upon the position of the node associated support of basis function associated with node
support of basis function associated with node
support of basis function associated with node
Figure 4.11: with the degree of freedom.
Chapter 5
General Properties of Finite Elements IT WOULD HAVE been observed that upto now we have not defined fi 52 nite elements in a precise manner. Various polygons like triangles, rectangles, etc. were loosely called finite elements. We rectify this omission and make precise the ideas expressed in the previous sections. Definition 5.1. A finite element is a triple (K, Σ, P) such that (i) K ⊂ Rn with a Lipschitz continuous boundary ∂K and Int K , φ. (ii) Σ is a finite set of linear forms over C ∞ (K). The set Σ is said to be the set of degrees of freedom of the finite element. (iii) P is a finite dimensional space of realvalued functions over K N and α , 1 ≤ i ≤ N such that Σ is Punisolvent: i.e. if Σ = {ϕi }i=1 i are any scalars, then there exists a unique function p ∈ P such that (5.1)
ϕi (p) = αi , 1 ≤ i ≤ N.
Condition (iii) of definition (5.1) is equivalent to the conditions that dim P = N = card Σ and that there exists a set of functions {p j }Nj=1 with ϕi (p j ) = δi j (1 ≤ i, j ≤ N), which forms a basis of P over R. Given any 53
54
5. General Properties of Finite Elements
p ∈ P we may write (5.2)
p=
N X
ϕi (p)pi .
i=1
53
P P Instead of (K, , P) one writes at times (K, K , PK ) for the finite element. In the various examples we cited in Sec. 4 our set of degrees of freedom for a finite element K (which was an nsimplex or hyperrectangle) had elements of the following type: Type 1: ϕ0i given by p 7→ p(a0i ). The points {a0i } were the vertices, the midpoints of sides, etc... 1 ). For instance, in the Hermite Type 2: ϕ1i,k given by p 7→ Dp(a1i )(ξi,k 1 = a −a , triangle of type (3) (cf. Example 4.6), we had a1i = ai , ξi,k i k where a1 , a2 , a3 were the vertices. 2 , ξ 2 ). For example, in the 18Type 3: ϕ2i,kl given by p 7→ D2 p(a2i )(ξi,k i,l 2 = e = ξ 2 , the unit vecdegreeoffreedom triangle, a2i = ai , ξi,k 1 i,1 ∂2 p (ai ) tor in the x1 direction so that we have D2 p(ai )(e1 , e1 ) = ∂x21 as a degree of freedom. (cf. Exercise 4.7).
In all these cases the points {ais } for s = 0, 1 and 2, are points of K and are called the nodes of the finite element. Definition 5.2. A finite element is called a Lagrange finite element if its degrees of freedom are only of Type 1. Otherwise it is called a Hermite finite element. (cf. Remark 4.2) Let (K, Σ, P) be a finite element and v : K → R be a “smooth” function on K. Then by virtue of the Punisolvency of Σ, there exists a unique element, say, πv ∈ P such that ϕi (πv) = ϕi (v) for all 1 ≤ i ≤ N, N . The function πv is called the Pinterpolate function where Σ = {ϕi }i=1 of v and the operator π : C ∞ (K) → P is called the Pinterpolation
5. General Properties of Finite Elements
55
operator. If {p j }Nj=1 is a basis for P satisfying ϕi (p j ) = δi j for 1 ≤ i, j ≤ N then we have the explicit expression π(·) =
(5.3)
N X
ϕi (·)pi .
i=1
Example 5.1. In the triangle of type (1) (see Example 4.1), P = P1 , Σ = {ϕi ; ϕi (p) = p(ai ), 1 ≤ i ≤ 3} and pi = λi , the barycentric coordinate 54 functions. Thus we also have πv =
(5.4)
3 X
v(ai )λi .
i=1
Exercise 5.1. Let K be a triangle with vertices a1 , a2 and a3 . Let ai j (i < j) be the midpoint of the side joining ai and a j . Define ΣK = {p 7→ p(ai j ), 1 ≤ i ≤ j ≤ 3}. Show that Σ is P1 unisolvent and that in general Vh 1 C 0 (Ω) for a triangulation made up of such finite elements. Exercise 5.2. Let K be a rectangle in R2 with vertices a1 , a2 , a3 , a4 . Let a5 , a6 , a7 , a8 be the midpoints of the sides as in Fig. 4.5. If Σ = {p 7→ p(ai ), 5 ≤ i ≤ 8}, show that Σ is not Q1 unisolvent. Let us now consider a family of finite elements of a given type. To be more specific, we will consider for instance a family of triangles of type (2) (see Example 4.2), but our subsequent descriptions extend to all types of finite elements in all dimensions. Pick, in particular, a triangle Kˆ with vertices {ˆa1 , aˆ 2 , aˆ 3 } from this family. Let the midpoints of the sides be {ˆa12 , aˆ 23 , aˆ 13 }. Set Pˆ = PKˆ = P2 and define accordingly the associated set of degrees of freedom for Kˆ as X X ˆ = {p 7→ p(ˆai ), 1 ≤ i ≤ 3; p 7→ p(ˆai j ), 1 ≤ i < j ≤ 3}. = Kˆ
ˆ P) ˆ Σ, ˆ as fixed in the In as much as we consider the finite element (K, sequel, it will be called the reference finite element of the family. Given any finite element K with vertices a1 , a2 , a3 in this family, there exists a unique invertible affine transformation of R2 i.e. of the
5. General Properties of Finite Elements
56
55
form F K ( xˆ) = BK xˆ + bK , where BK is an invertible 2 × 2 matrix and ˆ = K and F K (ˆai ) = ai , 1 ≤ i ≤ 3. It is bK ∈ R2 , such that F K (K) easily verified that F K (ˆai j ) = ai j for 1 ≤ i < j ≤ 3. Also, the space ˆ is precisely the space PK = P2 . Hence {p : K → R; p = pˆ ◦ F K−1 , pˆ ∈ P} the family {(K, Σ, P)} is equivalently defined by means of the following data: ˆ P), ˆ Σ, ˆ (i) A reference finite element (K, ˆ = K, ai = (ii) A family of affine mappings {F K } such that F K (K) F K (ˆai ), 1 ≤ i ≤ 3, ai j = F K (ˆai j ), 1 ≤ i < j ≤ 3,
and
ΣK = {p 7→ p(F K (ˆai )); p 7→ p(F K (ˆai j ))}, ˆ PK = {p : K → R; p = pˆ ◦ F K−1 , pˆ ∈ P}.
This special case leads to the following general definition. ˆ P) ˆ Σ, ˆ and (K, Σ, P) are affine Definition 5.3. Two finite elements (K, equivalent if there exists an affine transformation F on Rn such that (i) F( xˆ) = B xˆ + b, b ∈ Rn , B an invertible n × n matrix, ˆ (ii) K = F(K), (iii) ais = F(ˆais ), s = 0, 1, 2, 1 = Bξˆ 1 , ξ 2 = Bξˆ 2 , ξ 2 = Bξˆ 2 (iv) ξi,k i,k i,k i,k i,l i,l′
and ˆ (v) P{p : K 7→ R; p = pˆ ◦ F −1 , pˆ ∈ P}. This leads to the next definition. Definition 5.4. A family {(K, ΣK , PK )} of finite elements is called an affine family if all the finite elements (K, ΣK , PK ) are equivalent to a ˆ P). ˆ Σ, ˆ single reference finite element (K,
5. General Properties of Finite Elements 56
57
Let us see why the relations given by (iv) must be precisely of that form. We have by (v), p(x) = p( ˆ xˆ). This must be valid when we use the basis functions as well. We have: X X 1 ) pˆ 1i,k ( xˆ) pˆ ( xˆ) = pˆ (ˆa0i ) pˆ 0i ( xˆ) + D p(ˆ ˆ a1i )(ξˆi,k i
+
X
i,k
D
2
2 ˆ2 , ξi,l ) pˆ 2i,kl ( xˆ). pˆ (ˆa2i )(ξˆi,k
i,k,l
1 ) = Dp(a1 )Bξˆ 1 by a simple application of the chain Now D pˆ (ˆa1i )(ξˆi,k i i,k rule and therefore 1 1 D p(ˆ ˆ a1i )(ξˆi,k ) = Dp(a1i )ξi,k ,
by (iv).
By a similar treatment of the second derivative term, we get X X 1 )pi,k (x) pˆ ( xˆ) = p(a0i )p0i (x) + Dp(a1i )(ξi1 )(ξi,k i
+
X
i,k
D2 p(a2i )(ξik2 , ξil2 )pikl (x) = p(x).
i,k,l
Thus the relations (iv) and (v) are compatible. ˆ P) ˆ Σ, ˆ be affine equivalent with F K Theorem 5.1. Let (K, Σ, P) and (K, as the affine transformation. If v : K → R induces vˆ : Kˆ → R by ˆ (x = F K ( xˆ)), then πv vˆ ( xˆ) = v(x) for xˆ ∈ K, b = πˆ vˆ . N , Σ = {ϕ } N . By definition, Proof. Let Σˆ = {ϕˆ i }i=1 i i=1
ϕˆ i (πv) b = ϕi (πv) = ϕi (v), 1 ≤ i ≤ N. ϕˆ i (ˆπvˆ ) = ϕˆ i (ˆv) = ϕi (v), 1 ≤ i ≤ N.
Thus, ϕˆ i (ˆπvˆ ) = ϕˆ i (πv) b for 1 ≤ i ≤ N. Hence πˆ vˆ = πv b by uniqueness of the function πˆ vˆ . 57 Let us consider a polygonal domain Ω with a triangulation th . Suppose to each K ∈ th is associated a finite element, (K, ΣK , PK ), ΣK being the set of degrees of freedom, and PK the finite dimensional space such that ΣK is PK unisolvent. Then we have defined the interpolation operator πK . All these make sense locally i.e. at a particular finite element K. We now define the global counterparts of these terms. The comparison is given in the following table.
5. General Properties of Finite Elements
58
Table 5.1. Local definition 1.
Finite element K.
1.
Global definition S The set Ω = K K=tn
2. 3.
4.
5. 6.
7.
58
The boundary of K, ∂K. The space PK of functions K → R, which is finitedimensional. P N of The set K = {ϕi,K }i=1 degrees of freedom of K.
2. 3.
Basis functions of PK are N {pi,K }i=1 . The nodes of K are {a0i , a1i , a2i , . . .}.
5.
πK is the PK interpolation operator, defined by ϕi,K (πK v) = ϕi,K (v), for all P ϕi,K ∈ K .
Notice that, by definition, (5.5)
(πh v)K = πK (vK)
4.
6.
The boundary of Ω, ∂Ω = Γ. The space Vh of functions Ω → R, which is also finitedimensional. The set of degrees of freeP N , where dom h = {ϕi }i=1 ϕi (pK) = ϕi,K (pK). Basis function of Vh are {w j }. The nodes of th are by defS {Nodes of K} = inition, K∈th
7.
∪{b j } say. The Vh interpolation operator πh is defined by πh v ∈ Vh such that, ϕi,K (πh vK) = P ϕi,K (vK) for all ϕi,K ∈ K . for all
K ∈ th .
It is this property and the conclusion of theorem 5.1 that will be essential in our future error analysis. Definition 5.5. We say that a finite element of a given type is of class C 0 , resp. of class C 1 , if, whenever it is the generic finite element of a triangulation, the associated space Vh satisfies the inclusion Vh ⊂ C 0 (Ω), resp. Vh ⊂ C 1 (Ω). By extension, a triangulation is of class C 0 , resp. of class C 1 if it is made up of finite elements of class C 0 , resp. of class C 1 . Reference: A forthcoming book of Ciarlet and Raviart [5].
Chapter 6
Interpolation Theory in Sobolev Spaces WE OUTLINED THE internal approximation method in Sec. 3. We are 59 naturally interested in the convergence of the solutions uh ∈ Vh to the global solution u ∈ V. As a key step in this analysis we obtained the error estimate (cf. Theorem 3.1): (6.1)
u − uh  ≤ C inf u − vh . vh ∈Vh
To be more specific let us consider an example. Given Ω ⊂ R2 a polygon, consider the solution of the following problem, which is therefore posed in the space V = H01 (Ω): −∆u + au = f in Ω, (6.2) u = 0 on Γ. Let th be a triangulation of Ω by triangles of type (1), (2) or (3). Then uh ∈ Vh ⊂ H01 (Ω) and (6.1) reads as (6.3)
u − uh 1,Ω ≤ C inf u − vh 1,Ω . vh ∈Vh
We know ‘a priori’ that u ∈ H01 (Ω). Let us assume for the moment that u ∈ C 0 (Ω). (Such assumptions are made possible by the various 59
60
6. Interpolation Theory in Sobolev Spaces
regularity theorems. For instance, u ∈ H 2 (Ω) ⊂ C 0 (Ω) if f ∈ L2 (Ω) and Ω is a convex polygon). If u ∈ C 0 (Ω), then we may define the Vh interpolate of u, i.e., πh u by πh u(b j ) = u(b j ) for the nodes b j of the triangulation. Note also that πh uK = πK u (cf. (5.5)). Now from (6.3) we get, u − uh 1,Ω ≤ Cu − πh u1,Ω 21 X 2 u − πh u1,K = C K∈th
12 X = C u − πK u21,K K∈th
60
Thus the problem of estimating u − uh 1,Ω is reduced to the problem of estimating u − πK u1,K . This is one central problem in the finite element method and motivates the study of interpolation theory in Sobolev spaces. We consider more general types of Sobolev spaces for they are no more complicated for this purpose than those defined in Sec. 2. Definition 6.1. Let m ≥ 0 be an integer, and 1 ≤ p ≤ +∞. Then the Sobolev space W m,p (Ω) for Ω ⊂ Rn , open, is defined by W m,p (Ω) = {v ∈ L p (Ω); ∂α v ∈ L p (Ω) for all
α ≤ m}.
Remark 6.1. H m (Ω) = W m,2 (Ω). On the space W m,p (Ω) we define a norm  · m,p,Ω by 1/p Z X α p ∂ v dx (6.4) vm,p,Ω = Ω α≤m
and the seminorm  · m,p,Ω by Z (6.5) vm,p,Ω =
X
Ω α=m
1/p ∂α v p dx
6. Interpolation Theory in Sobolev Spaces
61
If k ≥ 1 is an integer, consider the space W k+1,p (Ω)/Pk . If v˙ stands for the equivalence class of v ∈ W k+1,p (Ω) we may define the analogues of (6.4) and (6.5) respectively by (6.6)
˙vk+1,p,Ω = inf v + pk+1,p,Ω p∈Pk
and
61
˙vk+1,p,Ω = vk+1,p,Ω .
(6.7)
These are obviously welldefined and  · k+1,p,Ω defines the quotient norm on the quotient space above. We then have the following key result, whose proof may be found in Neˇcas [20] for instance. Theorem 6.1. In W k+1,p (Ω)/Pk , the seminorm ˙vk+1,p,Ω is a norm equivalent to the quotient norm ˙vk+1,p,Ω , i.e., there exists a constant C = C(Ω) such that for all v˙ ∈ W k+1,p (Ω)/Pk (6.8)
˙vk+1,p,Ω ≤ ˙vk+1,p,Ω ≤ C˙vk+1,p,Ω .
Equivalently, we may state Theorem 6.2. There exists a constant C = C(Ω) such that for each v ∈ W k+1,p (Ω) (6.9)
inf. v + pk+1,p,Ω ≤ Cvk+1,p,Ω .
p∈Pk
(Note: This result holds if Ω has a continuous boundary and if it is bounded so that Pk ⊂ W k+1,p (Ω).) We now prove the following Theorem 6.3. Let W k+1,p (Ω) and W m,q (Ω) be such that W k+1,p (Ω) ֒→ W m,q (Ω) (continuous injection). Let π ∈ L (W k+1 , ρ(Ω), W m,q (Ω)), i.e. a continuous linear map, such that for each p ∈ Pk , πp = p. Then there exists C = C(Ω) such that for each v ∈ W k+1,p (Ω) v − πvm,q,Ω ≤ CI − πL (W k+1 ,P(Ω),W m,q (Ω)) vk+1,p,Ω 62
6. Interpolation Theory in Sobolev Spaces
62
Proof. For each v ∈ W k+1,p (Ω) and for each p ∈ Pk , we can write v − πv = (I − π)(v + p). Thus, v − πvm,q,Ω ≤ v − πvm,q,Ω
= I − πL (W k+1,p (Ω),W m,q (Ω)) v + pk+1,p,Ω ,
for all p ∈ PK . Hence, v − πvm,q,Ω ≤ I − πL (W k+1,p (Ω),W m,q (Ω)) inf v + pk+1,p,Ω p∈Pk
≤ CI − πL (W k+1,p (Ω),W m,q (Ω)) vk+1,p,Ω By theorem 6.2, this completes the proof.
ˆ of Rn are said to be affine equivDefinition 6.2. Two open subsets Ω, Ω alent if there exists an invertible affine map F mapping xˆ to B xˆ + b, B an ˆ = Ω. invertible (n × n) matrix and b ∈ Rn , such that F(Ω) ˆ are affine equivalent, then we have a bijection between their If Ω, Ω points given by xˆ ↔ x = F( xˆ). Also we have bijections between smooth ˆ defined by (v : Ω → R) ↔ (ˆv : Ω ˆ → R) where functions on Ω and Ω v(x) = vˆ ( xˆ). The following theorem gives estimates of vm,p,Ω and ˆvm,p,Ωˆ each in terms of the other. 63
ˆ ⊂ Rn be affine equivalent. Then there exist Theorem 6.4. Let Ω, Ω constants C, Cˆ such that for all v ∈ W m,p (Ω) (6.11)
ˆvm,p,Ωˆ ≤ CBm  det B−1/p vm,p,Ω
ˆ and for all vˆ ∈ W m,p (Ω) (6.12) Note:
vm,p,Ω ≤ CB−1 m  det B1/p ˆvm,p,Ωˆ .
6. Interpolation Theory in Sobolev Spaces
63
(i) It suffices to prove either (6.11) or (6.12). We get the other by ˆ We will prove the merely interchanging the roles of Ω and Ω. former. (ii) B is the usual norm of the linear transformation defined by B, Bx ˆ = Ω, F( xˆ) = B xˆ + b). . (Recall that F(Ω) viz. B = sup n x x∈R x,0
Proof. Let {e1 , . . . , en } be the standard basis for Rn . Let α be a multiindex with α = m. By choosing a suitable collection {e1α , . . . , emα } with appropriate number of repetitions from the basis, we may write, (∂α vˆ )( xˆ) = (Dm vˆ )( xˆ)(e1α , . . . , emα ), where Dm vˆ is the mth order Fr´echet derivative of vˆ and Dm vˆ ( xˆ) is consequently an mlinear form on Rn . Thus, ∂α vˆ ( xˆ) ≤ Dm vˆ ( xˆ) = sup Dm vˆ ( xˆ)(ξ1 , . . . , ξm ). ξi =1 1≤i≤m
Since this is true for all α = m, we get ˆvm,p,Ωˆ ≤ C1
(6.13)
Z
ˆ Ω
m
p
D vˆ ( xˆ) d xˆ
!1/p
≤ C2 ˆvm,p,Ωˆ .
The first inequality is a consequence of our preceding argument. The second follows by a straightforward argument. By composition of functions in differentiation: (6.14)
Dm vˆ ( xˆ)(ξ1 , . . . , ξm ) = Dm v(x)(Bξ1 , . . . , Bξm);
This gives
64
Dm vˆ ( xˆ) ≤ Dm v(x) Bm .
(6.15)
Hence the first inequality in (6.13) may be rewritten as Z ˆv p ˆ ≤ C1p Bmp Dm v(F( xˆ )) p d xˆ m,p,Ω
ˆ Ω
6. Interpolation Theory in Sobolev Spaces
64 =
p C1 Bmp  det B−1
Z
Ω
Dm v(x) p dx
≤ C p Bmp  det B−1 vm,p,Ω . by an inequality similar to the second inequality of (6.13). Raising to power 1/p on either side we get (6.11). This completes the proof. We now estimate the norms B and B−1  in terms of the ‘sizes’ of ˆ More precisely, if h, (resp. h) ˆ the supremum of the diameters Ω and Ω. ˆ we have the following: of all balls that can be inscribed in Ω, (resp. Ω), Theorem 6.5. B ≤ h/ρ, ˆ
and
ˆ B−1  ≤ h/ρ.
Proof. Again it suffices to establish one of these. Now, ! 1 B = sup Bξ . ˆ ξ=ρˆ ρ ˆ such that ξ = yˆ − zˆ. Then Let ξ ∈ Rn with ξ = ρ. ˆ Choose yˆ , zˆ ∈ Ω Bξ = Bˆy − Bˆz = y − z, where F(ˆy) = y, F(ˆz) = z. But y, z ∈ Ω and hence y − z ≤ h. Thus Bξ ≤ h. Hence B ≤ h/ρ, ˆ which completes the proof. We conclude this section with an important, often used, result. ˆ P) ˆ Σ, ˆ be a finite element. Let s(= 0, 1 or 2) be Theorem 6.6. Let (K, ˆ Assume that: the maximal order of derivatives occurring in Σ. 65
ˆ ֒→ C s (K) ˆ (i) W k+1,p (K) ˆ ֒→ W m,q (K) ˆ (ii) W k+1,p (K) ˆ (iii) Pk ⊂ Pˆ ⊂ W m,q (K)
ˆ P) ˆ Σ, ˆ such that for all affine Then there exists a constant C = C(K, equivalent finite elements (K, Σ, P) we have (6.16)
1
1
v − πK vm,q,K ≤ C(meas K) q − p
hk+1 K vk+1,p,K ρm K
for all v ∈ W k+1,p (K), where hK is the diameter of K and ρK is the supremum of diameters of all balls inscribed in K.
6. Interpolation Theory in Sobolev Spaces
65
ˆ for any polynomial p ∈ Pk we have πˆ p = p. We Proof. Since Pk ⊂ P, may write X X 1 πˆ vˆ = vˆ (ˆa0i ) pˆ 0i + (Dˆv(ˆa1i )(ξˆi,k )) pˆ 1i,k
(6.17)
i
i,k
X 2 ˆ2 = + (D2 vˆ (ˆa2i )(ξˆi,k , ξi,l )) pˆ 2i,k,l ,
.
i,k,l
all these sums being finite (the second and third may or may not be ˆ W m,q (K)). ˆ present). We claim that πˆ ∈ L (W k+1 , p(K), Since Pˆ ⊂ m,q m,q ˆ all the basis functions in (6.17) are in W (K). ˆ Thus, W (K) X ˆπvˆ m,q,Kˆ ≤ ˆv(ˆa0i )  pˆ 0i m,q,Kˆ i
(6.18)
+
X i,k
+
X i,k,l
1 Dˆv(ˆa1i )(ξi,k )  pˆ 1ik m,q,Kˆ 2 ˆ2 D2 vˆ (ˆa2i )(ξˆi,k , ξi,l )  pˆ 2ikl m,q,Kˆ
ˆ ֒→ C s (K) ˆ and all the numbers vˆ (ˆa0 ), etc. . . , are Since W k+l,p (K) i ˆ bounded by their essential supremum over K, ˆπvˆ m,q,Kˆ ≤ Cˆvk+1,p,Kˆ . 66
Hence the claim is valid. Now by virtue of (ii) and also our observation on preservation of polynomials, we may apply theorem 6.3 to πˆ . ˆ P) ˆ Σ, ˆ such that Hence there exists C = C(K, ˆv − πˆ vˆ m,q,Kˆ ≤ Cˆvk+1,p,Kˆ
for
ˆ vˆ ∈ W k+1,p (K).
Notice that πˆ vˆ = πd ˆ − πˆ vˆ = v −d πK v. Thus K v by Theorem 5.1. Thus v ˆ = K where F K ( xˆ) = BK xˆ + bK , we get if F K (K) (6.19)
m 1/q v − πK vm,q,K ≤ C1 B−1 ˆv − πˆ vˆ m,q,Kˆ , K   det B K 
by Theorem 6.4. Also by the same theorem (6.20)
ˆvk+1,p,Kˆ ≤ C2 BK k+1  det BK −1/p vk+l,p,K .
66
6. Interpolation Theory in Sobolev Spaces
Further  det BK  being the Jacobian of the transformation, we have hK meas K ˆ and BK  ≤  det BK  = , B−1 K  ≤ h/ρK by theorem 6.5. ρˆ meas Kˆ ˆ ρ, Since h, ˆ meas Kˆ are constants, combining (6.19), (6.20) and the preceding observation we complete the proof of the theorem. References: See Bramble and Hilbert [28], Bramble and Zl´qmal [2], ˇ sek Ciarlet and Raviart [7], Ciarlet and Wagschal [8], Strang [21], Zeniˇ [23, 24] and Zl´amal [25, 32].
Chapter 7
Applications to SecondOrder Problems Over Polygonal Domains WE APPLY THE results of the preceding section in studying the conver 67 gence of the finite element method, i.e. the convergence of the solutions uh of (Ph ) to the solution u of a problem (P) which corresponds to the choice V = H 1 (Ω) or H01 (Ω), which we saw in Sec. 2 led to secondorder problems. Let Ω be a polygonal domain throughout. Definition 7.1. A family (th ) of triangulations of Ω is regular is (i) for all th and for each K ∈ th , the finite elements (K, Σ, P) are ˆ P) ˆ Σ, ˆ called the all affine equivalent to a single finite element, (K, reference finite element of the family; (ii) there exists a constant σ such that for all th and for each K ∈ th we have (7.1)
hK ≤σ ρK
where hK , ρK are as in Theorem 6.6; 67
7. Applications to SecondOrder Problems...
68
(iii) for a given triangulation th if h = max hK ,
(7.2)
K∈th
then h → 0. Remark 7.1. The condition (ii) in definition 7.1 assures us that as h → 0 the triangles do not become “flat”; cf. Exercise 7.1.
68
Exercise 7.1. If n = 2 and the sets K are triangles, show that condition (ii) of definition 7.1 is valid if and only if there exists θ0 > 0 such that for all th and K ∈ th , θK ≥ θ0 > 0, θK being the smallest angle in K. Exercise 7.2. Consider the space Vh associate with th . Since Vh is finite dimensional all norms are equivalent and hence vh 0,∞,Ω ≤ Ch vh 0,Ω
for all
vh ∈ Vh ,
for some constant Ch , a priori dependent upon h, which we may evaluate as follows: If (th ) is a regular family of triangulations, show that there exists a constant C, independent of h, such that (7.3)
vh 0,∞,Ω ≤
C vh 0,Ω hn/2
for
vh ∈ Vh .
Also show that there exists a constant C such that (7.4)
vh 1,Ω ≤
C vh 0,Ω h
for all
vh ∈ Vh .
We now obtain an estimate for the error u − uh 1,Ω when the family of triangulations is regular, which also gives convergence. Theorem 7.1. Let (th ) be a regular family of triangulations on Ω of ˆ P). ˆ Σ, ˆ We class C 0 (i.e. Vh ⊂ C 0 (Ω)) with reference finite element (K, assume that there exists an integer k ≥ 1 such that ˆ (i) PK ⊂ Pˆ ⊂ H 1 (K) ˆ ֒→ C s (K) ˆ where s(= 0, 1, or 2) is the maximal order of (ii) H k+1 (K) P derivatives in ˆ .
7. Applications to SecondOrder Problems...
69
(iii) u ∈ H k+1 (Ω) (Regularity assumption).
Then there exists a constant C (independent of Vh ) such that u − uh 1,Ω ≤ Chk uk+1,Ω .
(7.5)
Proof. Since Vh ⊂ C 0 (Ω), PK ⊂ H 1 (K), we have Vh ⊂ V. By (ii) and 69 (iii) of the hypothesis we have that the Vh interpolate of u, viz. πh u is welldefined. Since πh u ∈ Vh , by our fundamental result (see Theorem 3.1 or relation (6.1)), it suffices to estimate u − πh u1,Ω . Now, ˆ ֒→ C s (K), ˆ H k+1 (K) ˆ ֒→ H 1 (K) ˆ (k ≥ 1), H k+1 (K) ˆ Pk ⊂ P ⊂ H 1 (K), and we may apply Theorem 6.6 with p = q = 2, m = 1 to get
(7.6)
u − πK u1,K ≤ Cuk+1,K
hk+1 K ρK
≤ Cuk+1,K hkK
(Since
hK ≤ σ). ρK
Similarly with m = 0 we get (7.7)
u − πK u0,K ≤ Cuk+1,K hk+1 K .
These together give (7.8)
u − πK u1,K ≤ ChkK uk+1,K .
Now since hK ≤ h, u − πh u1,Ω
21 X 2 u − πK u1,K = K∈th
21 X u2k+1,K ≤ C hk K∈th
K
= C h uk+1,Ω . This completes the proof.
7. Applications to SecondOrder Problems...
70
Example 7.1. Consider triangulations by triangles of type (1). Then 70 ˆ ֒→ C 0 (K). ˆ If u ∈ H 2 (Ω), k = 1, Pˆ = P1 and if n = 2 or 3, H 2 (K) Theorem 7.1 says that u − uh 1,Ω ≤ C hu2,Ω . We conclude the analysis of convergence in the norm  · 1,Ω with the following result. Theorem 7.2. Let (th ) be a regular family of triangulations of Ω, of ˆ Then (with the class C 0 . Let s = 0 or 1 and let P1 ⊂ Pˆ ⊂ H 1 (K). 1 1 assumption that u ∈ V = H (Ω) or H0 (Ω)) we have lim u − uh 1,Ω = 0
(7.9)
h→0
Proof. Let V = V ∩ W 2,∞ (Ω). Since s ≤ 1, W 2,∞ (·) ֒→ C s (·) and W 2,∞ (·) ֒→ H 1 (·). The second inclusion follows ‘a fortiori’ from the ˆ Thus we may apply Theorem first with s = 1. Also, P1 ⊂ Pˆ ⊂ H 1 (K). 6.6 with k = 1, p = ∞, m = 1, q = 2. Then for all v ∈ V , 1
v − πK v1,K ≤ C(meas K) 2 hv2,∞,K 1
≤ C(meas K) 2 hv2,∞,Ω . Summing over K, we get
v − πh v1,Ω
12 X meas K ≤ C hv2,∞,Ω K∈th
= C hv2,∞,Ω ,
since
P
(7.10)
K∈th
meas K = meas Ω, a constant. Thus, for all v ∈ V , lim v − πh v1,Ω = 0
h→0
71
Notice that V = V. Hence choose v0 ∈ V such that u−v0 1,Ω ≤ ǫ/2 where ǫ > 0 is any preassigned quantity. Then once v0 is chosen, by
7. Applications to SecondOrder Problems...
71
(7.10) choose h0 such that for all h ≤ h0 , v0 − πh v0 1,Ω ≤ ǫ/2. Now, by (6.1) u − uh 1,Ω ≤ Cu − πh v0 1,Ω
≤ C u − v0 1,Ω + v0 − πh v0 1,Ω
≤ Cǫ, for h ≤ h0 .
This gives (7.9) and completes the proof.
We now have, by Theorem 7.1, u − uh 0,Ω ≤ u − uh 1,Ω = 0(hk ). We now show, by another argument that u − uh 0,Ω = 0(hk+1 ), (at least in some cases) there by giving a more rapid convergence than expected. This is done by the AubinNitsche argument (also known as the duality argument). We describe this in an abstract setting. Let V be a normed space with norm denoted by  · . Let H be a Hilbert space with norm  ·  and inner product (·, ·) such that (7.11)
(i) V ֒→ H, (ii) V = H.
and
For secondorder problems: V = H 1 (Ω) or H01 (Ω) and H = L2 (Ω). Since H is a Hilbert space, we may identify it with its dual. Further since V is dense in H, we have that H may be identified with a subspace of V ′ , the dual of V. For, if g ∈ H, define g˜ ∈ V ′ by g˜ (v) = (g, v) · g˜ ∈ V ′ since ˜g(v) ≤ Cg v. If g˜ (v) = 0 for all v ∈ V, then (g, v) = 0 for all v ∈ H as well since V = H. Thus g = 0. This proves the identification. In the sequel we will set g = g˜ . 72 Recall that u and uh are the solutions of the problems: for all v ∈ V,
(P)
a(u, v) = f (v)
(Ph )
a(uh , vh ) = f (vh ) for all vh ∈ Vh ⊂ V,
and that the assumptions on (P) are as in the LaxMilgram lemma. Then we have the following theorem.
7. Applications to SecondOrder Problems...
72
Theorem 7.3. Let the spaces H and V satisfy (7.11). Then with our above mentioned notations, ) ( 1 inf ϕ − ϕh  , (7.12) u − uh  ≤ Mu − uh  sup g∈H g ϕh ∈Vh
where for each g ∈ H, ϕ ∈ V is the corresponding unique solution of the problem (P∗ )
a(v, ϕ) = (g, v)
for all v ∈ V,
and M the constant occurring in the inequality giving continuity of a(·, ·). Remark 7.2. Note that unlike in (P), we solve for the second argument of a(·, ·) in (P∗ ). This is called the adjoint problem of (P). The existence and uniqueness of the solution of (P∗ ) are proved in an identical manner. Note that if a(·, ·) is symmetric, then (P) is selfadjoint in the sense that (P) = (P∗ ). Proof. From the elementary theory of Hilbert spaces, we have u − uh  = sup
(7.13)
g∈H g,0
(g, u − uh ) . g
For a given g ∈ H, (7.14) 73
(g, u − uh ) = a(u − uh , ϕ)
Also if ϕh ∈ Vh we have, (7.15)
a(u − uh , ϕh ) = 0.
Thus (7.14) and (7.15) give (7.16)
(g, u − uh ) = a(u − uh , ϕ − ϕh ),
which gives us (g, u − uh ) ≤ Mu − uh  ϕ − ϕh ,
7. Applications to SecondOrder Problems... and hence
73
! ϕ − ϕh  u − uh  ≤ Mu − uh  sup . g g∈H g,0
by (7.13). Since this is true for any ϕh ∈ Vh we may take infimum over Vh to get (7.12), which completes the proof. For dimensions ≤ 3 and Lagrange finite elements we now show that u − uh 0,Ω = 0(hk+1 ). For this we need one more definition. Definition 7.2. Let V = H 1 (Ω) or H01 (Ω), H = L2 (Ω). The adjoint problem is said to be regular if the following hold: (i) for all g ∈ L2 (Ω), the solution ϕ of the adjoint problem for g belongs to H 2 (Ω) ∩ V; (ii) there exists a constant C such that for all g ∈ L2 (Ω) ϕ2,Ω ≤ Cg0,Ω ,
(7.17)
where ϕ is the solution of the adjoint problem for g. Theorem 7.4. Let (th ) be a regular family of triangulations on Ω with ˆ P). ˆ Σ, ˆ Let s = 0 and n ≤ 3. Suppose there reference finite element (K, ˆ Assume 74 exists an integer k ≥ 1 such that u ∈ H k+1 (Ω), Pk ⊂ Pˆ ⊂ H 1 (K). further that the adjoint problem is regular in the sense of Definition 7.2. Then there exists a constant C independent of h such that (7.18)
u − uh 0,Ω ≤ C hk+1 uk+1,Ω .
Proof. Since n ≤ 3, H 2 (·) ֒→ C 0 (·). Also, H 2 (·) ֒→ H 1 (·) and P1 ⊂ Pˆ ⊂ ˆ Thus for ϕ ∈ H 2 (Ω), by Theorem 7.1, H 1 (H). ϕ − πh ϕ1,Ω ≤ C hϕ2,Ω . Hence (7.19)
inf ϕ − ϕh 1,Ω ≤ C hϕ2,Ω .
ϕh ∈Vh
7. Applications to SecondOrder Problems...
74 By (7.12) and (7.19).
u − uh 0,Ω ≤ Mu − uh 1,Ω sup
g∈L2 (Ω)
! 1 Chϕ2,Ω . g0,Ω
By the regularity of (P∗ ), ϕ2,Ω ϕ2,Ω ≤ ≤ constant. g0,Ω g0,Ω
(7.20)
Thus, u − uh 0,Ω ≤ C hu − uh 1,Ω ≤ C h(hk uk+1,Ω )
(by theorem 7.1).
This gives (7.18) and completes the proof.
We finally give an estimate for the error in the L∞ norm. Theorem 7.5. Let (th ) be a regular family of triangulations on Ω ⊂ Rn , where n ≤ 3. Assume further that for all th and K ∈ th . (7.21) 75
0<τ≤
hK ≤, f rm[o]−−, h
τ being a constant.
ˆ ∩ L∞ (K). ˆ If (P∗ ) is regular, Let u ∈ H 2 (Ω) and P1 ⊂ Pˆ ⊂ H 1 (K) then there exists a constant C independent of h such that u − uh 0,∞,Ω ≤ C hu2,Ω ; if n = 2 √ (7.22) u − uh 0,∞,Ω ≤ C hu2,Ω if n = 3.
Proof. Assume n = 2. Now (7.23)
u − uh 0,∞,Ω ≤ u − πh u0,∞,Ω + πh u − uh 0,∞,Ω .
Note that since (uh − πh u) ∈ Vh , we may apply Exercise 7.1 to get (7.24)
uh − πh u0,∞,Ω ≤
C uh − πh u0,Ω . h
Thus, uh − πh u0,∞,Ω ≤
C uh − u0,Ω + u − πh u0,Ω h
7. Applications to SecondOrder Problems...
75
i Ch C1 h2 u2,Ω + C2 h2 u2,Ω h ≤ C hu2,Ω (by Theorem 7.4 and Theorem 6.6).
≤
ˆ Thus, Also H 2 (·) ֒→ C 0 (·); H 2 (·) ֒→ L∞ (·) and P1 ⊂ Pˆ ⊂ L∞ (K). by Theorem 6.6 with k = 1, p = 2, m = 0, q = ∞, 1
u − πK u0,∞,K ≤ C(meas K)− 2 h2 u2,K . Since n = 2, meas K ≥ Cρ2K ≥
Cτ2 2 C 2 h h ≥ σ2 K σ2
1
by (7.1), so that (meas K)− 2 ≤ C h−1 and therefore, u − πK u0,∞,K ≤ C hu2,K . Hence we obtain (7.22) for n = 2 since u − πh u0,∞,Ω = max u − πK u0,∞,K ≤ C hu2,Ω . K∈th
76
For n = 3, the only variation in the proof occurs in the fact that uh − πh u0,∞,Ω ≤
C uh − πh u0,Ω h3/2
as in Exercise 7.1 and that now meas K ≥ Cρ3K ≥ This completes the proof.
C Cτ3 3 3 · h ≥ ·h . K σ3 σ3
References: One may refer to Ciarlet and Raviart [6] for 0(h) convergence in the norm  · 0,∞,Ω for any n. See also Bramble and Thom´ee [1].
Chapter 8
Numerical Integration LET US START with a specific problem. Let Ω be a polygonal domain 77 in Rn . Consider the problem
(8.1)
n P ∂ ∂u − ∂xi ai j ∂x j = f in Ω, i, j=1 u = 0 on Γ = ∂Ω.
where the (ai j ) and f are functions over Ω which are smooth enough. Let us further assume that there exists α > 0 such that, for all ξ ∈ Rn , (8.2)
n X
i, j=1
ai j ξi ξ j ≥ α
n X
ξi2 .
i=1
It has been seen earlier (cf. Remark 2.2) that the above problem (8.1) is obtained from a problem (P) with a(·, ·) and f (·) being defined by (8.3)
R Pn ∂u ∂v a(u, v) = Ω i, j=1 ai j ∂x j ∂xi dx R f (v) = f v dx Ω
for u, v ∈ V = H01 (Ω). Approximating the solution by the finite element method, i.e. by constructing a regular family of triangulations (th ) with reference finite 77
8. Numerical Integration
78
ˆ P), ˆ Σ, ˆ we get the problems (Ph ), i.e., to find uh ∈ Vh such element (K, that (8.4)
for all vh ∈ Vh .
a(uh , vh ) = f (vh ),
M for V , then we may write If we choose a basis {wk }k=1 h
(8.5)
uh =
M X
uk wk .
k=1
Thus to solve (Ph ) we have to solve the linear system (8.6)
M X k=1
a(wk , wm )uk = f (wm ), 1 ≤ m ≤ M.
78
Notice that a(wk , wm ) = (8.7) f (wm ) =
n XZ X
K∈th
K i, j=1
K∈th
K
XZ
ai j
∂wk ∂wm dx ∂x j ∂xi
f wm dx.
Thus we have ended up with the computations of integrals over K ∈ th . These are, in general, difficult or impossible to evaluate exactly and one thus has to resort to numerical methods. We now study briefly how this may be done. ˆ = K, where, F K ( xˆ) = BK xˆ +bK , with det BK > Let us assume F K (K) 0. There is no loss in generality in the last assumption. Then if ϕ is a function over K, we have Z Z ϕ( ˆ xˆ)d xˆ (8.8) ϕ(x)dx = (det BK ) Kˆ
K
the functions ϕ and ϕˆ being in the usual correspondence. We then replace the expression in the righthand side by the following: (8.9)
Z
Kˆ
b ∼ ϕ( ˆ xˆ)dx
L X 1=1
ω ˆ 1 ϕ( ˆ bˆ 1 ).
8. Numerical Integration
79
In this section ∼ will denote the righthand side replacing the expression in the left hand side in similar relations). In (8.9) the quantities ω ˆ 1 are called the weights and the points bˆ 1 are called the nodes of the quadrature scheme. While in general we may assume ω ˆ1 ∈ R and bˆ 1 ∈ Rn , we will restrict ourselves to the most common case where ˆ 1 ≤ l ≤ L. ω ˆ 1 > 0 and bˆ 1 ∈ K, We may now define the error functional Eˆ by (8.10)
Eˆ (ϕ) ˆ =
Z
Kˆ
ϕ( ˆ xˆ)d xˆ −
L X
ω ˆ 1 ϕ( ˆ bˆ 1 ).
1=1
79
We will be interested in finding spaces of polynomials for which Eˆ (ϕ) ˆ = 0, i.e., again we need “polynomial invariance”, an idea already found in interpolation theory. The above quadrature scheme for Kˆ induces one on K as well since if we set ˆ 1, ω1,K = (det BK )ω (8.11) b1,k = F K (bˆ 1 ),
we then deduce the numerical quadrature scheme (8.12)
Z
K
ϕ(x)dx ∼
L X
ω1,K ϕ(b1,K ).
1=1
We shall therefore define the error functional Z L X ϕ(x)dx − (8.13) EK (ϕ) = ω1,K ϕ(b1,K ). K
l=1
Notice that the following relation holds: (8.14)
EK (ϕ) = (det BK )Eˆ (ϕ). ˆ
Example 8.1. Consider Kˆ to be an nsimplex in Rn . Let aˆ be its centroid. Let Z ˆ ϕ(ˆ ϕ( ˆ xˆ)d xˆ ∼ (meas K) ˆ a). Kˆ
8. Numerical Integration
80
Exercise 8.1. Show that E (ϕ) ˆ = 0 for ϕˆ ∈ P1 in Example 8.1. Example 8.2. Let Kˆ be a triangle in R2 . With the usual notations, set Z X 1 ˆ ϕ(ˆ ˆ ai j ). ϕ( ˆ xˆ)d xˆ ∼ (meas K) 3 Kˆ 1≤i< j≤3 80
Exercise 8.2. Show that Eˆ (ϕ) ˆ = 0 for ϕˆ ∈ P2 in Example 8.2. Example 8.3. Let Kˆ be as in Example 8.2. Let (cf. Fig. 8.1): Z 3 X X 1 . ˆ 3 ϕ(ˆ ˆ a ) + 8 ϕ(ˆ ˆ a ) + 27 ϕ(ˆ ˆ a ) ϕ( ˆ xˆ)d xˆ ∼ (meas K) i i j 60 ˆ K i=1 1≤i< j≤3
Figure 8.1: Exercise 8.3. Show that Eˆ (ϕ) ˆ = 0 for ϕˆ ∈ P3 , in Example 8.3. Let us now review the whole situation. We had the “original” approximation problem (Ph ): To find uh ∈ Vh such that a(uh , vh ) = f (vh ) for all vh ∈ Vh . This led to the solution of the linear system (8.6). By virtue of the quadrature scheme we arrive at a solution of a “modified” approximation problem (P∗h ): To solve the linear system (8.15)
M X k=1
ah (wk , wm )u∗k = fh (wm ), 1 ≤ m ≤ M,
8. Numerical Integration where
(8.16)
81
! n P PL P ∂wk ∂wm ah (wk , wm ) = ai j ∂x j ∂xi (b1,K ) l=1 ω1,K K∈th i, j=1 ! L P P ω1,K ( f wm )(b1,K ) . fh (wm ) = K∈th l=1
While uh was given by (8.5) we now obtain u∗h =
(8.17)
M X
u∗k wk .
k=1
Thus the problem (P∗h ) (not to be confused with any adjoint problem!) consists in finding u∗h ∈ Vh such that, for all wh ∈ Vh , 81 ah (u∗h , wh ) = fh (wh )
(8.18) where
(8.19)
! L n P P P ∂vh ∂wh ah (vh , wh ) = ai j ∂x j ∂xi (b1,K ) ω1,K K∈th l=1 i, j=1 L P P ω1,K ( f vh )(b1,K ), fh (vh ) = K∈th l=1
for vh , wh ∈ Vh .
Remark 8.1. The bilinear form ah (·, ·) : Vh × Vh → R and the linear form fh : Vh → R are not defined over V in general. For instance if V = H01 (Ω)(n = 2) in one of the examples, then as they require point values of the nodes, we see that they are not in general defined over V. Having obtained the approximate solution u∗h by numerical integration, we are naturally interested in its efficacy. Thus we require to know the error u−u∗h . We now carry out the error analysis, first in an abstract setting. Let us maintain our assumptions as in the LaxMilgram lemma and consider the problem (P). Then we have problems (P∗h ) to find u∗h ∈ Vh ⊂ V such that for all vh ∈ Vh , ah (u∗h , vh ) = fh (vh ) where fh ∈ Vh′ and ah (·, ·) is a bilinear form on Vh . Then we would like to answer the following questions:
8. Numerical Integration
82
(i) What are sufficient conditions such that (P∗h ) have unique solutions? (ii) Can we find an abstract error estimate for u − u∗h ? (iii) If u − uh  = 0(hk ), i.e., without numerical quadrature, under what conditions is this order of convergence preserved, i.e. when can we say u − u∗h  = 0(hk )? 82
The assumption of Vh ellipticity of the bilinear forms ah (·, ·) answers the first question (by the LaxMilgram lemma) and we will see in Theorem 8.2 under which assumptions it is valid. The following theorem answers the second question. Theorem 8.1. Let the bilinear forms ah (·, ·) be Vh elliptic uniformly with respect to h, i.e., there exists a constant e α > 0, independent of h, such that for all h and for all vh ∈ Vh , ah (vh , vh ) ≥ e αvh 2 .
(8.20)
Then the approximate problems (P∗h ) all have unique solutions u∗h , and further we have the estimate: u − u∗h  ≤
(8.21) ≤ C inf
vh ∈Vh
(
u − vh  + sup
wh ∈Vh
) !  f (wh ) − fh (wh ) a(vh , wh ) − ah (vh , wh ) + sup . wh  wh  wh ∈Vh
Remark 8.2. If a = ah , f = fh then we get our original estimate (3.1). Thus (8.21) generalizes our previous result. Remark 8.3. The terms involving a, ah and f , fh merely mean that if u∗h is to converge to u, then ah and fh must be “close to” a and f respectively. Their convergence to 0 with h may be viewed as “consistency conditions” which are so often found in Numerical Analysis. Proof. The existence and uniqueness of the u∗h are obvious by the LaxMilgram lemma applied to the Vh . Since, for all vh ∈ Vh , we have ah (u∗h , u∗h − vh ) = fh (u∗h − vh ),
8. Numerical Integration
83
a(u, u∗h − vh ) = f (u∗h − vh ), 83
we have the identity
(8.22)
ah (u∗h − vh , u∗h − vh ) = a(u − vh , u∗h − vh ) + {a(vh , u∗h − vh )
− ah (vh , u∗h − vh )} + { fh (u∗h − vh ) − f (u∗h − vh )}.
Hence by (8.20) we get e αu∗h − vh 2 ≤ Mu − vh  u∗h − vh 
+ a(vh , u∗h − vh ) − ah (vh , u∗h − vh )
+  fh (u∗h − vh ) − f (u∗h − vh ).
Thus, e αu∗h
a(vh , u∗h − vh ) − ah (vh , u∗h − vh ) − vh  ≤ Mu − vh  + u∗h − vh   fh (u∗h − vh ) − f (u∗h − vh ) + u∗h − vh  a(vh , wh ) − ah (vh , wh ) ≤ Mu − vh  + sup wh  wh ∈Vh  f (wh ) − fh (wh ) + sup wh  wh ∈Vh
since (u∗h − vh ) ∈ Vh . Hence, u − u∗h  ≤ u − vh  + u∗h − vh  M 1 a(vh , wh ) − ah (vh , wh ) ≤ 1+ u − vh  + sup e α e α wh ∈Vh wh   f (wh ) − fh (wh ) 1 . + sup e α wh ∈Vh wh 
Varying vh over Vh and taking the infimum, and replacing (1+ M/e α), 84 (1/e α) by a larger constant C, we get (8.21), which completes the proof.
8. Numerical Integration
84
The following theorem tells us when the uniform Vh ellipticity assumption of Theorem 8.1 is satisfied in the example we started with. Theorem 8.2. Let ah (·, ·) and fh (·) be as in (8.19). Assume further that L S (i) ω ˆ 1 > 0, 1 ≤ l ≤ L, (ii) Pˆ ⊂ Pk , and (iii) {bˆ l } contains a Pk′ −1 l=1
unisolvent subset. Then the ah (·, ·) are Vh elliptic uniformly with respect to h.
Proof. We must produce an α˜ > 0, free of h, such that ah (vh , vh ) ≥ αv ˜ h 21,Ω for all vh ∈ Vh . We have n L X XX ∂v ∂v h h (b ) ai j ah (vh , vh ) = ω1,K l,K ∂x ∂x j i i, j=1 K∈th l=1 (8.23) n !2 L XX X ∂pK ≥α ωl,K (bl.K ) ∂xi K∈t l=1 i=1 h
where pK = vh K . The inequality (8.23) is a result of the ellipticity condition (8.2) on the matrix (ai j ) and the fact that ωl,K > 0, since ω e1 > 0 and we assumed without loss in generality that det BK > 0. Now let pˆ K (b x) = pK (x), where x = BK xˆ + bK . Let BK = (bi j ), so that xj =
n X
b jl xˆl + bK, j .
l=1
Then
∂ pˆ K X ∂pK (x) ∂x j X ∂pK (x) = = b ji . ∂ xˆi ∂x j ∂ xˆi ∂x j j=1 j=1 n
Thus is ∂ pˆ K ( xˆ) ∂ pˆ K ( xˆ) ,..., Dˆ = ∂ xˆ1 ∂ xˆn
n
!
and
D=
! ∂pK (x) ∂pK (x) ,..., , ∂x1 ∂xn
ˆ 2 ≤ D2 BK 2 . Thus, we have Dˆ = DBK . Hence D !2 !2 n n X X ∂pK ∂ pˆ K ˆ −2 (8.24) (bl,K ) ≥ BK  (b1 ) . ∂xi ∂ xˆi i=1 i=1
8. Numerical Integration
85
85
Now suppose pˆ ∈ Pˆ is such that L X
(8.25)
l=1
!2 n X ∂ pˆ ˆ (b1 ) = 0. ω ˆ1 ∂ xˆi i=1
∂ pˆ ˆ (bl ) = 0 for all 1 ≤ l ≤ L and 1 ≤ i ≤ ∂ xˆi ∂ pˆ ∂ pˆ ∈ Pk′ −1 and hence = 0 by the Pk′ −1 n. Since Pˆ ⊂ Pk′ , we have ∂ xˆi ∂ xˆi ˆ 0 unisolvency. Thus pˆ ∈ P0 , and on the finite dimensional space P/P ˆ we have a norm defined by the (in practice we always have P0 ⊂ P) squareroot of the left hand side of (8.25). By the finite dimensionality this is equivalent to the norm defined by the  · 1,Kˆ norm on P. Hence we have a constant βˆ > 0 such that Then since ω ˆ 1 > 0, we have
(8.26)
L X l=1
!2 n X ∂ pˆ ˆ ˆ p (bl ) ≥ β ˆ 2l,Kˆ . ω ˆl ∂ x ˆ i i=1
We will apply this to pˆ K . We also have (8.27)
−2 −1 2  pˆ K 2l,K ≥ CB−1 K  (det B K ) pK l,K ,
by Theorem 6.4. Combining the inequalities (8.23), (8.24), (8.26) and (8.27), we get X −1 −2 −1 2 ˆ ah (vh , vh ) ≥ α (det BK )BK −2 βCB K  (det B K ) pK l,K K∈th
ˆ = αβC
X
K∈th
ˆ ≥ αβCγ
−2 2 (BK  B−1 K ) pK l,K
X
K∈th
pK 2l,K
2 ˆ = αβCγv ˜ h 2l,Ω h 1,Ω = αv −2 ≥ γ by Theorem 6.5. This proves the theorem. since (BK  B−1 K )
86
8. Numerical Integration
Let us now review our Examples 8.1 through 8.3 to see if the condi 86 tions of Theorem 8.2 are satisfied. Let n = 2 and consider Example 8.1. Clearly ω ˆ = meas Kˆ > 0. Also P ˆ P = P1 for triangles of type (1). Since K = {p(ˆa)} is P0 unisolvent, we have that for triangles of type (1) and the quadrature scheme of Example 8.1 the corresponding ah (·, ·) are Vh elliptic uniformly with respect to h. For triangles of type (2), Pˆ = P2 . The weights ω ˆ 1 are all > 0 in Example 8.2. Further we saw in Exercise 5.1 that {ˆai j }1≤i< j≤3 is P1 unisolvent. Hence Theorem 8.2 is valid for this quadrature scheme as well. For triangles of type (3), consider the quadrature scheme of Example 8.3. We have ω ˆ 1 > 0 and Pˆ = P3 . It was seen in Example 4.2 that the set {ai ; 1 ≤ i ≤ 3} ∪ {ai j ; 1 ≤ i < j ≤ 3} is P2 unisolvent. Hence the corresponding bilinear forms ah (·, ·) are Vh elliptic uniformly with respect to h. Exercise 8.4. Let (H,  · ) be a Hilbert space and V a subspace with norm  ·  such that V ֒→ H and V = H cf. Sec. 7. Then with the usual notations show that n1 inf (Mu − u∗h  ϕ − ϕh  + a(u∗h , ϕh ) − ah (u∗h , ϕh ) u − u∗h  ≤ sup g∈H g ϕh ∈Vh o + f (ϕh ) − fh (ϕh ))
where ϕ is the solution of adjoint problem for g.
We now turn our attention to the evaluation of the bound for u − u∗h  given by (8.21). For secondorder problems, for which the norm is  · 1,Ω , we will take as usual for vh ∈ Vh the element πh u ∈ Vh so that we now get the bound " a(πh u, wh ) − ah (πh u, wh ) ∗ u − uh 1,Ω ≤ C u − πh u1,Ω + sup wh 1,Ω wh ∈Vh # (8.28)  f (wh ) − fh (wh ) + sup . wh 1,Ω wh ∈Vh 87
Let us assume that we may apply Theorem 7.1, so that
8. Numerical Integration (8.29)
87
u − πh u1,Ω ≤ C hk uk+1,Ω .
In order to keep the same accuracy, we will therefore try to obtain estimates of the following form:
(8.30)
a(πh u, wh ) − ah (πh u, wh ) ≤ C(u)hk , sup wh 1,Ω wh ∈Vh  f (wh ) − fh (wh ) ≤ C( f )hk , wsup w  h 1,Ω ∈V h h
and these will in turn be obtained from “local” estimates. (cf. Theorem 8.4 and Exercise 8.5). As a preliminary step, we need two results which we prove now. The first of these is a historically important result in the interpolation theory in Sobolev spaces. Theorem 8.3 (BRAMBLEHILBERT LEMMA; cf. Bramble and Hilbert [27]). Let Ω ⊂ Rn be open with Lipschitz continuous boundary Γ. Let f ∈ W k+1,p (Ω)) which vanishes over Pk . Then there exists a constant C = C(Ω) such that, for all v ∈ W k+1,p (Ω), (8.31)
 f (v) ≤ C f ∗k+1,p,Ω vk+1,p,Ω .
Proof. For v ∈ W k+1,p (Ω) and all p ∈ Pk , we have f (v) = f (v + p), so that  f (v) =  f (v + p) ≤  f ∗k+1,p,Ω v + pk+1,p,Ω , and thus,
88
 f (v) ≤  f ∗k+1,p,Ω inf v + pk+1,p,Ω ≤
p∈Pk ∗ C f k+1,p,Ω vk+1,p,Ω ,
by Theorem 6.2, which completes the proof.
8. Numerical Integration
88
Lemma 8.1. Let ϕ ∈ W m,q (Ω), w ∈ W m,∞ (Ω). Then ϕw ∈ W m,q (Ω), and there exists a numerical constant C, independent of ϕ and w such that (8.32)
ϕwm,q,Ω ≤ C
M X j=0
ϕm− j,q,Ω w j,∞,Ω .
Proof. The result is an immediate consequence of the Leibniz formula: For any α = m, ∂α (ϕw) =
m X X
Cα,β ∂α−β (ϕ)∂β (w)
j=0 β= j
which yields (8.32).
We may now apply Lemma 8.1 and Theorem 8.3 to get the estimates (8.30). We do this in two stages (Theorem 8.4 and Exercise 8.5) in which, for the sake of simplicity, we present our results for the special case PK = P2 . Theorem 8.4. Let PK = P2 and consider a quadrature scheme such that ˆ ϕ) for all ϕˆ ∈ P2 , ξ( ˆ = 0. Then there exists a constant C, independent of K, such that for all ai j ∈ W 2,∞ (K) and for all p, p′ ∈ PK we have # " ∂p′ ∂p ∂p ∂p′ 1,K  0,K .  ≤ C h2K ai j 2,∞,K  (8.33) EK (ai j ) ∂x j ∂xi ∂x j ∂xi ∂p ∂p′ , ∈ P1 , it suffices to find a suitable esti∂x j ∂xi mate for EK (avw), for a ∈ W 2,∞ (K), v, w ∈ P1 . Further, since Proof. Since we have
(8.34) 89
EK (avw) = (det BK )Eˆ (ˆavˆ w), ˆ
ˆ and vˆ , we will first find an estimate for Eˆ (ˆavˆ w), ˆ with aˆ ∈ W 2,∞ (K) wˆ ∈ P1 . Let πˆ 0 wˆ be the orthogonal projection of wˆ onto the subspace P0 ˆ Then we may write in the sense of L2 (K). (8.35)
Eˆ (ˆavˆ w) ˆ = Eˆ (ˆavˆ πˆ 0 w) ˆ + Eˆ (ˆavˆ (wˆ − πˆ 0 w)). ˆ
8. Numerical Integration
89
(i) Estimate for Eˆ (ˆavˆ πˆ 0 w). ˆ ˆ → R defined by Consider the functional Eˆ : W 2,∞ (K) ˆ = ψˆ 7→ Eˆ (ψ)
Z
Kˆ
ˆ xˆ)d xˆ = ψ(
L X
ˆ bˆ 1 ). ω ˆ 1 ψ(
l=1
ˆ ψ ˆ ψ ˆ ≤ C ˆ 0,∞,Kˆ ≤ C ˆ 2,∞,Kˆ . Thus Eˆ is a continuous linear funcEˆ (ψ) 2,∞ ˆ Hence by Theorem 8.3, since Eˆ vanishes on P1 (⊂ tional on W (K). 1 P2 ) , we have a constant Cˆ such that (8.36)
ˆ ψ ˆ ≤ C ˆ 2,∞,Kˆ . Eˆ (ψ)
Thus ˆ avˆ πˆ 0 w Eˆ (ˆavˆ πˆ 0 w) ˆ ≤ Cˆ ˆ 2,∞,Kˆ ˆ avˆ 2,∞,Kˆ ˆπ0 w ≤ Cˆ ˆ 0,∞,Kˆ since πˆ 0 wˆ ∈ P0 is a constant function. By Lemma 8.1 (recall that vˆ ∈ P1 ), h i ˆ π0 w Eˆ (ˆavˆ πˆ 0 w) ˆ ≤ Cˆ ˆ 0,∞,Kˆ ˆa1,∞,Kˆ ˆv1,∞,Kˆ + ˆa2,∞,Kˆ ˆv0,∞,Kˆ . By the equivalence of the L2 and L∞ norms on P0 , and since the projection has norm less than that of the vector itself in any Hilbert space we have the chain of inequalities ˆ π0 w ˆ w ˆπ0 w ˆ 0,∞,Kˆ ≤ Cˆ ˆ 0,Kˆ ≤ C ˆ 0,Kˆ . 90
Similarly we may replace ˆv1,∞,Kˆ by ˆv1,Kˆ and ˆv0,∞,Kˆ by ˆv0,Kˆ since the L2 and L∞ norms are equivalent on P1 . Thus we get (8.37)
ˆ a1,∞,Kˆ ˆv1,Kˆ + ˆa2,∞,Kˆ ˆv0,Kˆ )w Eˆ (ˆavˆ πˆ 0 w) ˆ ≤ C(ˆ ˆ 0,Kˆ .
1 In this estimate, we do not use the “full” polynomial invariance of the quadrature scheme.
8. Numerical Integration
90 (ii) Estimate for Eˆ (ˆavˆ (wˆ − πˆ 0 w)). ˆ
ˆ Then Let wˆ ∈ P1 be fixed and let ϕˆ ∈ W 2,∞ (K). ˆ ϕ( Eˆ (ϕ( ˆ wˆ − πˆ 0 w)) ˆ ≤ C ˆ wˆ − πˆ 0 w) ˆ 0,∞,Kˆ ˆ ≤ Cϕ ˆ 0,∞,Kˆ wˆ − πˆ 0 w ˆ 0,∞,Kˆ
ˆ wˆ − πˆ 0 w ≤ C ˆ 0,∞,Kˆ ϕ ˆ 2,∞,Kˆ .
ˆ defined by ϕˆ 7→ Eˆ (ϕ( Thus the functional on W 2,∞ (K) ˆ wˆ − πˆ 0 w)) ˆ is ˆ w− continuous, linear with norm ≤ C ˆ πˆ 0 w ˆ 0,∞,Kˆ . Since for ϕˆ ∈ P1 , ϕ( ˆ wˆ − πˆ 0 w) ˆ ∈ P2 , we have that the functional vanishes on P1 . By Theorem 8.3, ˆ wˆ − πˆ 0 w Eˆ (ϕ( ˆ wˆ − πˆ 0 w)) ˆ ≤ C ˆ 0,∞,Kˆ ϕ ˆ 2,∞,Kˆ . Set ϕˆ = aˆ vˆ . Now, ˆ a2,∞,Kˆ ˆv0,∞,Kˆ + ˆa1,∞,Kˆ ˆv1,∞,Kˆ ). ˆavˆ 2,∞,Kˆ ≤ C(ˆ Again we may use the equivalence between the L∞ norms of vˆ and wˆ − πˆ 0 wˆ and the L2 norms of the same functions as in (i) since they belong to the finite dimensional space P1 . Also, by the triangle inequality, ˆ w wˆ − πˆ 0 w ˆ 0,Kˆ ≤ C ˆ 0,Kˆ . Thus we get (8.38)
ˆ a2,∞,Kˆ ˆv0,Kˆ + ˆa1,∞,Kˆ ˆv1,Kˆ )w Eˆ (ˆavˆ (wˆ − πˆ 0 w)) ˆ ≤ C(ˆ ˆ 0,Kˆ .
91
(iii) We can now complete the proof. Recall that EK (avw) = (det BK ) Eˆ (ˆavˆ w). ˆ Also, ˆam,∞,Kˆ ≤ C hm K am,∞,K
(8.39)
1
−2 ˆv2−m,Kˆ ≤ C h2−m K (det B K ) v2−m,K . 1
w ˆ 0,Kˆ ≤ C(det BK )− 2 w0,K ,
8. Numerical Integration
91
by Theorems 6.4 and 6.5. Combining (8.37), (8.38) and (8.39), we get EK (avw) ≤ C h2K (a1,∞,K v1,K + a2,∞,K v0,K )w0,K ≤ C h2K a2,∞,K v1,K w0,K .
Setting a = Ai j , v = the proof.
∂p′ ∂p ,w = we obtain (8.33), thus completing ∂x j ∂xi
We leave the second stage as an exercise: Exercise 8.5. Let PK = P2 and let the quadrature scheme be such that Eˆ (ϕ) ˆ = 0 for all ϕˆ ∈ P2 . Then show that for q such that W 2,q (K) ֒→ C 0 (K), there exists C independent of K such that for all f ∈ W 2,q (K) and all p ∈ PK , 1
1
EK ( f p) ≤ C h2K (det BK ) 2 − q  f 2,q,K p1,K . [Hint: If πˆ 1 is the orthogonal projection to P1 in the L2 sense then write Eˆ ( fˆpˆ ) = E ( fˆπˆ 1 pˆ ) + Eˆ ( fˆ( pˆ − πˆ 1 pˆ ))]. Remark 8.4. The inclusion W 2,q (K) ֒→ C 0 (K) is true if, for instance, 2 − nq > 0, by the Sobolev imbedding theorem. We now come to the final stage in the estimation of u − u∗h . Theorem 8.5. Let (th ) be a regular family of triangulations on Ω by 92 nsimplices of type (2). Let us assume that the Vh ellipticity is uniform with respect to h. Let Eˆ (ϕ) ˆ = 0 for all ϕˆ ∈ P2 . Then if u ∈ H 3 (Ω) ֒→ C 0 (Ω)(n ≤ 5), ai j ∈ W 2,∞ (Ω) and f ∈ W 2,q (Ω) for some q ≥ 2, we have the estimate (8.40)
u − u∗h 1,Ω ≤ C h2 [u3,Ω +  f 2,q,Ω ].
Proof. We estimate the various quantities in (8.28). We have: a(πh u, wh ) − ah (πh u, wh ) ≤
8. Numerical Integration
92
≤ ≤
n XX
K∈th i, j=1 n XX
K∈th i, j=1 2
≤h
n X
i, j=1
EK (ai j
∂(πh uK) ∂(wh K)  ∂x j ∂xi
C h2K ai j 2,∞,K 
∂(πh uK ) ∂(wh K) 1,K  0,K ∂x j ∂xi
12 X ∂(πh uK) 2 1,K ai j 2,∞,Ω  ∂x j K∈th
21 X ∂(wh K) 2  0,K ∂xi K∈th
(since hK ≤ h, and we may apply the CauchySchwarz inequality) ≤ C h2 πh u2,Ω wh 1,Ω Now, πh u2,Ω ≤ u2,Ω + u − πh u2,Ω ≤ Cu2,Ω ,
using Theorem 6.3 with P1 ⊂ PK = P2 . Therefore, for all wh ∈ Vh , we have a(πh u, wh ) − ah (πh u, wh ) ≤ Ch2 u2,Ω . wh 1,Ω
(8.41)
Similarly, we have  f (wh ) − fh (wh ) ≤ ≤ 93
Since q ≥ 2, X
K∈th
1 2
−
X
K∈th
X
K∈th
EK ( f wh K ) 1
1
C h2K (meas K) 2 − q  f 2,q,K wh 1,K .
1 ≥ 0 and by the general H¨older’s inequality, q 1
1
(meas K) 2 − q  f 2,q,K wh 1,K
21 − 1q X ≤ meas K K∈th
q1 X q  f  2,q,K K∈th
21 X 2 w  h 1,K K∈th
8. Numerical Integration
93
= C f 2,q,Ω wh 1,Ω . Hence we get, fot all wh ∈ Vh , (8.42)
 f (wh ) − fh (wh ) ≤ Ch2  f 2,q,Ω . wh 1,Ω
Combining (8.28), (8.29), (8.41) and (8.42) we get (8.40), thus completing the proof. Remark 8.5. The condition n ≤ 5 (needed for the continuous inclusion H 3 (Ω) ֒→ C 0 (Ω)) was already necessary for the definition of πh u. References: For a survey on numerical quadrature in general one may refer to Haber’s survey article [13]. For application to the finite element method, the Sec. 4.2 of the book by Strang and Fix [22] or Chapter 2 of the forthcoming book of Ciarlet and Raviart [5].
Chapter 9
The Obstacle Problem In Sec. 2 we cited the Obstacle Problem as an example of a nonlinear 94 abstract problem of Sec. 1. Let us recall a few facts about this to start with. Consider an elastic membrane (cf. Fig. 9.1) stretched over an open set Ω ⊂ R2 and fixed along the boundary Γ which is assumed to be Lipschitz continuous. Let a force of density Fdx act on the membrane. Let us assume the existence of an obstacle given by χ(x), for x ∈ Ω. Then vertical displacement given by u is the solution of the abstract problem where R P 2 ∂u ∂v dx a(u, v) = Ω (9.1) i=1 ∂xi ∂xi R f (v) = f v dx, f ∈ L2 (Ω) Ω
for u, v ∈ V = H01 (Ω), where f = F/t, t being the tension. The subset K is given by K = {v ∈ H01 (Ω); v ≥ χ a.e. inΩ}.
If v1 , v2 are in K and vi < χ in Ai with meas Ai = 0 for i = 1, 2, then λv1 + (1 − λ)v2 ≥ χ on (A1 ∪ A2 )c i.e. the complement of A1 ∪ A2 and meas (A1 ∪ A2 ) = 0. Thus K is convex. If v ∈ K, let vn ∈ K such that vn → v in H01 (Ω). Let vn ≥ χ in Acn , meas An = 0. Then all the vn are ≥ χ on (∪n An )c and meas (∪n An ) = 0. Hence v ≥ χ a.e. as well. Thus 95
9. The Obstacle Problem
96
v ∈ K and K is closed as well. We have the regularity assumption that χ ∈ H 2 (Ω). Of course it is
Figure 9.1: 95
The solution u satisfies (9.2)
J(u) = min J(v), v∈K
where J(v) = 21 a(v, v) − f (v) and is also characterized by the variational inequalities (cf. Theorem 1.1): (9.3)
a(u, v − u) ≥ f (v − u),
for all
v ∈ K.
We proposed as a problem to show that this problem is interpreted as the following classical problem (assuming u ∈ H01 (Ω) ∩ H 2 (Ω)). u ≥ χ in Ω, (9.4) −∆u = f where u > χ, u = 0 on Γ. We have a few regularity results which are listed below:
(i) If Ω is convex and Γ is a C 2 boundary then u ∈ H01 (Ω) ∩ H 2 (Ω). (ii) If f = 0 and Ω a convex polygon then also u ∈ H01 (Ω) ∩ H 2 (Ω). 96
(iii) The norm u2,Ω is bounded above by a function of  f 0,Ω and χ2,Ω in cases (i) and (ii).
9. The Obstacle Problem
97
Our aim in this section is to use the finite element method to approximate this problem and obtain error estimates. We list our assumptions now: Let Ω be a convex polygon, f ∈ L2 (Ω), χ ∈ H 2 (Ω) and let u ∈ 2 H (Ω) ∩ H01 (Ω). Remark 9.1. We cannot assume any more smoothness on u other than H 2 (Ω). For instance in the 1dimensional case if f = 0, and even if the function is very smooth the points of contact of u with χ will have discontinuous second derivatives in general; cf. Fig. 9.2.
Figure 9.2: With the above assumptions we proceed to the approximate problems, first in the abstract setting, as usual. We have the problems (Ph ) associated with the subspaces Vh ⊂ V = 1 H0 (Ω). We now choose closed convex subsets Kh ⊂ Vh . One has to bear in mind that, in general, Kh 1 K (we will see that this is the case in our approach, subsequently). We find uh ∈ Kh such that for all vh ∈ Kh , (9.5)
a(uh , vh − uh ) ≥ f (vh − uh ).
The existence and uniqueness of the uh follow from Theorem 1.1. Let H be a Hilbert space with norm  ·  and innerproduct (·, ·). Let 97 (V,  · ) be a subspace such that V ֒→ H, V = H. Then, as usual, if we identify H ′ and H, then H will be identified with a subspace of V ′ . (We
9. The Obstacle Problem
98
will take V = H01 (Ω) and H = L2 (Ω)). Also as in Sec. 1 (cf. proof of Theorem 1.2), for all u, v ∈ V we have (9.6)
a(u, v) = (Au)(v),
where A : V → V ′ is a linear map. We now pass on to an abstract error bound. Theorem 9.1 (FALK). Assume that f ∈ H, Au ∈ H. Then there exists a constant C, independent of Vh and Kh , such that (9.7)
"
# 21 u − uh  ≤ C inf (u − vh  + u − vh ) + inf uh − v . 2
vh ∈Kh
v∈K
(Note: The condition Au ∈ H = L2 (Ω) is satisfied if u ∈ H 2 (Ω) since Au = −∆u ∈ L2 (Ω).) Proof. Let α stand for the Vh ellipticity constant. Then
(9.8)
αu − uh 2 ≤ a(u − uh , u − uh )
= a(u, u) + a(uh , uh ) − a(u, uh ) − a(uh , u).
For any v ∈ K and any vh ∈ Kh , by (9.3) and (9.5), we have a(u, u) ≤ a(u, v) + f (u − v), (9.9) a(uh , uh ) ≤ a(uh , vh ) + f (uh − vh ). Substituting in (9.8) we get
αu − uh 2 ≤ a(u, v) + f (u − v) + a(uh , vh )
+ f (uh − vn ) − a(u, uh ) − a(uh , u)
= a(u, v − uh ) − f (v − uh ) + a(u, vh − u) − f (vh − u) + a(uh − u, vh − u)
= ( f − Au, u − vh ) + ( f − Au, uh − v) + a(uh − u, vh − u)
≤  f − Au u − vh  +  f − Au uh − v + Muh − u vh − u.
98
9. The Obstacle Problem
Notice that since (
r
99 α u − uh  − M
u − uh  u − vh  ≤
r
M u − vh )2 ≥ 0, we have α
1α M u − uh 2 + u − vh 2 , 2 M α
and hence αu − uh 2 ≤ C[u − vh  + uh − v] +
α M2 u − uh 2 + u − vh 2 2 2α
Or, u − uh 2 ≤ C[(u − vh  + u − vh 2 ) + uh − v]. Varying vh ∈ Kh and v ∈ K and extracting the square root after taking the infima we get (9.7). This completes the proof. Remark 9.2. If we have a linear problem then f = Au gives the solution and we get the original bound (3.1). Remark 9.3. From (9.8) we see that this estimate holds even if a(·, ·) is not symmetric. We now apply this to the specific membrane problem. Maintaining our assumptions on Ω, let th be a triangulation by triangles of type (1), and let Vh be the corresponding subspace of V = H01 (Ω). Remark 9.4. It is of no practical use if we go to more sophisticated finite elements, unlike the linear problem. Since u ∈ H 2 (Ω) is the maximum smoothness, we may atmost use our abstract estimate theorems only on the spaces P1 . One may be tempted to try for Kh those vh which are ≥ χ a.e. in Ω. However this is not of value from numerical and computational points of view for we do not easily know where exactly our piecewise linear 99 solution functions would touch χ. We set instead (9.10)
Kh = {vh ∈ Vh ; At all nodes b of th , vh (b) ≥ χ(b)}.
9. The Obstacle Problem
100 Remark 9.5.
Figure 9.3: As seen in Fig. 9.3, though for nodes b, vh (b) ≥ χ(b), it does not guarantee that vh ≥ χ a.e. Thus we see that Kh 1 K. Now the relation (9.10) is very easy to implement using the computer. We now have our main result on the error bound. Theorem 9.2 (FALK). There exists a constant C depending on  f 0,Ω and χ2,Ω such that for a regular family of triangulations (th ) as above we have u − uh 1,Ω ≤ C h.
(9.11)
Remark 9.6. The order of convergence is therefore the same as that for the linear problems when we use piecewise linear approximations. Proof. By Theorem 9.1, # 21 " 2 u − uh 1,Ω ≤ C inf u − vh 1,Ω + u − vh 0,Ω + inf uh − v0,Ω . vh ∈Kh
100
v∈K
(i) We first estimate the infimum over Kh . Note that if vh = πh u, then vh ∈ Vh . Also for all nodes b, vh (b) = πh u(b) = u(b) ≥ χ(b). Thus vh ∈ Kh as well. Thus, inf u − vh 21,Ω + u − vh 0,Ω ≤ u − πh u21,Ω + u − πh u0,Ω vh ∈Kh ≤ C h2 u22,Ω + u2,Ω .
9. The Obstacle Problem
101
(ii) For the infimum over K, consider v1 = max(uh , χ). Clearly v1 ≥ χ and v belongs to H 1 (Ω) because uh and χ also belong to H 1 (Ω) (this is a nontrivial result which we assume here). Hence v1 ∈ K, and thus, inf uh − v0,Ω ≤ uh − v1 0,Ω . v∈K
We have uh −
v1 20,Ω
=
Z
uh − χ2 dx, where Λh = {x; χ(x) ≥ uh (x)}.
Λh
If πh χ is the Vh interpolate of χ, then for all nodes b, uh (b) ≥ χ(b) = πh χ(b). Since both uh and πh χ are piecewise linear, we may now assert that uh ≥ πh χ everywhere. Thus uh − πh χ ≥ 0 on Ω. Thus for all x ∈ Λh , we have 0 < (χ − uh )(x) = (χ − uh )(x) ≤ (χ − πh χ)(x)
≤ (χ − πh χ)(x)
and for x ∈ Ω − Λh , (χ − uh )(x) = 0, so that uh − v1 0,Ω
21 Z ≤ χ − πh χ2 dx = χ − πh χ=0,Ω ≤ C h2 χ2,Ω . Ω
Hence
1
u − uh 1,Ω ≤ C(h2 ) 2 = Ch, where C depends on χ2,Ω and u2,Ω . However the regularity result (iii) 101 helps us to bound u2,Ω above by a constant C depending on  f 0,Ω and χ2,Ω which completes the proof of the theorem. References: Two important references are Falk [11] and [12]. For regularity results refer Brezis and Stampacchia [3] and Lewy and Stampacchia [16]. Another references is Mosco and Strang [19].
Chapter 10
Conforming Finite Element Method for the Plate Problem In Sec. 2 (cf. Example 2.4), as an example of fourthorder problem, we 102 described the plate problem. In abstract terms it is to find the solution of a(u, v) = f (v),
for all
v ∈ V,
where (10.2) K = V = H02 (Ω), Ω ⊂ R2 , ( )! 2 R ∂2 v ∂2 u ∂2 v ∂2 u ∂2 v a(u, v) = ∆u · ∆v + (1 − σ) 2 ∂ u dx, − − ∂x1 ∂x2 ∂x1 ∂x2 ∂x21 ∂x22 ∂x22 ∂x22 RΩ 2 f (v) = f v dx, f ∈ L (Ω). Ω
The problem was interpreted as the classical boundary value problem (10.3)
2 ∆ u = f in Ω, ∂u u = = 0 on Γ = ∂Ω. ∂ν 103
104
10. Conforming Finite Element Method for the Plate Problem
Remark 10.1. It was commented in Sec. 2 that the second term of the integrand in the definition of a(·, ·) does not contribute to the differential equation. Our method here will be equally applicable to both the cases. viz. with the second term present (the plate problem) or with that term absent (as it happens in Hydrodynamics). In our next section, on nonconforming methods, we will see that the second term is essential in order that we may apply that method. We assume that Ω is a polygonal domain in R2 . We saw in Sec. 3 that for fourthorder problems we need the inclusion Vh ⊂ C 1 (Ω). (cf. Exercise 3.1). Thus we need to use finite elements of class C 1 , such as the Argyris triangle, the BognerFoxSchmidt rectangle and so on (cf. Sec. 4). 103
When such finite elements can be imbedded in an affine family, then we have the approximation theory, for regular families of triangulations, available to us. We show that this is the case for the BognerFogSchmidt rectangle. However for the Argyris triangle or for the 18degreeoffreedom triangle such an imbedding is not possible and we have to modify the usual argument to obtain error estimates. The “minimal assumptions” for 0(h) convergence in the  · 2,Ω norm are that P2 ⊂ P and that u ∈ H 3 (Ω) ∩ H02 (Ω). We have that if Ω is a convex polygon and if f ∈ L2 (Ω), then u ∈ H 3 (Ω) ∩ H02 (Ω). This result is due to Knodrat´ev. We will go through the various examples of triangulations of class C 1 and study convergence in these cases. Example 10.1. The BognerFogSchmidt rectangle (cf. Exercise 4.9). Let PK = Q3 (dim PK = 16). We then have (cf. Fig. 10.1):
(10.4)
X K
) ( ∂p ∂2 p ∂p (ai ), (ai ), (ai ); 1 ≤ i ≤ 4 . = p(ai ), ∂x1 ∂x2 ∂x1 ∂x2
10. Conforming Finite Element Method for the Plate Problem
105
Figure 10.1: Equivalently, one may also use the PK unisolvent set
(10.5)
′ X K
= p(ai ), Dp(ai )(ai+1 − ai ), Dp(ai )(ai−1 − ai ), D2 p(ai )(ai+1 − ai , ai−1 − ai ); 1 ≤ i ≤ 4
(all indices being read modulo 4). 104 Recall that for an affine family of finite elements, the degrees of freedom p(a0i ), Dp(a1i )(ξik1 ), (D2 p(a2i )(ξik2 , ξil2 ) are such that (cf. Sec. 5): a0i = F(ˆa0i ), . . . , a2i = F(ˆa2i ), (10.6) ξ 1 = BK ξˆ 1 , . . . , ξ 2 = BK ξˆ2 i,k i,k i,1 i,1
for then πd ˆ vˆ which is essentially what we need for the abstract Kv = π error analysis. P In ′K note that (10.7)
ai+1 − ai = F(ˆai+1 ) − F(ˆai ) = BK (ˆai+1 − aˆ i ),
and so on. Thus it is clear that this rectangle can be imbedded in an affine family of finite elements. Now PK ⊂ Pˆ ⊂ Q3 for k = 3. By our abstract error analysis, we therefore have (10.8)
u − uh 2,Ω ≤ Ch2 u4,Ω .
assuming sufficient smoothness on u as usual.
106
105
10. Conforming Finite Element Method for the Plate Problem
A word about the boundary conditions. As in Exercise 3.1, we get ∂v = 0 on Γ, for v ∈ Vh . Thus in choosing our that Vh ⊂ H02 (Ω) if v = ∂ν basis functions we must assure ourselves that this condition is satisfied. This in turn depends on the values at the boundary nodes. Let b and c be two nodes on Γ such that the line joining them is parallel to (say) the x1 axis. Since we need v = 0 on this line, and since v will be a polynomial ∂v (b) = in x1 of degree ≤ 3 on this line we must have v(b) = v(c) = 0, ∂x1 ∂v ∂v ∂v ∂v = 0 on this line and since = (c) = 0. Also since we need ∂x1 ∂ν ∂ν ∂x2 ∂v ∂v (b) = (c) = 0 is a polynomial in x1 of degree ≤ 3, we need to set ∂x2 ∂x2 2 2 ∂ v ∂ v (b) = (c) = 0. Thus the degrees of freedom on all and ∂x1 ∂x2 ∂x, ∂x2 boundary nodes must be zero. The only “free” or “unknown” parameters are the degrees of freedom at the interior nodes. This takes care of the boundary conditions. Let us now turn to the Argyris triangle (cf. Example 4.7).
Figure 10.2: P We recall that PK = P5 , dim PK = 21, and K is given by (cf. Fig. 10.2) (10.9) X ∂p ∂3 p ∂p (ai ), . . . , 2 (ai ), 1 ≤ i ≤ 3; (ai j ), 1 ≤ i < j ≤ 3 . = p(ai ), ∂x1 ∂ν ∂x K
2
10. Conforming Finite Element Method for the Plate Problem
107
We may replace the first and second derivative values at the vertices by Dp(ai )(ai+1 − ai ), Dp(ai )(ai−1 − ai ), D2 p(ai )(ai+1 − ai , ai+1 − ai ), D2 p(ai )(ai+1 − ai , ai−1 − ai ), D2 p(ai )(ai−1 − ai , ai−1 − ai ) in order to get degrees of freedom for which the relations of the type (10.6) may be ∂p (ai j ), satisfied. However, one cannot replace the normal derivatives ∂ν 1 ≤ i < j ≤ 3, by such quantities since affine transformations do not preserve orthogonality. In order to estimate the errors we describe an “intermediary” finite element: Example 10.2. The Hermite Triangle of Type (5). Let PK = P5 (dim PK = 21). Define X n = p(ai ), Dp(ai )(ai+1 − ai ), Dp(ai )(ai−1 − ai ), K
D2 p(ai )(ai+1 − ai , ai+1 − ai ), D2 p(ai )(ai+1 − ai , ai−1 − ai ), D2 p(ai )(ai−1 − ai , ai−1 − ai ), 1 ≤ i ≤ 3;
o Dp(ai j )(ak − ai j ), 1 ≤ i < j ≤ 3, k , 1, , j .
That is to say, the only change compared to the Argyris triangle is that we have replaced the normal derivatives at ai j by the derivatives along the line joining ai j to ak , the opposite vertex. Symbolically we can represent such a triangle as in Fig. 10.3.
Figure 10.3: This element can be put in an affine family as is readily seen. If ΛK
106
108
10. Conforming Finite Element Method for the Plate Problem
is the associated interpolation operator, our error analysis yields (10.10)
v − ΛK vm,K
h6K ≤ C m v6,K , 0 ≤ m ≤ 6, ρK
for v ∈ H 6 (K). Remark 10.2. Though the Hermite triangle of type (5) yields an affine family, one cannot use it since Vh ⊂ C 0 (Ω) but, in general, Vh 1 C 1 (Ω) as is necessary for fourthorder problems. This is so because the adjacent triangles will not patch up, in general, in their derivatives along the medians; cf. 10.4.
Figure 10.4: 107
Again we show how to take care of the boundary conditions in the ∂v = 0 on Γ. Let us have two nodes Argyris triangle. We need again v = ∂ν ′ b, b , the vertices of a triangle lying on Γ with midpoint c. On this line v will be a polynomial of degree ≤ 5 in τ, an abscissa along this line. ∂v will be a polynomial in τ of degree ≤ 4 on this line. Hence for v = 0 ∂ν ∂v ∂v ′ ∂2 v on Γ we need to set, v(b) = v(b′ ) = 0, (b) = (b ) = 0, 2 (b) = ∂τ ∂τ ∂τ ∂v ∂v ′ ∂v ∂v ∂2 v ′ = 0 on Γ we set, (b) = (b ) = (c) = 0, (b ) = 0. For ∂ν ∂ν ∂ν ∂ν ∂τ2 ∂2 v ∂2 v ′ (b) = (b ) = 0. Thus the only free or unknown parameters ∂τ∂ν ∂τ∂ν ∂2 v are 2 at vertices on Γ and the degrees of freedom at all interior nodes. ∂ν
10. Conforming Finite Element Method for the Plate Problem
109
We now get an error estimate when we have triangulations of Argyris triangles. We use our usual terminology more loosely here1 . By a regular family of triangulations made up of Argyris triangles we mean hK ≤ σ, a that all th consist only of Argyris triangles and that for all K, ρK constant. We also assume that if h = max hK , then h → 0. K∈th
Theorem 10.1. For a regular family (th ) of triangulations made up of Argyris triangles (10.11)
v − πh vm,Ω ≤ C h6−m v6,Ω , 0 ≤ m ≤ 6. 108
−ν Proof. Let us denote the opposite vertex of ai j (i < j) by ak . Let → K −τ be the unit vector along the line be the unit outernormal at ai j and → K [ai , a j ], at ai j (cf. Fig. 10.5).
Figure 10.5: Let πK be the interpolation operator for the Argyris triangle K and let ΛK be that for the corresponding Hermite triangle of type (5). 1 Because we have to drop the assumption that all the finite elements are affine equivalent to a reference finite element.
10. Conforming Finite Element Method for the Plate Problem
110
Set δ = πK v − ΛK v. Then δ ∈ P5 . Now ∂δ ∂ ∂ (ai j ) = (πK v − ΛK v)(ai j ) = (v − ΛK v)(ai j ). ∂νK ∂νK ∂νK Also since πK v = ΛK v along any side of K (since the values of these polynomials of degree 5 as well as those of their first and second ∂δ = 0. derivatives agree at the endpoints), we have ∂τK Since Dδ(ai j )(ak − ai j ) =
∂δ −ν i + ∂δ (a )ha − a ,→ − (ai j )hak − ai j ,→ K ij k i j τ K i, ∂νK ∂τK
where h·, ·i is the Euclidean innerproduct, substituting for at ai j , we get
Dδ(ai j )(ak − ai j ) =
(10.12) 109
∂δ ∂δ and ∂νK ∂τK
∂ −ν i (v − ΛK v)(ai j )hak − ai j ,→ K ∂ν
Since δ ∈ P5 , using the unisolvency in the Hermite triangle we may express δ in terms of its basis functions. Since all degrees of freedom except those of the type Dδ(ai j )(ak − ai j ) are zero for δ, we have δ=
(10.13)
X
∂ −ν ip . (v − ΛK v)(ai j )hak − ai j ,→ K i jk ∂ν K 1≤i< j≤3 k,i,k, j
Now −ν i ≤ a − a  → −ν  ≤ h , hak − ai j ,→ K k ij K K
(10.14) and 
∂ (v − ΛK v)(ai j ) ≤ v − ΛK v1,∞,K ∂νK 6 1 h ≤ C(meas K)− 2 K v6,K . ρK
10. Conforming Finite Element Method for the Plate Problem
111
(Theorem 6.6 with m = 1, k = 5, p = 2, q = ∞). Also, meas K ≥ Cρ2K , and we have 
(10.15)
h6 ∂ (v − ΛK v)(ai j ) ≤ C 2K v6,K . ∂νK ρK
Finally, by Theorem 6.4 and 6.5 pi jk m,K ≤ C
(10.16)
hK  pˆ i jk m,Kˆ . ρm K
Combining (10.14), (10.15) and (10.16), we get δm,K ≤
(10.17)
X
1≤i< j≤3 k,i,k, j
≤C
h8K ρm+2 K

∂ −ν i p  (v − ΛK v)(ai j ) hak − ai j ,→ K i jk m,K ∂νK
v6,K ,
and hence, v − πK vm,K ≤ v − ΛK vm,K + δm,K h6K h2K ≤ C m 1 + 2 v6,K (Using (10.10) and (10.17)) ρK ρK hK ≤ σ). ≤ C h6−m v6,K (since hK ≤ h, ρK 110
This on summing over K gives (10.11), thus completing the proof. Exercise 10.1. Perform the same analysis for the 18degreeof freedom triangle. (cf. Exercise 4.7). For the interpolation theory of the HCTtriangle (cf. Exercise 4.8), the normal derivatives are handled as in the present case. However the arbitrariness of the interior point is an obstacle to be overcome. For a discussion of this, see Ciarlet [4].
112
10. Conforming Finite Element Method for the Plate Problem
Another finite element, similar in its principle to the HCTtriangle used in the conforming finite element method for the plate problem is the Fraeijs de Veubeke and Sander Quadrilateral. See Ciavaldini and N´ed´elec [9]. These are all essentially the finite elements used in the “conforming” methods to approximate the plate problem (we will define such methods at the beginning of Sec. 11).
Chapter 11
NonConforming Methods for the Plate Problem WE START WITH a brief classification of finite element methods. The 111 first class of methods are called conforming methods, which we have described upto now, except when we considered numerical integration. The second class consists of methods other than conforming. In the latter class we have the Nonconforming methods included: Given the abstract problem, the conforming methods deal with the finding of subspaces Vh ⊂ V and solving the problems (Ph )
ah (uh , vh ) = fh (vh ),
for all
vh ∈ Vh ,
where ah = a and fh = f for all h and uh ∈ Vh is the required solution. When we employ methods other than conforming we commit, in the terminology of G. Strang, “Variational Crimes”. (See Strang and Fix [22]). These may occur in the following ways: (i) When performing numerical integration, we may have ah and fh different from a and f respectively. However, Vh is a subspace of V; (ii) The boundary Γ of Ω may be curved. In this case triangles lying in the interior will be triangles of straight edges while those meeting the boundary will have curved edges like parabolas. These are 113
114
11. NonConforming Methods for the Plate Problem the socalled “isoparametric” finite elements. Hence if Ωh is the union of the finite elements of the triangulation th , then, in general, Ωh , Ω and consequently Vh 1 V (where Vh is a space of functions defined over Ωh ), ah , a, fh , f ; for a discussion of these, see Ciarlet and Raviart [30], [31].
(iii) When employing nonconforming methods (which will be dealt with subsequently) though Ωh = Ω, fh = f , we will have Vh 1 V and ah , a.
112
(iv) One may employ any combination of the above three.
Let us return to the plate problem. For a conforming method we need the inclusion Vh ⊂ H02 (Ω) which essentially results from the inclusion Vh ⊂ C 1 (Ω). Because of this necessity, when compared with secondorder problems, we either have the dimension of PK “large” (as in the case of the Argyris triangle) or that the structure of PK is complicated (as in the HCTtriangle). Also one would like to have just PK = P2 since u is only in H 3 (Ω) in most cases, but this is impossible by the ˇ sek result (cf. Remark 4.3) which stresses that at least polynomials Zeniˇ of degree 5 must be present in PK . Hence the desire to surmount these difficulties led to the devising of nonconforming methods, essentially developed by the Engineers. Since the root of all trouble is the inclusion Vh ⊂ H02 (Ω), we drop this condition. Thus we start with Vh ⊂ C 0 (Ω) and it is much easier from the computer programme view point. This of course, works only for a few finite elements, and we describe one of them. Example 11.1. The Adini’s rectangle; cf. Fig. 11.1.
11. NonConforming Methods for the Plate Problem
115
Figure 11.1:
The element K consists of a rectangle with vertices {ai , 1 ≤ i ≤ 4}; the space PK is given by PK = P3 ⊕ {x1 x32 } ⊕ {x31 x2 }, by which we mean polynomials of degree ≤ 4 whose only fourthdegree terms are those involving x1 x32 and x31 x2 . Thus P3 ⊂ PK . We have the set of degrees of freedom: X K
) ∂p ∂p (ai ), (ai ), 1 ≤ i ≤ 4 . = p(ai ), ∂x1 ∂x2 (
113
Of course this element can be used only for plates with sides parallel to the coordinate axes, such as rectangular plates. P Exercise 11.1. Show that in Example 11.1, K is PK unisolvent and that Adini’s rectangle is a finite element of class C 0 , and, in general, not of class C 1 . Thus we get ‘a priori’ that Vh ⊂ H 1 (Ω). For the boundary condition, we set all degrees of freedom on the boundary nodes as zero. This gives us that Vh ⊂ H01 (Ω). Thus the only ‘unknown’ or ‘free’ parameters are ∂vh is zero only the degrees of freedom at the interior nodes. Note that ∂ν at the boundary nodes, in general.
11. NonConforming Methods for the Plate Problem
116
In the abstract problem, we have a(·, ·) and f (·) given by
(11.1) 2 R ∂2 u ∂2 v ∂2 u ∂2 v ∂2 v ∆u∆v + (1 − σ) 2 ∂ u dx, a(u, v) = − − ∂x1 ∂x2 ∂x1 ∂x2 ∂x21 ∂x22 ∂x22 ∂x21 Ω R 2 f (v) = f vdx, f ∈ L (Ω). Ω
The second integral is defined over Vh as well. Thus for the discrete problem (Ph ) we may set fh = f . However while for u, v ∈ Vh the first integral is defined over each K ∈ th , we cannot define it over Ω, since we get Dirac measurelike terms along the boundary. To get over this, we now define !# ∂2 uh ∂2 uh ∂2 vh ∂2 uh ∂2 vh − dx − ∆uh ∆vh + (1 − σ) 2 ah (uh , vh ) = ∂x1 ∂x2 ∂x1 ∂x2 ∂x21 ∂x22 ∂x21 K∈th K (11.2)X Z " !# ∂2 uh ∂2 vh ∂2 uh ∂2 vh ∂2 uh ∂2 vh dx, σ∆uh ∆vh + (1 − σ) = + + 2 ∂x1 ∂x2 ∂x1 ∂x2 ∂x21 ∂x21 ∂x22 ∂x22 K∈t XZ "
h K
and we have the discrete problem (Ph ): To find uh ∈ Vh such that for all vh ∈ Vh (11.3)
ah (uh , vh ) = f (vh ).
114
We now prove the existence and uniqueness of the solution uh for (Ph ). We define on Vh the seminorm (11.4)
21 X vh 22,K 22,K . vh h = K∈th
Notice that this may be defined over V = H02 (Ω) as well and for v ∈ V, vh = v2,Ω . In the same way for u, v ∈ V, ah (u, v) = a(u, v). We now show that the seminorm  · h is indeed a norm on Vh . Let ∂vh vh ∈ Vh with vh h = 0. This gives that = constant over any K. But ∂x1 ∂vh at the common vertices given adjacent finite elements the value of ∂x1 ∂vh is constant over Ω. But this is zero on the coincide and hence ∂x1
11. NonConforming Methods for the Plate Problem
117
∂vh ∂vh = 0 on Ω. Similarly = 0 on Ω. Since ∂x1 ∂x2 Vh ⊂ C 0 (Ω) and vh = 0 on Γ, the above conditions give that vh = 0 over Ω. Thus (11.4) defines a norm on Vh . To show the existence and uniqueness of the solution of (Ph ), we show that ah (·, ·) is Vh elliptic. In fact we do more than this. We show that the ah (·, ·) are Vh elliptic uniformly with respect to h. Recall that from physical considerations, 0 < σ < 12 (see Sec. 2). Now XZ a(vh , vh ) = σ(∆vh )2 dx + (1 − σ)vh 2h K∈th K (11.5) boundary nodes. Hence
≥ (1 − σ)vh 2h . Remark 11.1. It was mentioned in passing in Sec. 10 that in order to apply nonconforming methods one needed the second term involving σ in the integral defining a(·, ·). The uniform Vh ellipticity could not be got in the Hydrodynamical case where this term is absent. We now proceed with the abstract error analysis. Theorem 11.1 (STRANG). Let ah (·, ·) be Vh elliptic uniformly with respect to h with e α > 0 so that for all vh ∈ Vh
(11.6)
ah (vh , vh ) ≥ e αvh h .
e such that for all uh , vh ∈ Vh Let in addition, there exist M
(11.7)
e h h vh h . a(uh , vh ) ≤ Mu
Assume that ah = a and  · h =  ·  on V. (These are needed to extend the definition of ah and  · h to V). Then there exists a constant C, independent of h, such that ) (  f (wh ) − ah (u, wh ) (11.8) u − uh h ≤ C inf u − vh h + sup vh ∈Vh wh h wh ∈Vh uh being the solution of (Ph ).
115
11. NonConforming Methods for the Plate Problem
118
Proof. For all vh ∈ Vh we have u − uh h ≤ u − vh h + uh − vh h .
(11.9)
Now for any vh ∈ Vh , f (uh − vh ) = ah (uh , uh − vh ), so that we may write e αuh − vh 2h ≤ ah (uh − vh , uh − vh ) and thus,
= ah (u − vh , uh − vh ) + f (uh − vh ) − ah (u, uh − vh ) e − vh h uh − vh h +  f (uh − vh ) − ah (u, uh − vh ), ≤ Mu e M u − vh h + e α e M ≤ u − vh h + e α
uh − vh h ≤
116
1  f (uh − vh ) − ah (u, uh − vh ) e α uh − vh h 1  f (wh ) − ah (u, wh ) sup . e α wh ∈Vh wh h
Substituting in (11.9) and varying vh ∈ Vh and taking the infimum we get (11.8). This completes the proof. Remark 11.2. In case the method is conforming, then ah (u, wh ) = a (u, wh ) = f (wh ) and the second term disappears in (11.8), leaving us with the original bound (3.1). We note that for the ah (·, ·) defined for the plate problem by (11.2), the conditions of Theorem 11.1 are satisfied. The condition (11.6) is embodied in (11.5). The condition (11.7) follows from the similar (continuity) condition on a(·, ·) and an application of the CauchySchwarz inequality. Exercise 11.2. Let (H) be a Hilbert space with innerproduct (·, ·) and norm  · . Let (V,  · ) be a subspace such that V ֒→ H and V = H. Let Vh ⊂ H. Define Eh (u, v) = ( f, v) − ah (u, v),
for all
u, v ∈ Vh ∪ V.
11. NonConforming Methods for the Plate Problem
119
Then show that 1 e − uh h sup inf ϕ − ϕh h u − uh  ≤ Mu g∈H g ϕh ∈Vh # " 1 inf (Eh (u, ϕ − ϕh ) + Eh (ϕ, u − uh )) + sup g∈H g ϕh ∈Vh where for all v ∈ V, a(v, ϕ) = g(v). We now go on with the error analysis and study the order of convergence. We assume that u ∈ H 3 (Ω) ∩ H02 (Ω) which is quite realistic from the regularity results. Now, since πh u ∈ Vh , we have that inf u − vh h ≤ u − πh uh .
vh ∈Vh
117
Again applying error bounds for each K and then summing over all K we get u − πh uh ≤ Chu3,Ω . Thus inf u − vh h ≤ Chu3,Ω .
(11.10)
vh ∈Vh
Our aim is to get a similar estimate for the second term in (11.8). In fact we show that (11.11)
sup wh ∈Vh
 f (wh ) − ah (u, wh ) ≤ Chu3,Ω . wh h
This entails more work. We define Eh (u, wh ) = f (wh ) − ah (u, wh )
(11.12)
for u ∈ H 3 (Ω) ∩ H02 (Ω), wh ∈ Vh . Since wh ∈ H01 (Ω), there exists a sequence {wnh } in D(Ω) converging to wh in H01 (Ω). Hence, Z Ω
f wnh dx =
Z " Ω
∆u∆wnh + (1 − σ) 2
∂2 u ∂2 wnh ∂2 u ∂2 wnh ∂2 u ∂2 wnh − 2 − 2 ∂x1 ∂x2 ∂x1 ∂x2 ∂x1 ∂x22 ∂x2 ∂x21
!#
dx
11. NonConforming Methods for the Plate Problem
120 =−
Z
(grad ∆u)(grad wnh )dx,
Ω
by Green’s formula. The term involving (1 − σ), by Lemma 2.2, can be converted to an integral over Γ. All integrals over Γ vanish since wnh ∈ D(Ω). Since both sides of the above relation are continuous linear functionals on H01 (Ω), we can pass to the limit to obtain (11.13)
f (wh ) =
Z
f wh dx = −
Ω
Z
(grad ∆u)(grad wh )dx
Ω
for all wh ∈ Vh . Now, ah (u, wh ) =
XZ
K∈th
K
" #! ∂2 u ∂2 wh ∂2 u ∂2 wh ∂2 u ∂2 wh − 2 ∆u∆wh + (1 − σ) 2 dx − ∂x1 ∂x2 ∂x1 ∂x2 ∂x1 ∂x22 ∂x22 ∂x21
X Z − grad ∆u grad wh dx = K∈th
(11.14)
+
I
K
∂wK ∆u h dγ + (1 − σ) ∂νK
∂K
118
I
∂K
! ∂2 u ∂wh ∂2 u ∂wh dγ , − 2 + ∂ν ∂ν ∂τ ∂τ ∂τK K K K K
by Green’s formula (2.15) and Lemma 2.2 again. Notice however that by standard orientation arguments and the continuity of wh over Ω,
(11.15)
XI
K∈th
∂K
∂2 u ∂wh dγ = 0. ∂νK ∂τK ∂τK
Using (11.13), (11.14) and (11.15), we substitute in (11.12) to get
(11.16)
X I ∂2 u ∂wh −∆u + (1 − σ) Eh (u, wh ) = dγ. 2 ∂ν ∂τ K K K∈t h ∂K
11. NonConforming Methods for the Plate Problem
121
Figure 11.2: Splitting the boundary into four parts as in Fig. 11.2, we get (11.17)
Eh (u, wh ) =
X
K∈th
∆1,K
! !! ∂wh ∂wh u, + ∆2,K u, , ∂x1 ∂x2
where for j = 1, 2, we define ∆ j,K (11.18)
! Z ∂wh = u, ∂x j
K ′j
!! ∂2 u ∂wh ∂wh −∆u + (1 − σ) − Λ dγ K ∂x j ∂τ2K ∂x j
!! Z ∂wh ∂2 u ∂wh − ΛK dγ, − −∆u + (1 − σ) 2 ∂x j ∂τK ∂x j ′′ Kj
ΛK being the Q1 interpolation operator associated with the values at the 119 four vertices. Note that (11.17) is the same as (11.15). This is evident provided we show that the contribution of the terms involving ΛK is zero. This is so because ! wh = 0 on the boundary Γ and on the common boundaries, ∂wh is linear and equal in value for both adjacent finite elements ΛK ∂x j since it agrees at the vertices, but occurs with opposite signs as is obvi′ ′′ ous from Fig. 11.3 (K11 = K12 ).
11. NonConforming Methods for the Plate Problem
122
Figure 11.3: We also record that
(11.19)
∂wh ∈ ∂ j PK where ∂x j K (
) ∂p ∂ j PK = : K → R; p ∈ PK . ∂x j
We now prove a result analogous to the BrambleHilbert lemma which will help us to estimate that ∆ j,K ’s and hence Eh (u, wh ).
120
Theorem 11.2 (BILINEAR LEMMA). Let Ω ⊂ Rn be open with Lipschitz continuous boundary Γ. Let W be a subspace of W l+1,q (Ω) such that P1 ⊂ W. Let b be a continuous bilinear form over W k+1,p (Ω) × W such that b(p, w) = 0 for all p ∈ Pk , w ∈ W b(v, q) = 0 for all v ∈ W k+1,p (Ω), q ∈ P1 . Then there exists a constant C = C(Ω) such that for all v ∈ W k+1,p (Ω), w ∈ W,
(11.21)
b(v, w) ≤ Cb vk+1,p,Ω wl+1,q,Ω .
Proof. For a given w ∈ W, b(·, w) : v 7→ b(v, w) is a continuous linear form on W k+1,p (Ω) vanishing on Pk . Hence by the BrambleHilbert lemma, (11.22)
b(v, w) ≤ Cb(·, w)∗k+1,p,Ω vk+1,q,Ω .
11. NonConforming Methods for the Plate Problem
123
However for all q ∈ P1 , b(v, w) = b(v, w + q) so that b(v, w) ≤ b vk+1,p,Ω w + ql+1,q,Ω and hence b(v, w) ≤ b vk+1,p,Ω inf w + ql+1,q,Ω q∈P1
≤ Cb vk+1,p,Ω wl+1,q,Ω , by the theorem 6.2, so that b(·, w)∗k+1,p,Ω ≤ Cb wl+1,q,Ω . and substituting in (11.22), we get (11.21), which completes the proof. We may now prove the theorem on our error estimate and order of convergence. Theorem 11.3. For a regular family (th ) of triangulations made up of Adini’s rectangles (11.23)
u − uh h ≤ C hu3,Ω .
Proof. Let us first estimate ∆ j,K  for j = 1, 2. Set ∂2 u ϕ = −∆u + (1 − σ) ∈ H 1 (K), since u ∈ H 3 (Ω), 2 ∂τ K ∂wh ∈ ∂1 P K . v = ∂x1 Define
(11.24)
δ1,K (ϕ, v) =
Z
K1′
ϕ(v − ΛK v)dγ −
Z
K1′′
ϕ(v − ΛK V)dγ,
for v ∈ ∂1 P, ϕ ∈ H 1 (K). If h2 is the length of K1′ (and K1′′ ) and h1 is that of K2′ (and K2′′ ) we have by a simple change of variable (11.25)
δ1,K (ϕ, v) = h2 δ1,Kˆ (ϕ, ˆ vˆ ),
121
11. NonConforming Methods for the Plate Problem
124
where Kˆ is the reference finite element. Since P0 ⊂ Q1 which is preˆ δ1,Kˆ (ϕ, served by ΛK we have that for all vˆ ∈ P0 , ϕˆ ∈ H 1 (K), ˆ vˆ ) = 0. d Now let ϕˆ ∈ P0 and vˆ ∈ ∂1 P. We wish to show that δ1,Kˆ (ϕ, ˆ vˆ ) = 0: ˆ We may take for K, the unit square. Since ϕˆ ∈ P0 , its value on Kˆ is a constant, say, b0 . Now let vˆ ∈ ∂d ˆ is of the form 1 P. Then v (11.26) vˆ = a0 + a1 x1 + a2 x2 + a3 x21 + a4 x1 x2 + a5 x22 + a6 x21 x2 + a7 x32 . Taking the values at the four vertices we get ΛK vˆ = a0 + (a1 + a3 )x1 + (a2 + a5 + a7 )x2 + (a4 + a6 )x1 x2 .
(11.27)
Now K1′ is the line x1 = 1 and K1′′ is the line x1 = 0. Thus = −(a5 + a7 )x2 + a5 x22 + a7 x32 , vˆ − ΛK vˆ x1 =0 vˆ − ΛK vˆ = −(a5 + a7 )x2 + a5 x22 + a7 x32 . x =1 1
Hence, Z
′
ϕ(ˆ ˆ v − ΛK vˆ )dγ =
Z1 0
(11.28)K1 =
Z
K1′′
122
b0 (−(a5 + a7 )x2 + a5 x22 + a7 x32 )dx2 ϕ(ˆ ˆ v − ΛK vˆ )dγ.
Thus δ1,Kˆ (ϕ, ˆ vˆ ) = 0 for ϕˆ ∈ P0 , vˆ ∈ ∂d 1 P. Note further that the bilinear form δ1,Kˆ is continuous, for δ1,Kˆ (ϕ, ˆ vˆ ) ≤ Cϕ ˆ L2 (∂K) vL2 (∂K) ˆ ˆ ˆ . ≤ Cϕ ˆ 1,Kˆ ˆv1,Kˆ , by the Trace theorem (cf. Th. 2.3). Thus we may apply the bilinear lemma to the bilinear form δ1,Kˆ with l = k = 0 to get (11.29)
δ1,Kˆ (ϕ, ˆ vˆ ) ≤ Cϕ ˆ 1,Kˆ ˆv1,Kˆ .
11. NonConforming Methods for the Plate Problem
125
We also have the relations 1 ˆ 1,Kˆ ≤ CBK   det BK − 2 ϕ1,K , ϕ (11.30) ˆv ˆ ≤ CBK   det BK − 12 v1,K . 1, K
Now BK  ≤ C hK and  det BK  = meas K/ meas Kˆ ≥ Cρ2K , and thus C hK 1 BK   det BK − 2 ≤ ≤ C. Also h2 ≤ hK , so that ρK (11.31) where
∆1,K (ϕ, v) ≤ C hK u3,K wh 2,K ∂2 u ϕ = −∆u + (1 − σ) ∂τ2K , v = ∂wh , ∂x1
and similarly, (11.32)
∆2,K (ϕ, v) ≤ C hK u3,K wh 2,K .
These inequalities lead us to the estimate (11.33)
 f (wh ) − ah (u, wh ) ≤ Chu3,Ω wh h ,
for a regular family of triangulations made up of Adini’s rectangles. 123 Thus varying wh over Vh and taking the supremum, we get the estimate (11.11). Using (11.10) and (11.11) and substituting in (11.8) we get the required estimate as given in (11.23). This completes the proof. Remark 11.3. By the Duality Argument, Lesaint and Lascaux [15] have proved that u − uh 1,Ω ≤ Ch2 u4,Ω assuming u ∈ H 4 (Ω). They have also got an improved 0(h2 ) convergence order in the  · h norm, when all the rectangles are equal  a “superconvergence” result. We close this section with a brief description of other types of finite elements used in nonconforming methods.
126
11. NonConforming Methods for the Plate Problem
Example 11.2. The Zienkiewicz triangle (cf. Exercise 4.6) cf. Fig. 11.4.
Figure 11.4:
We get Vh ⊂ C 0 (Ω) only and hence the method is nonconforming. It does not always yield convergence. The method works if the sides of all triangles are parallel to three directions only, as in Fig. 11.5.
Figure 11.5:
124
This is not so if the number of directions is four, as in Fig. 11.6. (The Union Jack Problem).
11. NonConforming Methods for the Plate Problem
127
Figure 11.6: Example 11.3. Morley’s Triangle (cf. Fig. 11.7).
Figure 11.7: Here PK = P2 . We always get convergence for regular families, of course. In fact if u ∈ H 4 (Ω), then (11.34)
u − uh h ≤ Chu4,Ω .
What is astonishing is that this finite element is not even of class C 0 . Example 11.4. Fraeijs de Veubeke triangle. This finite element is again a triangle. Apart from the values of the polynomials at the vertices and the midpoints of the sides, we also take the average normal derivative along the sides. Here the space PK , which we will not describe, satisfies the inclusion P2 ⊂ P K ⊂ P3 ,
11. NonConforming Methods for the Plate Problem
128 and X K
Z ∂p = dγ, 1 ≤ i ≤ 3 p(a ), 1 ≤ i ≤ 3; p(a ), 1 ≤ i < j ≤ 3; . i ij ∂ν K′ i
125
The finite element is shown symbolically in Fig. 11.8.
Figure 11.8: Here also the finite element is not of class C 0 in general, but the method always yields convergence. References: For general reference on nonconforming methods, see Strang and Fix [22], for the bilinear lemma, see Ciarlet [29]. For a detailed study of the Zienkiewicz triangle, Moreley’s triangle and Fraeijs de Veubeke triangle, see Lascaux and Lesaint [15]. For a nonconforming method with penalty, see Babuska and Zl´amal [26].
Bibliography ´ V. [1] BRAMBLE, J.H: THOMEE, Interior maximum norm estimates for some simple finite element methods. Revue FranC ¸ aise d’Automatique, Informatique, Recherche Operationelle (RAIRO), R2, (1974), pp.518. ´ [2] BRAMBLE, J.H; ZLAMAL, M. Triangular elements in the finite element methods, Mathematics of Computation, 24, (1970), pp. 809820. [3] BREZIS, H; STAMPACCHIA, G. Sur la r´egularit´e de la solution d’in´equations elliptiques, Bull. Soc. Math. France, 96(1968), pp. 153180. [4] CIARLET, P.G. Sur l’´el´ement de Clough et Tocher, RAIRO, R2, (1974), pp. 1927. [5] CIARLET, P.G; RAVIART, P.A. ´ ements Finis pour les Probl`ems aux Limites La M´ethodes des El´ Elliptiques, (to appear). [6] CIARLET, P.G; RAVIART, P.A. Maximum principle and uniform convergence for the finite element method, Computer Methods in Applied Mechanics and Engineering, 2, (1973), pp. 1731. 129
126
BIBLIOGRAPHY
130 [7] CIARLET, P.G; RAVIART, P.A.
General Lagrange and Hermite interpolation in Rn with applications to finite element methods, Arch. Rational Mech. Anal. 46, (1972), pp. 177199. [8] CIARLET, P.G; WAGSCHAL, C. Multipoint Taylor formulas and applications to the finite element method, Numer. Math., 17, (1971), pp. 84100. [9] CIAVALDINI, J.F; NEDELEC, J.C. Sur l’´el´ement de Fraeijs de Veubeke et Sander, RAIRO, R2, (1974), pp. 2945. [10] DUVAUT, G; LIONS, J.L. Les In´equations en M´ecanique et en Physique, Dunod, Paris (1972). [11] FALK, R.S. Error estimates for the approximation of a class of variational inequalities, Math. Comput. (to appear). 127
[12] FALK, R.S. Approximation of an elliptic boundary value problem with unilateral constraints, RAIRO, (to appear). [13] HABER, S. Numerical evaluation of multiple integrals, SIAM Review, 12, (1970), pp. 481526. [14] LANDAU, L; LIFCHITZ, E. Th´eorie de l’Elasticit´e, Mir, Moscow (1967). [15] LASCAUX, P; LESAINT, P. Some nonconforming finite elements for the plate bending problem, RAIRO, (to appear).
BIBLIOGRAPHY
131
[16] LEWY, H; STAMPACCHIA, G. On the regularity of the solution of a variational inequality, Comm. Pure Appl. Math., 22, (1969), pp. 153188. [17] LIONS, J.L; MAGENES, E. Probl`emes aux Limites non Homogenes et Applications. Vol. 1, Dunod, Paris, (1968). [18] LIONS, J.L; STAMPACCHIA, G. Variational inequalities, Comm. Pure Appl. Math., 20, (1967). pp. 493519. [19] MOSCO, U; STRANG, G. Onesided approximation and variational inequalities, Bull. Amer. Math. Soc. 80 (1974), pp. 308312. ˇ J. [20] NECAS, Les M´ethodes Directes en Theorie des Equations Elliptiques, Masson, Paris, (1967). [21] STRANG, G. Approximation in the finite element method, Number. Math., 19(1972), pp. 8198. [22] STRANG, G; FIX, G.J. An Analysis of the Finite Element Method, PrenticeHall, Englewood Cliffs, (1973). ˇ ˇ [23] ZENI SEK, A. Interpolation polynomials on the triangle, Numer. Math., 15, (1970), pp. 283296. ˇ ˇ [24] ZENI SEK, A. Polynomial approximation on tetrahedrons in the finite element method, J. Approximation Theory, 7, (1973), pp. 334351.
BIBLIOGRAPHY
132 ´ [25] ZLAMAL, M.
On the finite element method, Numer. Math., 12, (1968), pp. 394409. ˇ ´ [26] BABUSKA, I.; ZLAMAL, M. Nonconforming elements in the finite element method with penalty, SIAM J. Numer. Anal., 10(1973), pp. 863875. [27] BRAMBLE, J.H.; HILBERT, S.R. Estimation of linear functionals on Sobolev spaces with application to Fourier transforms and spline interpolation, SIAM J. Numer. Anal., 7, (1970), pp. 113124. [28] BRAMBLE, J.H.; HILBERT, S.R. Bounds for a class of linear functionals with applications to Hermite interpolation, Numer. Math., 16, (1971), pp. 362369. [29] CIARLET, P.G. Conforming and nonconforming finite element methods for solving the plate problem, in Conference on the Numerical Solution of Differential Equations (G.A. WATSON, Editor), pp. 2131, Lecture Notes in Mathematics, Vol. 363, SpringerVerlag, Berlin (1974). [30] CIARLET, P.G.; RAVIART, P.A. Interpolation theory over curved elements, with applications to finite element methods, Comp. Methods in Appl. Mech. and Engineering, 1, (1972), pp. 217249. [31] CIARLET, P.G.; RAVIART, P.A. The combined effect of curved boundaries and numerical integration in isoparametric finite element methods, in The Mathematical Foundations of the Finite Element Method with Applications to Partial Differential Equations (A.K. AZIZ, Editor), pp. 409474, Academic Press, New York, (1972).
128
BIBLIOGRAPHY
133
´ [32] ZLAMAL, M. A finite element procedure of the second order of accuracy, Numer. Math., 16, (1970), pp. 394402.
BIBLIOGRAPHY
134 129
Additional References Books (Additional References) BABUSKA, I; AZIZ, A.K. First part of ‘The Mathematical Foundations of the Finite Element Method with applications to Partial Differential Equatios’ (Symposium, University of Maryland, Baltimore, June2630, 1972) (Edited by A.K. AZIZ) Academic Press, New York, (1972). NORRIE, D.H; DE VRIES, G. The Finite Element Method  Fundamentals and Applications, Academic Press, New York, (1973). ODEN, J.T. Finite Elements of Nonlinear Continue, McGraw Hill, New York, (1972). ZIENKIEWICZ, O.C. The finite Element Method in Engineering Science, McGraw Hill, London, (1971).
130
Journals (where articles dealing with the finite element method may be regularly found). Revue Franc¸aise d’Automatique, Informatique, Recherche Operationelle (RAIRO), “Red” Series Computer Methods in Applied Mechanics and Engineering. International Journal for Numerical Methods in Engineering. Mathematics of Computation. SIAM Journal on Numerical Analysis Numerische Mathematik. Proceedings of Conferences totally or partially devoted to the finite element method: The Mathematics of Finite Elements and Applications, (Conference, Brunel University, 1820 April, 1972), Edited by J.R. WHITEMAN, Academic Press, London, 1973. Conference on the Numerical Solution of Differential Equations, (Conference, Dundee, 0306 July 1973), Edited by G.A. WATSON, Lecture Notes in Mathematics, Vol. 363, SpringerVerlag, New York, 1974. The Mathematical Foundations of the Finite Element Method with
BIBLIOGRAPHY
135
Applications to Partial Differential Equations—, (Symposium, University of Maryland, Baltimore, 2630 June. 1972), Edited by A.K. AZIA, Academic Press. New York, 1972. Computing Methods in Applied Science and Engineering, Vol. 1 and 2 (International Symposium IRIA LABORIA, Versailles, 1721 December, 1973) Edited by R. GLOWINSKI and J.L. LIONS, Lecture notes in Computer Science, Vol., 10 and 11, SpringerVerlag, New York, 1974. Topics in Numerical Analysis (Royal Irish Academy Conference on Numerical Analysis, Dublin, 1972), Edited by J.J.H. MILLER, Academic Press, New York, 1973. Mathematical Aspects of Finite Elements in Partial Differential Equations (Symposium, 0103 April, 1974, Mathematics Research Center, The University of Wisconsin, Madison), Edited by C. de BOOR, Academic Press, New York, 1974. Proceedings of the Second Conference on the Mathematics of Finite Elements and Applications, (Brunel University, April 0710, 1975, to appear.)