INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING Int. J. Numer. Meth. Engng. 47, 427– 475 (2000) Guaranteed computable bounds for the exact error in the nite element solution—Part II: bounds for the energy norm of the error in two dimensions† T. Strouboulis1 , I. Babuska2;∗ and S. K. Gangaraj1 1 Department of Aerospace Engineering; Texas A&M University; College Station; TX 77843-3141; U.S.A. 2 Texas Institute for Computational and Applied Mathematics; The University of Texas at Austin; Austin; TX 78712; U.S.A. SUMMARY This paper addresses the computation of guaranteed upper and lower bounds for the energy norm of the exact error in the nite element solution. These bounds are constructed in terms of the solutions of local residual problems with equilibrated residual loads and are rather sharp, even for coarse meshes. The sharpness of the bounds can be further improved by employing a few iterations of a relatively inexpensive iterative scheme. The main result is that the bounds are guaranteed for the energy norm of the exact error, unlike the bounds which have been proposed in [13, 14] which are guaranteed only for the energy norm of the error with respect to an enriched (truth-mesh) nite element solution. Copyright ? 2000 John Wiley & Sons, Ltd. KEY WORDS: a posteriori error estimation; nite element method; guaranteed error bounds; upper and lower error bounds 1. INTRODUCTION In this paper we address the problem of computation of upper and lower bounds EuUFE , and EuLFE , for the energy norm of the exact error |||uEX − uFE ||| , namely EuLFE 6|||uEX − uFE ||| 6EuUFE (1) where uEX and uFE are the exact and nite element solutions of the Laplacian in two dimensions, respectively. The main tools for computing the bounds are solutions of local residual problems of the form: ∗ Correspondence to: I. Babuska, TICAM, The University of Texas at Austin, TAY 2-400, Austin, TX 78712, U.S.A. E-mail: babuska@brahma.ticam.utexas.edu † Dedicated to the memory of Dr. Richard H. Gallagher Contract=grant Contract=grant Contract=grant Contract=grant Contract=grant sponsor: sponsor: sponsor: sponsor: sponsor: U.S. Army Research Oce; contract=grant number: DAAL03-G-028 National Science Foundation; contract=grant number: MSS-9025110 U.S. Oce of Naval Research; contract=grant numbers: N00014-96-1-0021 and N00014-96-1-1015 U.S. Oce of Naval Research; contract=grant number: N00014-90-J-1030 National Science Foundation; contract=grant numbers: DMS-91-20877 and DMS-95-01841 CCC 0029-5981/2000/020427–49$17.50 Copyright ? 2000 John Wiley & Sons, Ltd. Received 17 February 1999 428 T. STROUBOULIS, I. BABUSKA AND S. K. GANGARAJ ! Find ê Q !; uFE ∈ Q! such that EQ ! B! ê Q !; uFE ; v = R!; uFE (v) ∀v ∈ Q! (2) where ! is the typical subdomain in which the residual problem is posed, Q! is the solution EQ space, B! is the bilinear form over !, and R!; uFE is an equilibrated local residuum. For the one-dimensional version of this paper, see [1] and for previous work on the equilibrated residual problems, see [2–11]. Letting Q! = U! , the energy space in the subdomain, we get the exact ! indicators |||ê U !; uFE |||! , and we can use them to get the guaranteed upper bound rP ! 2 |||uEX − uFE ||| 6E uUFE = |||ê U (3) !; uFE |||! ! def p where ||| · ||| = B (·; ·), is the energy norm over . Here E uUFE is the exact upper bound or upper estimator, which is not computable and is often approximated by employing Q! = S! ⊂ U! , where S! is a discrete subdomain space obtained by the h, p or hp restriction of the nite element space in the subdomain, and we have rP EuSFE = |||ê S!;! uFE |||2! 6E uUFE (4) ! We will call EuSFE! the computed estimate. It can be proved that EuSFE! , is a guaranteed upper estimate for the error with respect to an enriched or truth mesh solution uFE e S |||uFE e − uFE ||| 6EuFE (5) Note however that in many cases we have EuSFE 6|||uEX − uFE ||| 6EuUFE (6) and hence EuSFE cannot be relied upon as an upper bound for |||uEX − uFE ||| , which is the quantity S of interest and not |||uFE e − uFE ||| . In [1] and in this paper we show that EuFE can be corrected by realizing that r P S! 2 ! (7) E uU = (EuS )2 + |||ê U !; u − ê !; u |||! FE FE FE ! FE and by constructing a computable upper estimate S! ! |||ê U !; uFE − ê !; uFE |||! 6Eeˆ S!;!u (8) FE to get the computable upper bound |||uEX − uFE ||| 6EuUFE s = (EuSFE )2 + P ! Eeˆ S!;!u 2 (9) FE A lower bound for the energy norm can be obtained from the observation that, for any function ’ in the energy space U R(’) R(v) 6 sup = |||R|||−1; = |||uEX − uFE ||| |||’||| v∈U |||v||| Copyright ? 2000 John Wiley & Sons, Ltd. (10) Int. J. Numer. Meth. Engng. 47, 427–475 (2000) GUARANTEED COMPUTABLE BOUNDS FOR THE EXACT ERROR—PART II 429 where R is the global residuum. However, unless ’ is close to the exact error uEX −uFE , this bound could be too small and hence useless. To construct an admissible ’, we rst construct the global Q Q! indicator ê Q uFE patching together the subdomain error indicator functions, namely ê uFE |! = ê !; uFE . Q Q However ê uFE 6∈ U, is not admissible because ê uFE is discontinuous at the vertices and the edges across the subdomains. Nevertheless, we can construct a gap function eˆ Qu , such that êQ uFE + FE eˆ Qu ∈ U, is admissible and we have the lower bound FE EQ uFE = R ê Q uFE + eˆ Q u FE |||ê Q uFE + eˆ Q ||| u 6|||uEX − uFE ||| (11) FE Another lower bound may be obtained by exploiting the orthogonality properties of the error and the error indicators to obtain the exact lower estimator r P = (E uUFE )2 − |||eˆ Uu |||2! 6|||uEX − uFE ||| (12) EU uFE ! FE and its computable version r P E SuFE = (EuSFE )2 − |||eˆ Su |||2! 6E U uFE 6|||uEX − uFE ||| (13) The lower estimate (12) was obtained from the identity r P |||uEX − uFE ||| = (EuSFE )2 − |||eˆ Su |||2! + |||qEX |||2 (14) where qEX ∈ U, is the exact equilibrium correction which is obtained from P ∀v ∈ U B (qEX ; v) = − B! eˆ Su ; v (15) ! FE ! ! FE FE Employing now nite element approximation qFE of qEX , which is computed from the same space as uFE , we have |||qEX |||2 = |||qFE |||2 + |||qEX − qFE |||2 (16) Note that qFE can be computed by utilizing the already available factorized global stiness matrix. We can now repeat the above steps (9) and (13) to construct computable bounds for |||qEX −qFE ||| , namely E SqFE 6|||qEX − qFE ||| 6EqUFE Using (16) and (17) into (14) we get the enhanced bounds r P ẼuU;FE(1) = (EuSFE )2 − |||eˆ Su |||2! + |||qFE |||2 + (EqUFE )2 ! and ẼuL;FE(1) = FE r P (EuSFE )2 − |||eˆ Su |||2! + |||qFE |||2 + (E SqFE )2 Copyright ? 2000 John Wiley & Sons, Ltd. ! FE (17) (18) (19) Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 430 T. STROUBOULIS, I. BABUSKA AND S. K. GANGARAJ Figure 1. The polygonal domain for the model problem: (a) The domain , its boundary @ which consists of the Dirichlet parts D1 ; D2 ; D3 (shown with thick gray line), and the Neumann part N (shown with thin line); (b) the coarse regular mesh of curvilinear quadrilaterals which is used for the computation of the nite element solution Figure 2. Examples of partition of the mesh into patches: (a) The patches coincide with the elements; (b) the patches consist of the vertex patches for the re-entrant corners, and the elements that remain after the vertex patches for the re-entrant corners have been formed Repeating the above procedure in a recursive way we obtain an iteration for improving the bounds. Each iteration employs the factorization of the global stiness matrix which was obtained during the computation of the nite element solution and the factorization of the local stiness matrices which are obtained during the solution of the local residual problems. At the end of each iteration we can compute the reliability interval Ru(m) = FE m) ẼuU;FE(m) − ẼuL;( FE |||uFE ||| (20) As we will see below, a few iterations are sucient to produce dramatic improvement in the accuracy of the bounds and the size of the reliability interval. Following this Introduction, in Section 2 we give the description of the boundary value problem, and the local residual problems. In Section 3 we give the construction of the computable guaranteed upper bound, in Section 4 we give the computable guaranteed lower bound and in Section 5 we give the recursion for improving the bounds. 2. THE MODEL PROBLEM, ITS FINITE ELEMENT APPROXIMATION, THE RESIDUAL EQUATION AND THE EQUILIBRATED LOCAL RESIDUAL PROBLEMS We will employ the mixed boundary-value problem for the Laplacian as our model problem. Let be the polygonal domain shown in Figure 1(a) with boundary ≡ @ ; = D ∪ N ; D ∩ N = ∅, with D = D1 ∪ D2 ∪ D3 , as shown in Figure 1(a). The model problem reads: Find uEX which satises − uEX = 0 uEX | 1 D ∪ @uEX @n N Copyright ? 2000 John Wiley & Sons, Ltd. 3 D = 0; =0 in (21) uEX | D2 = 1 (22) (23) N Int. J. Numer. Meth. Engng. 47, 427–475 (2000) GUARANTEED COMPUTABLE BOUNDS FOR THE EXACT ERROR—PART II The variational formulation of the general mixed boundary-value problem reads n o def def p Find uEX ∈ U ; u = u|kukU = B (u; u)¡∞; u| D = u such that Z Z Z def def B (uEX ; v) = ∇uEX · ∇v = L (v) = fv + gv ∀v ∈ U ; 0 431 (24) N In the example employed throughout the paper we used zero loads f ≡ 0; g ≡ 0, zero Dirichlet 1D = 0; u| 3D = 0, and constant non-zero Dirichlet boundary boundary conditions on 1D and 3D , u| 2 1 2 2D = 1, where , D , D , 3D , and N , are shown in Figure 1(a). condition on D , u| def Let Th = {j }nelem j = 1 , denote a nite element mesh, which is a subdivision of into nonnelem S = overlapping curvilinear quadrilaterals, namely j ; j ∩ k = ∅ for j 6= k. Here we will j=1 employ the coarse regular mesh shown in Figure 1(b) and the meshes which are obtained from the nested subdivision this mesh. By regular mesh we mean a mesh for which the intersection of the boundary of any two elements, @j ∩ @k ; j 6= k is either empty, or a vertex, or an entire edge for both the elements j , and k . Given a mesh Th we will construct the corresponding set of mapped bi-p nite element (trial) solutions n o p def (25) Sh;p u ( ) = u ∈ C 0 ( ) | u| ◦ F ∈ Ŝ ∀ ∈ Th ; u| D = u ⊆ U ; u and the corresponding space of test-functions Sh;p 0 ( ) ⊆ U ; 0 . Here F : ˆ 7→ is the transformation def which maps the master square ˆ = (−1; 1)2 into the curvilinear quadrilateral (in the results of p this paper F was constructed using the blending function approach given in [12]), Ŝ is the bi-p element space over ˆ (here we employed p = 2). We will denote the nite element solution by uS hp and we will obtain it by solving the following discrete problem: Find uS hp ∈ Sh;p u ( ) such that (26) B uS hp ; v = L (v) ∀vh ∈ Sh;p 0 ( ) In general we have uS hp 6= uEX . This can be veried by computing the residuals corresponding to uS hp , namely, the interior residuals def r;uS p = f + uS hp | (27) h and the jumps in the normal derivative T (∇uS hp |k − ∇uS hp |j )(x)·n (x); x ∈ = @j @k ; j 6= k def x∈⊆ N (28) J; uS p (x) = 2(g(x) − ∇uS hp (x)·n N (x)); h 0; x∈⊆ D T where n is the unit normal on = @j @k , which points into k , and n N is the exterior unit normal on N . We will denote the exact error in uS hp by def p = uEX − u p eSEX Sh h Copyright ? 2000 John Wiley & Sons, Ltd. (29) Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 432 T. STROUBOULIS, I. BABUSKA AND S. K. GANGARAJ We will also employ the error relative to an enriched nite element solution uS p+kn , namely h=2 e p+k Sh=2 n p Sh def = uS p+kn − uS hp (30) h=2 p+k Here uSp+k is the nite element solution from Sh=2 n ; u ( ), the set of trial solutions of degree p+k on h=2n the mesh Th=2 which is obtained from T , using n uniform nested subdivisions. The enriched nite n h S p+kn element solution uS p+kn is often called the truth-mesh solution, and the error e Sh=2 p , the truth-mesh h=2 h error; see for example [13,14]. S p+kn h=2 . Substituting uEX = eSEX Let us now give the variational problems satised by eSEX p , and e p p +uS p , S h h h h into the variational statement B (uEX ; v) = L (v), for all v ∈ U ; 0 , we get that the exact error eSEX p , h is the solution of the variational problem: Find eSEX p ∈ U ; 0 such that h def p ; v) = R ; u p (v) = L (v) − B (u p ; v) B (eSEX Sh S h ∀v ∈ U ; 0 h (31) This variational problem is called the residual equation, andp+k the linear functional R ;uS p : S n h U ; 0 7→ R, the global residuum. Similarly, letting uS p+kn = uS hp + e Sh=2 p , into the variational statement h=2 for uS p+kn , namely h h=2 p+k ∀v ∈ Sh=2 n ; 0 ( ) B (uS p+kn ; v) = L (v) h=2 (32) S p+kn is the solution of the variational problem: it can be seen that e Sh=2 p p+k Sh=2 n Find e S p h h p+k ∈ Sh=2 n ; 0 ( ) such that S p+kn ; v = R ; uS p (v) B e Sh=2 p Noting that B (uS hp ; v) = P ∈Th h p+k ∀v ∈ Sh=2 n ; 0 ( ) h B (uS hp ; v), and integrating by parts the term B (uS hp ; v) we get ! Z P 1 r; uS p v + J; uS p v R ; uS p (v) = h h h h ∈Th ∈Th ⊆ @ 2 n o def p def B (w; w)¡∞; w|@∩ D = 0 7→ R, given by Here R;UNEQ uS p : U; 0 = w| ||w||U = h Z Z P 1 def (v) = r v + J; uS p v R;UNEQ ; uS p uS p h h h ⊆ @ 2 P (33) R;UNEQ uS p (v| ) = P Z (34) (35) is the unequilibrated element residuum, namely the linear functional which corresponds to the loads r; uS p and 12 J; uS p , ⊆ @. We call R;UNEQ uS p unequilibrated because, in general, we have h h h Z R;UNEQ u p (1) = S h Copyright ? 2000 John Wiley & Sons, Ltd. r; uS p + h P 1 ⊆ @ 2 Z J; uS p 6= 0 (36) h Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 433 GUARANTEED COMPUTABLE BOUNDS FOR THE EXACT ERROR—PART II Next we construct a new set of element residua Z Z Z P P def EQ; p UNEQ J; u p v = R (v) + sg(; ); uS p v R; u p (v) = r; uS p v + S h S ⊆ @ h ⊆ @ h (37) h where J; u p = 12 J; uS p + sg(; ); uS p S h h (38) h and sg(; ) = ± 1, such that sg(; ) + sg(; ∗ ) = 0, for every interior edge = @ ∩ @∗ , and sg(; ) = 0, for ⊆ @ , and ; uS p ∈ L2 (), is the edge correction for the edge , which will be h specied below. Note that J; u ∗ S p h + J; u p = J; uS p S ∀ = @ ∩ @∗ h h and hence we have an exact splitting of the global residuum, namely P EQ; p R; u p (v| ) ∀v ∈ U ; 0 R ; uS p (v) = S ∈Th h (39) h Following Ladeveze [4] the edge corrections ; uS p will be constructed such that the orthogonality h condition p p R;EQ; u p (v) = 0 ∀v ∈ Sh; 0 () = def S h n p v v | ◦ F ∈ Ŝ ; v|@∩ o =0 D (40) is satised in each element (see also [4–10] and Appendix I for the construction). Note that for each element with @ ∩ D = ∅, we have 1 ∈ Sh;p 0 (), and (40) implies the equilibrium condition Z Z P (1) = r + J; u p = 0 (41) R;EQ;p p ; u u p S S h ⊆ @ h S h We will now introduce the local residual problems which will be used in the a posteriori error def estimation. Let h = {!j }npatches denote a partition of the mesh Th into non-overlapping patches j=1 of elements = npatches S j=1 ! j ; ! j = nel(! S j) k=1 ! k j T !j !k = ∅; j 6= k (42) ! where nel(!j ), is the number of elements in the patch !j , and k j , is the kth element in !j . Figures 2(a) and 2(b) show two examples of partition of the mesh Th of Figure 1(b) into two sets of patches. Figure 2(a) shows a partition for which the patches are the same as the elements. Figure 2(b) shows a partition which employs the vertex patches for the re-entrant corners (each of these patches consists of the elements connected to the vertex of each re-entrant corner), and that rest of the patches coincide with the elements that remain after the vertex patches have been formed. The a posteriori error estimates will be based on the following patch n o residual problems: def p def EX Find ê !j ; u p ∈ U!j ;0 = v ||v||U!j = B!j (v; v)¡ ∞; v|@!j ∩ D = 0 such that S h def P EQ; p p R;EQ; B!j ê EX !j ; u p ; v = R!j ; u p (v) = u p (v| ) S h Copyright ? 2000 John Wiley & Sons, Ltd. S h ∈Th ⊆ !j S ∀v ∈ U!j ;0 (43) h Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 434 T. STROUBOULIS, I. BABUSKA AND S. K. GANGARAJ In the case that @!j ∩ D = , (43) is a mixed boundary value problem and it has a unique solution. In the cases that @!j ∩ D = ∅, (43) is a Neumann problem; in this case however R!EQ;p j ; uS p satises h P EQ;p EQ;p EX R; uS p (1) = 0, and ê !j ; u p exists and is unique up to the consistency condition R!j ; uS p (1) = h S h ∈Th ⊆ !j h a constant, which will be xed by imposing the condition Z ê EX !j ; u p = 0 S !j (44) h We will refer to (43) as the patch residual problem in !j ; in the case that !j coincides with an element, we will call (43) the element residual problem. The functions ê EX !j ; u p will be called the S h patch residual error indicator functions, and the quantities EX EX !j ;u p = ||ê !j ; u p ||U!j S S h (45) h the exact patch error indicators. Summing the squares of the indicators, we get the exact patch residual estimator s 2 r P P def EX 2 E h ; u p = = ||ê EX (46) EX !j ; u p ||U! !j ;u p S S !j ∈ h h S !j ∈ h h h j EX EX Below we will prove that ||eSEX p ||U 6E h ; u p , which means that E h ; u S h S h estimator for ||eSEX p ||U , the energy norm of the exact error. p h is a guaranteed upper h EX The exact indicators EX !j ; u p , and the exact estimator E h ; u S ! Th=2j n S h p h are not computable and must be approximated. Let denote a nite element mesh in !j which is obtained by employing n uniform nested subdivisions of the restriction of Th in !j , and n o p+k def ! p+k v ∈ C 0 (!j ) v |!j ◦ F!j ∈ Ŝ ∀!j ∈ Th=2j n ; v|@!j ∩ D = 0 (47) Sh=2 n ;0 (!j ) = p+k the restriction of the test space Sh=2 n ; 0 ( ) in !j . We will approximate the exact error indicator S p+kn h=2 function ê EX !j ; u p , by the computed error indicator function ê !j ; uS p , which is the solution of the S h h following discrete problem: S p+kn p+k Find ê !h=2 j ; uS p ∈ Sh=2n ; 0 (!j ) such that h S p+kn EQ;p B!j ê !h=2 j ; uS p ; v = R!j ; u p (v) S h In the cases that @!j ∩ h p+k ∀v ∈ Sh=2 n ;0 (!j ) (48) S p+kn D = ∅, we will x ê !h=2 j ; uS p by imposing the additional condition h Z p+k Sh=2 n !j Copyright ? 2000 John Wiley & Sons, Ltd. ê !j ; uS p = 0 (49) h Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 435 GUARANTEED COMPUTABLE BOUNDS FOR THE EXACT ERROR—PART II We will also employ the computed patch residual error indicators S p+kn S p+kn h=2 !h=2 j ; uS p = ||ê !j ; uS p ||U!j h (50) h and the computed patch residual estimator s S p+kn def = E h=2 h ;uS p P !j ∈ h h S p+kn S p+kn 2 !h=2 j ; uS p s P = !j ∈ h h S p+kn S p+kn 2 ||ê !h=2 j ; uS p ||U! h S p+kn ||U 6E h=2 , and hence E h=2 Below we will prove that ||e Sh=2 p h ;u p h ;u S h for ||e p+k Sh=2 n p Sh h S (51) j is a guaranteed upper estimator p h ||U , the energy norm of the truth-mesh error. We will also show that in many cases S p+kn S p+kn EX p ||U 6E ||e Sh=2 ||U 6E h=2 6||eSEX p h ; u h; u p S h S h h (52) p h S p+kn , to be a guaranteed upper estimator for ||eSEX and hence we cannot rely E h=2 p ||U . h ;u p S h h 3. THE UPPER BOUND FOR THE ENERGY NORM OF THE EXACT ERROR EX We will now prove the main results about the estimates E h; u p S h omit the subscript uS hp from the various quantities. S p+kn and E h=2 . Below we will often h; u p S h Theorem 3.1 (Upper bounds). Consider the setting and the notations given in Section 1. We have EX p ||U 6E ||eSEX h ; u p h S h (53) and S p+kn S p+kn ||e Sh=2 ||U 6E h=2 p h; u S h p h (54) Further; we have S p+kn EX !h=2 j ;uS p 6!j ; u h S p h (55) and hence S p+kn E h=2 6E EX h; u h; u p S Copyright ? 2000 John Wiley & Sons, Ltd. h S p h (56) Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 436 T. STROUBOULIS, I. BABUSKA AND S. K. GANGARAJ Proof. Employing the denition of the energy norm using duality, the residual equation and the exact splitting of the global residuum, we get B eSEX p ;v p ||U = sup ||eSEX h P R (v) = sup 6 sup v∈U ; 0 ||v||U v∈U ; 0 h ||v||U v∈U ; 0 !j ∈ h |R!EQ;p (v|!j )| j (57) ||v||U It follows that v EQ;p 2 u (v| ) R !j !j ||v||U! u P P (v) R!EQ;p j j EX t 6 sup ||eS p ||U 6 sup h ||v||U!j ||v||U ||v||2U! v∈U ; 0 !j ∈ h !j ∈ h v∈U!j ;0 (58) j where we have used the Cauchy–Schwarz inequality. From the denition of ê EX !j , as the solution EQ;p of B!j (ê EX !j ; v) = R!j (v), for all v ∈ U!j ;0 , we get def ||L(U!j ;0 ;R) = sup ||R!EQ;p j v∈U!j ;0 B!j (êEX R!EQ;p (v) !j ; v) j = sup = ||êEX !j ||U!j ||v||U!j ||v||U!j v∈U!j ;0 and substituting into (58) we get (53). S p+kn S p+kn ||U 6E h=2 The proof of the bound ||e Sh=2 p h ;u h S p h (59) follows exactly the same steps (57)–(59) where the p+k p+k supremums are taken over Sh=2 n ;0 ( ) and Sh=2n ;0 (!j ) instead of U ; 0 and U!j ;0 respectively. In summary, we have v 2 u EQ;p p+k u P (v) R!EQ;p Sh=2 R (v) n j t 6 supv∈S p+kn (!j ) (60) ||e S p ||U = sup h=2 ;0 h ||v||U ||v||2U! !j ∈ h v∈S p+k ( ) h=2n ;0 j S p+kn we get and using the denition of ê !h=2 j S p+kn B!j ê !h=2 ;v S p+kn j sup = ||ê !h=2 ||U!j j ||v||U!j v∈S p+k (!j ) R!EQ;p (v) j = sup ||v|| p+k U!j v∈S (!j ) h=2n ;0 (61) h=2n ;0 S p+kn p+k EX Let us now prove that !h=2 j ;uS p 6!j ;u p . Noting that Sh=2n ;0 (!j ) ⊂ U!j ;0 , and subtracting (48) from S h h (43) we get the orthogonality B!j ê EX !j ; u S p h S p+kn − ê !h=2 j ; uS p ; v = 0 h p+k ∀v ∈ Sh=2 n ;0 (!j ) (62) and hence S p+kn h=2 2 2 EX ||ê EX !j ; u p ||U!j = ||ê !j ; uS p ||U!j + ||ê !j ; u S h Copyright ? 2000 John Wiley & Sons, Ltd. h S S p+kn p h 2 − ê !h=2 j ; uS p ||U! j (63) h Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 437 GUARANTEED COMPUTABLE BOUNDS FOR THE EXACT ERROR—PART II p+k p+k p+k Sh=2 Sh=2 Sh=2 n n n p+k where we have used that B!j ê EX !j ; u p − ê !j ; uS p ; ê !j ; uS p = 0, because ê !j ; uS p ∈ Sh=2n ;0 (!j ). This can S h h h h also be written as S p+kn !h=2 j ;uS p s EX !j ;u = h 2 p S h S p+kn h=2 2 − ||ê EX !j ; u p − ê !j ; uS p ||U! S (64) j h h S p+kn EX and we have !h=2 j ;uS p 6!j ;u p . Further, summing (63) over all the patches in h , we get S h P !j ∈ h h 2 ||ê EX !j ; u p ||U! = S j h P P S p+kn !j ∈ h 2 ||ê !h=2 j ; uS p ||U! + j h !j ∈ h S p+kn ||ê EX !j ; u p S h 2 − ê !h=2 j ; uS p ||U! h j (65) which can also be written as S p+kn E h=2 h ;uS p s 2 EX E h ;u = S h p h − P S p+kn !j ∈ h h=2 2 ||ê EX !j ;u p − ê !j ;uS p ||U! S j h h (66) S p+kn EX 6E . and we get that E h=2 h ;u p h ;u p S S h h S p+kn h=2 It is rather easy to prove that the estimate ||eSEX p ||U 6E h ;u S h p h does not hold in general. Counterexample 3.1. Consider the case of the example problem where Th is the mesh shown in Figure 1(b), and let k = 1; n = 3. For this example S3 EX p ||U ¡ ||e p ||U ≈ 0·894||eSEX ETh=8 S h ;u p S S p+kn Hence we may have E h=2 h ;u S for ||eSEX p ||U . h h h S p+kn p h h=2 ¡ ||eSEX p ||U , and E h ;u S h p h cannot be relied to be an upper estimator h We will now get a computable upper bound for ||eSEX p ||U by employing the identity h s EX = E h ;u p S h S p+kn E h=2 h ;u p S h 2 + P S p+kn !j ∈ h and constructing an upper estimate for ||ê EX !j ; u p S h h=2 2 ||ê EX !j ;u p − ê !j ;uS p ||U! S h h j (67) S p+kn − ê !h=2 j ; uS p ||U!j , in each patch !j . We will rst h consider the case that h = Th , in which the patches are the elements. We have the following result: Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 438 T. STROUBOULIS, I. BABUSKA AND S. K. GANGARAJ S p+k h Theorem 3.2 (Computable explicit upper estimator for ||ê EX ; u p −ê ; uS p ||U ). Consider the S h h setting and the notations given in Section 2. Then S p+k ||ê EX ; u p S h − ê ;huS p ||U 6EEXPL p+k S ; eˆ;uh h (68) p h S where def ; eˆ;uh P S p+k h EEXPL p+k = C1 (p + k)||r; u p + ê ; u p ||L2 () + S S S h h p h S The constants C1 and C2 are given by s max |J | C1 (q)6 Ĉ1 (q); K ⊂ @ C2 (p + k)||J; u s C2 (q)6 p S h − @ Shp+k ê ; uS p ||L2 () @n h max |J | Ĉ2 (q) K (69) (70) where K = min(|J | min (J−T J−1 )) (71) J (resp. Jj ) is the Jacobian of the transformation of the element (resp. edge); and min (J−T J−1 ) is the minimum eigenvalue of the matrix J−T J−1 ; and Ĉ1 (q) = sup v∈V ˆ ˆ with V qˆ = q ||v̂||L2 ( ˆ) ; ||v̂||U ˆ Ĉ2 (q) = sup q v∈V ˆ ˆ ||v̂||L2 ()ˆ ||v̂||U ˆ Z q v ∈ U ˆ B ˆ (v; w) = 0 ∀w ∈ Ŝ ; v = 0 (72) def (73) ˆ Proof. Letting the patches !j coincide with the elements of the mesh, we get the exact element error indicator functions ê EX ; u p , as the solutions of the variational problems: S Find ê EX ; uS p h h ∈ U; 0 such that EQ;p B ê EX ; u p ; v = R; u p (v) S S h ∀v ∈ U; 0 (74) h S p+k Similarly, the computed element error indicator functions ê ;huS p , are the solutions of the discrete h problems: Find ê Shp+k ; uS p h ∈ Sh;p+k 0 () such that S p+k B ê ;huS p ; v = R;EQ;p u p (v) S h h ∀v ∈ Sh;p+k 0 () (75) Subtracting (75) from (74) for v ∈ Sh;p+k 0 (), and using the notation S p+k def e S p+k = ê EX ; u eˆ;uh Copyright ? 2000 John Wiley & Sons, Ltd. S p h p h − ê ;huS p (76) h S Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 439 GUARANTEED COMPUTABLE BOUNDS FOR THE EXACT ERROR—PART II S p+k for the exact error in ê ;huS p , we get h B e S p+k ; v = REQ;p eˆ;uh p+k S ;uS p ;eˆ h h p h S S p+k def h (v) = R;EQ;p u p (v) − B (ê; uS p ; v) S ∀v ∈ U; 0 h h (77) and the orthogonality condition B e S p+k ; v = REQ;p eˆ;uh S p+k ;uS p ;eˆ h p S h (v) = 0 ∀v ∈ Sh;p+k 0 () (78) h Further, for each element with @ ∩ D = ∅, from (44) and (49), we have Z Z p+k Z S EX e S p+k = ê ; u p − ê ;huS p = 0 eˆ;uh S p h S h (79) h S p+k It follows that e S p+k , the error in ê ;huS p , is the solution of the variational problem: eˆ;uh h p h S p+k Find e S p+k ∈ V such that ; 0 eˆ;uh p h S B e S p+k ; v = REQ;p p+k S ;uS p ;eˆ h h eˆ;uh p S h where p+k def = V ; 0 p+k ∀v ∈ V ; 0 (v) (80) Z p+k v ∈ U B (v; w) = 0 ∀w ∈ S0 (); v = 0 if @ ∩ D =∅ (81) From the duality, we get REQ;p S p+k (v) ||e S p+k ||U = eˆ;uh p S h sup p+k v∈V ; 0 ;uS p ;eˆ h h ||v||U = ||REQ;p S p+k ;uS p ;eˆ h ||L(U; 0 ;R) (82) h S p+k Further, integrating by parts the term B (ê;huS p ; v) in (77) and using (37) we get h Z REQ;p p+k S ;uS p ;eˆ h h (v) = r; uS p + ê h Shp+k ; uS p h v+ P ⊂ @ Z J; u S p h − @ Shp+k ê; uS p @n h v (83) From the triangle inequality in (83) we have Z Z @ Shp+k Shp+k P EQ;p r; uS p + ê; uS p v+ J; u p − ê; uS p v: |R p+k (v)|6 S S h @n h h h h ;uS p ;eˆ ⊂ @ h and employing the Cauchy–Schwarz inequality we get P @ Shp+k S p+k ê; uS p ||v||L2 () |REQ;p S p+k (v)|6r; uS p + ê h 2 ||v||L2 () + J; uS p − @n L () L2 () h h h ;uS p ;eˆ h ⊂ @ (84) h Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 440 T. STROUBOULIS, I. BABUSKA AND S. K. GANGARAJ and from (82) we get ||e S p+k ||U 6 sup eˆ;uh p+k v∈V ; 0 p h S S p+k r; uS p + ê h h P ||v||L2 () ||v||L2 () @ Shp+k + − ê J ; uS p p ; u S @n L2 () ||v||U L2 () ||v||U h h ⊂ @ (85) and letting C1 (p + k) = sup p+k v∈V ; 0 ||v||L2 () ; ||v||U C2 (p + k) = sup p+k v∈V ; 0 ||v||L2 () ||v||U (86) and we obtain S p+k ||e S p+k ||U 6C1 (p + k)r; uS p + ê;huS p eˆ;uh h h p h S L2 () + P ⊂ @ C2 (p + k)J; u p S h − @ Shp+k ê; uS p @n L2 () h (87) The constants C1 , and C2 , are given in terms of the corresponding constants Ĉ1 , and Ĉ2 , in the master element , ˆ as shown below. Employing that q ||v||L2 () 6 max |J | || v̂||L2 ( ˆ) ; q ||v||L2 () 6 max |J | || v̂||L2 ()ˆ (88) where v ∈ U , J = @x=@x̂, is the Jacobian matrix and |J |, is the determinant of the Jacobian matrix, and p ||v||U ¿ K || v̂||U ˆ (89) K = min(|J |min (J−T J−1 )) (90) where and (88) and (89) are proven in Appendix II; we get s C1 (q)6 max |J | Ĉ1 (q); K Copyright ? 2000 John Wiley & Sons, Ltd. s C2 (q)6 max |J | Ĉ2 (q) K (91) Int. J. Numer. Meth. Engng. 47, 427–475 (2000) GUARANTEED COMPUTABLE BOUNDS FOR THE EXACT ERROR—PART II 441 where Ĉ1 (q) = sup q v∈V ˆ ˆ ||v̂||L2 ( ˆ) ; ||v̂||U ˆ Ĉ2 (q) = sup q v∈V ˆ ˆ ||v̂||L2 ()ˆ ||v̂||U ˆ (92) Let us now give some details about how the constants Ĉ1 (q), and Ĉ2 (q) are computed. Computation of Ĉ1 and Ĉ2 : The computation of Ĉ1 (q), and Ĉ2 (q) will be done by noting that p p Ĉ2 (q) = lim Ĉ2 (q) Ĉ1 (q) = lim Ĉ1 (q); p→∞ p→∞ (93) where p Ĉ1 (q) = sup v̂∈Vq=ˆ p and Vˆ q=p = def || v̂||L2 ( ˆ) ; || v̂||U ˆ p Ĉ2 (q) = sup v̂∈Vq=ˆ p || v̂||L2 ()ˆ || v̂||U ˆ Z p q v ∈ Ŝ B ˆ (v; w) = 0 ∀w ∈ Ŝ ; v=0 (94) (95) ˆ p p In order to compute Ĉ1 (q), and Ĉ2 (q) we need to construct a basis for Vˆ q=p. Let i ; i = 1; 2; : : : ; (p + 1)2 , be the bi-p element shape functions and let Z i+1 ˆ ˜ Z 1 ; i = 1; 2; : : : ; (p + 1)2 − 1 (96) i = i+1 − 1 ˆ We then have Z ˆ ˜i = 0 (97) Group the above functions ˜i , into lower order Li = ˜i ; i = 1; 2; : : : ; (q + 1)2 − 1 (98) j = 1; 2; : : : ; (p + 1)2 − (q + 1)2 (99) and higher order ˜ H j = i+(q+1)2 −1 ; and compute the matrices D and E with entries L Dij = − B ˆ (H i ; j ); Ejk = B ˆ (Lj ; Lk ) (100) def for i = 1; 2; : : : ; N q=p = (p + 1)2 − (q + 1)2 , and j; k = 1; 2; : : : ; (q + 1)2 − 1, and determine F = E−1 . Now let ! 2 2 (q+1) P (q+1) P H ˜ Dij Fjk Lk ; i = 1; 2; : : : ; N q=p; (101) i = i + k=1 Copyright ? 2000 John Wiley & Sons, Ltd. j=1 Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 442 T. STROUBOULIS, I. BABUSKA AND S. K. GANGARAJ and note that the ˜ i , are orthogonal to any Lm , for m = 1; 2; : : : ; (q + 1)2 − 1, because ! 2 2 (q+1) P (q+1) P L H L ˜ Dij Fjk Bˆ (L ; Lm ) Bˆ (i ; m ) = Bˆ (i ; m ) + j=1 k=1 = −Dim + = −Dim + 2 (q+1) P 2 (q+1) P k=1 j=1 2 (q+1) P j=1 Dij Fjk 2 (q+1) P Di;j k=1 | k ! Ekm ! Fjk Ekm {z =0 } jm R where jm is the Kronecker delta. Further, we also have ˆ ˜ i = 0, from (97). Hence ˜ i ; i = 1; 2; : : : ; N q=p; form the basis of Vˆq=p. Therefore, any v ∈ Vˆq=p, can be expressed in the form v̂ = q= p NP i=1 ai ˜ i (102) and we have ||v̂||2L2 () ˆ =a A T q= p p Aq= ij = a; and ||v̂||2Uˆ = aT Bq=pa; Bijq=p = Z ˆ ˆ ˜i ˜j (103) @˜i @˜j @˜i @˜j + @x̂1 @x̂1 @x̂2 @x̂2 Similarly, we have ||v̂||2L2 ()ˆ = aT Cq=pa; Z Cijq=p = (104) Z ˆ ˜i ˜j (105) It follows that p aT Aq=pa (Cˆ1 )2 = max T q=p ; q= p a B a a∈RN and hence p Cˆ1 = r max kq=p; k=1; :::; N q=p p aT Cq=pa (Cˆ2 )2 = max T q=p q= p a B a a∈RN p Cˆ2 = r max k=1; :::; N q=p kq=p (106) (107) where kq=p and kq=p are the eigenvalues of the generalized eigenvalue problems Aq=py = kq=pBq=py; Cq=py = kq=pBq=py (108) p p We computed Cˆ1 ; and Cˆ2 using the Harwell Library routine EA12AD (see [15]). For each q, p p we employed a suciently high p to obtain convergence of Cˆ1 (q); and Cˆ1 (q) to 12 signicant digits. Table I shows the converged values of Cˆ1 (q) and Cˆ2 (q) for q = 0; 1; 2; : : : ; 8, and the values of p for which the convergence was achieved. Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 427–475 (2000) GUARANTEED COMPUTABLE BOUNDS FOR THE EXACT ERROR—PART II 443 Table I. Guaranteed upper bound for the energy norm of the error: The constants Cˆ1 (q) and Cˆ2 (q) for q = 1; 2; : : : ; 8, and the corresponding values of p which was used to compute them q Cˆ1 (q) p Cˆ2 (q) p 0 1 2 3 4 5 6 7 8 0·636619772367 0·630301952842 0·609347491519 0·568383672094 0·530296680623 0·500349587345 0·464059359874 0·430324985456 0·400678456146 8 9 12 14 16 18 20 22 24 0·816496580927 0·810315743719 0·781457077283 0·756306743339 0·690556011657 0·621895745689 0·558985734593 0·498567556786 0·426746908655 11 14 15 17 19 21 24 26 28 ˆ ˆ Note that the values R of C1 (0) and C2 (0) can also be computed analytically. Any function v̂ ∈ Uˆ with zero average ˆ v̂ = 0 can be expanded in a cosine series in the form ∞ P v̂ = m; n=0 m+n6=0 am; n cos m n (1 + x̂1 ) cos (1 + x̂2 ) 2 2 (109) We then have ||v̂||2Uˆ = 2 4 ∞ P m; n=0 m+n6=0 (m2 + n2 )a2m; n ; ||v̂||2L2 () ˆ = ∞ P m; n=0 m+n6=0 a2m; n (110) and hence ∞ P 2 Cˆ1 (0) = sup 0 v∈ ˆ V̂ˆ ||v̂||2L2 () ˆ ||v̂||2Uˆ = 4 sup 2 am; n m; n=0 m+n6=0 ∞ P m; n=0 m+n6=0 a2m; n (111) (m2 + n2 )a2m; n The above supremum is achieved for the cases m = 0 and n = 1, or m = 1 and n = 0, and we get 2 Cˆ1 (0) = ≈ 0·636619772367 (112) which agrees with the value given in Table I. Further, in order to compute Cˆ2 (0), let us consider the edge ˆ1 which joins the vertices (−1; −1) and (1; −1) in the master element. We have ∞ 2 ∞ ∞ 2 ∞ P P P P 1 2 nam; n am; n = ||v̂||L2 (ˆ1 ) = m=1 n=1 m=1 n=1 n ∞ ∞ ∞ ∞ P ∞ P P 1 P 2 2 2 P = 6 n a n2 a2m; n (113) m; n 2 6 m=1 n=1 m=1 n=1 n n=1 Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 444 T. STROUBOULIS, I. BABUSKA AND S. K. GANGARAJ ∞ P where we used Cauchy–Schwarz inequality, and the identity 1=n2 = 2 =6. Therefore, we have n=1 ∞ P ||v̂||2L2 (ˆ1 ) 2 Cˆ2 (0) = sup ||v̂||2Uˆ 0 v∈ ˆ V̂ˆ m; n=0 m+n6=0 ∞ P 2 2 6 sup 3 am; n m; n=0 m+n6=0 n2 a2m; n (114) (m + n2 )a2m; n The above supremum is achieved for the case m = 0, and for all n and we get r 2 ˆ ≈ 0·816496580927 C2 (0) = 3 which agrees with the value given in Table I. Above we have constructed the guaranteed upper estimate v 2 2 u P U; Shp+k def u S p+k EX EX EEXPL ||eS p ||U 6ETh ; u p 6ETh ; u p = t EThh ; u p + p+k S S h S h S h ∈Th h ; eˆ;uh (115) (116) p S h This estimate is valid for the special case that h ≡ Th , and that the solutions of the element S p+k residual problems ê ;EXu p , are approximated by the p method, namely ê ;huS p . S h h Let us now consider the more general case where the element residual problems are solved by S p+kn , is computed. As before, we have the h method, namely the case for which the estimate ETh=2 h; u p S s S p+kn ETh=2 h; u = ETEX h ;u p S 2 p S h h + P ∈Th h S p+kn ||ê EX ; u − ê ;h=2uS p ||2U p S h (117) h and ||ê EX ; u S P S p+kn p h − ê ;h=2uS p ||2U = h 0 ∈Th=2 n ||ê EX ; u S S p+kn p h − ê ;h=2uS p ||2U0 (118) h where Th=2 n is the mesh obtained from n nested subdivisions of the element . S p+kn h=2 We will now construct an upper estimate for ||ê EX ; u p − ê; uS p ||U , exactly as we constructed the S h h upper estimate for ||uEX − uS hp ||U in (116) above. We have v u p+k+‘ 2 p+k p+k+‘ P S Sh=2 u S n n n h=2 EX EX 2 h=2 + || ê − ê ||ê EX p+k 0 p+k = t E 0 ; u p − ê ; uS p ||U 6E S n ;eˆ ;eˆ ||U0 S S h h n ; eˆ h=2 Th=2 n ; eˆ h=2 Th=2 n 0 ∈Th=2 n (119) Using the explicit upper estimate S p+k+‘ n h=2 EXPL ||êEX 0 ; eˆ − ê 0 p+k+‘ ; eˆ ||U0 6E S n (120) 0 ; eˆ;uh=2p S h Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 427–475 (2000) GUARANTEED COMPUTABLE BOUNDS FOR THE EXACT ERROR—PART II S p+k+‘ n S p+kn h=2 where we have used the subscript ê in êEX 0 ; eˆ , and ê0 ; eˆ (119) we get S p+kn h=2 EX ||êEX ; u p − ê; uS p ||U 6E S h h p+k h=2n S ; eˆ Th=2 n 445 indicate ê h=2 , and employing (120) into 6E p+k+‘ U; Sh=2 n S (121) p+k n ; eˆ h=2 Th=2 n where E p+k+‘ U; Sh=2 n S p+k n ; eˆ h=2 Th=2 n v u p+k+‘ 2 P EXPL 2 S E S p+k+‘ =u t E h=2n S p+kn + def ; eˆ h=2 Th=2 n n 0 ∈Th=2 n (122) 0 ; eˆ;uh=2p S h The guaranteed upper bound for the case that h 6= Th follows similar steps as depicted in Figure 3. We have thus constructed the guaranteed upper estimator for ||e SEX p ||U . h Theorem 3.3 (Guaranteed upper bound for ||eSEX p ||U ). Consider the above setting and the notah tions given above. We have p+k p+k+‘ U; Sh=2 n =Sh=2n def EX p ||U 6E 6E = ||eSEX p h ; uS h ; uS p h h h v u p+k 2 2 P U; S p+k+‘ u Sh=2n h=2n p+k t E ; u p + E S n h ! S Th=2j n ;eˆ!h=2 j !j ∈ h h (123) Let us now give an example which illustrates the above constructions. Example 1 (The upper bound for the energy norm of the exact error). We considered the model problem and the mesh Th shown in Figure 1 and biquadratic elements p = 2. We computed the energy norm of the exact error ||eSEX 2 ||U , by employing an overkill mesh Th;ovk , shown in h Figure 4 which is constructed such that the relative error ||eSEX 2 ||U =||uS 2 ||U , is 0·01 per cent accurate. h h This mesh was obtained by rening the original mesh Th , uniformly three times then by rening the elements adjacent to the re-entrant corners ve more times. The overkill error was computed using the GFEM (see [16]) space of degree p = 8 and by including the rst term of the asymptotic expansion of the solution of the Laplacian at each reentrant corner. , by solving the element residual problems using We also computed the exact upper bound ETEX h the rened meshes shown in Figure 5 and p = 8; these meshes are the restriction of the mesh shown in Figure 4 which was used in the overkill. We have ||eSEX 2 ||U h ||uS 2 ||U = 29·17 per cent6 h ETEX h = 48·22 per cent ||uS 2 ||U h Next, we obtained the element error indicators by employing the p and the h methods to solve the element residual problems. Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 446 T. STROUBOULIS, I. BABUSKA AND S. K. GANGARAJ Figure 3. Guaranteed upper bound for the energy norm of the error: Schematic depiction of the steps involved in the construction of the guaranteed upper bound for ||eEX p ||U : (a) The mesh Th which is employed in the computation of the S h ! nite element solution uS p ; (b) The partition of Th into the patches !j and the meshes Th=2j n in the computation of the h S p+k n error indicator functions ê !h=2 j; u S p h ! ! ; (c) The mesh Th=21 n employed in the patch !1 ; (d) The elements 0 ∈ Th=21 n which are used for the solution of the element residual problems S p+k (a) Element residual problems solved using the p method: We computed EThh , for p + k = 3; S p+k 4; : : : ; 8. Figure 6 shows the graphs of the eectivity indices of the estimators EThh =||eSEX 2 ||U , S p+k S p+k h S p+k h ||U¡EThh ¡||eSEX and ||eSh2 ||U =||eSEX 2 ||U , for p + k = 3; 4; : : : ; 8. Note that, we have ||e 2 2 ||U . S h h h S3 h n (b) Element residual problems solved using the h method: We computed EThh=2 , for n = 0; 1; 2; 3. S3 S3 n n h=2 ||U = Figure 7 shows the graphs of the eectivity indices EThh=2 =||eSEX 2 ||U , and the ratios ||e 2 S h S3 h n h=2 ||eSEX ||U ¡ 2 ||U , for n = 0; 1; 2; 3. Note that as expected from Theorem 3.1, we have, ||e 2 S S3 h h n EThh=2 ¡||eSEX 2 ||U . h Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 447 GUARANTEED COMPUTABLE BOUNDS FOR THE EXACT ERROR—PART II Figure 4. The overkill mesh Th;ovk used for the computation of ||eEX p ||U : This mesh was obtained by rening the original S h mesh Th uniformly three times, then by rening the elements adjacent to the re-entrant corners ve more times Figure 5. The overkill meshes employed for the computation of the exact element error indicators EX !j : These meshes are the restrictions of the overkill mesh of Figure 4 in the elements of the original coarse mesh Let us also construct the computable guaranteed upper bounds by employing the solutions of the element residual problems obtained by the p and the h methods. (a) Element residual problems solved using the p method: Figure 8 compares the graphs of the U; Shp+k eectivity indices ETh as expected we have EX EX =||eSEX 2 ||U , and ETh =||eS 2 ||U , for p + k = 3; 4; : : : ; 8. We note that h h U; Shp+k EX ||eSEX 2 ||U ¡ET ¡ET h h h (b) Element residual problems solved using the h method: In this case we employed k = 1 and ‘ = 1 and computed the upper bound E p+k p+k+‘ U; Sh=2 n =Sh=2n Figure 9 compares the graphs of the eectivity indices ETEX =||eSEX 2 ||U , h h 3 U; Sh=2 n , which we will denote by ETh . EX ||eSEX 2 ||U =||eS 2 ||U h h and U; S 3 n ETh h=2 =||eSEX 2 ||U , for n = 0; 1; 2; 3 and as expected, we have h 3 U; Sh=2 n EX ||eSEX 2 ||U ¡ET ¡ET h h h Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 448 T. STROUBOULIS, I. BABUSKA AND S. K. GANGARAJ Figure 6. The upper bound for energy norm of the error: Approximate solution of the element residual problems using S p+k h the p method. Graphs of the eectivity index of the exact estimator ETEX=||eEX 2 ||U , the computed estimator ET h S h h p+k h S2 h and the ratio of the energy norm of the truth-mesh error to the energy norm of the exact error ||e S p+k p + k = 3; 4; : : : ; 8. Note that ET h h S p+k S =||eEX 2 ||U , ||U =||eEX || , S2 U S h versus h is an upper estimator for ||e h2 ||U , but it underestimates the energy norm of the exact S h error ||eEX 2 ||U S h S p+kn p+k U; Sh=2 n Let us also mention that the bounds E EX , E h=2 , and E h h h S p+kn E h=2 (Xint ), h depend on the particular choice of U; S p+kn the edge corrections, i.e. we have E EX = E EX (Xint ), and E h h=2 (Xint ) where Xint is the h h vector of edge-corrections for the interior edges. In Appendix I we give the Ladeveze’s method for constructing the edge corrections which are unique up to a set of free constants, one for each interior vertex of the mesh. The constants can be chosen equal to zero and we call this choice as the Bank–Weiser equilibration [10] or the constants can be chosen such that they minimize a weighted norm at each vertex (see [4–6]) and we call this choice as Ladeveze equilibration. Another factor which aects these bounds is the order of equilibration. As we have seen in Appendix I, we can construct the edge corrections Xint , such that the orthogonality condition p ∀v ∈ Sh;0 () EQ;p R;u p (v) = 0 S h (124) holds for each element. However, a weaker condition is sucient for the existence of the estimate, namely the consistency condition EQ;0 R;u p (1) = 0 S (125) h In general, we can employ a qth-order equilibration for which EQ; q R;u p (v) = 0 S h Copyright ? 2000 John Wiley & Sons, Ltd. q ∀v ∈ Sh;0 () (126) Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 449 GUARANTEED COMPUTABLE BOUNDS FOR THE EXACT ERROR—PART II Figure 7. The upper bound for energy norm of the error: Approximate solution of the element residual problems using S3 n h=2 the h method. Graphs of the eectivity index of the exact estimator ETEX=||eEX =||eEX 2 ||U , the computed estimator ET 2 ||U , h S 3 Sh=2 n and the ratio of the energy norm of the truth-mesh error to the energy norm of the exact error ||e p+k h S2 h n = 0; 1; 2; 3; 4. Note that we have ||e S S h h S2 p+k S ||U ¡ET h h h h ||U =||eEX || , S2 U for h ¡||eEX 2 ||U S h with 06q6p. We will express the dependence of the bounds on the order and the type of BW(q) Lad(q) (resp. E EX; ) is the equilibration by introducing additional indices, for example E EX; h h exact upper estimator based on Bank–Weiser (resp. Ladeveze) equilibration of order q. Example 2 (The inuence of the choice of equilibration on the upper bound). Figure 10(a) S p+k ; EQ(q) EQ(q) h (resp. Figure 10(b)) shows the graphs of ETEX; =||eSEX =||eSEX resp. 2 ||U , and ETh 2 ||U , h h h U; Shp+k ;EQ(q) =||eSEX for EQ = BW, and EQ = Lad, and for q = 1; 2. From these results it can ETh 2 ||U h EQ(q) be seen that for the example problem considered here ETEX; is practically independent of h S p+k ; EQ(q) U; Shp+k ; EQ(q) , and the upper bounds ETh the choice of equilibration and EThh signicantly by the choice of the equilibration. Let us now summarize the conclusions in this section: are not inuenced 1. We have the guaranteed upper estimator EX p ||U 6E ||eSEX h = r P h !j ⊂ h 2 ||êEX ;u p ||U! S h j (127) given by 2. However, the computable version of E EX h S p+kn E h=2 h Copyright ? 2000 John Wiley & Sons, Ltd. s = P !j ⊂ h S p+kn ||ê;uh=2S p ||2U! h j (128) Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 450 T. STROUBOULIS, I. BABUSKA AND S. K. GANGARAJ Figure 8. The upper bound for energy norm of the error: Approximate solutions of the element residual problems using p+k h U; S the p method. Graphs of the eectivity indices ET S p+k U; S ET h h → p+k h U; S =||eEX 2 ||U for p + k = 3; 4; : : : ; 8. Note that ET h ETEX , h h h ¿ETEX and h as p + k → ∞ Figure 9. The upper bound for energy norm of the error: Approximate solutions of the element residual problems 3 U; Sh=2 n using the h method. Graphs of the eectivity indices ET 3 U; Sh=2 n ET h → 3 U; Sh=2 n =||eEX 2 ||U for n = 0; 1; 2; 3. Note that ET S h ETEX , h as h h=2n h ¿ETEX and h →0 is not a guaranteed upper bound for the energy norm of the exact error ||eSEX p ||U and we h could have S p+kn EX p ||U 6E ¡||eSEX E h=2 h h h Copyright ? 2000 John Wiley & Sons, Ltd. (129) Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 451 GUARANTEED COMPUTABLE BOUNDS FOR THE EXACT ERROR—PART II Figure 10. The inuence of the choice of equilibration on the upper bound: Eect of the type and the order of EX; EQ(q) the equilibration. (a) Graphs of ET h p+k U; S ET h h ; EQ(q) =||eEX || S2 U S p+k h =||eEX 2 ||U , and ET S h h ; EQ(q) EX; EQ(q) =||eEX 2 ||U ; (b) graphs of ET S h h =||eEX 2 ||U and S h for EQ = BW, and EQ = Lad, and for q = 1; 2. Note that for the considered example, the eect h of the order and the type of equilibration is practically negligible 3. Employing the identity v u p+k 2 p+k u S n P Sh=2 n EX EX 2 + || ê − ê E h = t E h=2 ;uS p ||U! ;u p h !j ⊂ h Copyright ? 2000 John Wiley & Sons, Ltd. S h h j (130) Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 452 T. STROUBOULIS, I. BABUSKA AND S. K. GANGARAJ and constructing a computable upper estimate for S p+kn p+k+‘ U; Sh=2 n h Th=2j n ;eˆ!h=2 j h=2 ||êEX ;u p − ê;uS p ||U!j 6E S h S ! (131) p+k n we get the guaranteed computable upper estimator v 2 u p+k 2 p+k p+k+‘ u Sh=2n P U; S U; Sh=2 n h=2n EX t p ||U 6E ||eSEX E E 6E = + p+k h h h S h !j ⊂ h ! n (132) Th=2j n ;eˆ!h=2 j 4. There is no signicant dierence between the values of the estimators for various types and orders of equilibration. 4. THE LOWER BOUND FOR THE ENERGY NORM OF THE EXACT ERROR The residual error indicator functions êEX !j ; u p , can be put together into a function dened over , S which we will denote by êEX uS p h ∈ L ( ), namely 2 h def EX êEX u p |!j = ê!j ; u p ; S S h !j ∈ h (133) h Note that êEX u p is smooth in the interior of elements but it may have jump discontinuities at the S h vertices and across the interior edges of the mesh. We will denote these discontinuities by eˆEX u p G S h def EX EX = <êEX u p = = êu p (x + 0n ) − êu p (x − 0n ) S S h h S (134) h eˆEX u p S and we will call G h , the gap of êEX u p on the edge . Here , is an interior edge and n , is the unit S h eˆEX u p S normal assigned to . For the edges on the boundary ⊂ @ , we let G h ≡ 0. Let ETh denote the set of all edges in the mesh and let us dene the broken energy space with prescribed gaps on all the edges, namely ( def UBroken = Ex G eˆu p eˆEX u p S 2 v ∈ L ( ) v| ∈ U;0 ; <v= = − G h ; ∈ ETh ) (135) T s h h Broken EX eu p and we have êEX u p ∈ U eˆ Ex . Note also that the exact error eS p does not have gaps, namely G S ≡ EX S h 0, for all ∈ ETh and G h p Tus h h p ∈ U ; 0 eSEX h Copyright ? 2000 John Wiley & Sons, Ltd. h (136) Int. J. Numer. Meth. Engng. 47, 427–475 (2000) GUARANTEED COMPUTABLE BOUNDS FOR THE EXACT ERROR—PART II 453 The dierence between the exact error and the exact error indicator def EX EX p eEX ˆEX = eS − êu p u p S h S h (137) h ∈ UBroken and eEX satises the following orthogonality condition: Then we have eEX eˆ Ex ˆEX ˆEX G u p u p p Tus h h S h S h Lemma 4.1. Consider the above setting and the notations. Then P ∈Th B eEX ; v =0 EX ˆ ∀v ∈ U ; 0 u p (138) S h Proof. For any v ∈ U ; 0 , we have P P p ;v B eEX ; v = B eSEX − B!j êEX !j ; u p ; v ˆEX u p ∈Th h S S !j ∈ h h = R ; uS p (v) − h P !j ∈ h h R!EQ;p j ; uS p (v) = 0 (139) h From the above orthogonality it follows that eEX has the minimal energy norm as shown in ˆEX u p S h the following result. Lemma 4.2. Consider the above setting and notations. Then we have rP ∈Th rP ||eEX ||2 ¡ ˆEX U u p S h ∈Th ||eˆEX ||2U u ∀eˆEX ∈ UBroken Ex u G eˆ p h p h S S p Tus h h (140) Proof. Note that ∈ U ; 0 = eEX ˆEX − eˆEX u u p S h (141) p h S and hence P ∈ Th ||eˆEX ||2U = u p S h P ∈ Th ||eEX − ||2U = ˆEX u p S h P ∈ Th where we have used the orthogonality condition (4:6) that and omitting ||||2U , (140) follows. ||eEX ||2 + ||||2U ˆEX U u p (142) S h P ∈ Th B eEX ; v = 0, for all v ∈ U ; 0 , ˆEX u p S h Lemma 4.3. Consider the above setting and the notations. Then p ||U = ||eSEX h Copyright ? 2000 John Wiley & Sons, Ltd. r P EX )2 − (E ||eEX ||2 ˆEX U h ∈ Th u p (143) S h Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 454 T. STROUBOULIS, I. BABUSKA AND S. K. GANGARAJ EX Proof. We have eEX = eSEX p − êu p , and hence ˆEX u p S h S h P ∈ Th h 2 ||eEX ||2 = ||eSEX p || U−2 ˆEX U u p h S h P ∈Th P EX 2 + B êEX ||êEX u p ; eS p u p ||U S h h S ∈Th (144) h and noting that P P EX EX EX EX EX 2 = B êEX ; e REQ;p (eSEX p p ) = R (e p ) = B (e p ; e p ) = ||e p || u p U S S S S S S ∈Th h h h ∈ Th h h h (145) h we get (143). Combining (140) and (143) we get the following guaranteed lower bound for ||eSEX p ||U . h Theorem 4.1. Consider the above setting and the notations. Then r P def L EX p ||U EX ) = ( (E )2 − ||eˆEX ||2U 6 ||eSEX ∀eˆEX ∈ UBroken E Ex p e ˆ ; u h h; u p u u u G eˆ S h S p h S h p h ∈Th h S p Tus h h p h S (146) and L L eˆEX 6 E E h h; u u p S h p S h EX p eEX ˆEX = ||eS ||U u p (147) h S h Let us now give one particular method of construction of the gap function eˆEX . Let XTh denote u p h S the set of vertices in the mesh. For X ∈ XTh , note that êEX uS p (X ) is well dened because h EX 1+ (!j ) êEX u p |!j = ê!j ; u p ∈ H S where ¿0, and êEX !j ; u S p h S h (148) h is continuous at the vertices of !j . Figure 11(b) depicts the error indicator functions along the common edges of the shaded elements shown in Figure 11(a) and illustrates eˆEX u p the gap at the vertex node X . Therefore, we dene the vertex gap GX; , at the vertex X of , S h def = êEX GX; eˆEX u p | (X ) − u p S h S h nX 1 P êEX | X (X ); nX j=1 uS hp j X ∈ XTh (149) X is the set of elements which are connected to X . Let e(1) be given by piecewise where {Xj }nj=1 ˆEX u p S h mapped bilinear functions | e(1) ˆEX u p S h =− 4 P k=1 eˆEX u p GX ; ’Xk S h k (150) where ’Xk ∈ S01 (), is the mapped bilinear basis function corresponding to the vertex Xk . Figure , along the common edges of the shaded elements shown in 11(c) depicts the gap function e(1) ˆEX u p S h Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 455 GUARANTEED COMPUTABLE BOUNDS FOR THE EXACT ERROR—PART II Figure 11. The lower bound for energy norm of the error: Construction of the gap functions for the lower bound. (a) The mesh indicating the four shaded elements for which the error indicator functions and the gap functions are (1) shown; (b) the error indicator functions along the edges AX and BX ; (c) the gap function EX , along the common edges eˆu (2) of the shaded elements; (d) the gap function eˆEX u p; S h 1, 2 p h S along the common edges of the shaded elements (1) Figure 11(a). The function êEX u p + eˆEX is continuous at the vertices but has gaps across the edges. S u p h S h Let ( | e(2) ˆEX u p S h def ∈ U; <eˆEX +(1) = = u eˆEX p h S u p (1) v ∈ U;0 v| = − 12 <êEX uS p + eˆEX = ; ⊂ @ h S u p ) (151) S h h then we have = e(1) + e(2) ∈ UBroken e(1)+(2) Ex ˆEX ˆEX ˆEX G eˆ u p u p u p h h h S Copyright ? 2000 John Wiley & Sons, Ltd. S S p Tus h h (152) Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 456 T. STROUBOULIS, I. BABUSKA AND S. K. GANGARAJ L and we can now compute the lower bound E h; u S L E h; u S e(1)+(2) = ˆEX ;OPT p h u p S h (2) eˆEX u p S h e(1)+(2) . We can also optimize this bound by ˆEX p h u p S h L max (1)+(2) E eˆEX h ; uS p u p h ∈ U EX S (1) h ; <eˆu + = eˆEX p h S (153) u p S h 1=2 (), This maximum can also be expressed in terms of a set of unknown edge-functions ∈ H0;0 (2) ∗ by letting eˆ EX ∈ U; <eˆEX +(1) = , be such that for every = @ ∩ @ , we have u u ;E eˆ EX p h p h S S u p S h (1) | = − <êEX (2) u + eˆEX = ; eˆEX ; p h u p S (1) (2) | ∗ = − (1 − )<êEX u + eˆEX = eˆEX ; u p S p h u p S h u p S S h (154) S h h and ; v = − B (1) ; v ∀v ∈ U;0 B (2) eˆEX ; eˆEX u p u p h h S (155) S and we have L E h; u p S h L (1)+(2) = max E h; u eˆEX ;OPT u p p S h ∈ C 0 () ∈E Th S h + (2) (1) eˆEX eˆEX ; u p u p h h S (156) S (1) This optimum value is, in general, less than ||eSEX p ||U . Note that we could have also dened EX eˆ u p h S h in terms of the unknown weights at the vertices of the mesh and then the optimal value would (2) , along the common edges coincide with ||eSEX p ||U . Figure 11(d) depicts the gap functions EX eˆ ;1=2 u p h S h of the shaded elements shown in Figure 11(a). Remark 4.1. The solution of the Dirichlet problem (155) exists and is unique. This is veried by (1) (2) 1=2 noting that <êEX u +eˆEX = ∈ H0;0 (); ∈ @ and for ∈ Th . Therefore, eˆEX ;1=2 can be obtained by the p h S u p u p h h S S superposition of the solutions of four Dirichlet boundary value problems for which homogeneous (1) boundary conditions are imposed on three edges and <êEX u +eˆEX = are imposed on the fourth edge. p h S u p S h Let us note that the lower bound L ((1)+(2) ) E eˆEX h ; uS p u p h ê EX !j ; u S p h is based on the exact error indicator functions S h and it is not computable. Let us now consider the case that the computed error indicator Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 457 GUARANTEED COMPUTABLE BOUNDS FOR THE EXACT ERROR—PART II S p+kn Broken functions ê !h=2 j ; uS p are employed. We also note that for S p+k ∈ U S p+k , we have eˆ;uh h h eˆ;u p p h S S GT h h S p+kn |!j = ê EX S p+k |!j − eˆEX !j ; u u eˆ;uh S p h S p h S − ê !h=2 j ; uS p p h (157) h and hence || S p+k ||2U!j = || S p+k ||2U!j + ||ê EX !j ; u eˆ;uh eˆ;uh p h S S p+kn S p h p h 2 − ê !h=2 j ; uS p ||U! j where we have used the orthogonality (3:10) that B!j ê EX !j ; u S p+k Sh=2 n ;0 (!j ) and that (158) h S p+k S p+k |!j ∈ Sh=2 n ;0 (!j ). eˆ;uh p p h S p+kn − ê !h=2 j ; uS p ; v = 0, for all v ∈ h Using (158) in the result in Theorem 4.1, we get the S h computable lower bound L E h ; uS p h v u P u S p+kn 2 p ||U || S p+k ||2U 6 ||eSEX S p+k = t E h=2h ; u p − eˆ;uh S p h S eˆ;uh ∈ Th h (159) h p h S Let us now compute the lower bound for the example problem. Example 3. (The lower bound for the energy norm of the error). Figure 12 (resp. Figure. 13) L;Shp+k 3 L; Sh=2 n L EX =||eSEX 2 ||U ) with ETh =||eS 2 ||U for p + k = 3; 4; : : : ; 8 (resp. h h p+k L;Sh=2 n L to indicate the lower bound E S p+k : We n = 0; 1; 2; 3). Here used the notation ETh h; u p compares ETh =||eSEX 2 ||U (resp. ETh h S h eˆ;uh p h S note that we have L; Shp+k E Th ¡ ETLh ; 3 L; Sh=2 n ETh ¡ETLh i.e. the lower bounds obtained from the approximate solutions of the element residual problems are smaller than ETLh . Figures. 12 and 13 show that the lower bounds can be very small relative to the energy norm of the exact error and hence it is necessary to compute enhanced lower bounds. Let us now summarize the conclusions in this section: 1. We have the guaranteed lower bound L E h; u S p h eˆEX u p h S def r EX 2 P p ||U = ||eˆEX ||2U 6 ||eSEX E h ; u p − u S h ∈Th p h S (160) h ∈ UBroken is the gap function which is constructed such that êEX ∈ U ; 0 . where eˆEX u p + eˆEX eˆ Ex u u p h S G p Tus h h Copyright ? 2000 John Wiley & Sons, Ltd. S h p h S Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 458 T. STROUBOULIS, I. BABUSKA AND S. K. GANGARAJ Figure 12. The lower bound for energy norm of the error: Approximate solutions of the element residual problems by p+k h L; S the p method. Graphs of the eectivity indices of the lower bounds ET have L;S 3 ET h h h =||eEX 2 ||U for p + k = 2; 3; : : : ; 8. Note that we S h = 0 and for p + k65 the lower bound for the energy norm of the exact error is practically useless Figure 13. The lower bound for energy norm of the error: Approximate solutions of the element residual problems by 3 L; Sh=2 n the h method. Graphs of the eectivity indices of the lower bounds ET h L; S 3 ET h h =||eEX 2 ||U for n = 0; 1; 2; 3. Note that we have S h = 0 and for n62 the lower bound for the energy norm of the exact error is practically useless Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 459 GUARANTEED COMPUTABLE BOUNDS FOR THE EXACT ERROR—PART II 2. The computable version of the lower bound is given by L E h; u S S p+k p h eˆ;uh v u S p+k P u 2 n p ||U || S p+k ||2U 6 ||eSEX = t E h=2h ; u p − S p h S eˆ;uh ∈Th h (161) h p h S S p+kn such that êuh=2 where S p+k ∈ UBroken p + S p+k ∈ U ; 0 . eˆ Ex S eˆ;uh G p Tus h h p h S eˆ;uh h p h S 3. The computed lower bound can be very small relative to ||eSEX p ||U , in some cases and we h L L . have E h ; u p S p+k 6 E h ; u p eˆEX u S eˆ;uh h S p h S p h S h 5. ITERATIVE IMPROVEMENT OF THE BOUNDS FOR THE ENERGY NORM OF THE EXACT ERROR In this section we give the construction of the improved upper and lower bounds. Let us rst consider the case that the element residual problems are solved exactly. Let def EX ẽEX u p | = êu p + eˆEX u S S h (162) p h h S denote enhanced error indicator function and let def EX EX qeEX ˜EX = eS p − ẽu p u p S h S h (163) h be its dierence with the exact error eSEX p . For any v ∈ U ; 0 , we have h EX EX B qeEX ˜EX ; v = B eS p ; v − B ẽ ; v u p S h h P p ;v − B êEX ;v = B eSEX u p + eˆEX u h ∈Th p ;v − = B eSEX | h P ∈Th {z S p h h S REQ;p (v) − } P ∈Th B eˆEX ;v u p S h (164) 0 is the solution of the variational problem: It follows that qeEX ˜EX u p S h Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 460 T. STROUBOULIS, I. BABUSKA AND S. K. GANGARAJ Find qeEX ∈ U ; 0 such that ˜EX u p S h P B eˆEX ; v ∀v ∈ U ; 0 B qeEX ˜EX ; v = − u u p p h ∈Th S h (165) S We will call qeEX ˜EX the exact equilibrium correction. We now have: u p S h Lemma 5.1. Consider the above setting and the notations. Then 2 EX 2 p || ||eSEX U = (E h ) − h P ∈Th ||eˆEX ||2U + ||qeEX ||2 ˜EX U u u p p h S (166) S h EX EX Proof. Substituting (162) into (163) and using eˆEX EX = eS p − êu p , we get u p S h S h h EX eˆEX = eˆEX EX − qe˜EX u p h S u p u p h h (167) S S and therefore P ∈Th ||eˆEX ||2U = u p h S P ∈Th 2 ||eˆEX EX ||U + u p S h P ∈Th ||qeEX ||2 ˜EX U u p (168) S h EX = 0, from Lemma 4.1 in which we proved the orthogonality where we have used B eˆEX EX ; qe˜EX u u P that ∈Th B eˆEX EX ; v = 0, for all v ∈ U ; 0 . Using (168) in Lemma 4.3 where we showed that u r P EX )2 − 2 (E ||eˆEX ||eSEX p ||U = EX ||U , we get the proof of (166). h p h p h S S p h S h ∈Th u p S h Sp h Let qe˜EX be the nite element approximation of qeEX , namely, the solution of the following ˜EX u p u p h h S S problem: Find q Shp e˜EX u p ∈ Shp such that S h P Shp ;v = − B eˆEX ;v B qe˜EX u u p p h ∈Th S h S ∀v ∈ Shp (169) By subtracting (169) from (165) we get the orthogonality condition Shp B qeEX ˜EX − qe˜EX ; v = 0 u p u p h h S Copyright ? 2000 John Wiley & Sons, Ltd. ∀v ∈ Shp (170) S Int. J. Numer. Meth. Engng. 47, 427–475 (2000) GUARANTEED COMPUTABLE BOUNDS FOR THE EXACT ERROR—PART II 461 and it follows that Sp Sp 2 2 EX 2 h h ||qeEX ˜EX ||U = ||qe˜EX ||U + ||qe˜EX − qe˜EX ||U u p u p u p u p h h h h S S S (171) S Let Sp def h ũ(0) = uShp + ẽEX u p + qe˜EX Sp S h h (172) u p S h denote the enhanced nite element solution. From (163) and (172) we have Sp (0) h qeEX ˜EX − qe˜EX = uEX − ũS p u p u p h h S (173) h S satises the orthogonality condition and ũ(0) Sp h ;v =0 B uEX − ũ(0) Sp h ∀v ∈ Shp (174) We then have L EX E 6||uEX − ũ(0) || 6E Sp U ; ũ(0) ; ũ(0) h L resp. E ; ũ(0) EX where E ; ũ(0) h h p S h p S h p h h h S (175) p h S is the upper bound (resp. lower bound) which is computed by substituting ũ(0) in place of uShp in Section 3 (resp. Section 4). Using this in (171), we get Sp h 2 2 Shp Shp L EX 2 EX ||qe˜EX ||2U + E 6||q || 6||q ||2 + E EX (0) U e˜ e˜EX U ; ũ ; ũ(0) h u p S p h S h u p u p h h S h S (176) p h S Substituting (176) into (166), we obtain (1) U; (1) p ẼL; h ; u p 6||uEX − uSh ||U 6Ẽ h ; u p S S h (177) h where s (1) def ẼU; h ; uS p = h EX E h; u S 2 p h Copyright ? 2000 John Wiley & Sons, Ltd. − P ∈Th Sp EX h ||eˆEX ||2U + ||qe˜EX ||2U + E u ; ũ(0) p S h u p S h h p h 2 (178) S Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 462 T. STROUBOULIS, I. BABUSKA AND S. K. GANGARAJ is the enhanced upper bound and s 2 2 P Shp L; (1) def 2 2 L EX EX E − || || + ||q || + E Ẽ h ; u p = (0) EX eˆu U U e˜ h; u p ; ũ S S h h u p p h ∈Th S h S h (179) p h S is the enhanced lower bound. We note that 2 2 2 2 P (1) (1) EX L − ẼL; − E = ||eˆEX ||2U ẼU; (0) h ; u p h ; u p = E ;ũ(0) u ;ũ S S h h h h p h S p h (180) p h ∈Th S S is the gap function corresponding to ũ(0) . It follows that the dierence between the where eˆEX u p Sp h S h upper and the lower bounds can be decreased iteratively by decreasing the gaps in the enhanced (1) U; (1) nite element solution. We also note that the enhanced bounds ẼL; h ; u p and Ẽ h ; u p , are based on S h S h the exact solutions of the residual problems êEX u p and hence they are not computable. Following S h similar steps as in [1] we get the computable enhanced bounds and using the symbol ũ(m) , to denote ũ(m) we obtain following recursive relation: Sp h U; S p+kn def Ẽ h ; ũh=2(m) = v u S p+k 2 U; S p+k 2 P u Shp n n 2 2 t E h=2 − || + ||q || + E h ; ũh=2(m+1) p+k || (m) p+k U U S n S h ;ũ eˆ h=2 (m) ∈Th n (181) e˜ h=2 (m) u˜ u˜ for m = 0; 1; : : : ; M − 1 where S p+kn = ũ(m) + ẽũh=2 ũ(m+1) (m) + q Sp Sp h h Shp S (182) p+k n e˜ h=2 (m) u˜ and for m = M , we let U; S p+kn U Ẽ h ;ũ(M ) = E h ;ũh=2 (M ) (183) Similarly a recursive relation can be derived for the lower bound, v u S p+k 2 P L; S p+k 2 p+k L; Sh=2 Shp n def u n n 2 2 Ẽ h ; ũ(m) = t E h=2 − || + ||q || + Ẽ h ;ũh=2(m+1) p+k || (m) p+k U U S S h ;ũ n ∈Th eˆ h=2 (m) u˜ n (184) e˜ h=2 (m) u˜ for m = 0; 1; 2; : : : ; M − 1 and for m = M , we let v u p+k 2 P u S n L L − || S p+kn ||2U Ẽ h ; ũ(M ) = Ẽ h ;ũ(M ) = t E h=2 (M ) h ;ũ ∈Th eˆ h=2 (M ) (185) u˜ Let us now employ these constructions for the example problem. Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 463 GUARANTEED COMPUTABLE BOUNDS FOR THE EXACT ERROR—PART II Figure 14. Iterative improvement of the bounds for the energy norm of the error: Exact solutions of the element residU; (m) L; (m) ual problems. Graphs of the eectivity indices ẼTh =||eEX =||eEX 2 ||U , and ẼTh 2 ||U , for m = 0; 1; 2; 3. Note that both S S h h the bounds practically converge in three iterations to the energy norm of the exact error Example 4 (Iterative improvement of the upper and the lower bounds). Let us rst consider the case that the element residual problems are solved exactly. Figure 14 shows the graphs of (m) L; (m) =||eSEX =||eSEX the eectivity indices ẼU; 2 ||U , and ẼTh 2 ||U , for m = 0; 1; 2; 3. We note that both the Th h h bounds practically converge in three iterations to the energy norm of the exact error. Next we will report the computable improved bounds. As before, we will denote the estimate S p+k obtained by the p-method (resp. h-method) of solution of the element residual problems by EThh S3 n resp. ETh=2 ; recall that these estimates are not upper bounds for the energy norm of the exh U; Sh8 ; (m) act error. Let us denote the corresponding enhanced guaranteed upper bounds by ẼTh U; S 3 ; (m) ẼTh h=8 , respectively and the lower bounds by L; S 8 ; (m) ẼTh h , and Figure 15 compares the graphs of the eectivity indices of L; Sh8 ; (m) and ẼTh h S8 S 8 ; (m) h EX h =||eSEX 2 ||U with ||e 2 ||U =||eS 2 ||U . We note that after three iterations ẼTh S h h h derestimates and U; S 8 ; (m) ẼTh h S 3 ; (m) S3 always un- is the guaranteed upper bound for the energy norm of the exact error. Figure 16 compares the graphs of the eectivity indices ẼTh=8 h L; S 3 ; (m) h L; Sh8 ; (m) , and ẼTh S 8 ; (m) practically converge to the truth-mesh error estimate, for m = 0; 1; 2; 3. Further, ẼThh ||eSEX p ||U h , and L; S 3 ; (m) ẼTh h=8 , respectively. S 8 ; (m) U; Sh8 ; (m) ẼThh =||eSEX =||eSEX 2 ||U , ẼTh 2 ||U , 3 U; Sh=8 ; (m) =||eSEX 2 ||U , ẼTh h =||eSEX 2 ||U , h n h=2 ||U =||eSEX and ẼTh h=8 =||eSEX 2 ||U with ||e 2 2 ||U . As before, the estimates converge to the energy Sh h h norm of the exact error after three iterations. Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 464 T. STROUBOULIS, I. BABUSKA AND S. K. GANGARAJ Figure 15. Iterative improvement of the bounds for the energy norm of the error: Approximate solutions of the elS 8 ; (m) ement residual problems using the p method. Graphs of the eectivity indices ẼTh h and L; S 8 ; (m) ẼT h =||eEX || , h S2 U S 8 ; (m) ẼTh , h for m = 0; 1; 2; 3. Note that after three iterations h S8 U; Sh8 ; (m) verged to the truth-mesh error estimate ||e h2 ||U , and ẼT S h h and U; Sh8 ; (m) =||eEX 2 ||U , ẼT S h L; S 8 ; (m) ẼT h , h h =||eEX 2 ||U , S h have practically con- , is a guaranteed upper bound for the energy norm of the exact error From Figures 14 and 15 it can be seen that the enhanced bounds converge to constant values which are dierent from ||eSEX p ||U . The converged values of the bounds depend on values of p + k h S p+kn and n used for computing the error indicator functions êuh=2 p . S h U; Shp+k ; (m) Figure 17 shows the graphs of the eectivity indices ẼTh ||eSEX 2 ||U , for p + h U; S 3 n ; (m) =||eSEX ẼTh h=2 2 ||U , h L; Shp+k ; (m) =||eSEX 2 ||U and ẼTh h = k = 3; 4; 6; 8; while Figure 18 shows the graphs of the eectivity indices L; S n ; (m) and E˜Th h=2 =||eSEX 2 ||U , for n = 1; 2; 3. It can be seen that 3 h p+k U; S lim E˜Th h ; (m) p+k→∞ p+k = L; S lim E˜Th h p+k→∞ ; (m) p ||U = ||eSEX h (186) and U; S n ; (m) L; S n ; (m) p ||U E˜Th h=2 E˜Th h=2 = lim = ||eSEX lim n n 3 h=2 →0 3 h=2 →0 h (187) The relative dierence between the upper bound and the lower bound will be called the reliability interval S p+kn ; (m) def = RThh=2 p+k U; Sh=2 n ; (m) E˜Th L; S p+k ; (m) h=2n − E˜Th ||uS 2 ||U (188) h Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 465 GUARANTEED COMPUTABLE BOUNDS FOR THE EXACT ERROR—PART II Figure 16. Iterative improvement of the bounds for the energy norm of the error: Approximate solutions of the eleS 3 ; (m) ment residual problems using the h method. Graphs of the eectivity indices E˜T h=8 h 3 L; Sh=8 ; (m) and E˜T h =||eEX || , S2 U h 3 ˜U; Sh=8 ; (m) =||eEX =||eEX 2 ||U , E 2 ||U , S h S3 3 L; Sh=8 h h for m = 0; 1; 2; 3. Note that after three iterations E˜T h=8 , and E˜T to the truth-mesh error estimate ||e 3 Sh=2 n s2 h 3 U; Sh=8 ||U , and E˜T h S h , practically converge is the guaranteed upper bound for the energy norm of the exact error The reliability intervals for the bounds shown in Figures 17 and 18 are shown in Figures 19 and 20 respectively. From Figure 19 (resp. Figure 20) we see that the reliability interval is less than 0·1 for p + k ¿ 6 (resp. n ¿ 2). From the above examples we conclude that 1. The improved upper and lower bounds converge to the energy norm of the exact error only in the case that the element residual problems are solved exactly. 2. In the case that the approximate solutions of the element residual problems are employed, the improved estimates converge to the the truth-mesh error estimate and not the energy norm of the exact error. 3. The guaranteed upper bound for the energy norm of the exact error is obtained by taking into account the explicit estimate correction term as shown above. 4. The improved upper and lower bounds converge after a few iterations, i.e. the reliability interval converges to a constant value. The reliability interval can be decreased by solving the element residual problems accurately. Let us now give examples to show the inuence of the approximate solutions of the residual problems on the reliability interval. Example 5 (Bounds for the energy norm of the error using the patch residual method). Let us group the elements which are connected to the re-entrant corners into patches, !1 ; !2 ; : : : ; !Npat (see Figure 21(a)) and solve the patch-residual problems. We will denote the enhanced upper and lower bounds which are computed by solving patch residual problems shown in Figure 21 by Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 466 T. STROUBOULIS, I. BABUSKA AND S. K. GANGARAJ Figure 17. Iterative improvement of the bounds for the energy norm of the error: Approximate solutions of the element p+k ; (m) h U; S residual problems using the p method. Graphs of the eectivity indices E˜T h p+k ; (m) h L; S ˜ =||eEX 2 ||U , and ET S for p + k = 3; 4; 6; 8, and m = 0; 1; 2; 3 h h =||eEX 2 ||U , S h Figure 18. Iterative improvement of the bounds for the energy norm of the error: Approximate solutions of the element 3 U; Sh=2 n ; (m) residual problems using the h method. Graphs of the eectivity indices E˜T h n = 1; 2; 3, and m = 0; 1; 2; 3 3 L; Sh=2 n ; (m) ˜ =||eEX 2 ||U , and ET S h =||eEX 2 ||U , for S h h E U; Sh=2n ; (m) and EL; Sh=2n ; (m) respectively. Figure 22 shows the graphs of the eectivity indices of the 3 U; S 3 n ; (m) ˜L; Sh=2n ; (m) =||eEX =||eEX improved upper bounds E˜ h=2 2 ||U and the improved lower bounds E 2 ||U 3 3 Th Sh Th Sh for n = 0; 1; 2; 3. By comparing Figure 22 with Figure 18 we note that the reliability interval obtained by the patch residual method is smaller than the reliability interval obtained by the Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 427–475 (2000) GUARANTEED COMPUTABLE BOUNDS FOR THE EXACT ERROR—PART II 467 Figure 19. Iterative improvement of the bounds for the energy norm of the error: Approximate solutions of the element S p+k ; (m) , for p + k = 3; 4; 6; 8, and m = 0; 1; 2; 3. residual problems using the p method. Graphs of the reliability intervals RT h h Note that the reliability interval is less than 0·1 for p + k ¿ 6 Figure 20. Iterative improvement of the bounds for the energy norm of the error: Approximate solutions of the element S 3 n ; (m) , for n = 1; 2; 3, and m = 0; 1; 2; 3. residual problems using the h method. Graphs of the reliability intervals RT h=2 h Note that the reliability interval is less than 0·1 for n ¿ 2 element residual method, i.e. S3 n ; (m) R h=2 h S3 n ; (m) ¡RThh=2 (189) The above example shows that reliability interval can be improved by solving the patch residual problems in the elements which are connected to the re-entrant corners. Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 468 T. STROUBOULIS, I. BABUSKA AND S. K. GANGARAJ Figure 21. Bounds for the energy norm of the error using the patch residual method: Approximate solution of the patch residual problems. The meshes of curvilinear quadrilaterals of degree p = 3 constructed within the elements to solve the p+k ; (m) h=2n U; S patch residual problems and obtain the estimate E h : (a) n = 0; (b) n = 1; (c) n = 2 Figure 22. Iterative improvement of the bounds for the energy norm of the error: Approximate solutions of the patch residual problems as shown in Figures 21(a)–(c). Graphs of the eectivity indices of the improved upper bounds 3 U; Sh=2 n ; (m) E˜ 3 L; Sh=2 n ; (m) ˜ =||eEX 2 ||U and the improved lower bounds E S h h h =||eEX 2 ||U for n = 0; 1; 2; 3, and m = 0; 1; 2; 3 S h Let us now consider dierent methods of constructions of the edge corrections which were mentioned in Section 1. Example 6 (Comparison of the bounds based on dierent methods of construction of the edge corrections). In this example we show that type of splitting of the jumps for the construction of the equilibrated residuals in (37), practically, does not inuence the converged values of the reliability intervals. We considered Ladeveze and Bank–Weiser methods of construction of the edge corrections as in Example 2 and computed the values of the upper and the lower bounds. Figure 23 compares the graphs of the eectivity indices of the improved upper bounds 3 U; S 3 ; (m) ˜L; Sh ; (m) =||eEX E˜ h =||eEX 2 ||U and the improved lower bounds E 2 ||U for the above types of Th Sh Copyright ? 2000 John Wiley & Sons, Ltd. Th Sh Int. J. Numer. Meth. Engng. 47, 427–475 (2000) GUARANTEED COMPUTABLE BOUNDS FOR THE EXACT ERROR—PART II 469 Figure 23. Comparison of the bounds based on dierent methods of constructions of the edge corrections: Approximate solutions of the element residual problems as shown in Figure 2(a). Graphs of the eectivity indices of the improved 3 3 U; S ; (m) ˜L; Sh ; (m) =||eEX upper bounds E˜T h =||eEX 2 ||U and the improved lower bounds ET 2 ||U for q = 1; 2 for Ladeveze and the h S S h h h Bank–Weiser methods of construction of the edge corrections splittings. We note that for m¿3 the eectivity indices of the bounds obtained from the four types of the edge corrections are practically close. In this section we have shown the iterative improvement of the bounds for the energy norm of the error and we note the following: 1. If the local residual problems are solved exactly, the upper and the lower bounds converge to the energy norm of the exact error. If the local residual problems are solved approximately, the lower bound converges to the energy norm of the error with respect to the enhanced solution (truth-mesh error). The guaranteed upper bound is obtained by taking into account the error in the approximate solutions of the residual problems. 2. The converged value of the reliability interval decreases if the local residual problems are solved accurately. We also note that the converged value of the reliability interval obtained using the patch residual problems is smaller than the reliability interval obtained using the element residual problems. 3. The converged values of the upper and the lower bounds are practically independent of the method of construction of the equilibrated element residuals. 6. CONCLUSIONS In this paper we addressed the problem of a posteriori estimation of the bounds for the energy norm of the exact error in the energy norm. The major conclusions of this paper are: Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 470 T. STROUBOULIS, I. BABUSKA AND S. K. GANGARAJ 1. The upper bound for the global energy norm of the error is obtained by employing the exact solutions of the local residual problems. 2. In order to construct guaranteed upper bound based on the computed error indicators, we must construct an upper estimate for the element energy norm of the error approximate solutions of the local residual problems. We constructed an explicit estimate which is suciently accurate for computing the guaranteed upper bound. 3. A lower bound for the global energy norm of the error can be obtained by constructing the ‘gap function’ for the solutions of element residual problems. 4. The improved upper and lower bounds for the global energy norm of the error can be constructed by employing the ‘equilibrium corrections’. The improved upper and lower bounds, the global energy norm of the error, can be improved by employing an iterative procedure which employs the factorized global stiness matrix and the solutions to the element residual problems. The objective of this paper was to illustrate the concepts and the procedure to compute the bounds for the energy norm of the exact error in the nite element solution given in Reference [1], for two-dimensional heat-conduction problems. In forthcoming papers, we will address to construction of bounds for the outputs (temperature, heat-ux at a point, ux-intensity factors, etc.) computed either from the nite element solution or a recovered solution obtained by local averaging (e.g. the ZZ-SPR). We will also address the case of the elasticity problem. APPENDIX I I.1. Construction of the edge corrections for the equilibrated element residuals Let us describe the detailed construction of the edge corrections ; uS p by the Ladeveze method h [2–5] and the Bank–Weiser method [9]. For simplicity we will omit the subscript uS hp . We let = p+1 P i=1 (i) (i) (190) where (i) ; i = 1; 2, are the mapped linear functions on and (i+1) ; i = 2; : : : ; p. are the mapped polynomial functions of degree i. We will construct ; ⊂ @, such that REQ (v) = 0, for all v ∈ p (). This is achieved by determining (i) such that Sh;0 REQ (; j ) = 0 (191) where ; j ; j = 1; 2; : : : ; p + 1, are the mapped element shape functions. Below will employ a special form of (i) such that they satisfy the following orthogonality condition Z 0 if j¡i (i) ; j = (192) 1 if j = i where i = 3; : : : ; p + 1, and j = 1; 2; : : : ; p + 1. This enables us to determine (i) ; i = 1; 2, by solving a small system of equations and (i) ; i = 3; : : : ; p + 1, can be determined explicitly, as shown below. Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 427–475 (2000) GUARANTEED COMPUTABLE BOUNDS FOR THE EXACT ERROR—PART II 471 Figure 24. Construction of the edge corrections for the equilibrated element residuals: An example which depicts the enumeration of the elements jX , and the edges kX in the patch associated with the vertex X nel(X ) Let X denote an interior vertex node, jX j=1 , a set of elements connected to X , and X nedge(X ) the set of edges in the mesh patch connected X (see for example Figure 24). Let k k=1 1 X ∈ Sh be the hat function associated with the vertex X . In the rst step we construct by employing the orthogonality condition that R (v) = 0, for v ∈ Sh;p 0 ( ): R (X ) = nel(X P) j=1 REQ;p (X |jX ) = 0 X and we get the following set of equations: Z P ( ) + sg(kX ; jX )kX X = 0; RUNEQ X X j kX ⊂ @jX j = 1; 2; : : : ; nel(X ) and using (190) together with the orthogonalities in (192) we get Z P (1) (2) (2) UNEQ sg(kX ; jX ) (1) R X (X ) + X = 0; X X + X X j ⊂@jX k kX (193) j k k k j = 1; 2; : : : ; nel(X ) (194) (195) (2) Further, we will also construct (1) X and X such that the following orthogonalities are satised: k k Z Z kX (1) X X = 1; k kX (2) X X = 0 k (196) Therefore, we can now write the system of equations (195) as follows: MX XX = bX Copyright ? 2000 John Wiley & Sons, Ltd. (197) Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 472 T. STROUBOULIS, I. BABUSKA AND S. K. GANGARAJ where b X = n onel(!X ) −RUNEQ (X ) X j j=1 ; M ∈ Rnel; nedge and kX = (1) for k = 1; 2; : : : ; nedge. In the X k example shown in Figure 24 the above system of equations is (X ) −RUNEQ X 1X 1 1 0 0 −1 UNEQ X ( ) −R −1 X X 1 0 0 2 2 = 0 −1 1 0 3X (X ) −RUNEQ X 3 0 0 −1 1 X −RUNEQ ( ) 4 4X (198) X We note that the solution to (197) exists but it is not unique (see [2–5] for details). The Bank–Weiser equilibration involves choosing a particular solution X nedge = 0; X nedge−j =− nedge P k=nedge−j+1 RUNEQ (X ); X k j = 1; 2; : : : ; nedge − 1 (199) The Ladeveze-type of splitting involves choosing a particular solution for (197), which minimizes the ‘2 norm, s X 2 nedge P k X (200) ||X ||‘2 ;{| X |−2 }nedge = X k k=1 | k=1 k | After the linear edge corrections have been determined, the higher-order corrections can be computed from Z i−1 P (j) (j) ; i ; i = 2; : : : ; p + 1; ⊆ @ (201) (i) = − RUNEQ (; i ) − j=1 Note that the higher-order corrections can be determined explicitly and are dened uniquely for each edge. APPENDIX II II.1. Bounds for the norms of the functions in an element def Lemma I.1. Let ˆ = (−1; 1) × (−1; 1) denote the master element and let x = x( x̂1 ; x̂2 ) denote the transformation of the master element co-ordinates ( x̂1 ; x̂2 ) to x and let @x1 @x1 @ x̂2 def @ x̂ 1 (202) J = @x2 @x2 @ x̂1 @ x̂2 denote the Jacobian matrix of the transformation. Then we have r max |J | ||v̂||L2 (ˆ) ||v||L2 () 6 Copyright ? 2000 John Wiley & Sons, Ltd. (203) Int. J. Numer. Meth. Engng. 47, 427–475 (2000) GUARANTEED COMPUTABLE BOUNDS FOR THE EXACT ERROR—PART II and 473 p ||v||U ¿ K ||v̂||Uˆ (204) K = min(|J |min ((J )−T (J )−1 )) (205) where def |J |¿0 is the determinant of the the Jacobian matrix in (202) and min ((J )−T (J )−1 ) is the minimum eigenvalue of (J )−T (J )−1 . Proof. The proof of (203) is obtained by noting that Z Z v(x1 ; x2 ) dx1 dx2 = v̂( x̂1 ; x̂2 )|J | d x̂1 d x̂2 (206) ˆ and hence Z ||v||2L2 () = Z v2 = ˆ Z v̂2 |J |6 max |J | v̂2 (207) ˆ Further 1 @x2 @ x̂1 1 @x1 @ x̂1 = ; =− @x1 |J | @ x̂2 @x2 |J | @ x̂2 1 @x2 @ x̂2 1 @x1 @ x̂2 =− ; = @x1 |J | @ x̂1 @x2 |J | @ x̂1 and we obtain @v @x @ x̂1 @x1 def 1 = ∇v = @v @ x̂1 @x2 @x2 T @v̂ @v̂ def ; and where ∇x̂ v̂ = @ x̂1 @ x̂2 @ x̂2 @v̂ @x1 @ x̂1 = (J )−1 ∇x̂ v̂ @v̂ @ x̂2 @x2 @ x̂2 (208) @ x̂2 @x1 @ x̂1 @x1 (209) @ x̂2 1 @x2 (J )−1 = |J | @ x̂1 − @x2 − is the inverse of the Jacobian matrix (202). By employing (206) and (207), we get Z Z 2 T ||v||U = (∇v) (∇v) = ((J )−1 ∇x̂ v̂)T ((J )−1 ∇x̂ v̂)|J | Z = ˆ ˆ (∇x̂ v̂)T (J )−T (J )−1 (∇x̂ v̂)|J | Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 427–475 (2000) 474 T. STROUBOULIS, I. BABUSKA AND S. K. GANGARAJ Z ¿ ˆ min ((J )−T (J )−1 )|J |(∇x̂ v̂)T (∇x̂ v̂) ¿ K Z ˆ (∇x̂ v̂)T (∇x̂ v̂) (210) where min ((J )−T (J )−1 ) is the minimum eigenvalue of the matrix (J )−T (J )−1 2 2 @x1 @x2 @x2 @x2 @x1 @x1 − + − @ x̂2 @ x̂2 @ x̂1 @ x̂2 @ x̂1 @ x̂2 1 = 2 2 |J |2 @x @x @x @x @x1 @x2 1 1 2 2 − + − @ x̂1 @ x̂2 @ x̂1 @ x̂2 @ x̂1 @ x̂1 (211) ACKNOWLEDGEMENTS The work of T. Strouboulis and S. K. Gangaraj was supported by the U.S. Army Research Oce under Grant DAAL03-G-028, by the National Science Foundation under Grant MSS-9025110 and by the U.S. Oce of Naval Research under Grants N00014-96-1-0021 and Grant N00014-96-11015. The support of Dr. Kailasam Iyer of the U.S. Army Research Oce and Dr. Dick Lau of the Oce of Naval Research is greatly appreciated. The work of I. Babuska, was supported by the U.S. Oce of Naval Research under Grant N00014-90-J-1030 and by the National Science Foundation under Grants DMS-91-20877 and DMS-95-01841. REFERENCES 1. Babuska I, Strouboulis T, Gangaraj SK. Guaranteed computable bounds for the exact error in the nite element solution—Part I: one dimensional model problem. Computer Methods in Applied Mechanics and Engineering 1999; 176:51–79. 2. Babuska I, Strouboulis T, Upadhyay CS, Gangaraj SK, Copps K. Validation of a posteriori error estimators by numerical approach. International Journal for Numerical Methods in Engineering 1994; 37:1073–1123. 3. Babuska I, Strouboulis T, Upadhyay CS, Gangaraj SK. A model study of element residual estimators for linear elliptic problems: the quality of the estimators in the interior of meshes of triangles and quadrilaterals. Computers and Structures 1995; 57:1009–1028. 4. Ladeveze P. Comparaison de modeles de milieux continus, These de Doctorat d’Etat es Science, Mathematiques, L’Universite P. et M. Curie, Paris, 1975. 5. Ladeveze P, Leguillon D. Error estimate procedure in the nite element method and applications. SIAM Journal on Numerical Analysis 1983; 20(3):485–509. 6. Ladeveze P, Maunder EAW. A general method for recovering equilibrating element tractions. Computer Methods in Applied Mechanics and Engineering 1996; 137:111–151. 7. Ainsworth M, Oden JT. A-posteriori error estimation in nite element analysis. Computer Methods in Applied Mechanics and Engineering: Computational Mechanics Advances 1997; 142:1–88. 8. Ainsworth M, Oden JT. A-posteriori error estimators for second order elliptic systems, Part 2. An optimal order process for calculating self-equilibrating uxes. Computational Mathematics and Applications 1993; 26:75–87. 9. Stein E, Ahamad R. An equilibrium method for stress calculation using nite element displacement models. Computer Methods in Applied Mechanics and Engineering 1977; 10:175–198. 10. Bank RE, Weiser A. Some a-posteriori error estimators for elliptic partial dierential equations. Mathematics of Computation 1985; 44:283–301. 11. Babuska I, Strouboulis T. The Finite Element Method and its Reliability. Oxford University Press: Oxford, 2000, to appear. Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 427–475 (2000) GUARANTEED COMPUTABLE BOUNDS FOR THE EXACT ERROR—PART II 475 12. Szabo BA. Babuska I. Finite Element Analysis. Wiley: New York, 1991. 13. Paraschivoiu M, Peraire J, Patera AT. A-posteriori nite element bounds for linear functional outputs of elliptic partial dierential equations. Computer Methods in Applied Mechanics and Engineering 1997; 150:289–312. 14. Paraschivoiu M, Patera AT. A hierarchical duality approach to bounds for the bounds of partial dierential equations. Computer Methods in Applied Mechanics and Engineering 1998; 158:389–407. 15. Harwell Subroutine Library, AEA Technology, Oxfordshire, England. 16. Babuska I, Strouboulis T, Copps K. The design and analysis of the generalized nite element method. Computer Methods in Applied Mechanics and Engineering, 1999, in press. Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 427–475 (2000)

1/--страниц