- Email: [email protected]
A posteriori error estimates for fourth-order elliptic problems
Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559 www.elsevier.com/locate/cma
A posteriori error estimates for fourth-order elliptic problems Slimane Adjerid Department of Mathematics and Interdisciplinary Center for Applied Mathematics, Virginia Polytechnic Institute and State University, Blacksburg, VA 24061, USA Received 20 February 2001; accepted 19 November 2001
Abstract We extend the dichotomy principle of Babuska and Yu [Math. Numerica Sinica 13 (1991) 89; Math. Numerica Sinica 13 (1991) 307] and Adjerid et al. [Math. Mod. Meth. Appl. S. 9 (1999) 261; SIAM J. Sci. Comput. 21 (1999) 728] for estimating the finite element discretization error to fourth-order elliptic problems. We show how to construct a posteriori error estimates from jumps of the third partial derivatives of the finite element solution when the finite element space consists of piecewise polynomials of odd-degree and from the interior residuals for even-degree approximations on meshes of square elements. These estimates are shown to converge to the true error under mesh refinement. We also show that these a posteriori error estimates are asymptotically correct for more general finite element spaces. Computational results from several examples show that the error estimates are accurate and efficient on rectangular meshes. 2002 Elsevier Science B.V. All rights reserved.
1. Introduction A posteriori error estimators and/or indicators are an integral part of adaptive finite element methods for partial differential equations. They are used for second-order elliptic [9–11,26–28] and parabolic [3,4,7,8] problems. More recently a posteriori error estimates have been developed for hyperbolic problems [5,6,24]. Local error estimates are used to control adaptive enrichment by refining regions having large errors and coarsening regions with small errors. On the other hand, global a posteriori error estimates are used to control solution accuracy. Ideal a posteriori error estimates should be computationally simple and provide guaranteed asymptotical exactness in the sense that the estimates should converge to the true errors from above as the finite element space is enriched. Even–odd error estimates have been used successfully by Yu [27,28] for elliptic problems and Adjerid et al. [3,4] for parabolic problems. They are asymptotically correct for tensor product [3,27,28] and hierarchical [4] finite element spaces under mesh refinement. Most of the existing error estimators are applicable to second-order problems. However, many important physical situations, such as elasticity [23], fluid mechanics [18] and thin film growth [21,29] are modeled by boundary-value problems of fourth order. The first a posteriori error estimate for fourthorder problems appeared in Verfurth [26]. The author is not aware of any other attempt to compute efficient a posteriori error estimates for higher-order problems. In this paper we show how to extend the even–odd error estimation procedures to linear fourth-order boundary value problems. We use C 1 finite 0045-7825/02/$ - see front matter 2002 Elsevier Science B.V. All rights reserved. PII: S 0 0 4 5 - 7 8 2 5 ( 0 1 ) 0 0 4 1 2 - 1
2540
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
element approximations on a square mesh and compute error estimates using jumps of the third partial derivatives of odd-degree finite element solutions and by solving local residual problems for even-degree solutions. The even–odd error estimates are both very simple and efficient to compute. Furthermore, the evendegree error estimators do not involve any neighbor information, which makes them ideal for parallel computations. Moreover, numerical results suggest that even-degree estimators perform better than odddegree estimators. Both even- and odd-degree estimators are extended to more general hierarchical finite element spaces. More recent approaches to a posteriori error estimation include guaranteed estimates [12,13], weighted and goal-oriented estimates [15,20] and bound estimates for functionals [22]. Here we will show that our estimates are asymptotically correct in the H 2 norm. In Section 3, we prove several lemmas. We establish asymptotic exactness of these error estimates on square elements for both bi p even-degree Section 4 and odd-degree Section 5. These error estimates are further extended to the more economical hierarchical approximations in Section 6. We test our error estimators on three examples in Section 7 and present numerical results which show that the actual theory holds for square- and rectangular-element meshes. Finally, numerical results [4] suggest that the even–odd error estimators are applicable to more general meshes as well as problems with singularities and we expect this to continue to hold for fourth-order problems.
2. Problem formulation Let us consider the linear scalar fourth-order partial differential equation x ¼ ðx1 ; x2 Þ 2 X ¼ ½0; 12 ;
Lu ¼ f ; with Lu ¼
X X
ð1Þjbj Db ðaa;b ðxÞDa uÞ;
ð2:1aÞ
ð2:1bÞ
jaj 6 2 jbj 6 2
where jaj ¼ a1 þ a2 , jbj ¼ b1 þ b2 and a2 Da ¼ oa1 x1 ox2 ;
ð2:1cÞ
subject to the Dirichlet boundary conditions uðxÞ ¼ 0;
ou ðxÞ ¼ 0; ov
x 2 oX;
ð2:1dÞ
where ou=ov is the outward normal derivative of u. The functions aa;b ðxÞ, jaj, jbj 6 2, are smooth functions and L is a positive definite operator on the Sobolev space H02 . The Galerkin form of (2.1a)–(2.1d) consists of determining u 2 H02 satisfying Aðv; uÞ ¼ ðv; f Þ;
8 v 2 H02 ;
ð2:2aÞ
2
where the strain energy and L inner products, respectively, are # Z Z "X X Aðv; uÞ ¼ aa;b ; Da uDb v dx1 dx2 X
and ðv; uÞ ¼
Z Z uv dx1 dx2 : X
ð2:2bÞ
jaj 6 2 jbj 6 2
ð2:2cÞ
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
Sobolev spaces H s , s P 0 are equipped with the usual inner product, norm and seminorm X ðDa v; Da uÞ; ðv; uÞs ¼
2541
ð2:3aÞ
jaj 6 s 2
kuks ¼ ðu; uÞs ; X 2 jujs ¼ ðDa u; Da uÞ;
ð2:3bÞ ð2:3cÞ
jaj¼s
where a 0 subscript on H 2 implies that functions also satisfy (2.1d). Finite element solutions of (2.2a) are obtained by approximating H02 by a finite dimensional subspace S0N ;p and determining U 2 S0N ;p such that AðV ; U Þ ¼ ðV ; f Þ;
8 V 2 S0N ;p :
ð2:4Þ
Let us partition X into a uniform mesh of square elements Di , i ¼ 1; 2; . . . ; N and define S N ;p as S N ;p ¼ fw 2 H 2 j wðxÞ 2 Qp ðDi Þ; x 2 Di ; i ¼ 1; 2; . . . ; N g;
ð2:5Þ
where Qp ðDi Þ is the space of bi-polynomial functions that are products of univariate polynomials of degree p in x1 and x2 on Di . Again with the subscript 0 S0N ;p ¼ S N ;p \ H02 . In the following lemmas we describe standard interpolation and a priori discretization error estimates for conforming finite element solutions of (2.1a)–(2.1d) that will be used in the subsequent analysis. Lemma 1. Let W 2 S0N ;p be an interpolant of w 2 H02 \ H pþ1 that is exact when w 2 Qp ðXÞ. Then, there exists a constant C > 0 such that kW wks 6 Chpþ1s kwkpþ1 ;
s ¼ 0; 1; 2;
ð2:6aÞ
where
pffiffiffiffi h ¼ 1= N : H02
ð2:6bÞ pþ1
Further, let u 2 \H and U 2 S a constant C > 0 such that ku U ks 6 Chpþ1s kukpþ1 ;
Proof. Cf. Oden and Carey [19].
N ;p
be solutions of (2.1a)–(2.1d) and (2.4), respectively. Then, there exists
s ¼ 0; 1; 2:
ð2:6cÞ
We also recall the trace theorem. Lemma 2. Let X be a bounded open set of Rn with a Lipschitz boundary oX. Then there exists a constant, C > 0, such that kvk0;oX 6 Cðh1=2 kvk0;X þ h1=2 kvk1;X Þ;
Proof. Cf. Grisvard [17].
8 v 2 H 1:
ð2:7Þ
b 2 Pp ð½1; 12 Þ. Then there exist conLemma 3. Let W 2 Pp ðDÞ and its mapping to the canonical element W stants C1 ; C2 > 0 such that
2542
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559 2 8 W^ 2 Pp ð½1; 1 Þ; s P 0
C1 kW^ ksþ1 < kW^ ks < C2 kW^ ksþ1 ;
ð2:8Þ
and C1 kW ksþ1 < hkW ks < C2 kW ksþ1 ;
8 W 2 Pp ðDÞ; s P 0;
ð2:9Þ
where Pp ðDÞ is the space of polynomials of degree p on an element D:
Proof. Eq. (2.8) is obtained from the equivalence of all norms on finite dimensional spaces. Similarly, 2 Eq. (2.9) is straightforward using the mapping of the square element ½z h=2; z þ h=2 to the canonical 2 element ½1; 1 . 2
Let us construct a hierarchical basis for S N ;p on the canonical element ½1; 1 as a tensor product of the one-dimensional hierarchical basis consisting of the four standard cubic Hermite polynomial shape functions Nlm , l ¼ 1; 1; m ¼ 0; 1, complemented by the interior shape functions N0i , i P 4, such that 0 dN1 ð1Þ ¼ 0; dx
0 N1 ð1Þ ¼ 1;
0 N1 ð1Þ ¼ 0;
1 N1 ð1Þ ¼ 0;
1 dN1 ð1Þ ¼ 1; dx
1 dN1 ð1Þ ¼ 0; dx
dN10 ð1Þ ¼ 0; dx
N10 ð1Þ ¼ 0;
N10 ð1Þ ¼ 1;
N11 ð1Þ ¼ 0;
dN11 ð1Þ ¼ 0; dx
dN11 ð1Þ ¼ 1: dx
ð2:10aÞ ð2:10bÞ ð2:10cÞ ð2:10dÞ
Higher-order interior shape functions are obtained by integrating Lobatto polynomials on ½1; 1 to obtain Z n N0i ðnÞ ¼ Pi1 ðxÞ dx; i ¼ 4; 5; . . . ; p; ð2:11Þ 1
where Pi is the Lobatto polynomial of degree i Pi ðxÞ Pi2 ðxÞ Pi ðxÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; 2ð2i 1Þ
ð2:12Þ
with Pi being the Legendre polynomial of degree i. Using the orthogonality of Legendre polynomials yields 1 Piþ1 Pi1 ð2i 3ÞPiþ1 2ð2i 1ÞPi1 þ ð2i þ 1ÞPi3 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi N0iþ1 ¼ pffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi ¼ : ð2:13Þ 2i þ 1 2i 1 2i 3 2ð2i 1Þð4i2 4i 3Þ One can easily verify that the shape functions N0i , i P 4, satisfy R1 (i) 1 ðd2 N0i =dn2 Þðd2 N0j =dn2 Þ dn ¼ dij , i; j P 4, where dij is the Kronecker symbol. (ii) N0i has i 4 simple roots in ð1; 1Þ and two double roots at 1. 2 (iii) N0i ¼ ðn2 1Þ N~ i4 ðnÞ, where N~ 2l has 2l non-zero roots in ð1; 1Þ symmetrically disposed with respect to the origin, fnk ; k ¼ 1; 2; . . . ; lg. On the other-hand N~ 2lþ1 has 2l+1 roots in ð1; 1Þ, f0g [ fnk ; k ¼ 1; 2; . . . ; lg.
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
2543
The two-dimensional basis on ½1; 1 ½1; 1 is defined by • 16 vertex modes: Nijk;m ðn; gÞ ¼ Nik ðnÞNjm ðgÞ;
i; j ¼ 1; 1; k; m ¼ 0; 1:
ð2:14aÞ
• 16ðp 3Þ edges modes: j;n Ni;1 ¼ Nij ðnÞN0n ðgÞ;
i ¼ 1; 1; j ¼ 0; 1; n ¼ 4; . . . ; p;
ð2:14bÞ
j;n ¼ Nij ðgÞN0n ðnÞ; Ni;2
i ¼ 1; 1; j ¼ 0; 1; n ¼ 4; . . . ; p;
ð2:14cÞ
j;n ¼ Nij ðnÞN0n ðgÞ; Ni;3
i ¼ 1; 1; j ¼ 0; 1; n ¼ 4; . . . ; p;
ð2:14dÞ
j;n ¼ Nij ðnÞN0n ðgÞ; Ni;4
i ¼ 1; 1; j ¼ 0; 1; n ¼ 4; . . . ; p;
ð2:14eÞ
j; k ¼ 4; . . . ; p:
ð2:14fÞ
2
• ðp 3Þ interior modes: N0j;k ¼ N0j ðnÞN0k ðgÞ; 3. Preliminary considerations Let p~½z be the univariate operator that interpolates functions in H02 ðz h=2; z þ h=2Þ and their first derivatives such that p~uðz h=2Þ ¼ uðz h=2Þ;
ð3:1aÞ
d~ pu du ðz h=2Þ ¼ ðz h=2Þ; dx dx
ð3:1bÞ
p~uðz nj h=2Þ ¼ uðz nj h=2Þ;
j ¼ 1; 2; . . . ; s; for p 3 ¼ 2s
ð3:1cÞ
p~uðz nj h=2Þ ¼ uðz nj h=2Þ;
j ¼ 0; 1; 2; . . . ; s; for ðp 3Þ ¼ 2s þ 1;
ð3:1dÞ
or where nj are the simple roots of wpþ1 ðz; zÞ ¼ z
pþ1
N0pþ1 .
Thus, for p P 3, the function
pþ1
p~½zz
ð3:2Þ
N0pþ1
have identical roots with first derivatives which both vanish at z h=2. Hence, wpþ1 ðz; zÞ, and w0pþ1 ðz; zÞ, and w00pþ1 ðz; zÞ are, respectively, proportional to N0pþ1 , Lobatto and Legendre polynomials on 0 ½z h=2; z þ h=2, with ð Þ denoting ordinary differentiation. We use p~ to define a two-dimensional interpolation operator pi on element i satisfying pi uðxÞ ¼ p~½x1;i ~ p½x2;i uðxÞ 2 Qp ðDi Þ, x 2 Di . Since the mesh is uniform, we will omit the elemental index i and the dependence of p~½xj;i and wðxj;i ; xi Þ on the coordinates of the cell center ðx1;i ; x2;i Þ, i ¼ 1; 2; . . . ; N , whenever confusion is unlikely. The functions wpþ1 ðxj Þ, j ¼ 1; 2, provide a basis for the dominant contribution to the spatial discretization error on element D for both odd- and even-order finite element approximations. Indeed, we shall show that estimates EðxÞ of eðxÞ have the form EðxÞ ¼ b1 wpþ1 ðx1 Þ þ b2 wpþ1 ðx2 Þ;
x 2 D:
In order to estimate the finite element discretization error, we need to prove the following series of lemmas.
2544
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
Lemma 4. Let pu be an interpolant of u 2 H pþ2 , then uðxÞ puðxÞ ¼ /ðxÞ þ cðxÞ;
x 2 D;
ð3:3aÞ
where /ðxÞ ¼ b1 wpþ1 ðx1 Þ þ b2 wpþ1 ðx2 Þ;
ð3:3bÞ
kDa /k0;D 6 Chpþ1jaj kukpþ1;D ;
ð3:3cÞ
kcks;D 6 Chpþ2s kukpþ2;D ;
jaj 6 2;
s ¼ 0; 1; . . . ; p þ 1;
ð3:3dÞ
and koxj cks;D 6 Chpþ1s kukpþ2;D ; Z D
o2xj wpþ1 ðxj Þxkj dx1 dx2 ¼ 0;
s ¼ 0; 1; j ¼ 1; 2; k ¼ 0; 1; . . . ; p 2; p P 3:
ð3:3eÞ ð3:3fÞ
Remark. Local Sobolev norms are defined like their global counterparts (2.3b) and (2.3c) with X replaced by D. Proof. Let pu 2 Qp and V ¼ w þ b1 xpþ1 þ b2 x2pþ1 such that w 2 Qp , pV ¼ pu and V interpolates u at two 1 other points. Then the interpolation error can be written as u pu ¼ u V þ V pu: ð3:4Þ Using the fact that pu ¼ pV we have pþ1 p~½x2 x2pþ1 Þ: u pu ¼ u V þ V pV ¼ u V þ b1 ðx1pþ1 p~½x1 xpþ1 1 Þ þ b2 ðx2
ð3:5Þ
u pu ¼ c þ /;
ð3:6Þ
Thus where c ¼ u V . Using standard interpolation error estimates (2.6a) we obtain (3.3c)–(3.3e). Finally, we note that wpþ1 and N0pþ1 are polynomials of degree p þ 1 with identical roots and therefore are proportional, wpþ1 / N0pþ1 , p P 3. This leads to the orthogonality result (3.3f). Lemma 5. Let Pu 2 S0N ;p be an interpolant of u 2 H pþ2 for x 2 X that agrees with pu when x 2 D, then jAðW ; u PuÞj 6 Chp kukpþ2 kW k2 ;
8 W 2 S0N ;p :
ð3:7Þ
Proof. Let us consider first Z að2;0Þ;ð2;0Þ ðxÞ o2x1 W o2x1 ðu puÞ dx1 dx2 D Z Z ¼ að2;0Þ;ð2;0Þ o2x1 W o2x1 ðu puÞ dx1 dx2 þ að2;0Þ;ð2;0Þ ðxÞ að2;0Þ;ð2;0Þ o2x1 W o2x1 ðu puÞ dx1 dx2 ; ð3:8Þ D
D
being the center of the element D. xÞ, with x where að2;0Þ;ð2;0Þ ¼ að2;0Þ;ð2;0Þ ð Using the continuity of aa;b and applying the Schwarz inequality with the estimates (2.6a) we have
Z
2 2
að2;0Þ;ð2;0Þ ox1 W ox1 ðu puÞ dx1 dx2
< Chku Puk2;D kW k2;D < Chp kukpþ2;D kW k2;D : að2;0Þ;ð2;0Þ ðxÞ
D
ð3:9Þ
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
2545
Using (3.3a) and (3.3b), the first term in the right-hand side of (3.8) becomes Z Z o2x1 W o2x1 ðu puÞ dx1 dx2 ¼ o2x1 W ðo2x1 / þ o2x1 cÞ dx1 dx2 :
ð3:10Þ
Using (3.3b) and the fact that w00pþ1 is proportional to the Legendre polynomial Pp1 we have Z o2x1 W o2x1 / dx1 dx2 ¼ 0; 8 W 2 S N ;p :
ð3:11Þ
D
D
D
Applying the Schwarz inequality and using (3.3e) leads to
Z
o2 W o2 c dx1 dx2 6 Chp kuk kW k ; 8 W 2 S N ;p : pþ2 2 x1 x1
ð3:12Þ
D
Following the same steps we prove jðo2x2 W ; o2x2 ðu puÞÞj 6 Chp kukpþ2;D kW k2D ;
8 W 2 S N ;p :
Applying Green’s theorem we obtain Z Z Z o2x2 W o2x1 ðu puÞ dx1 dx2 ¼ ox2 W o2x1 ðu puÞdx1 C 1
D
Z
ð3:13Þ
Cþ 1
ox2 W o2x1 ðu puÞ dx1
ox2 W o2x1 ox2 ðu puÞ dx1 dx2 :
ð3:14Þ
Using (3.3a) we have Z Z 2 ox2 W ox1 ox2 ðu puÞ dx1 dx2 ¼ ox2 W o2x1 ox2 ð/ þ cÞ dx1 dx2 :
ð3:15Þ
Using the fact that o2x1 ox2 / ¼ 0 we have Z Z ox2 W o2x1 ox2 ðu puÞ dx1 dx2 ¼ ox2 W o2x1 ox2 c dx1 dx2 :
ð3:16Þ
D
D
D
D
D
Integration by parts yields Z Z Z 2 2 2 ox2 W ox1 ox2 c dx1 dx2 ¼ ox2 W ox1 c dx1 dx2 þ D
D
Integrating the boundary terms by parts we obtain Z Z Z ox2 W o2x1 ox2 c dx1 dx2 ¼ o2x2 W o2x1 c dx1 dx2 D
D
C 1
C 1
ox2 W
o2x1 c dx1
Z Cþ 1
ox1 ox2 W ox1 c dx1 þ
o2x2 W o2x1 c dx1 :
Z Cþ 1
ox1 o2x2 W ox1 c dx1 :
Applying the Schwarz inequality and using (3.3e), we immediately have
Z
o2 W o2 c dx1 dx2 6 kW k kck 6 Chp kuk kW k : 2 2 pþ2 2 x2 x1
ð3:17Þ
ð3:18Þ
ð3:19Þ
D
The boundary terms in (3.18) are estimated by applying the Schwarz inequality with (2.7) and (3.3d) to obtain
Z
ox1 ox2 W ox1 c dx1 6 Cðh1=2 kW k2 þ h1=2 kW k3 Þðh1=2 kck1 þ h1=2 kck2 Þ
C 1
6 Chp kukpþ2;D ðkW k2;D þ hkW k3;D Þ 6 Chp kukpþ2;D kW k2;D :
ð3:20Þ
2546
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
We further integrate by parts the boundary terms in (3.14) and use (3.1a)–(3.1d) to obtain Z Z ox2 W o2x1 ðu puÞ dx1 ¼ ox1 ox2 W ox1 ðu puÞ dx1 : C 1
ð3:21Þ
C 1
Since ox1 ox2 W and ox1 ðu puÞ are continuous along the horizontal edge C shared by two neighboring elements Di and Dj , the contribution of such terms to (3.7) can be written as Z D ð aDa;bi aa;bj Þ ox1 ox2 W ox1 ðu puÞ dx1 : ð3:22Þ C
Applying the Schwarz inequality and using (2.6a), (2.7), (2.8) and the smoothness of the coefficients of L we obtain
Z
D
Dj i
ð aa;b Þ ox1 ox2 W ox1 ðu puÞ dx1
< CðkW k2 þ hkW k3 Þðku puk1 þ hku puk2 Þ
aa;b C
< CkW k2;Di \Dj ku puk1;Di \Dj < Chp kW k2;Di \Dj kukpþ2;Di \Dj : Following similar steps we prove
Z
a b p
a D ðu puÞD W dx dx a;b 1 2 < Ch kukpþ2;X kW k2;X ;
8 a; b:
ð3:23Þ
ð3:24Þ
X
Combining (3.8)–(3.24) completes the proof of (3.7).
In the following lemma we state a superconvergence result. Lemma 6. Let U and u 2 H pþ2 be solution of (2.4) and (2.2a), respectively and Pu 2 S N ;p interpolate u as described in Lemma 5, then kU Puk2 < Chp kukpþ2 :
ð3:25Þ
Proof. Subtracting AðW ; puÞ from the orthogonality condition 8 W 2 S N ;p
ð3:26Þ
AðW ; U puÞ ¼ AðW ; u puÞ:
ð3:27Þ
AðW ; U uÞ ¼ 0; leads to
Using Lemma 5 with W ¼ U pu and the ellipticity of the bilinear form A we complete the proof. In the following theorem we show that the finite element error can be divided into a leading component / plus a higher-order component h. Lemma 7. Under the conditions of Lemma 6, there exists 1 such that kek22 ¼ ku Puk22 þ 1
ð3:28aÞ
j1 j 6 CðuÞh2p1 :
ð3:28bÞ
with
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
2547
Proof. Adding and subtracting Pu to e yields (3.28a) with 2
1 ¼ 2ðu Pu; Pu U Þ2 þ kPu U k2 : Applying the Schwarz inequality and using (2.6a) and (3.25) yields (3.28b).
ð3:29Þ
Theorem 1. If u and U are as defined in Lemma 6, then e ¼ u U ¼ U þ H;
ð3:30aÞ
where the restrictions of U and H to D are / of (3.3b) and h ¼ c þ pu U :
ð3:30bÞ
Furthermore, kUk22 6 Ch2p2 kuk2pþ1 ;
ð3:30cÞ
jHk22 6 Ch2p kuk2pþ1 :
ð3:30dÞ
Proof. Adding and subtracting pu to e and using (3.3a)–(3.3d) and (3.25) we prove the theorem.
4. Even-degree a posteriori error estimation Error estimates of even-degree approximations are computed using interior residuals by solving local Petrov–Galerkin problems for the error. In order to show this we substitute u¼eþU
ð4:1Þ
in (2.2a) to obtain Aðv; eÞ ¼ gðvÞ;
8 v 2 H02 ;
ð4:2Þ
where gðvÞ ¼ ðv; f Þ Aðv; U Þ:
ð4:3Þ
The error is approximated by E ¼ b1 wpþ1 ðx1 Þ þ b2 wpþ1 ðx2 Þ;
ð4:4Þ
on element D and the test function v is selected as vj ðxÞ ¼ ðxj xj;i Þdðx1 Þdðx2 Þ;
ð4:5Þ
where dðzÞ ¼
wpþ1 ðz; zÞ : z z
ð4:6Þ
Since wpþ1 ðxj Þ and vj ðxÞ; j ¼ 1; 2 vanish on oD, the error estimate satisfies the local Dirichlet problems AD ðvj ; EÞ ¼ gD ðvj Þ;
j ¼ 1; 2;
ð4:7Þ
where the subscript D denotes that inner products are restricted to D. Using (3.30a)–(3.30d), the finite element problem (4.7) becomes AD ðv; /Þ ¼ gD ðvÞ AD ðv; hÞ;
8 v 2 H02 :
ð4:8Þ
2548
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
Which in turn can be written as D ðv; /Þ ¼ gD ðvÞ AD ðv; hÞ FD ðv; /Þ; A where D ðv; /Þ ¼ A
Z Z "X X D
and FD ðv; /Þ ¼
ð4:9Þ
# a
b
aa;b ð xÞD /D v dx1 dx2
ð4:10Þ
jaj¼2 jbj¼2
Z Z "X X D
þ
8 v 2 H02 ;
# a
ðaa;b ðxÞ aa;b ð xÞÞD /D v dx1 dx2
jaj¼2 jbj¼2
Z Z "X X D
b
# a
b
aa;b ðxÞD /D v dx1 dx2 :
ð4:11Þ
jaj<2 jbj<2
Using symmetry properties of wpþ1 ðzÞ and vj ðxÞ we obtain the uncoupled constant-coefficient problem on element D Z Z D dðx2 Þðw00pþ1 ðx1 ÞÞ2 dx1 dx2 ¼ gðv1 Þ Aðv1 ; hÞ F ðv1 ; /Þ; ð4:12aÞ að2;0Þ;ð2;0Þ b1 D
aDð0;2Þ;ð0;2Þ b2
Z Z
2
dðx1 Þðw00pþ1 ðx2 ÞÞ dx1 dx2 ¼ gðv2 Þ Aðv2 ; hÞ F ðv2 ; /Þ:
ð4:12bÞ
D
The finite element problem (4.12a) and (4.12b) may be further simplified by neglecting the higher order terms to obtain the following system for bj , j ¼ 1; 2, Z Z 2 b1 dðx2 Þðw00pþ1 ðx1 ÞÞ dx1 dx2 ¼ gðv1 Þ; ð4:13aÞ aDð2;0Þ;ð2;0Þ D
aDð0;2Þ;ð0;2Þ b2
Z Z
dðx1 Þðw00pþ1 ðx2 ÞÞ2 dx1 dx2 ¼ gðv2 Þ:
ð4:13bÞ
D
The error estimate E in (4.4) with (4.13a) and (4.13b) is shown to be asymptotically correct in the next theorem. However, prior to this we need the following properties of wpþ1 , d and vj , j ¼ 1; 2: Lemma 8. Let p > 3 be an even integer, then there exist C > 0 and c > 0 such that 2
ch2ðpsÞþ4 < kwpþ1 ks;D < Ch2ðpsÞþ4 ; Z Z
s ¼ 0; 1; 2;
dðxj Þdx1 dx2 ¼ Chpþ2 ;
ð4:14aÞ ð4:14bÞ
D
Z Z
2
dðxk Þw00 ðxj Þ dx1 dx2 ¼ Ch3p ;
j; k ¼ 1; 2;
ð4:14cÞ
D 2
ch4p < kvj k2;D < Ch4p : Proof. A direct computation reveals the results.
ð4:14dÞ
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
2549
In the following lemma we establish bounds for the difference between bj and its approximation bj , j ¼ 1; 2. Lemma 9. Let p > 3 be an even integer and u 2 H02 \ H pþ2 . If the coefficients of L in (2.1b) are smooth functions then N X
ðb2ji b2ji Þ 6
i¼1
C kuk2pþ2 ; h
j ¼ 1; 2:
ð4:15Þ
Proof. Subtracting (4.13a) and (4.13b) from (4.12a) and (4.12b) we obtain 2
ðbji bji Þ2 ¼
ðAðvj ; hÞ þ F ðvj ; /ÞÞ : RR 2 2 dðx2 Þw00pþ1 ðx1 Þ dx1 dx2 Þ Di
i ð aDð2;0Þ;ð2;0Þ
ð4:16Þ
Applying the Schwarz inequality and using (4.14a)–(4.14d), (3.30d) and the continuity of aa;b ðxÞ we have jAðvj ; hÞj2 6 Ckvj k22;Di khk22;Di 6 Ch6p kuk2pþ2;Di :
ð4:17Þ
Using the continuity of aa;b (x) with (3.3c) we also prove jF ðvj ; /Þj2 6 Ch2 kvj k22;Di k/k22;Di 6 Ch6p kuk2pþ2;Di :
ð4:18Þ
Combining (4.16)–(4.18) with (4.14a)–(4.14d) yields ðbji bji Þ2 < Ckuk2pþ2;Di :
ð4:19Þ
A summation over all elements leads to N X
2
2
ðbji bji Þ < Ckukpþ2 :
ð4:20Þ
i¼1
Using (3.3b) and (4.14a)–(4.14d) we obtain 2
b2ji ¼
ko2xj /k0;Di kw00pþ1 ðxj Þk20;Di
<
C kuk2pþ2;Di : h2
ð4:21Þ
We may solve (4.13a) and (4.13b) for bj to obtain 2
b2ji ¼
i ð aDð2;0Þ;ð2;0Þ
gðvj Þ
RR
2
Di
dðx2 Þw00pþ1 ðx1 Þ dx1 dx2 Þ
2
:
ð4:22Þ
If we apply the Schwarz inequality, use (4.2) and (4.14a)–(4.14d), sum over all elements and use (2.6c) we obtain N X
2
b2ji 6 Ch2p kek2 6
i¼1
C 2 kukpþ2 : h2
ð4:23Þ
Combining (4.21) and (4.23) N X i¼1
½b2ji þ b2ji 6
C 2 kukpþ2 : h2
ð4:24Þ
2550
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
Applying the Schwarz inequality #1=2 " #1=2 " N N N X X X 2 2 2 2 2 ½bji bji 6 C ðbji bji Þ ðbji þ bji Þ i¼1
i¼1
ð4:25Þ
i¼1
while using (4.20) and (4.24) leads to (4.15).
Now we will state and prove the main result of this section. Theorem 2. Let u 2 H02 \ H pþ2 and U 2 S N ;p be solutions of (2.2a) and (2.4), respectively. If p > 3 is an even integer, then kek22 ¼
N X 2 X i¼1
b2ji kw00pþ1 ðxj Þk20;Di þ ;
ð4:26aÞ
j¼1
where 2
jj 6 Ch2p1 kukpþ2 :
ð4:26bÞ
Proof. Consider kek22 ¼
N h i X ko2x1 ek20;Di þ ko2x2 ek20;Di þ kox1 ox2 ek20;Di þ kek21
ð4:27Þ
i¼1
and use (3.30a)–(3.30d) and (3.3b) to obtain (Z Z " ! N 2 X X 2 2 2 00 kek2 ¼ bji jwpþ1 ðxj Þj þ jo2x1 hj2 þ 2o2x1 / o2x1 h þ jo2x2 hj2 Di
i¼1
þ
j¼1
2o2x2 / o2x2 h
# 2
þ jox1 ox2 hj
) dx1 dx2
2
þ kek1 :
ð4:28Þ
2
Adding and subtracting b2ji jw00pþ1 ðxj Þj , to the above integrand we obtain (4.26a) with (Z Z " ! N 2 X X 2 2 00 2 ðbji bji Þjwpþ1 ðxj Þj þ jo2x1 hj2 þ 2o2x1 / o2x1 h þ jo2x2 hj2 ¼ i¼1
Di
j¼1
#
þ 2o2x2 / o2x2 h þ jox1 ox2 hj2 dx1 dx2
) þ kek21 :
ð4:29Þ
Applying the Schwarz and triangle inequalities and using the estimates (2.6c), (3.3c), (3.3d), (4.14a)–(4.14d) and (4.15) yields (4.26b). 5. Odd-degree a posteriori estimator If u is a smooth solution of (2.1a)–(2.1d) on D, then error estimates E of odd-degree approximations may be computed in terms of jumps in derivatives o3xj U , j ¼ 1; 2 at the vertices of D. Using e ¼ u U we compute the jumps in derivatives at pk ¼ ðp1k ; p2k Þ, k ¼ 1; 2; 3; 4, of D as
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
½o3xj eðpk Þj ¼ ½o3xj U ðpk Þj b1 ½o3xj wpþ1 ðp1k Þj þ b2 ½o3xj wpþ1 ðp2k Þj ;
2551
j ¼ 1; 2; k ¼ 1; 2; 3; 4; x 2 D; ð5:1Þ
where ½qðpÞj denotes the jump in q at point p in the xj direction. Since jumps in the finite element solution o3xj U and o3xj wpþ1 , j ¼ 1; 2, are known, (5.1) yields a linear algebraic system for b1 and b2 at each vertex pk . On a rectangular mesh the set of systems (5.1) have a unique solution. More general meshes are treated in [4]. kEk22;D is the average of local errors based on the pairs of bj , j ¼ 1; 2. In the next lemma we state some useful properties of Legendre polynomials. Lemma 10. Let Pp be the Legendre polynomial of degree p, then Z 1 2 ; Pp2 ðtÞ dt ¼ 2p þ1 1 0 0 Ppþ1 ¼ ð2p þ 1ÞPp þ Pp1 ;
Pp ðþ1Þ ¼ 1;
pþ1
pðp þ 1Þ ; 2
ð5:2bÞ
p > 1;
Pp ð1Þ ¼ ð1Þp ;
Pp0 ð1Þ ¼ ð1Þ
ð5:2aÞ
ð5:2cÞ
p P 0;
ð5:2dÞ
p > 0:
Proof. See Abramowitz and Stegun [1].
Theorem 3. Let u 2 H02 \ H pþ2 , Pu be as defined in Lemma 5, and p P 3 be an odd integer, then ku Puk22 ¼
N X 2 X 4 X
h4 64p2 ðp 1Þ2 ð2p 1Þ
i¼1
j¼1
½o3xj Puðpk Þ2j þ 2 ;
ð5:3aÞ
k¼1
where 2
j2 j 6 Ch2p1 kukpþ2 :
ð5:3bÞ
Proof. The proof will be established in two steps. In the first step we assume u ¼ uj ¼ xjpþ1 , j ¼ 1; 2 and show that ko2xj ðuj Puj Þk2 ¼
h4
4 X
64p2 ðp 1Þ2 ð2p 1Þ
k¼1
½o3xj Puj 2j :
ð5:4Þ
Consider two neighboring elements Di and Dik and compute the jump in rj ðxj Þ ¼ o3xj wpþ1 ¼ w000 pþ1 ðxj Þ in the jth direction at common vertices. On a square-element mesh we have ½o3xj ðuj Puj Þj ¼ 2rj ðh=2Þ;
ð5:5Þ P4
ko2xj ðuj
2
Puj Þk ¼
2 kw00pþ1 ðxj Þk
k¼1
2
½o3xj Puj j
16ðrj ðh=2ÞÞ
2
¼ Ap
4 X k¼1
2
½o3xj Puj j ;
ð5:6aÞ
2552
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
where 2
Ap ¼
kw00pþ1 kD 16ðrj ðh=2ÞÞ2
ð5:6bÞ
:
Transforming the integrals from a physical element D to the canonical element, applying the chain rule and using the fact that w00pþ1 / Pp1 to obtain n n2 dn n d n d w ðzÞ ¼ ð2=hÞ w ðnÞ ¼ cð2=hÞ Pp1 ðnÞ; n pþ1 pþ1 dzn dn dnn2 Thus, Z 1 h4 h4 2 Ap ¼ : Pp1 ðtÞ dt ¼ 2 2 27 p2 ðp 1Þ 1 64p2 ðp 1Þ ð2p 1Þ
n ¼ 2; 3:
ð5:7Þ
ð5:8Þ
This proves (5.4). Now, we assume w to be a p-degree polynomial with respect to x1 and x2 . Thus, ½oxj ðw PwÞðpk Þ ¼ 0;
w Pw ¼ 0; b1 x1pþ1
b2 x2pþ1S
Let v ¼ w þ þ neighbors of D with XD ¼ Pv0 ¼ Pu and
j ¼ 1; 2:
ð5:9Þ
and Di , i ¼ 1; 2; 3; 4 denote, respectively, the right, top, left, and bottom Dk [ D. We define a function v such that vjD ¼ v0 and vjDk ¼ vk such that
4 k¼1
v0 ðqj Þ ¼ uðqj Þ;
ð5:10Þ
where qj , j ¼ 1; 2 are two other points in D, Pvk ¼ Pu;
for x 2 Dk ;
k ¼ 1; 2; 3; 4;
ð5:11Þ
and the following continuity conditions: 3 þ o3x1 Pv1 ðpþ k Þ ¼ ox1 Pv0 ðpk Þ;
k ¼ 1; 2;
ð5:12aÞ
3 þ o3x2 Pv2 ðpþ k Þ ¼ ox2 Pv0 ðpk Þ;
k ¼ 2; 3;
ð5:12bÞ
3 o3x1 Pv3 ðp k Þ ¼ ox1 Pv0 ðpk Þ;
k ¼ 3; 4;
ð5:12cÞ
3 o3x2 Pv4 ðp k Þ ¼ ox2 Pv0 ðpk Þ;
k ¼ 4; 1:
ð5:12dÞ
The use of the standard a priori estimate (2.6a) leads to kv0 kpþ1 6 kukpþ1
ð5:13Þ
and ku v0 ks;D 6 Chpþ2s kukpþ2;D ;
s ¼ 0; 1; 2:
ð5:14Þ
Using the continuity conditions (5.12a)–(5.12d) ½o3xj Puðpk Þj ¼ ½o3xj Pv0 ðpk Þj ;
j ¼ 1; 2; k ¼ 1; 2; 3; 4:
ð5:15Þ
For instance, when j ¼ k ¼ 1, 3 ½o3x1 Puðp1 Þ1 ¼ o3x1 Pv1 ðpþ 1 Þ ox1 Pv0 ðp1 Þ:
Using (5.12a)–(5.12d) we have
ð5:16Þ
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
2553
½o3x1 Puðp1 Þ1 ¼ ½o3x1 Pv0 ðp1 Þ1 :
ð5:17Þ
Adding and subtracting v0 to u Pu we obtain 2
2
2
ku Puk2;D ¼ ku v0 k2;D þ 2ðu v0 ; v0 PuÞ2;D þ kv0 Puk2;D :
ð5:18Þ
Using (5.10) we can write 2
2
ku Puk2;D ¼ kv0 Pv0 k2;D þ D ;
ð5:19aÞ
where 2
D ¼ ku v0 k2;D þ 2ðu v0 ; v0 PuÞ2;D :
ð5:19bÞ
Combining (5.19a), (5.4), (5.15) we obtain 2
ku Puk2;D ¼ Ap
4 X 2 X k¼1
2
½o3xj Puðpk Þj þ D :
ð5:20Þ
j¼1
Finally, summing over all elements and using (2.6a) we establish (5.3a) and (5.3b) where 2
2 ¼ ku v0 k2 þ 2ðu v0 ; v0 PuÞ2
ð5:21Þ
and 2
2 6 Ch2p1 kukpþ2 :
ð5:22Þ
Now we are ready to prove the main theorem of this section. Theorem 4. Let u 2 H02 \ H pþ2 , U 2 S N ;p be solutions of (2.2a) and (2.4), respectively. If p P 3 is an odd integer, then N X 2 X 4 X h4 2 2 ku U k2 ¼ ½o3x U ðpk Þj þ ; ð5:23aÞ 2 64p2 ðp 1Þ ð2p 1Þ i¼1 j¼1 k¼1 j where 2
jj 6 Ch2p1 kukpþ2 :
ð5:23bÞ
Proof. Adding and subtracting o3xj Puðpk Þ leads to 2
2
2
½o3xj U ðpk Þj ¼ ½o3xj Puðpk Þj þ 2½o3xj Puðpk Þj ½o3xj ðU Puðpk ÞÞj þ ½o3xj ðU Puðpk ÞÞj :
ð5:24Þ
Summing over j, k, multiplying by Ap from (5.8) and summing over all elements to obtain Ap
N X 2 X 4 X i¼1
j¼1
2
½o3xj U ðpk Þj ¼ Ap
k¼1
N X 2 X 4 X i¼1
j¼1
2
½o3xj Puðpk Þj þ 3 ;
ð5:25aÞ
k¼1
where 3 ¼ 2Ap
N X 2 X 4 X i¼1
j¼1
k¼1
½o3xj Puðpk Þj ½o3xj ðU Puðpk ÞÞj þ Ap
N X 2 X 4 X i¼1
j¼1
½o3xj ðU Puðpk ÞÞ2j :
ð5:25bÞ
k¼1
Let D0 be the canonical element 1 6 n1 , n2 6 1 and use norm equivalence on the finite dimensional space Qp ðD0 Þ to show that
2554
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
max ðjo3n1 vj þ jo3n2 vjÞ 6 Ckvk2;D0 ;
ðn1 ;n2 Þ2D0
8 v 2 Qp ðD0 Þ:
ð5:26Þ
A subsequent linear mapping of D0 to an element Di , i ¼ 1; 2; . . . ; N yields max ðjo3x1 wj þ jo3x2 wjÞ 6
ðx1 ;x2 Þ2Di
C kwk2;Di ; h2
8 w 2 Qp ðDi Þ:
ð5:27Þ
S4 Let Di ¼ n¼0 Di;n , be the union of an element Di;0 and its four neighbors Di;n , n ¼ 1; 2; 3; 4 having one common edge with it. Using (5.27) ½o3xj ðU PuÞðpk Þj 6 C max fjo3x1 ðU PuÞj þ jo3x2 ðU PuÞjg ðx1 ;x2 Þ2Di
6
C kU Puk2;D ; i h2
j ¼ 1; 2; k ¼ 1; 2; 3; 4:
ð5:28Þ
Using (5.28) and (5.8) we obtain Ap
N X 4 X 2 X i¼1
k¼1
½o3xj ðU PuÞðpk Þ2j 6 C
j¼1
N X
kU Puk22;D 6 CkU Puk22 : i
ð5:29Þ
i¼1
Similarly, using (5.3a) and (5.3b) we show Ap
N X 2 X 4 X i¼1
j¼1
½o3xj Puðpk Þ2j 6 Cku Puk22 :
ð5:30Þ
k¼1
Applying the Schwarz inequality we have
N X 2 X 4 X
½o3xj Puðpk Þj ½o3xj ðU PuÞðpk Þj
2Ap
i¼1 j¼1 k¼1 " #1=2 " #1=2 N X 2 X 4 N X 2 X 4 X X 2 2 3 3 6 Ap ½oxj Puðpk Þj Ap ½oxj ðU PuÞðpk Þj : i¼1
j¼1
k¼1
i¼1
j¼1
ð5:31Þ
k¼1
Combining (5.25b), (5.29), (5.30) and (5.31) yields j3 j 6 CkU Puk22 þ ku Puk2 kU Puk2 :
ð5:32Þ
The estimates (2.6a) and (3.25) imply that j3 j 6 Ch2p1 kuk2pþ2 :
ð5:33Þ
Combining (3.28a) and (3.28b), (5.3a) and (5.3b), (5.25a) and (5.33) yields (5.23a) and (5.23b) with ¼ 1 þ 2 3 ; which completes the proof.
ð5:34Þ
6. Error estimators for other finite element spaces The convergence rate is determined by the highest degree polynomial that can be interpolated exactly and the solution smoothness. Thus, piecewise tensor product spaces contain many higher-order terms that do not increase the convergence rate. Other hierarchical bases [2,16,25] lead to better conditioned systems and
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
2555
more efficient methods. Hierarchical bases with fewer degrees of freedom for the same order of accuracy are only available for second-order elliptic and parabolic problems. We will construct C 1 hierarchical shape functions and finite element spaces for fourth-order elliptic problems for which the even–odd a posteriori error procedures of Theorems 2 and 4 are still applicable and which provide optimal rates of convergence with fewer degrees than the classical tensor-product shape functions. The terms x1pþ1 and x2pþ1 are the only monomials missing from a bi p polynomial approximation for it to contain a complete ðp þ 1Þ-degree polynomial. These are the only monomial terms needed to make the error estimates of Sections 4 and 5 asymptotically correct when the other monomial terms of degree p þ 1 are present in the solution space S N ;p . Theorem 5. Under the conditions of Theorems 2 and 4, let Qp ðDi Þ be the restriction of S N ;p to Di (cf. (2.5)) and let Mp ðDi Þ be a space of complete polynomials of degree p on Di . If Qp satisfies Mp Qp Mpþ1 ;
Mpþ1 Qp [ fx1pþ1 ; x2pþ1 g;
ð6:1Þ
then the error estimates (4.26a) and (4.26b) or (5.23a) and (5.23b) apply when p is even or odd, respectively. Proof. The proof goes along the same lines as that of Theorems 4 and 2.
The conditions (6.1) may be used to construct such finite element spaces for which the error estimates (4.26a) and (4.26b) or (5.23a) and (5.23b) are asymptotically correct. However, for efficiency reasons we will use the smallest finite element spaces that satisfy (6.1) Q3 ¼ fNij0;0 ; Nij0;1 ; Nij1;0 ; Nij1;1 ; i; j ¼ 1; 1g;
ð6:2aÞ
j;4 j;4 j;4 j;4 Q4 ¼ Q3 [ fNi;1 ; Ni;2 ; Ni;3 ; Ni;4 ; i ¼ 1; 1; j ¼ 0; 1g;
ð6:2bÞ
j;5 j;5 j;5 j;5 Q5 ¼ Q4 [ fNi;1 ; Ni;2 ; Ni;3 ; Ni;4 ; i ¼ 1; 1; j ¼ 0; 1g;
ð6:2cÞ
j;6 j;6 j;6 j;6 Q6 ¼ Q5 [ fNi;1 ; Ni;2 ; Ni;3 ; Ni;4 ; i ¼ 1; 1; j ¼ 0; 1g [ N04;4 :
ð6:2dÞ
As described above, the set Q3 is the standard bi-cubic C 1 approximation. However, Q4 , Q5 , and Q6 are smaller than bi p, p ¼ 4; 5; 6 spaces with only one interior mode in Q6 . Furthermore, we note that with increasing p, the set Qp contains little more than half the terms in the bi p approximation space. 7. Computational examples We illustrate the performance of the odd–even error estimators for both bi p and hierarchical approximations by solving a fourth-order boundary-value problem using square and rectangular meshes. We will show that the effectivity index g ¼ kEk2 =kek2 approaches unity under mesh refinement. Example 1. We solve the fourth-order boundary-value problem
ð7:1Þ
2556
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
o4 uðx; yÞ o4 uðx; yÞ o4 uðx; yÞ þ2 2 2 þ þ uðx; yÞ ¼ f ðx; yÞ; 4 ox ox oy oy 4
ðx; yÞ 2 X
ð7:2aÞ
with the boundary conditions uðx; yÞ ¼ 0;
ouðx; yÞ ¼ 0; ov
ðx; yÞ 2 oX:
ð7:2bÞ
We select X ¼ ½0; 1 ½0; 1 and the right hand side f ðx; yÞ such that the exact solution is given by 2
2
uðx; yÞ ¼ sinðpxÞ sinðpyÞ :
ð7:2cÞ
We solve (7.2a)–(7.2c) on uniform square-elements having 16, 64, 256, and 1024 elements using bi p approximations of orders 3–6. Errors and global effectivity indices in the H 2 norm are presented in Table 1. The a posteriori estimators behave as predicted by the theory and converge to the true errors under h-refinement. Results of Table 1, as reported in [4], indicate that even-degree estimators perform better that odd-degree estimators. Example 2. We solve (7.2a)–(7.2c) on X ¼ ½0; 2 ½0; 1 on rectangular meshes having 16, 64, 256, and 1024 elements using hierarchical bases defined in (6.2a)–(6.2d) with degrees from 3 to 6. We present the true errors and effectivity indices in the H 2 norm in Table 2. Although we did not do the analysis for rectangular-element meshes, the results in Table 2 show that estimators converge to the true error. The global effectivity indices converge to unity under mesh refinement for the smaller sets Qp of (6.2a)–(6.2d). Again, even-degree estimators perform better than odd-degree estimators. Example 3. Consider the fourth-order boundary-value problem (7.2a) on ½0; 12 where Dirichlet boundary conditions and f ðx; yÞ are selected such that the exact solution is 0:2
2
2
uðx; yÞ ¼ ðx2 þ y 2 Þ x2 ðx 1Þ y 2 ðy 1Þ :
ð7:3Þ
Table 1 Errors and effectivity indices for Example 1 on N-element uniform meshes with piecewise bi p polynomial approximations on square mesh N 16 64 256 1024
p¼3
p¼4
p¼5
p¼6
kek2
g
kek2
g
kek2
g
kek2
g
0.1085(þ1) 0.2761(þ0) 0.6936(1) 0.1736(1)
0.6177 0.9084 0.9876 1.0021
0.1440(00) 0.1833(1) 0.2302(2) 0.2881(3)
0.9941 0.9991 1.0003 1.0006
0.1428(1) 0.9075(3) 0.5696(4) 0.3564(5)
0.5565 0.8841 0.9772 0.9955
0.1130(2) 0.3584(4) 0.1124(5) 0.3516(7)
1.0030 1.0009 1.0003 1.0002
Table 2 Errors and effectivity indices for Example 2 on N-element uniform meshes with piecewise hierarchical polynomial approximations on a rectangular mesh N 16 64 256 1024
p¼3
p¼4
p¼5
p¼6
kek2
g
kek2
g
kek2
g
kek2
g
0.1869(þ1) 0.1119(þ1) 0.2847(1) 0.7150(1)
1.7401 0.7683 0.9474 0.9975
0.1491(þ1) 0.1451(þ0) 0.1847(1) 0.2320(2)
1.0391 1.0085 1.0028 1.0017
0.5120(1) 0.1433(1) 0.9094(3) 0.5707(4)
1.5351 0.6843 0.9163 0.9823
0.4908(1) 0.1133(2) 0.3586(4) 0.1125(5)
1.0987 1.027 1.007 1.002
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
2557
Table 3 Errors and effectivity indices for Example 3 on N-element meshes with piecewise bi p polynomial approximations N
p¼5
p¼6
kek2
g
kek2
g
Uniform meshes 16 64 256
0.5168(4) 0.5151(5) 0.4913(6)
0.4920 0.5773 0.5859
0.1185(4) 0.1182(5) 0.1190(6)
1.6493 1.6167 1.5512
Graded meshes 16 64 256
0.1863(4) 0.1238(5) 0.7850(7)
0.5373 0.8637 0.9630
0.6827(5) 0.2541(6) 0.8167(8)
1.1848 1.0178 1.0037
This solution exhibits a singular behavior near the origin which will limit the rate of convergence in h. On uniform meshes the singularity pollutes the solution and the error estimate globally [14]. Since our error estimates do not account for pollution effects, they will perform poorly. If we eliminate pollution errors by solving the problem on highly graded meshes that resolve the singularity and recover optimal converge rates, we would obtain reasonable accuracy. We begin by solving (7.2a), (7.2b), (7.3) using uniform meshes having 16, 64 and 256 elements with bi 5 and bi 6 polynomial approximations. Effectivity indices shown in Table 3 exhibit large deviations from unity. In order to improve the accuracy of the error estimations, we construct highly graded rectangular meshes using the vertices ðxi ; yj Þ b b pffiffiffiffi i j xi ¼ pffiffiffiffi ; yj ¼ pffiffiffiffi ; i; j ¼ 0; 1; . . . ; N ð7:4Þ N N and solve the problem (7.2a), (7.2b), (7.3) with N ¼ 16, 64, 256. We used bi 5 and bi 6 polynomial approximations with b ¼ 2 and b ¼ 3, respectively. Errors and global effectivity indices are displayed in Table 3. Effectivity indices for graded meshes are much closer to unity than those for uniform meshes. This example shows that the present theory applies to more general situations than those stated previously.
8. Conclusion We have developed asymptotically exact a posteriori error estimates for two-dimensional fourth-order boundary value problems based on the even–odd dichotomy of Babuska and Yu with tensor product bi p finite element approximations. Although, both even- and odd-degree a posteriori error estimators are computationally simple, the even-degree estimator does not involve neighbors communication which is more efficient for parallel implementation and appears to converge to the true error faster as the mesh is h-refined. Finally, we showed that both even- and odd-degree estimators converge to the true error using a more general and efficient hierarchical finite element approximations. Computational results show that the even–odd estimators are asymptotically correct under more general conditions than those indicated by the present theory. For instance, results from Example 3 indicate that they apply on graded quadrilateral meshes in the presence of a singularity. Extension of the analysis to hexahedral meshes and higher-order elliptic equations is straight forward. Adjerid et al. [4] showed that even–odd estimators converged to the true error for second-order problems with both Dirichlet and Neumann boundary conditions on some
2558
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
structured and moderately distorted quadrilateral meshes. We expect this to continue to hold for fourthorder problems. Convergence analyses for nonlinear problems and more general meshes including adaptively refined and/or unstructured quadrilateral, triangular and tetrahedral meshes are less obvious and will be investigated in the future. Likewise, nothing is known about the convergence of even–odd error estimates on curvilinear elements where a theory is needed.
Acknowledgements This Research was supported by the National Science Foundation (grant no. ASC 9720227).
References [1] M. Abramowitz, I.A. Stegun (Eds.), Handbook of Mathematical Functions, Dover, New York, 1965. [2] S. Adjerid, M. Aiffa, J.E. Flaherty, Hierarchical finite element bases for triangular and tetrahedral elements, Comput. Meth. Appl. Mech. Engrg. 190 (2000) 2925–2941. [3] S. Adjerid, I. Babuska, J.E. Flaherty, A posteriori error estimation for the finite element method-of-lines solution of parabolic problems, Math. Model. Meth. Appl. Sci. 9 (1999) 261–286. [4] S. Adjerid, B. Belguendouz, J.E. Flaherty, A posteriori finite element error estimation for diffusive systems, SIAM J. Sci. Comput. 21 (1999) 728–746. [5] S. Adjerid, J.E. Flaherty, L. Krivodonova, A posteriori error estimation for discontinuous Galerkin solutions of hyperbolic problems, Comput. Meth. Appl. Mech. Engrg. 191 (2002) 1097–1112. [6] S. Adjerid, J.E. Flaherty, L. Krivodonova, A posteriori discontinuous Galerkin error estimation for two-dimensional hyperbolic problems, 2000, in preparation. [7] S. Adjerid, J.E. Flaherty, P.K. Moore, Y. Wang, High-order adaptive methods for parabolic systems, Physica D 60 (1992) 94–111. [8] S. Adjerid, J.E. Flaherty, Y. Wang, A posteriori error estimation with finite element methods of lines for one-dimensional parabolic systems, Numer. Math. 65 (1993) 1–21. [9] M. Ainsworth, J.T. Oden, A procedure for a posteriori error estimation for h-p finite element method, Comput. Meth. Appl. Mech. Engrg. 101 (1992) 72–96. [10] M. Ainsworth, J.T. Oden, A unified approach to a posteriori error estimation using element residual methods, Numer. Math. 65 (1993) 23–50. [11] M. Ainsworth, J.T. Oden, A Posteriori Error Estimation in Finite Element Analysis, Wiley, New York, 2000. [12] I. Babuska, T. Strouboulis, S.K. Gangaraj, Guaranteed computable bounds for the exact error in the finite element solution. Part I: one-dimensional model problem, Comput. Meth. Appl. Mech. Engrg. 176 (1999) 51–79. [13] I. Babuska, T. Strouboulis, S.K. Gangaraj, Guaranteed computable bounds for the exact error in the finite element solution. Part II: bounds for the energy norm of the error in two-dimensions, Int. J. Numer. Meth. Eng. 47 (2000) 427–475. [14] I. Babuska, T. Strouboulis, A. Mathur, C.S. Upadhyay, Pollution-error in the h-version of the finite-element method and the local quality of a-posteriori error estimators, Finite Elem. Anal. Des. 17 (1994) 273–321. [15] R. Becker, R. Rannacher, A feedback Approach to error control in finite element methods: basic analysis and examples, EAST– WEST J. Numer. Math. 4 (1996) 237–264. [16] P. Carnevali, R.V. Morric, Y. Tsuji, B. Taylor, New basis functions and computational procedures for p-version finite element analysis, Int. J. Numer. Meth. Eng. 36 (1993) 3759–3779. [17] P. Grisvard, Elliptic Problems in Nonsmooth Domains, Pitman, Boston, 1985. [18] L.G. Leal, Laminar Flow and Convective Transport Processes, Butterworth-Heinemann, Boston, 1992. [19] J.T. Oden, G.F. Carey, Finite Elements, Mathematical Aspects, Prentice-Hall, Englewood Cliffs, 1983. [20] J.T. Oden, S. Prudhomme, Goal-oriented error estimation and adaptivity for the finite element method, Comput. Math. Appl. 41 (2001) 735–756. [21] M. Ortiz, E.A. Repetto, H. Si, A Continuum model for kinetic roughening and coarsening in thin films, J. Mech. Phys. Solids 47 (1999) 697–730. [22] M. Paraschivoiu, J. Preraie, Y. Maday, A.T. Patera, Fast bounds for outputs of partial differential equations, in: J. Borggaard, J. Burns, E. Cliff, S. Schreck (Eds.), Computational Methods for Optimal Design and Control, Birkhauser, Boston, 1998. [23] L.A. Segel, Mathematics Applied to Continuum Mechanics, Dover, New York, 1987.
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
2559
[24] E. S€ uli, P. Houston, Finite element methods for hyperbolic problems: a posteriori error analysis and adaptivity, in: I.S. Duff, G.A. Watson (Eds.), The State of the Art in Numerical Analysis, Clarendon Press, Oxford, 1997, pp. 441–471. [25] B. Szabo, I. Babuska, Finite Element Analysis, John Wiley, New York, 1991. [26] R. Verfurth, A Review of a Posteriori Error Estimation and Adaptive Mesh Refinement Techniques, Teubner-Wiley, New York, 1996. [27] D. Yu, Asymptotically exact a-posteriori error estimator for elements of bi-even degree, Math. Numerica Sinica 13 (1991) 89–101. [28] D. Yu, Asymptotically exact a-posteriori error estimator for elements of bi-odd degree, Math. Numerica Sinica 13 (1991) 307–314. [29] A. Zangwill, Some causes and a consequence of epitaxial roughening, J. Cryst. Growth 163 (1996) 8–21.
A posteriori error estimates for fourth-order elliptic problems Slimane Adjerid Department of Mathematics and Interdisciplinary Center for Applied Mathematics, Virginia Polytechnic Institute and State University, Blacksburg, VA 24061, USA Received 20 February 2001; accepted 19 November 2001
Abstract We extend the dichotomy principle of Babuska and Yu [Math. Numerica Sinica 13 (1991) 89; Math. Numerica Sinica 13 (1991) 307] and Adjerid et al. [Math. Mod. Meth. Appl. S. 9 (1999) 261; SIAM J. Sci. Comput. 21 (1999) 728] for estimating the finite element discretization error to fourth-order elliptic problems. We show how to construct a posteriori error estimates from jumps of the third partial derivatives of the finite element solution when the finite element space consists of piecewise polynomials of odd-degree and from the interior residuals for even-degree approximations on meshes of square elements. These estimates are shown to converge to the true error under mesh refinement. We also show that these a posteriori error estimates are asymptotically correct for more general finite element spaces. Computational results from several examples show that the error estimates are accurate and efficient on rectangular meshes. 2002 Elsevier Science B.V. All rights reserved.
1. Introduction A posteriori error estimators and/or indicators are an integral part of adaptive finite element methods for partial differential equations. They are used for second-order elliptic [9–11,26–28] and parabolic [3,4,7,8] problems. More recently a posteriori error estimates have been developed for hyperbolic problems [5,6,24]. Local error estimates are used to control adaptive enrichment by refining regions having large errors and coarsening regions with small errors. On the other hand, global a posteriori error estimates are used to control solution accuracy. Ideal a posteriori error estimates should be computationally simple and provide guaranteed asymptotical exactness in the sense that the estimates should converge to the true errors from above as the finite element space is enriched. Even–odd error estimates have been used successfully by Yu [27,28] for elliptic problems and Adjerid et al. [3,4] for parabolic problems. They are asymptotically correct for tensor product [3,27,28] and hierarchical [4] finite element spaces under mesh refinement. Most of the existing error estimators are applicable to second-order problems. However, many important physical situations, such as elasticity [23], fluid mechanics [18] and thin film growth [21,29] are modeled by boundary-value problems of fourth order. The first a posteriori error estimate for fourthorder problems appeared in Verfurth [26]. The author is not aware of any other attempt to compute efficient a posteriori error estimates for higher-order problems. In this paper we show how to extend the even–odd error estimation procedures to linear fourth-order boundary value problems. We use C 1 finite 0045-7825/02/$ - see front matter 2002 Elsevier Science B.V. All rights reserved. PII: S 0 0 4 5 - 7 8 2 5 ( 0 1 ) 0 0 4 1 2 - 1
2540
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
element approximations on a square mesh and compute error estimates using jumps of the third partial derivatives of odd-degree finite element solutions and by solving local residual problems for even-degree solutions. The even–odd error estimates are both very simple and efficient to compute. Furthermore, the evendegree error estimators do not involve any neighbor information, which makes them ideal for parallel computations. Moreover, numerical results suggest that even-degree estimators perform better than odddegree estimators. Both even- and odd-degree estimators are extended to more general hierarchical finite element spaces. More recent approaches to a posteriori error estimation include guaranteed estimates [12,13], weighted and goal-oriented estimates [15,20] and bound estimates for functionals [22]. Here we will show that our estimates are asymptotically correct in the H 2 norm. In Section 3, we prove several lemmas. We establish asymptotic exactness of these error estimates on square elements for both bi p even-degree Section 4 and odd-degree Section 5. These error estimates are further extended to the more economical hierarchical approximations in Section 6. We test our error estimators on three examples in Section 7 and present numerical results which show that the actual theory holds for square- and rectangular-element meshes. Finally, numerical results [4] suggest that the even–odd error estimators are applicable to more general meshes as well as problems with singularities and we expect this to continue to hold for fourth-order problems.
2. Problem formulation Let us consider the linear scalar fourth-order partial differential equation x ¼ ðx1 ; x2 Þ 2 X ¼ ½0; 12 ;
Lu ¼ f ; with Lu ¼
X X
ð1Þjbj Db ðaa;b ðxÞDa uÞ;
ð2:1aÞ
ð2:1bÞ
jaj 6 2 jbj 6 2
where jaj ¼ a1 þ a2 , jbj ¼ b1 þ b2 and a2 Da ¼ oa1 x1 ox2 ;
ð2:1cÞ
subject to the Dirichlet boundary conditions uðxÞ ¼ 0;
ou ðxÞ ¼ 0; ov
x 2 oX;
ð2:1dÞ
where ou=ov is the outward normal derivative of u. The functions aa;b ðxÞ, jaj, jbj 6 2, are smooth functions and L is a positive definite operator on the Sobolev space H02 . The Galerkin form of (2.1a)–(2.1d) consists of determining u 2 H02 satisfying Aðv; uÞ ¼ ðv; f Þ;
8 v 2 H02 ;
ð2:2aÞ
2
where the strain energy and L inner products, respectively, are # Z Z "X X Aðv; uÞ ¼ aa;b ; Da uDb v dx1 dx2 X
and ðv; uÞ ¼
Z Z uv dx1 dx2 : X
ð2:2bÞ
jaj 6 2 jbj 6 2
ð2:2cÞ
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
Sobolev spaces H s , s P 0 are equipped with the usual inner product, norm and seminorm X ðDa v; Da uÞ; ðv; uÞs ¼
2541
ð2:3aÞ
jaj 6 s 2
kuks ¼ ðu; uÞs ; X 2 jujs ¼ ðDa u; Da uÞ;
ð2:3bÞ ð2:3cÞ
jaj¼s
where a 0 subscript on H 2 implies that functions also satisfy (2.1d). Finite element solutions of (2.2a) are obtained by approximating H02 by a finite dimensional subspace S0N ;p and determining U 2 S0N ;p such that AðV ; U Þ ¼ ðV ; f Þ;
8 V 2 S0N ;p :
ð2:4Þ
Let us partition X into a uniform mesh of square elements Di , i ¼ 1; 2; . . . ; N and define S N ;p as S N ;p ¼ fw 2 H 2 j wðxÞ 2 Qp ðDi Þ; x 2 Di ; i ¼ 1; 2; . . . ; N g;
ð2:5Þ
where Qp ðDi Þ is the space of bi-polynomial functions that are products of univariate polynomials of degree p in x1 and x2 on Di . Again with the subscript 0 S0N ;p ¼ S N ;p \ H02 . In the following lemmas we describe standard interpolation and a priori discretization error estimates for conforming finite element solutions of (2.1a)–(2.1d) that will be used in the subsequent analysis. Lemma 1. Let W 2 S0N ;p be an interpolant of w 2 H02 \ H pþ1 that is exact when w 2 Qp ðXÞ. Then, there exists a constant C > 0 such that kW wks 6 Chpþ1s kwkpþ1 ;
s ¼ 0; 1; 2;
ð2:6aÞ
where
pffiffiffiffi h ¼ 1= N : H02
ð2:6bÞ pþ1
Further, let u 2 \H and U 2 S a constant C > 0 such that ku U ks 6 Chpþ1s kukpþ1 ;
Proof. Cf. Oden and Carey [19].
N ;p
be solutions of (2.1a)–(2.1d) and (2.4), respectively. Then, there exists
s ¼ 0; 1; 2:
ð2:6cÞ
We also recall the trace theorem. Lemma 2. Let X be a bounded open set of Rn with a Lipschitz boundary oX. Then there exists a constant, C > 0, such that kvk0;oX 6 Cðh1=2 kvk0;X þ h1=2 kvk1;X Þ;
Proof. Cf. Grisvard [17].
8 v 2 H 1:
ð2:7Þ
b 2 Pp ð½1; 12 Þ. Then there exist conLemma 3. Let W 2 Pp ðDÞ and its mapping to the canonical element W stants C1 ; C2 > 0 such that
2542
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559 2 8 W^ 2 Pp ð½1; 1 Þ; s P 0
C1 kW^ ksþ1 < kW^ ks < C2 kW^ ksþ1 ;
ð2:8Þ
and C1 kW ksþ1 < hkW ks < C2 kW ksþ1 ;
8 W 2 Pp ðDÞ; s P 0;
ð2:9Þ
where Pp ðDÞ is the space of polynomials of degree p on an element D:
Proof. Eq. (2.8) is obtained from the equivalence of all norms on finite dimensional spaces. Similarly, 2 Eq. (2.9) is straightforward using the mapping of the square element ½z h=2; z þ h=2 to the canonical 2 element ½1; 1 . 2
Let us construct a hierarchical basis for S N ;p on the canonical element ½1; 1 as a tensor product of the one-dimensional hierarchical basis consisting of the four standard cubic Hermite polynomial shape functions Nlm , l ¼ 1; 1; m ¼ 0; 1, complemented by the interior shape functions N0i , i P 4, such that 0 dN1 ð1Þ ¼ 0; dx
0 N1 ð1Þ ¼ 1;
0 N1 ð1Þ ¼ 0;
1 N1 ð1Þ ¼ 0;
1 dN1 ð1Þ ¼ 1; dx
1 dN1 ð1Þ ¼ 0; dx
dN10 ð1Þ ¼ 0; dx
N10 ð1Þ ¼ 0;
N10 ð1Þ ¼ 1;
N11 ð1Þ ¼ 0;
dN11 ð1Þ ¼ 0; dx
dN11 ð1Þ ¼ 1: dx
ð2:10aÞ ð2:10bÞ ð2:10cÞ ð2:10dÞ
Higher-order interior shape functions are obtained by integrating Lobatto polynomials on ½1; 1 to obtain Z n N0i ðnÞ ¼ Pi1 ðxÞ dx; i ¼ 4; 5; . . . ; p; ð2:11Þ 1
where Pi is the Lobatto polynomial of degree i Pi ðxÞ Pi2 ðxÞ Pi ðxÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; 2ð2i 1Þ
ð2:12Þ
with Pi being the Legendre polynomial of degree i. Using the orthogonality of Legendre polynomials yields 1 Piþ1 Pi1 ð2i 3ÞPiþ1 2ð2i 1ÞPi1 þ ð2i þ 1ÞPi3 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi N0iþ1 ¼ pffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi ¼ : ð2:13Þ 2i þ 1 2i 1 2i 3 2ð2i 1Þð4i2 4i 3Þ One can easily verify that the shape functions N0i , i P 4, satisfy R1 (i) 1 ðd2 N0i =dn2 Þðd2 N0j =dn2 Þ dn ¼ dij , i; j P 4, where dij is the Kronecker symbol. (ii) N0i has i 4 simple roots in ð1; 1Þ and two double roots at 1. 2 (iii) N0i ¼ ðn2 1Þ N~ i4 ðnÞ, where N~ 2l has 2l non-zero roots in ð1; 1Þ symmetrically disposed with respect to the origin, fnk ; k ¼ 1; 2; . . . ; lg. On the other-hand N~ 2lþ1 has 2l+1 roots in ð1; 1Þ, f0g [ fnk ; k ¼ 1; 2; . . . ; lg.
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
2543
The two-dimensional basis on ½1; 1 ½1; 1 is defined by • 16 vertex modes: Nijk;m ðn; gÞ ¼ Nik ðnÞNjm ðgÞ;
i; j ¼ 1; 1; k; m ¼ 0; 1:
ð2:14aÞ
• 16ðp 3Þ edges modes: j;n Ni;1 ¼ Nij ðnÞN0n ðgÞ;
i ¼ 1; 1; j ¼ 0; 1; n ¼ 4; . . . ; p;
ð2:14bÞ
j;n ¼ Nij ðgÞN0n ðnÞ; Ni;2
i ¼ 1; 1; j ¼ 0; 1; n ¼ 4; . . . ; p;
ð2:14cÞ
j;n ¼ Nij ðnÞN0n ðgÞ; Ni;3
i ¼ 1; 1; j ¼ 0; 1; n ¼ 4; . . . ; p;
ð2:14dÞ
j;n ¼ Nij ðnÞN0n ðgÞ; Ni;4
i ¼ 1; 1; j ¼ 0; 1; n ¼ 4; . . . ; p;
ð2:14eÞ
j; k ¼ 4; . . . ; p:
ð2:14fÞ
2
• ðp 3Þ interior modes: N0j;k ¼ N0j ðnÞN0k ðgÞ; 3. Preliminary considerations Let p~½z be the univariate operator that interpolates functions in H02 ðz h=2; z þ h=2Þ and their first derivatives such that p~uðz h=2Þ ¼ uðz h=2Þ;
ð3:1aÞ
d~ pu du ðz h=2Þ ¼ ðz h=2Þ; dx dx
ð3:1bÞ
p~uðz nj h=2Þ ¼ uðz nj h=2Þ;
j ¼ 1; 2; . . . ; s; for p 3 ¼ 2s
ð3:1cÞ
p~uðz nj h=2Þ ¼ uðz nj h=2Þ;
j ¼ 0; 1; 2; . . . ; s; for ðp 3Þ ¼ 2s þ 1;
ð3:1dÞ
or where nj are the simple roots of wpþ1 ðz; zÞ ¼ z
pþ1
N0pþ1 .
Thus, for p P 3, the function
pþ1
p~½zz
ð3:2Þ
N0pþ1
have identical roots with first derivatives which both vanish at z h=2. Hence, wpþ1 ðz; zÞ, and w0pþ1 ðz; zÞ, and w00pþ1 ðz; zÞ are, respectively, proportional to N0pþ1 , Lobatto and Legendre polynomials on 0 ½z h=2; z þ h=2, with ð Þ denoting ordinary differentiation. We use p~ to define a two-dimensional interpolation operator pi on element i satisfying pi uðxÞ ¼ p~½x1;i ~ p½x2;i uðxÞ 2 Qp ðDi Þ, x 2 Di . Since the mesh is uniform, we will omit the elemental index i and the dependence of p~½xj;i and wðxj;i ; xi Þ on the coordinates of the cell center ðx1;i ; x2;i Þ, i ¼ 1; 2; . . . ; N , whenever confusion is unlikely. The functions wpþ1 ðxj Þ, j ¼ 1; 2, provide a basis for the dominant contribution to the spatial discretization error on element D for both odd- and even-order finite element approximations. Indeed, we shall show that estimates EðxÞ of eðxÞ have the form EðxÞ ¼ b1 wpþ1 ðx1 Þ þ b2 wpþ1 ðx2 Þ;
x 2 D:
In order to estimate the finite element discretization error, we need to prove the following series of lemmas.
2544
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
Lemma 4. Let pu be an interpolant of u 2 H pþ2 , then uðxÞ puðxÞ ¼ /ðxÞ þ cðxÞ;
x 2 D;
ð3:3aÞ
where /ðxÞ ¼ b1 wpþ1 ðx1 Þ þ b2 wpþ1 ðx2 Þ;
ð3:3bÞ
kDa /k0;D 6 Chpþ1jaj kukpþ1;D ;
ð3:3cÞ
kcks;D 6 Chpþ2s kukpþ2;D ;
jaj 6 2;
s ¼ 0; 1; . . . ; p þ 1;
ð3:3dÞ
and koxj cks;D 6 Chpþ1s kukpþ2;D ; Z D
o2xj wpþ1 ðxj Þxkj dx1 dx2 ¼ 0;
s ¼ 0; 1; j ¼ 1; 2; k ¼ 0; 1; . . . ; p 2; p P 3:
ð3:3eÞ ð3:3fÞ
Remark. Local Sobolev norms are defined like their global counterparts (2.3b) and (2.3c) with X replaced by D. Proof. Let pu 2 Qp and V ¼ w þ b1 xpþ1 þ b2 x2pþ1 such that w 2 Qp , pV ¼ pu and V interpolates u at two 1 other points. Then the interpolation error can be written as u pu ¼ u V þ V pu: ð3:4Þ Using the fact that pu ¼ pV we have pþ1 p~½x2 x2pþ1 Þ: u pu ¼ u V þ V pV ¼ u V þ b1 ðx1pþ1 p~½x1 xpþ1 1 Þ þ b2 ðx2
ð3:5Þ
u pu ¼ c þ /;
ð3:6Þ
Thus where c ¼ u V . Using standard interpolation error estimates (2.6a) we obtain (3.3c)–(3.3e). Finally, we note that wpþ1 and N0pþ1 are polynomials of degree p þ 1 with identical roots and therefore are proportional, wpþ1 / N0pþ1 , p P 3. This leads to the orthogonality result (3.3f). Lemma 5. Let Pu 2 S0N ;p be an interpolant of u 2 H pþ2 for x 2 X that agrees with pu when x 2 D, then jAðW ; u PuÞj 6 Chp kukpþ2 kW k2 ;
8 W 2 S0N ;p :
ð3:7Þ
Proof. Let us consider first Z að2;0Þ;ð2;0Þ ðxÞ o2x1 W o2x1 ðu puÞ dx1 dx2 D Z Z ¼ að2;0Þ;ð2;0Þ o2x1 W o2x1 ðu puÞ dx1 dx2 þ að2;0Þ;ð2;0Þ ðxÞ að2;0Þ;ð2;0Þ o2x1 W o2x1 ðu puÞ dx1 dx2 ; ð3:8Þ D
D
being the center of the element D. xÞ, with x where að2;0Þ;ð2;0Þ ¼ að2;0Þ;ð2;0Þ ð Using the continuity of aa;b and applying the Schwarz inequality with the estimates (2.6a) we have
Z
2 2
að2;0Þ;ð2;0Þ ox1 W ox1 ðu puÞ dx1 dx2
< Chku Puk2;D kW k2;D < Chp kukpþ2;D kW k2;D : að2;0Þ;ð2;0Þ ðxÞ
D
ð3:9Þ
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
2545
Using (3.3a) and (3.3b), the first term in the right-hand side of (3.8) becomes Z Z o2x1 W o2x1 ðu puÞ dx1 dx2 ¼ o2x1 W ðo2x1 / þ o2x1 cÞ dx1 dx2 :
ð3:10Þ
Using (3.3b) and the fact that w00pþ1 is proportional to the Legendre polynomial Pp1 we have Z o2x1 W o2x1 / dx1 dx2 ¼ 0; 8 W 2 S N ;p :
ð3:11Þ
D
D
D
Applying the Schwarz inequality and using (3.3e) leads to
Z
o2 W o2 c dx1 dx2 6 Chp kuk kW k ; 8 W 2 S N ;p : pþ2 2 x1 x1
ð3:12Þ
D
Following the same steps we prove jðo2x2 W ; o2x2 ðu puÞÞj 6 Chp kukpþ2;D kW k2D ;
8 W 2 S N ;p :
Applying Green’s theorem we obtain Z Z Z o2x2 W o2x1 ðu puÞ dx1 dx2 ¼ ox2 W o2x1 ðu puÞdx1 C 1
D
Z
ð3:13Þ
Cþ 1
ox2 W o2x1 ðu puÞ dx1
ox2 W o2x1 ox2 ðu puÞ dx1 dx2 :
ð3:14Þ
Using (3.3a) we have Z Z 2 ox2 W ox1 ox2 ðu puÞ dx1 dx2 ¼ ox2 W o2x1 ox2 ð/ þ cÞ dx1 dx2 :
ð3:15Þ
Using the fact that o2x1 ox2 / ¼ 0 we have Z Z ox2 W o2x1 ox2 ðu puÞ dx1 dx2 ¼ ox2 W o2x1 ox2 c dx1 dx2 :
ð3:16Þ
D
D
D
D
D
Integration by parts yields Z Z Z 2 2 2 ox2 W ox1 ox2 c dx1 dx2 ¼ ox2 W ox1 c dx1 dx2 þ D
D
Integrating the boundary terms by parts we obtain Z Z Z ox2 W o2x1 ox2 c dx1 dx2 ¼ o2x2 W o2x1 c dx1 dx2 D
D
C 1
C 1
ox2 W
o2x1 c dx1
Z Cþ 1
ox1 ox2 W ox1 c dx1 þ
o2x2 W o2x1 c dx1 :
Z Cþ 1
ox1 o2x2 W ox1 c dx1 :
Applying the Schwarz inequality and using (3.3e), we immediately have
Z
o2 W o2 c dx1 dx2 6 kW k kck 6 Chp kuk kW k : 2 2 pþ2 2 x2 x1
ð3:17Þ
ð3:18Þ
ð3:19Þ
D
The boundary terms in (3.18) are estimated by applying the Schwarz inequality with (2.7) and (3.3d) to obtain
Z
ox1 ox2 W ox1 c dx1 6 Cðh1=2 kW k2 þ h1=2 kW k3 Þðh1=2 kck1 þ h1=2 kck2 Þ
C 1
6 Chp kukpþ2;D ðkW k2;D þ hkW k3;D Þ 6 Chp kukpþ2;D kW k2;D :
ð3:20Þ
2546
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
We further integrate by parts the boundary terms in (3.14) and use (3.1a)–(3.1d) to obtain Z Z ox2 W o2x1 ðu puÞ dx1 ¼ ox1 ox2 W ox1 ðu puÞ dx1 : C 1
ð3:21Þ
C 1
Since ox1 ox2 W and ox1 ðu puÞ are continuous along the horizontal edge C shared by two neighboring elements Di and Dj , the contribution of such terms to (3.7) can be written as Z D ð aDa;bi aa;bj Þ ox1 ox2 W ox1 ðu puÞ dx1 : ð3:22Þ C
Applying the Schwarz inequality and using (2.6a), (2.7), (2.8) and the smoothness of the coefficients of L we obtain
Z
D
Dj i
ð aa;b Þ ox1 ox2 W ox1 ðu puÞ dx1
< CðkW k2 þ hkW k3 Þðku puk1 þ hku puk2 Þ
aa;b C
< CkW k2;Di \Dj ku puk1;Di \Dj < Chp kW k2;Di \Dj kukpþ2;Di \Dj : Following similar steps we prove
Z
a b p
a D ðu puÞD W dx dx a;b 1 2 < Ch kukpþ2;X kW k2;X ;
8 a; b:
ð3:23Þ
ð3:24Þ
X
Combining (3.8)–(3.24) completes the proof of (3.7).
In the following lemma we state a superconvergence result. Lemma 6. Let U and u 2 H pþ2 be solution of (2.4) and (2.2a), respectively and Pu 2 S N ;p interpolate u as described in Lemma 5, then kU Puk2 < Chp kukpþ2 :
ð3:25Þ
Proof. Subtracting AðW ; puÞ from the orthogonality condition 8 W 2 S N ;p
ð3:26Þ
AðW ; U puÞ ¼ AðW ; u puÞ:
ð3:27Þ
AðW ; U uÞ ¼ 0; leads to
Using Lemma 5 with W ¼ U pu and the ellipticity of the bilinear form A we complete the proof. In the following theorem we show that the finite element error can be divided into a leading component / plus a higher-order component h. Lemma 7. Under the conditions of Lemma 6, there exists 1 such that kek22 ¼ ku Puk22 þ 1
ð3:28aÞ
j1 j 6 CðuÞh2p1 :
ð3:28bÞ
with
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
2547
Proof. Adding and subtracting Pu to e yields (3.28a) with 2
1 ¼ 2ðu Pu; Pu U Þ2 þ kPu U k2 : Applying the Schwarz inequality and using (2.6a) and (3.25) yields (3.28b).
ð3:29Þ
Theorem 1. If u and U are as defined in Lemma 6, then e ¼ u U ¼ U þ H;
ð3:30aÞ
where the restrictions of U and H to D are / of (3.3b) and h ¼ c þ pu U :
ð3:30bÞ
Furthermore, kUk22 6 Ch2p2 kuk2pþ1 ;
ð3:30cÞ
jHk22 6 Ch2p kuk2pþ1 :
ð3:30dÞ
Proof. Adding and subtracting pu to e and using (3.3a)–(3.3d) and (3.25) we prove the theorem.
4. Even-degree a posteriori error estimation Error estimates of even-degree approximations are computed using interior residuals by solving local Petrov–Galerkin problems for the error. In order to show this we substitute u¼eþU
ð4:1Þ
in (2.2a) to obtain Aðv; eÞ ¼ gðvÞ;
8 v 2 H02 ;
ð4:2Þ
where gðvÞ ¼ ðv; f Þ Aðv; U Þ:
ð4:3Þ
The error is approximated by E ¼ b1 wpþ1 ðx1 Þ þ b2 wpþ1 ðx2 Þ;
ð4:4Þ
on element D and the test function v is selected as vj ðxÞ ¼ ðxj xj;i Þdðx1 Þdðx2 Þ;
ð4:5Þ
where dðzÞ ¼
wpþ1 ðz; zÞ : z z
ð4:6Þ
Since wpþ1 ðxj Þ and vj ðxÞ; j ¼ 1; 2 vanish on oD, the error estimate satisfies the local Dirichlet problems AD ðvj ; EÞ ¼ gD ðvj Þ;
j ¼ 1; 2;
ð4:7Þ
where the subscript D denotes that inner products are restricted to D. Using (3.30a)–(3.30d), the finite element problem (4.7) becomes AD ðv; /Þ ¼ gD ðvÞ AD ðv; hÞ;
8 v 2 H02 :
ð4:8Þ
2548
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
Which in turn can be written as D ðv; /Þ ¼ gD ðvÞ AD ðv; hÞ FD ðv; /Þ; A where D ðv; /Þ ¼ A
Z Z "X X D
and FD ðv; /Þ ¼
ð4:9Þ
# a
b
aa;b ð xÞD /D v dx1 dx2
ð4:10Þ
jaj¼2 jbj¼2
Z Z "X X D
þ
8 v 2 H02 ;
# a
ðaa;b ðxÞ aa;b ð xÞÞD /D v dx1 dx2
jaj¼2 jbj¼2
Z Z "X X D
b
# a
b
aa;b ðxÞD /D v dx1 dx2 :
ð4:11Þ
jaj<2 jbj<2
Using symmetry properties of wpþ1 ðzÞ and vj ðxÞ we obtain the uncoupled constant-coefficient problem on element D Z Z D dðx2 Þðw00pþ1 ðx1 ÞÞ2 dx1 dx2 ¼ gðv1 Þ Aðv1 ; hÞ F ðv1 ; /Þ; ð4:12aÞ að2;0Þ;ð2;0Þ b1 D
aDð0;2Þ;ð0;2Þ b2
Z Z
2
dðx1 Þðw00pþ1 ðx2 ÞÞ dx1 dx2 ¼ gðv2 Þ Aðv2 ; hÞ F ðv2 ; /Þ:
ð4:12bÞ
D
The finite element problem (4.12a) and (4.12b) may be further simplified by neglecting the higher order terms to obtain the following system for bj , j ¼ 1; 2, Z Z 2 b1 dðx2 Þðw00pþ1 ðx1 ÞÞ dx1 dx2 ¼ gðv1 Þ; ð4:13aÞ aDð2;0Þ;ð2;0Þ D
aDð0;2Þ;ð0;2Þ b2
Z Z
dðx1 Þðw00pþ1 ðx2 ÞÞ2 dx1 dx2 ¼ gðv2 Þ:
ð4:13bÞ
D
The error estimate E in (4.4) with (4.13a) and (4.13b) is shown to be asymptotically correct in the next theorem. However, prior to this we need the following properties of wpþ1 , d and vj , j ¼ 1; 2: Lemma 8. Let p > 3 be an even integer, then there exist C > 0 and c > 0 such that 2
ch2ðpsÞþ4 < kwpþ1 ks;D < Ch2ðpsÞþ4 ; Z Z
s ¼ 0; 1; 2;
dðxj Þdx1 dx2 ¼ Chpþ2 ;
ð4:14aÞ ð4:14bÞ
D
Z Z
2
dðxk Þw00 ðxj Þ dx1 dx2 ¼ Ch3p ;
j; k ¼ 1; 2;
ð4:14cÞ
D 2
ch4p < kvj k2;D < Ch4p : Proof. A direct computation reveals the results.
ð4:14dÞ
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
2549
In the following lemma we establish bounds for the difference between bj and its approximation bj , j ¼ 1; 2. Lemma 9. Let p > 3 be an even integer and u 2 H02 \ H pþ2 . If the coefficients of L in (2.1b) are smooth functions then N X
ðb2ji b2ji Þ 6
i¼1
C kuk2pþ2 ; h
j ¼ 1; 2:
ð4:15Þ
Proof. Subtracting (4.13a) and (4.13b) from (4.12a) and (4.12b) we obtain 2
ðbji bji Þ2 ¼
ðAðvj ; hÞ þ F ðvj ; /ÞÞ : RR 2 2 dðx2 Þw00pþ1 ðx1 Þ dx1 dx2 Þ Di
i ð aDð2;0Þ;ð2;0Þ
ð4:16Þ
Applying the Schwarz inequality and using (4.14a)–(4.14d), (3.30d) and the continuity of aa;b ðxÞ we have jAðvj ; hÞj2 6 Ckvj k22;Di khk22;Di 6 Ch6p kuk2pþ2;Di :
ð4:17Þ
Using the continuity of aa;b (x) with (3.3c) we also prove jF ðvj ; /Þj2 6 Ch2 kvj k22;Di k/k22;Di 6 Ch6p kuk2pþ2;Di :
ð4:18Þ
Combining (4.16)–(4.18) with (4.14a)–(4.14d) yields ðbji bji Þ2 < Ckuk2pþ2;Di :
ð4:19Þ
A summation over all elements leads to N X
2
2
ðbji bji Þ < Ckukpþ2 :
ð4:20Þ
i¼1
Using (3.3b) and (4.14a)–(4.14d) we obtain 2
b2ji ¼
ko2xj /k0;Di kw00pþ1 ðxj Þk20;Di
<
C kuk2pþ2;Di : h2
ð4:21Þ
We may solve (4.13a) and (4.13b) for bj to obtain 2
b2ji ¼
i ð aDð2;0Þ;ð2;0Þ
gðvj Þ
RR
2
Di
dðx2 Þw00pþ1 ðx1 Þ dx1 dx2 Þ
2
:
ð4:22Þ
If we apply the Schwarz inequality, use (4.2) and (4.14a)–(4.14d), sum over all elements and use (2.6c) we obtain N X
2
b2ji 6 Ch2p kek2 6
i¼1
C 2 kukpþ2 : h2
ð4:23Þ
Combining (4.21) and (4.23) N X i¼1
½b2ji þ b2ji 6
C 2 kukpþ2 : h2
ð4:24Þ
2550
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
Applying the Schwarz inequality #1=2 " #1=2 " N N N X X X 2 2 2 2 2 ½bji bji 6 C ðbji bji Þ ðbji þ bji Þ i¼1
i¼1
ð4:25Þ
i¼1
while using (4.20) and (4.24) leads to (4.15).
Now we will state and prove the main result of this section. Theorem 2. Let u 2 H02 \ H pþ2 and U 2 S N ;p be solutions of (2.2a) and (2.4), respectively. If p > 3 is an even integer, then kek22 ¼
N X 2 X i¼1
b2ji kw00pþ1 ðxj Þk20;Di þ ;
ð4:26aÞ
j¼1
where 2
jj 6 Ch2p1 kukpþ2 :
ð4:26bÞ
Proof. Consider kek22 ¼
N h i X ko2x1 ek20;Di þ ko2x2 ek20;Di þ kox1 ox2 ek20;Di þ kek21
ð4:27Þ
i¼1
and use (3.30a)–(3.30d) and (3.3b) to obtain (Z Z " ! N 2 X X 2 2 2 00 kek2 ¼ bji jwpþ1 ðxj Þj þ jo2x1 hj2 þ 2o2x1 / o2x1 h þ jo2x2 hj2 Di
i¼1
þ
j¼1
2o2x2 / o2x2 h
# 2
þ jox1 ox2 hj
) dx1 dx2
2
þ kek1 :
ð4:28Þ
2
Adding and subtracting b2ji jw00pþ1 ðxj Þj , to the above integrand we obtain (4.26a) with (Z Z " ! N 2 X X 2 2 00 2 ðbji bji Þjwpþ1 ðxj Þj þ jo2x1 hj2 þ 2o2x1 / o2x1 h þ jo2x2 hj2 ¼ i¼1
Di
j¼1
#
þ 2o2x2 / o2x2 h þ jox1 ox2 hj2 dx1 dx2
) þ kek21 :
ð4:29Þ
Applying the Schwarz and triangle inequalities and using the estimates (2.6c), (3.3c), (3.3d), (4.14a)–(4.14d) and (4.15) yields (4.26b). 5. Odd-degree a posteriori estimator If u is a smooth solution of (2.1a)–(2.1d) on D, then error estimates E of odd-degree approximations may be computed in terms of jumps in derivatives o3xj U , j ¼ 1; 2 at the vertices of D. Using e ¼ u U we compute the jumps in derivatives at pk ¼ ðp1k ; p2k Þ, k ¼ 1; 2; 3; 4, of D as
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
½o3xj eðpk Þj ¼ ½o3xj U ðpk Þj b1 ½o3xj wpþ1 ðp1k Þj þ b2 ½o3xj wpþ1 ðp2k Þj ;
2551
j ¼ 1; 2; k ¼ 1; 2; 3; 4; x 2 D; ð5:1Þ
where ½qðpÞj denotes the jump in q at point p in the xj direction. Since jumps in the finite element solution o3xj U and o3xj wpþ1 , j ¼ 1; 2, are known, (5.1) yields a linear algebraic system for b1 and b2 at each vertex pk . On a rectangular mesh the set of systems (5.1) have a unique solution. More general meshes are treated in [4]. kEk22;D is the average of local errors based on the pairs of bj , j ¼ 1; 2. In the next lemma we state some useful properties of Legendre polynomials. Lemma 10. Let Pp be the Legendre polynomial of degree p, then Z 1 2 ; Pp2 ðtÞ dt ¼ 2p þ1 1 0 0 Ppþ1 ¼ ð2p þ 1ÞPp þ Pp1 ;
Pp ðþ1Þ ¼ 1;
pþ1
pðp þ 1Þ ; 2
ð5:2bÞ
p > 1;
Pp ð1Þ ¼ ð1Þp ;
Pp0 ð1Þ ¼ ð1Þ
ð5:2aÞ
ð5:2cÞ
p P 0;
ð5:2dÞ
p > 0:
Proof. See Abramowitz and Stegun [1].
Theorem 3. Let u 2 H02 \ H pþ2 , Pu be as defined in Lemma 5, and p P 3 be an odd integer, then ku Puk22 ¼
N X 2 X 4 X
h4 64p2 ðp 1Þ2 ð2p 1Þ
i¼1
j¼1
½o3xj Puðpk Þ2j þ 2 ;
ð5:3aÞ
k¼1
where 2
j2 j 6 Ch2p1 kukpþ2 :
ð5:3bÞ
Proof. The proof will be established in two steps. In the first step we assume u ¼ uj ¼ xjpþ1 , j ¼ 1; 2 and show that ko2xj ðuj Puj Þk2 ¼
h4
4 X
64p2 ðp 1Þ2 ð2p 1Þ
k¼1
½o3xj Puj 2j :
ð5:4Þ
Consider two neighboring elements Di and Dik and compute the jump in rj ðxj Þ ¼ o3xj wpþ1 ¼ w000 pþ1 ðxj Þ in the jth direction at common vertices. On a square-element mesh we have ½o3xj ðuj Puj Þj ¼ 2rj ðh=2Þ;
ð5:5Þ P4
ko2xj ðuj
2
Puj Þk ¼
2 kw00pþ1 ðxj Þk
k¼1
2
½o3xj Puj j
16ðrj ðh=2ÞÞ
2
¼ Ap
4 X k¼1
2
½o3xj Puj j ;
ð5:6aÞ
2552
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
where 2
Ap ¼
kw00pþ1 kD 16ðrj ðh=2ÞÞ2
ð5:6bÞ
:
Transforming the integrals from a physical element D to the canonical element, applying the chain rule and using the fact that w00pþ1 / Pp1 to obtain n n2 dn n d n d w ðzÞ ¼ ð2=hÞ w ðnÞ ¼ cð2=hÞ Pp1 ðnÞ; n pþ1 pþ1 dzn dn dnn2 Thus, Z 1 h4 h4 2 Ap ¼ : Pp1 ðtÞ dt ¼ 2 2 27 p2 ðp 1Þ 1 64p2 ðp 1Þ ð2p 1Þ
n ¼ 2; 3:
ð5:7Þ
ð5:8Þ
This proves (5.4). Now, we assume w to be a p-degree polynomial with respect to x1 and x2 . Thus, ½oxj ðw PwÞðpk Þ ¼ 0;
w Pw ¼ 0; b1 x1pþ1
b2 x2pþ1S
Let v ¼ w þ þ neighbors of D with XD ¼ Pv0 ¼ Pu and
j ¼ 1; 2:
ð5:9Þ
and Di , i ¼ 1; 2; 3; 4 denote, respectively, the right, top, left, and bottom Dk [ D. We define a function v such that vjD ¼ v0 and vjDk ¼ vk such that
4 k¼1
v0 ðqj Þ ¼ uðqj Þ;
ð5:10Þ
where qj , j ¼ 1; 2 are two other points in D, Pvk ¼ Pu;
for x 2 Dk ;
k ¼ 1; 2; 3; 4;
ð5:11Þ
and the following continuity conditions: 3 þ o3x1 Pv1 ðpþ k Þ ¼ ox1 Pv0 ðpk Þ;
k ¼ 1; 2;
ð5:12aÞ
3 þ o3x2 Pv2 ðpþ k Þ ¼ ox2 Pv0 ðpk Þ;
k ¼ 2; 3;
ð5:12bÞ
3 o3x1 Pv3 ðp k Þ ¼ ox1 Pv0 ðpk Þ;
k ¼ 3; 4;
ð5:12cÞ
3 o3x2 Pv4 ðp k Þ ¼ ox2 Pv0 ðpk Þ;
k ¼ 4; 1:
ð5:12dÞ
The use of the standard a priori estimate (2.6a) leads to kv0 kpþ1 6 kukpþ1
ð5:13Þ
and ku v0 ks;D 6 Chpþ2s kukpþ2;D ;
s ¼ 0; 1; 2:
ð5:14Þ
Using the continuity conditions (5.12a)–(5.12d) ½o3xj Puðpk Þj ¼ ½o3xj Pv0 ðpk Þj ;
j ¼ 1; 2; k ¼ 1; 2; 3; 4:
ð5:15Þ
For instance, when j ¼ k ¼ 1, 3 ½o3x1 Puðp1 Þ1 ¼ o3x1 Pv1 ðpþ 1 Þ ox1 Pv0 ðp1 Þ:
Using (5.12a)–(5.12d) we have
ð5:16Þ
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
2553
½o3x1 Puðp1 Þ1 ¼ ½o3x1 Pv0 ðp1 Þ1 :
ð5:17Þ
Adding and subtracting v0 to u Pu we obtain 2
2
2
ku Puk2;D ¼ ku v0 k2;D þ 2ðu v0 ; v0 PuÞ2;D þ kv0 Puk2;D :
ð5:18Þ
Using (5.10) we can write 2
2
ku Puk2;D ¼ kv0 Pv0 k2;D þ D ;
ð5:19aÞ
where 2
D ¼ ku v0 k2;D þ 2ðu v0 ; v0 PuÞ2;D :
ð5:19bÞ
Combining (5.19a), (5.4), (5.15) we obtain 2
ku Puk2;D ¼ Ap
4 X 2 X k¼1
2
½o3xj Puðpk Þj þ D :
ð5:20Þ
j¼1
Finally, summing over all elements and using (2.6a) we establish (5.3a) and (5.3b) where 2
2 ¼ ku v0 k2 þ 2ðu v0 ; v0 PuÞ2
ð5:21Þ
and 2
2 6 Ch2p1 kukpþ2 :
ð5:22Þ
Now we are ready to prove the main theorem of this section. Theorem 4. Let u 2 H02 \ H pþ2 , U 2 S N ;p be solutions of (2.2a) and (2.4), respectively. If p P 3 is an odd integer, then N X 2 X 4 X h4 2 2 ku U k2 ¼ ½o3x U ðpk Þj þ ; ð5:23aÞ 2 64p2 ðp 1Þ ð2p 1Þ i¼1 j¼1 k¼1 j where 2
jj 6 Ch2p1 kukpþ2 :
ð5:23bÞ
Proof. Adding and subtracting o3xj Puðpk Þ leads to 2
2
2
½o3xj U ðpk Þj ¼ ½o3xj Puðpk Þj þ 2½o3xj Puðpk Þj ½o3xj ðU Puðpk ÞÞj þ ½o3xj ðU Puðpk ÞÞj :
ð5:24Þ
Summing over j, k, multiplying by Ap from (5.8) and summing over all elements to obtain Ap
N X 2 X 4 X i¼1
j¼1
2
½o3xj U ðpk Þj ¼ Ap
k¼1
N X 2 X 4 X i¼1
j¼1
2
½o3xj Puðpk Þj þ 3 ;
ð5:25aÞ
k¼1
where 3 ¼ 2Ap
N X 2 X 4 X i¼1
j¼1
k¼1
½o3xj Puðpk Þj ½o3xj ðU Puðpk ÞÞj þ Ap
N X 2 X 4 X i¼1
j¼1
½o3xj ðU Puðpk ÞÞ2j :
ð5:25bÞ
k¼1
Let D0 be the canonical element 1 6 n1 , n2 6 1 and use norm equivalence on the finite dimensional space Qp ðD0 Þ to show that
2554
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
max ðjo3n1 vj þ jo3n2 vjÞ 6 Ckvk2;D0 ;
ðn1 ;n2 Þ2D0
8 v 2 Qp ðD0 Þ:
ð5:26Þ
A subsequent linear mapping of D0 to an element Di , i ¼ 1; 2; . . . ; N yields max ðjo3x1 wj þ jo3x2 wjÞ 6
ðx1 ;x2 Þ2Di
C kwk2;Di ; h2
8 w 2 Qp ðDi Þ:
ð5:27Þ
S4 Let Di ¼ n¼0 Di;n , be the union of an element Di;0 and its four neighbors Di;n , n ¼ 1; 2; 3; 4 having one common edge with it. Using (5.27) ½o3xj ðU PuÞðpk Þj 6 C max fjo3x1 ðU PuÞj þ jo3x2 ðU PuÞjg ðx1 ;x2 Þ2Di
6
C kU Puk2;D ; i h2
j ¼ 1; 2; k ¼ 1; 2; 3; 4:
ð5:28Þ
Using (5.28) and (5.8) we obtain Ap
N X 4 X 2 X i¼1
k¼1
½o3xj ðU PuÞðpk Þ2j 6 C
j¼1
N X
kU Puk22;D 6 CkU Puk22 : i
ð5:29Þ
i¼1
Similarly, using (5.3a) and (5.3b) we show Ap
N X 2 X 4 X i¼1
j¼1
½o3xj Puðpk Þ2j 6 Cku Puk22 :
ð5:30Þ
k¼1
Applying the Schwarz inequality we have
N X 2 X 4 X
½o3xj Puðpk Þj ½o3xj ðU PuÞðpk Þj
2Ap
i¼1 j¼1 k¼1 " #1=2 " #1=2 N X 2 X 4 N X 2 X 4 X X 2 2 3 3 6 Ap ½oxj Puðpk Þj Ap ½oxj ðU PuÞðpk Þj : i¼1
j¼1
k¼1
i¼1
j¼1
ð5:31Þ
k¼1
Combining (5.25b), (5.29), (5.30) and (5.31) yields j3 j 6 CkU Puk22 þ ku Puk2 kU Puk2 :
ð5:32Þ
The estimates (2.6a) and (3.25) imply that j3 j 6 Ch2p1 kuk2pþ2 :
ð5:33Þ
Combining (3.28a) and (3.28b), (5.3a) and (5.3b), (5.25a) and (5.33) yields (5.23a) and (5.23b) with ¼ 1 þ 2 3 ; which completes the proof.
ð5:34Þ
6. Error estimators for other finite element spaces The convergence rate is determined by the highest degree polynomial that can be interpolated exactly and the solution smoothness. Thus, piecewise tensor product spaces contain many higher-order terms that do not increase the convergence rate. Other hierarchical bases [2,16,25] lead to better conditioned systems and
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
2555
more efficient methods. Hierarchical bases with fewer degrees of freedom for the same order of accuracy are only available for second-order elliptic and parabolic problems. We will construct C 1 hierarchical shape functions and finite element spaces for fourth-order elliptic problems for which the even–odd a posteriori error procedures of Theorems 2 and 4 are still applicable and which provide optimal rates of convergence with fewer degrees than the classical tensor-product shape functions. The terms x1pþ1 and x2pþ1 are the only monomials missing from a bi p polynomial approximation for it to contain a complete ðp þ 1Þ-degree polynomial. These are the only monomial terms needed to make the error estimates of Sections 4 and 5 asymptotically correct when the other monomial terms of degree p þ 1 are present in the solution space S N ;p . Theorem 5. Under the conditions of Theorems 2 and 4, let Qp ðDi Þ be the restriction of S N ;p to Di (cf. (2.5)) and let Mp ðDi Þ be a space of complete polynomials of degree p on Di . If Qp satisfies Mp Qp Mpþ1 ;
Mpþ1 Qp [ fx1pþ1 ; x2pþ1 g;
ð6:1Þ
then the error estimates (4.26a) and (4.26b) or (5.23a) and (5.23b) apply when p is even or odd, respectively. Proof. The proof goes along the same lines as that of Theorems 4 and 2.
The conditions (6.1) may be used to construct such finite element spaces for which the error estimates (4.26a) and (4.26b) or (5.23a) and (5.23b) are asymptotically correct. However, for efficiency reasons we will use the smallest finite element spaces that satisfy (6.1) Q3 ¼ fNij0;0 ; Nij0;1 ; Nij1;0 ; Nij1;1 ; i; j ¼ 1; 1g;
ð6:2aÞ
j;4 j;4 j;4 j;4 Q4 ¼ Q3 [ fNi;1 ; Ni;2 ; Ni;3 ; Ni;4 ; i ¼ 1; 1; j ¼ 0; 1g;
ð6:2bÞ
j;5 j;5 j;5 j;5 Q5 ¼ Q4 [ fNi;1 ; Ni;2 ; Ni;3 ; Ni;4 ; i ¼ 1; 1; j ¼ 0; 1g;
ð6:2cÞ
j;6 j;6 j;6 j;6 Q6 ¼ Q5 [ fNi;1 ; Ni;2 ; Ni;3 ; Ni;4 ; i ¼ 1; 1; j ¼ 0; 1g [ N04;4 :
ð6:2dÞ
As described above, the set Q3 is the standard bi-cubic C 1 approximation. However, Q4 , Q5 , and Q6 are smaller than bi p, p ¼ 4; 5; 6 spaces with only one interior mode in Q6 . Furthermore, we note that with increasing p, the set Qp contains little more than half the terms in the bi p approximation space. 7. Computational examples We illustrate the performance of the odd–even error estimators for both bi p and hierarchical approximations by solving a fourth-order boundary-value problem using square and rectangular meshes. We will show that the effectivity index g ¼ kEk2 =kek2 approaches unity under mesh refinement. Example 1. We solve the fourth-order boundary-value problem
ð7:1Þ
2556
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
o4 uðx; yÞ o4 uðx; yÞ o4 uðx; yÞ þ2 2 2 þ þ uðx; yÞ ¼ f ðx; yÞ; 4 ox ox oy oy 4
ðx; yÞ 2 X
ð7:2aÞ
with the boundary conditions uðx; yÞ ¼ 0;
ouðx; yÞ ¼ 0; ov
ðx; yÞ 2 oX:
ð7:2bÞ
We select X ¼ ½0; 1 ½0; 1 and the right hand side f ðx; yÞ such that the exact solution is given by 2
2
uðx; yÞ ¼ sinðpxÞ sinðpyÞ :
ð7:2cÞ
We solve (7.2a)–(7.2c) on uniform square-elements having 16, 64, 256, and 1024 elements using bi p approximations of orders 3–6. Errors and global effectivity indices in the H 2 norm are presented in Table 1. The a posteriori estimators behave as predicted by the theory and converge to the true errors under h-refinement. Results of Table 1, as reported in [4], indicate that even-degree estimators perform better that odd-degree estimators. Example 2. We solve (7.2a)–(7.2c) on X ¼ ½0; 2 ½0; 1 on rectangular meshes having 16, 64, 256, and 1024 elements using hierarchical bases defined in (6.2a)–(6.2d) with degrees from 3 to 6. We present the true errors and effectivity indices in the H 2 norm in Table 2. Although we did not do the analysis for rectangular-element meshes, the results in Table 2 show that estimators converge to the true error. The global effectivity indices converge to unity under mesh refinement for the smaller sets Qp of (6.2a)–(6.2d). Again, even-degree estimators perform better than odd-degree estimators. Example 3. Consider the fourth-order boundary-value problem (7.2a) on ½0; 12 where Dirichlet boundary conditions and f ðx; yÞ are selected such that the exact solution is 0:2
2
2
uðx; yÞ ¼ ðx2 þ y 2 Þ x2 ðx 1Þ y 2 ðy 1Þ :
ð7:3Þ
Table 1 Errors and effectivity indices for Example 1 on N-element uniform meshes with piecewise bi p polynomial approximations on square mesh N 16 64 256 1024
p¼3
p¼4
p¼5
p¼6
kek2
g
kek2
g
kek2
g
kek2
g
0.1085(þ1) 0.2761(þ0) 0.6936(1) 0.1736(1)
0.6177 0.9084 0.9876 1.0021
0.1440(00) 0.1833(1) 0.2302(2) 0.2881(3)
0.9941 0.9991 1.0003 1.0006
0.1428(1) 0.9075(3) 0.5696(4) 0.3564(5)
0.5565 0.8841 0.9772 0.9955
0.1130(2) 0.3584(4) 0.1124(5) 0.3516(7)
1.0030 1.0009 1.0003 1.0002
Table 2 Errors and effectivity indices for Example 2 on N-element uniform meshes with piecewise hierarchical polynomial approximations on a rectangular mesh N 16 64 256 1024
p¼3
p¼4
p¼5
p¼6
kek2
g
kek2
g
kek2
g
kek2
g
0.1869(þ1) 0.1119(þ1) 0.2847(1) 0.7150(1)
1.7401 0.7683 0.9474 0.9975
0.1491(þ1) 0.1451(þ0) 0.1847(1) 0.2320(2)
1.0391 1.0085 1.0028 1.0017
0.5120(1) 0.1433(1) 0.9094(3) 0.5707(4)
1.5351 0.6843 0.9163 0.9823
0.4908(1) 0.1133(2) 0.3586(4) 0.1125(5)
1.0987 1.027 1.007 1.002
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
2557
Table 3 Errors and effectivity indices for Example 3 on N-element meshes with piecewise bi p polynomial approximations N
p¼5
p¼6
kek2
g
kek2
g
Uniform meshes 16 64 256
0.5168(4) 0.5151(5) 0.4913(6)
0.4920 0.5773 0.5859
0.1185(4) 0.1182(5) 0.1190(6)
1.6493 1.6167 1.5512
Graded meshes 16 64 256
0.1863(4) 0.1238(5) 0.7850(7)
0.5373 0.8637 0.9630
0.6827(5) 0.2541(6) 0.8167(8)
1.1848 1.0178 1.0037
This solution exhibits a singular behavior near the origin which will limit the rate of convergence in h. On uniform meshes the singularity pollutes the solution and the error estimate globally [14]. Since our error estimates do not account for pollution effects, they will perform poorly. If we eliminate pollution errors by solving the problem on highly graded meshes that resolve the singularity and recover optimal converge rates, we would obtain reasonable accuracy. We begin by solving (7.2a), (7.2b), (7.3) using uniform meshes having 16, 64 and 256 elements with bi 5 and bi 6 polynomial approximations. Effectivity indices shown in Table 3 exhibit large deviations from unity. In order to improve the accuracy of the error estimations, we construct highly graded rectangular meshes using the vertices ðxi ; yj Þ b b pffiffiffiffi i j xi ¼ pffiffiffiffi ; yj ¼ pffiffiffiffi ; i; j ¼ 0; 1; . . . ; N ð7:4Þ N N and solve the problem (7.2a), (7.2b), (7.3) with N ¼ 16, 64, 256. We used bi 5 and bi 6 polynomial approximations with b ¼ 2 and b ¼ 3, respectively. Errors and global effectivity indices are displayed in Table 3. Effectivity indices for graded meshes are much closer to unity than those for uniform meshes. This example shows that the present theory applies to more general situations than those stated previously.
8. Conclusion We have developed asymptotically exact a posteriori error estimates for two-dimensional fourth-order boundary value problems based on the even–odd dichotomy of Babuska and Yu with tensor product bi p finite element approximations. Although, both even- and odd-degree a posteriori error estimators are computationally simple, the even-degree estimator does not involve neighbors communication which is more efficient for parallel implementation and appears to converge to the true error faster as the mesh is h-refined. Finally, we showed that both even- and odd-degree estimators converge to the true error using a more general and efficient hierarchical finite element approximations. Computational results show that the even–odd estimators are asymptotically correct under more general conditions than those indicated by the present theory. For instance, results from Example 3 indicate that they apply on graded quadrilateral meshes in the presence of a singularity. Extension of the analysis to hexahedral meshes and higher-order elliptic equations is straight forward. Adjerid et al. [4] showed that even–odd estimators converged to the true error for second-order problems with both Dirichlet and Neumann boundary conditions on some
2558
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
structured and moderately distorted quadrilateral meshes. We expect this to continue to hold for fourthorder problems. Convergence analyses for nonlinear problems and more general meshes including adaptively refined and/or unstructured quadrilateral, triangular and tetrahedral meshes are less obvious and will be investigated in the future. Likewise, nothing is known about the convergence of even–odd error estimates on curvilinear elements where a theory is needed.
Acknowledgements This Research was supported by the National Science Foundation (grant no. ASC 9720227).
References [1] M. Abramowitz, I.A. Stegun (Eds.), Handbook of Mathematical Functions, Dover, New York, 1965. [2] S. Adjerid, M. Aiffa, J.E. Flaherty, Hierarchical finite element bases for triangular and tetrahedral elements, Comput. Meth. Appl. Mech. Engrg. 190 (2000) 2925–2941. [3] S. Adjerid, I. Babuska, J.E. Flaherty, A posteriori error estimation for the finite element method-of-lines solution of parabolic problems, Math. Model. Meth. Appl. Sci. 9 (1999) 261–286. [4] S. Adjerid, B. Belguendouz, J.E. Flaherty, A posteriori finite element error estimation for diffusive systems, SIAM J. Sci. Comput. 21 (1999) 728–746. [5] S. Adjerid, J.E. Flaherty, L. Krivodonova, A posteriori error estimation for discontinuous Galerkin solutions of hyperbolic problems, Comput. Meth. Appl. Mech. Engrg. 191 (2002) 1097–1112. [6] S. Adjerid, J.E. Flaherty, L. Krivodonova, A posteriori discontinuous Galerkin error estimation for two-dimensional hyperbolic problems, 2000, in preparation. [7] S. Adjerid, J.E. Flaherty, P.K. Moore, Y. Wang, High-order adaptive methods for parabolic systems, Physica D 60 (1992) 94–111. [8] S. Adjerid, J.E. Flaherty, Y. Wang, A posteriori error estimation with finite element methods of lines for one-dimensional parabolic systems, Numer. Math. 65 (1993) 1–21. [9] M. Ainsworth, J.T. Oden, A procedure for a posteriori error estimation for h-p finite element method, Comput. Meth. Appl. Mech. Engrg. 101 (1992) 72–96. [10] M. Ainsworth, J.T. Oden, A unified approach to a posteriori error estimation using element residual methods, Numer. Math. 65 (1993) 23–50. [11] M. Ainsworth, J.T. Oden, A Posteriori Error Estimation in Finite Element Analysis, Wiley, New York, 2000. [12] I. Babuska, T. Strouboulis, S.K. Gangaraj, Guaranteed computable bounds for the exact error in the finite element solution. Part I: one-dimensional model problem, Comput. Meth. Appl. Mech. Engrg. 176 (1999) 51–79. [13] I. Babuska, T. Strouboulis, S.K. Gangaraj, Guaranteed computable bounds for the exact error in the finite element solution. Part II: bounds for the energy norm of the error in two-dimensions, Int. J. Numer. Meth. Eng. 47 (2000) 427–475. [14] I. Babuska, T. Strouboulis, A. Mathur, C.S. Upadhyay, Pollution-error in the h-version of the finite-element method and the local quality of a-posteriori error estimators, Finite Elem. Anal. Des. 17 (1994) 273–321. [15] R. Becker, R. Rannacher, A feedback Approach to error control in finite element methods: basic analysis and examples, EAST– WEST J. Numer. Math. 4 (1996) 237–264. [16] P. Carnevali, R.V. Morric, Y. Tsuji, B. Taylor, New basis functions and computational procedures for p-version finite element analysis, Int. J. Numer. Meth. Eng. 36 (1993) 3759–3779. [17] P. Grisvard, Elliptic Problems in Nonsmooth Domains, Pitman, Boston, 1985. [18] L.G. Leal, Laminar Flow and Convective Transport Processes, Butterworth-Heinemann, Boston, 1992. [19] J.T. Oden, G.F. Carey, Finite Elements, Mathematical Aspects, Prentice-Hall, Englewood Cliffs, 1983. [20] J.T. Oden, S. Prudhomme, Goal-oriented error estimation and adaptivity for the finite element method, Comput. Math. Appl. 41 (2001) 735–756. [21] M. Ortiz, E.A. Repetto, H. Si, A Continuum model for kinetic roughening and coarsening in thin films, J. Mech. Phys. Solids 47 (1999) 697–730. [22] M. Paraschivoiu, J. Preraie, Y. Maday, A.T. Patera, Fast bounds for outputs of partial differential equations, in: J. Borggaard, J. Burns, E. Cliff, S. Schreck (Eds.), Computational Methods for Optimal Design and Control, Birkhauser, Boston, 1998. [23] L.A. Segel, Mathematics Applied to Continuum Mechanics, Dover, New York, 1987.
S. Adjerid / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2539–2559
2559
[24] E. S€ uli, P. Houston, Finite element methods for hyperbolic problems: a posteriori error analysis and adaptivity, in: I.S. Duff, G.A. Watson (Eds.), The State of the Art in Numerical Analysis, Clarendon Press, Oxford, 1997, pp. 441–471. [25] B. Szabo, I. Babuska, Finite Element Analysis, John Wiley, New York, 1991. [26] R. Verfurth, A Review of a Posteriori Error Estimation and Adaptive Mesh Refinement Techniques, Teubner-Wiley, New York, 1996. [27] D. Yu, Asymptotically exact a-posteriori error estimator for elements of bi-even degree, Math. Numerica Sinica 13 (1991) 89–101. [28] D. Yu, Asymptotically exact a-posteriori error estimator for elements of bi-odd degree, Math. Numerica Sinica 13 (1991) 307–314. [29] A. Zangwill, Some causes and a consequence of epitaxial roughening, J. Cryst. Growth 163 (1996) 8–21.