- Email: [email protected]
Quality assessment of the a-posteriori error estimation in finite elements
Finite Elements in Analysis and Design 11 (1992) 285-306
285
Elsevier FINEL 250
Quality assessment of the a-posteriori error estimation in finite elements Ivo Babu~ka, Lothar Plank a and Rodolfo Rodriquez a,b a Institute for Physical Science and Technology and Department of Mathematics, University of Maryland, College Park, MD 20740, USA b Departamento de Matemdtica, Facuhad de Ciencias Exactas, Universidad Nacional de La Plata, C.C. 172, 1900 La Plata, Argentina Received October 1991 Revised March 1992 Abstract. The paper addresses results of numerical experimentation with an a-posteriori error estimator. The
paper concentrates on the problem of the selection of benchmarks which is directed to a numerical verification of basic theoretical features of the estimator.
Introduction
During recent years much interest has been focused on the design of a-posteriori error estimates in the finite element method, their experimental verification and their use for adaptive procedures. We refer for example to [1-24] and further citations there. The a-posteriori error estimators are often derived by purely heuristic reasoning and without mathematical analysis. If a mathematical analysis is made, it adresses most often only asymptotic properties and special meshes. Usually the estimators are numerically analyzed (verified) on a set of examples (benchmarks) which are selected at best by heuristic reasoning. Often these benchmark computations are used for the design of a correction factor which is used to improve the quality of the estimators for this set of benchmarks. This pragmatic approach to test and compare the estimators by benchmark computations has serious shortcomings. The quality of an estimator is usually very sensitive to the structure of the solution and the meshes used. The performance of an estimator can be very different. It may depend on whether the meshes are nearly translation invariant or general, or whether the meshes are general or such that the error indicators are nearly equal on all elements (i.e. the meshes are equilibrated), or whether the solution is smooth, has point singularities or is generally unsmooth, etc. Some of these aspects are adressed in [5]. Hence the benchmark computations could motivate misleading conclusions. Each benchmark computation should be a representative of a more or less precisely defined class of computational problems and the conclusions used for this class of problems. In this paper we will address various aspects of the relation between theoretical understanding of an estimator and the benchmark selection. We will adress a certain estimator which is used in practice, and discuss its properties. We do not intend to compare it with any other error estimator. Correspondence to: Ivo Babu~ka, Institute for Physical Science and Technology and Department of Mathematics, University of Maryland, College Park, MD 20740, USA. 0168-874X/92/$05.00 © 1992 - Elsevier Science Publishers B.V. All rights reserved
286
1. Babu~ka et al. / Quality assessment of the a-posteriori error estimation
Model problem ket /2 ~ E2 be a b o u n d e d polygonal domain and let 1" be its boundary. A s s u m e further that .0 = U i~_ t~i, where gl i is a polygonal domain (with b o u n d a r y a/2 i) and that /2i c~ S'lj = 0 for i 4:j. We will be interested in solving the following problem:
-div(aVu)=f 0u a--=g an
in/2,
on I~:,
u=0
on I)~,
l't~Ul~.=l',
(1)
where n is the o u t w a r d normal to /2. We assume that f ~ H ~(/2), g ~ H~(1)~), a >/a 0 >~ 0 is constant on /2i, i = 1. . . . . m. We will assume that either m e a s ( F o) ~ 0 or, if F D = 9i. f a r d x d y + frg ds = 0 and u ( A ) = 0, A ~ . (Because of the assumption about f and g. u ~ H~(/2), a > 1, and so u ( A ) = 0 makes sense.) We will consider the weak solution of (1). To this end we d e n o t e by Hk(/2) the usual Sobolev space and let H,l{u
= 0 on
F u r t h e r we let
B ( u , /')
=-~f
SU 0t! 0U ~)l') ff - - - - + - - - dx dy ax ax
ay a y ,
(2)
be the bilinear form defined on H(}(/2) × H~I(~) and
X(t!) = f~ft' d x d y +
frgL,
ds
(3)
be a linear functional on H(}(/2). T h e n the weak solution u of (1) in H(I exists, is unique and
B(u,c)=F(t')
Vc~H,I.
(4)
If 1"0 = ~, we consider H I ( / 2 ) instead of H01(/2) and assume u ( A ) = 0, which is well-defined, as we have in fact higher regularity for u. Let us consider a family ~ of triangular meshes M on /2 satisfying a uniform minimal angle condition. We d e n o t e by 3 - the closed triangles of the mesh M, U 3 - = ~ . W e will assume that if J - ~ M then i n t e r i o r ( Y ) c / 2 t for some l, i.e. 3 - lies in a domain where a is a constant. To every mesh M we associate a p a r a m e t e r t(M)>1 O, for example t ( M ) = N where N is the n u m b e r of the nodes that do not belong to 1"D or t ( M ) = h where h = m a x ( d ( J ) ) and d ( 3 - ) denotes the diameter of ~ . W e will always assume that h = m a x ( d ( 3 j ) ) ~ 0 as t ( M ) ~ O. F u r t h e r let S ( M ) be the set of all continuous functions o n / 2 whose restriction on a g ~ M is a linear function. Obviously S( M ) c Hl(g2). Let So(M) = H~(12) C~S ( M ) and S A ( M ) = {u S ( M ) , u ( A ) = 0}. T h e n the finite element solution u s ~ So(M) satisfies
B(us, c)=B(u,u)
Vc~S0(M),
(5)
w h e r e u is the exact solution of (1). If 1"D = ~, we replace So(M) by SA(M). We will consider a family of finite element solutions u s associated to ~t" and assume that II u - US(M)tl E = B(u US(M), U -- US(M))1/2 ~ 0 as t ( M ) ~ O. By (5) we associate to any solution u its finite element solution u s = us(M, u) and an error estimator ~ ( u , M ) that will be defined in the next section. In fact ~ d e p e n d s on u s and M but because of (5) we write ~ ( u , M). This estimator approximates the energy n o r m of the error, i.e. ~ ( u , M ) = II e [[ E = (B(e, e)) 1/2, where e = u - u s . -
287
I. Babugka et al. / Quality assessment o f the a-posteriori error estimation
Error estimator
Let M ~atv, ~ / , ~ t ~ M be two triangles with common side 3"l,i as shown in Fig. 1. By nt, i we denote the outer normal from ~ t to 5~; we will also write ~ t = ° ~ i n and ~ / = 3-out, as ~in is the triangle for which the normal nt,/ is the outward one. Let us now define for every side 3"l,i
[[a~. ](x) = (aVus) l~root .nl, i - (aVus) l:~, .nl,i,
(6)
the jump of a(aUs/Ont, i) across the side 3"t,i. This value is independent of the choice of nt, ~, i.e. [[a(OUs/anl,i)](x) = ~a(bUs/OUi,t)]l(x) and so we will write ~a(OUs/On)](x). For each triangle 3- ~ M let
R : = - [div(aVu) I ~-+ div(aVus) l : ] + f I : .
(7)
We may write for any v ~ H d the residual equation
B(e,v)=
~ f e:vdxdy+ 5r~M
g-a
vds+ ~ J ~ a T l l v d s ,
"~
yEG
3'
(a)
Otl
where G stands for the set of all interior edges of M and n for the outward normal to 12. (If
Fo = 0 eqn. (8) is the same, but v ~ Hi(12).) For any 3 - ~ M we define 1
i----(fR:
11:=
dx d y
(9)
where t3-1 is the area of 3". (1I: is the L 2 ( 3 -) projection of the interior residual onto the set of constants). Further for each 3' ~ FN, where 3' is a side of a triangle ,Y, we define
1
aUs I
l-l, = T i f,( g - -Tnn l dS,
(10)
where 13" I is the length of 3'. ( H , is the L2(y) projection of the boundary residual onto the set of constants.) For each edge y of the mesh let
{
[[2(~Us/On)]]
j~=
H~,
if y ~ G,
if 3'c_rN,
(11)
if y c F o ,
and for each triangle 3 - ~ M we define the error indicator r / : : n : :=
13-1
"'S~dxdy+7
~
13/I a--i--~ds
,
(12)
Fig. 1. Scheme of neighbouring elements.
L Babu~ka et aL / Quality assessment of the a-posteriori error estimation
288
which approximates the energy error on ~5r. Further we define the error estimator ~ ( u , M ) :=
r/.~-
•,
(13)
" 17 ~ M
where K is a suitable constant which will be defined later. It is convenient to write g*(u,
M ) :=
Y[ r/.~
(14)
.'f ~ M
for the choice K = 1. Let us note that the error indicator is meaningful for more general f and g than we have assumed. Nevertheless our generality is sufficient for practical purposes. Obviously the set of problems we are discussing is precisely described.
Quality of the error estimator
In the previous section, we defined the error estimator g~(u, M ) which approximates IIe IIE. In the following we will discuss the relation between g'(u, M ) and I1e II E, especially the effectivity index of the estimator ~: = ~ ( u , M ) / l l e l l t . : . Before addressing this problem we introduce some notation which we will need in the following. Let S P ( M ) be the space of all continuous function on J2 whose restriction on 3~M is a polynomial of degree p. Let S P ( M ) = S p h i l l ( s O ) and u~ ~ S g ( M ) be the finite element solution of problem (1), i.e. B(u}',v)=B(u,v)
Vv~Sg(M).
(If F D = fl, we use S P ( M ) instead of S g ( M ) . ) Obviously Uls' is the p-version solution of problem (1). For more about the h-p-version we refer to [6]. We have for any k > 1:
II U~s'
-
u
II ~ ~ C p - ¢ ~ - j) II u It Hkh min{p'k- I}.
(15)
with C dependent on k, but not on p and h. Remark. Although in (15) the maximal element side is used, an analogous element by element estimate holds. Now we are able to state the basic theorems about N * ( u , M). (See [1] for the detailed proof and the next section for an outline of the main ideas.) Theorem 1. For any integer p there is a constant K f depending only on the minimal angle ce.i o f the triangles used and on p, such that
Il e HF. ~ K f ~ * ( u, M ) + Hu - u~ ll E + O( h 3/2)
(16)
and the constant K f grows logarithmically with p. Theorem 2. that
There is a constant K l > 0, depending only on %7 and on the jumps of a, such
~ * ( u , M ) < ~ K l llell E + O( h3/2).
(17)
L Babugka et al. / Quality assessment o f the a-posteriori error estimation
289
Table 1 Constants C~ and C~- for various angles a :
a
C:
7.5 15.0 22.5 30.0 37.5 45.0 52.5 60.0
0.051 0.072 0.087 0.099 0.108 0.115 0.119 0.121
C.~
t
p=2
p=4
p=6
p=8
2.390 1.682 1.363 1.169 1.035 0.939 0.875 0.850
3.306 2.341 1.918 1.670 1.508 1.400 1.334 1.309
3.660 2.609 2.156 1.895 1.727 1.615 1.547 1.522
3.988 2.839 2.343 2.058 1.876 1.757 1.684 1.657
The constants K~ and K~' are computable. We have K~ = max C~-,
(18)
3~M
K 1 = ( max C~r)Q,,
(19)
Qa =
(20)
where
max
:ng-'*O
(al:/al~,).
For 2 ~
(21)
3.45 sin-l/2( ot:/2) <~C~.~ 5.85 sin-1/2( a : / 2 ) ,
(22)
and
where t~9. is the minimal angle of 3 . In Table 1 we give the values of C~- and C~-. Let us note that for h small only very few elements have a side on the interface 092; we can in fact use Qa = 1. Remark. We have dealt only with the Laplace equation. In this case we use the minimal angle % r of the triangle. For a general differential equation we have to replace it by ot~which is the minimal angle of the triangle after a local transformation of the general equation to Laplace's equation (this is always possible). The t e r m O ( h 3/2) in (16) and (17) depends only on II f tl nl(n) and II g II nkr), but not on u. For any u the term II u - uff II E can be made arbitrarily small by (15) using sufficiently large p. We note that the smoothness of u is governed by the smoothness of f and g, transition of boundary conditions and smoothness of F. It is unessential that the estimates (16) and (17) have an asymptotic character. For example if f = const on every 3 - and g = const on every triangle side of FN, the term O(h 3/2) disappears. Hence as t ( M ) ~ O, the term O(h 3/2) is negligible in general. Theorems 1 and 2 indicate that the error estimator ~ * ( u , M ) can be either smaller or bigger than the true error. Further the error estimator loses accuracy when a 3 ~ 0 and for unsmooth u the reliability of the estimator decreases because a higher p is needed in the estimator. The coefficient x, which scales the error estimator (~'(u, M ) = K ~ ' * ( u , M)), can be calculated in different ways. We can use the geometrical average of the bounds in (18), (19).
290
1. Babugka et al. / Quality assessment of the a-posteriori error estimation
Another possibility is to choose K such that the error estimator is asymptotically exact for uniform meshes of equilateral triangles and f = 0. In this case we get x = ( 8 ~ ) - 1 / 2 = 0.26864 which will be used for all examples. Theorems 1 and 2 indicate the expected theoretical properties of the estimator well. The structure of Theorems 1 and 2 clearly suggests the set of benchmark experiments that should be performed.
Outline of the proof of Theorems 1 and 2 For simplicity assume that f and g are constant on each g resp l ~ F N. (The general case is discussed in [1].) With ep := u~ - u S E S(p, we have I l e l l E < ][ep]lE+ ]]U--Uf][E
(23)
II epll2: = B ( ep, ep) = B(uJs' - u, % ) + B( u - u s, % ) = B ( e, % ) .
(24)
and
For any continuous function v let v I denote the linear interpolant of v on 3 . Using (8), we have:
y"fJ'(e"-4) dsl
l,ep[,~=B(e,ep)=.,ic. ~ [fY :x(ep-elp) ~,2M
rlJ
~-T~(ep-elp)
1 +~,c~ 3
<
Y'. J~M
dxdy
.......
Va l j ~-jl
/ca..
r" ( e t , - e I ) ds
2) )/ 2
~JC~-II ep II t.:,J,
where
1 -~lfT(v-
CP~-=
sup
Vl) d x d y I 2 + ' ~
'
"~e,,\P,,
1
I~-
f/ (U - UI) dsl 2
) 1/2
Ica~
f lV v l
2 dx dy
Here, Pp and P0 denote the space of all polynomials of degree p or 0 on 5r, respectively. Hence ]l epll E < ~ K ~ g ' * ( u , M ) , with K f = sup C~-.
J~M
Using (23) we get (16) without the term O(h3/2). This is related to the simplification assumption about f and g made in this section. Let us now address inequality (17). For any triangle Y as shown in Fig. 2, let ¢0 be the cubic bubble function vanishing on a3- with ¢0(Q0) = 1, where Q0 is the barycenter of 3-. Let
L Babu~ka et al. / Quality assessment o f the a-posteriori error estimation
291
Fig. 2. Notation for triangular elements.
further ~ ° ( 3 - ) be the space of quadratic functions vanishing at the vertices of .Y- and let {~°t}3=1 be the canonical basis of this space (i.e. q~i(Qj)= ~ij, i = 1, 2, 3, j = 1. . . . . 6) and let V: be the space spanned by {~ot}~=0. Set w ~ = E3=oWiq~i~ V: with coefficients w0, Wl, w2, (-03 chosen such that /i2
dx d y = [.Y-I J[ y "': dx dy, a ftJtoo: ds = •
where
j2 li 11i [ f 7 - ' T d s ,
(25)
i = 1, 2, 3,
(26)
li[aJt i
[a]t, is the average of a over the triangles sharing the edge li, i.e. / =
[a]ti
| /
x
I 2(al~'i"+al:°u') ~a[:
if li = a~in ~ &Y-°ut'
(27)
if lic&Y'NF.
The coefficients w 0. . . . . o)3 satisfying (25) and (26) are OJi
t°O =
2[a]l i '
a-----~--:
V
2 i= 1 [ a ] l i
li
).
Since for any interior edge l the value wv) are the same on both triangles 3- sharing the edge, the function w composed from w : on every triangle is continuous and hence oo~ H~I'2). Using (23), (24) we get: ( g ~ * ( u , M ) ) 2=
]~ r/2 3-~M
17z " ' ~ d x d y + ~ 1 E I11 ,: a leo:
j2
]
=E
13-If
=E
[a]t /j~o f J l f ~ o : dx dy+ ~1 ,~.~a~aljl: , : ds ]
,~~M
ds
~-~M
-<
where [a]t Co =
max
1~&fr,3-~M
a I ~-
Since oJ ~ H i we get from (8)
(g'*(u, M))2<~QB(e, w) < C~ Ilell Elleo lIE.
(28)
I. Babu~ka et al. / Quality assessment o f the a-posteriori error estimation
292
3 Now let for r = {r i}i=o
C9=
i l Vt.,r l I12
sup
"
[Irll.~O
dx dy
3
r 2
,
1
o + 5 E
ri-
i=1
where v~ E E~ is the function satisfying f~v r d x dy = I J l r o , T h e r e f o r e w z is the function v r corresponding to
I lilJz, [a]li
IJ-l/Tj
r° -
- - - ,
r i --
alj
,
fliUr ds = I li[r i, i = 1, 2, 3.
i=1,2,3
and so
I1~o~ IIu2 . , z - a [ . ~ f lIVVrtl 2 d x d y
2 ,
~ c;~,l~( =c.~( IJ
a]: 2
211'2j
+ Z121 3
1{ 2
i=l
[a]6
( a~ll~, )
l
01:7
a
f.,2p,2 ~< t _ . j ~ a
2
7"/aw,
where C~ =
alj
max ,c~J,~M [a],"
Hence II o)I[ e ~< { max " ,Y-~M
C~-}C'~'*(u, M)
and so, using (28), we obtain ~*(u,
M ) ~< Qa( max C ) ) II etl ~, J~M
where
Qa = C.C~
a[:7 = max - .gem a [ ~,
and we get (17).
Experimental results In this section we will experimentally verify the reliability of various conclusions that can be m a d e from T h e o r e m s 1 and 2 for the given estimator.
Influence of mesh topology T h e o r e m 1 was proven for general meshes without any restrictions on the mesh topology. H e n c e , it can be expected that particular m e s h e s with different topologies could influence the
I. Babu~ka et aL / Quality assessment of the a-posteriori error estimation
// //
W
293
//
mesh ratio 1/1
//// //// mesh ratio 1/2
Fig. 3. Regular three-direction mesh.
effectivity index in the same order of magnitude as indicated by the theoretical upper and lower bounds. In addition the influence of the angle has to be visible (and could also depend on the topology). To check the validity of this conclusion for a smooth solution we consider the following problem: For 12 = (0, 1) x (0, 1), let u be the solution of -Au =f
in 12,
u =0
on ~e.
(29)
We choose f such that u = sin wx sin -try.
Three-direction meshes Let us consider the three-directional meshes which are obtained by uniform refinement of a basic mesh of the type 1/n, n = 1, 2, 4, 8 as shown in Fig. 3. In Table 2 we report the effectivity index for these meshes. We also report the u p p e r and lower bounds for the effectivity index from Theorems 1, 2, when p = 2 was used.
Criss-cross meshes Let us consider now the "criss-cross" meshes as shown in Fig. 4. The effectivity indices for these meshes are given in Table 3.
Diamond-shaped meshes Finally, we consider the "diamond-shaped" meshes as shown in Fig. 5. The effectivity indices are given in Table 4.
Table 2 Effectivity index for problem (29) (u = sinrrx sin~ry) and 3-direction meshes Number of elements
Mesh
8 16 32 64 128 256 512 1024 2048 4096
0.852
Upper bound Lower bound
2.34 0.29
1/1
1/2
1/4
1/8
0.993 1.070
1.289 1.183
1.147
1.804 1.503
1.252 1.180
2.068 1.571
1.277 1.192
2.145 1.594
1.287 2.87 0.21
2.169 3.87 0.15
5.38 0.12
1. Babu~ka et al. / Quality assessment of the a-posteriori error estimation
294
/ \
,
\/\
\/\/ /\/\
mesh ratio 1/1
, mesh ratio 1/2
Fig. 4. Regular criss-cross mesh.
Table 3 Effectivity index for problem (29) (u = sin'rrx sin'rry) and "criss-cross'" meshes Number of elements
Mesh
8 16 32 64 128 256 512 1024 2048 4096
1.35(I
Upper bound Lower bound
1/1
1/2
1/4
I/8
{L939 1.1136
1.1}49 1.053
1.077
1.349 1.179
1.{185 1.088
1.518 1.214
1.093 1.I19(I
1.562 1.223
1.I195 2.34 I).29
2.87 {}.21
1.574 3.87 0.15
5.38 I). 12
Comparing the tables we can make the following conclusions, which are in agreement with Theorems 1, 2: (1) For h ~ 0 the effectivity index stays in the theoretical bounds. In addition we see that in our concrete case the effectivity index converges to a limiting value. This is not surprising, because the mesh is translation invariant and the solution can be periodically extended. In this case various results about superconvergence can be used [21]. (2) The effectivity index depends on the topology. The error estimator is not asymptotically exact. (3) The dependence of the effectivity index on the angle is clearly visible.
mesh ratio 1/1
mesh ratio 1/2 Fig. 5. Regular "diamond-shaped" mesh.
L Babu~ka et al. / Quality assessment of the a-posteriori error estimation
295
Table 4 Effectivity index for problem (29) (u = sinrrx sinlry) and "diamond-shaped" meshes Number of elements
Mesh
6 12 20 24 40 48 72 80 144 160 272 288 544 576 1056 1088 2112 2176 4160 4224 8448
1.465
Upper bound Lower bound
2.34 0.29
1/1
1/2
1/4
1/8
0.750 1.017 0.523 1.026 0.385 1.107 1.111 1.188 0.959 1.131 1.558 1.234 1.708 1.138 1.740 1.246 2.274 1.139 1.793 2.515 2.87 0.21
3.87 0.15
5.38 0.12
(4) The theoretical bounds for the effectivity index are relatively pessimistic. This is because these bounds are independent of the topology and the solution. Better bounds in our case could be obtained by exploiting the superconvergence effect [5].
Influence of the solution structure Theorems 1 and 2 have been shown for a general solution. We expect that the effectivity index may depend on the solution structure. To check this conclusion consider the solution of the following problem on g2 = (0, 1) × (0, 1): -Au=f au --=0
by
ing2, for y = 0 a n d
u=0 y=l
for x = 0 and x = l, (30)
with f such that u = sin ~ x . We will use meshes as in Fig. 3 with the orientation shown there (mesh 1 / 2 and 1 / 4 ) and with the opposite orientation ( 2 / 1 and 4 / 1 ) . From Table 5 we see in fact that the effectivity index depends on the relation between the mesh and the structure of the solution (especially the relation between the second derivatives and the mesh orientation). Hence in general we have to deal with uncertainties in the effectivity index. (No correction factor could improve the situation.) We note that our factor K = 0.2684... could be changed, but in general (for all funotions) the range in the effectivity index could not be improved. Again, we see that Theorems 1 and 2 reliably predict the outcome of the numerical experiments.
296
I. Babugka et al. / Quality assessment 0["the a-posteriori error estimation
Table 5 Effectivity index for problem (30) (u = sin'n'x) and 3-direction meshes Number of elements
Mesh
8 16 32 64 128 256 512 1024 2048 4096
1.013
Upper bound Lower bound
2.34 0.29
1/1
l/2
1/4
4/1
2/1
1.564
[).717
1.1117
2.258
0.5117
1.598
/).783
1.131
2.273
11.553
1.608
il.8011
1.137
2.278
().566
1.611
11.804
1.139
2.279
11.569
1,612
0.806
2.87 0.21
3.87 0.t5
2.87 0.21
3.~7 1/.15
Influence o f the smoothness o f the solution T h e o r e m s 1 a n d 2 c o n t a i n h i g h e r - o r d e r t e r m s w h i c h we n e g l e c t e d in t h e r e p o r t e d t h e o r e t i c a l b o u n d s in T a b l e s 2 - 5 w h e n c o m p u t i n g t h e b o u n d s . T h e s e t e r m s i n f l u e n c e t h e c o n v e r g e n c e o f t h e e f f e c t i v i t y i n d e x as h --* 0. T o c h e c k t h e e f f e c t o f t h e s e t e r m s w e will c o n s i d e r o n c e m o r e p r o b l e m (29), u s i n g n o w t h e s o l u t i o n
u ( x , y ) = sin a'rrx sin a w y .
(31)
( W e use r e g u l a r t h r e e - d i r e c t i o n m e s h e s , as in Fig. 3 w i t h m e s h r a t i o 1 / 1 . ) We see from Table 6 that the insufficient smoothness influences the effectivity index and that the effectivity index can behave erratically. Nevertheless the bounds--where higher order terms are neglected--are still valid. W e s e e in g e n e r a l t h a t t h e e f f e c t i v i t y i n d e x decreases with decreasing smoothness of the solution. The lower bound goes down with i n c r e a s i n g p , b u t t h e u p p e r b o u n d is n o t i n f l u e n c e d . H e n c e w e c a n e x p e c t t h a t t h e e f f e c t i v i t y i n d e x will go d o w n ( a s s u m i n g h e u r i s t i c a l l y t h a t t h e e f f e c t i v i t y i n d e x m o v e s as t h e m i d d l e o f t h e interval). T h e s a m e e f f e c t will b e s e e n later. O n c e m o r e w e s e e g o o d a g r e e m e n t w i t h t h e c o n c l u s i o n s s t e m m i n g f r o m T h e o r e m 1 a n d t h a t w e c a n n o t h o p e for t h e e x i s t e n c e o f a universal correction factor.
Table 6 Effectivity index for problem (31) (u = sinawx sina'rry) and 3-direction meshes Number of elements
a= 1
a=2
a=4
8 32 128 512 2048 8192
0.852 1.071/ 1.147 1.180 1.192
1/.833 0.876 1.080 1.161 1.187
0.438 0.932 0.887 1.090 1.168
Upper bound Lower bound
2.34 0.29
2.34 0.29
2.34 0.29
a = 1l
0.379 0.5311 0.7811 1.022 1.148 2.34 0.29
L Babu~ka et al. / Quality assessment o f the a-posteriori error estimation
297
Table 7 Effectively index for problem (32) (u(x, y) = sinarx sinh'rry) and 3-direction meshes Number of elements
Mesh 1/1
1/2
1/4
1/8
0.841
8
16 32 64 128 256 512 1024 2048 4096
0.958 0.978
1.269 1.102
1.040
1.795 1.457
1.166 1.065
2.042 1.534
1.190 1.071
2.136 1.560
1.199
Upper bound Lower bound
2.34 0.29
2.87 0.21
2.166 3.87 0.15
5.38 0.12
Harmonic solutions The error estimator contains terms related to the residual and the jumps, respectively. To isolate the influence of the jump term, we consider the following problem (again i n / 2 = (0, 1) x (0, 1)): -Au=0
inO,
u=0
forx=O,x=landy=O,
0u
-- =g On
for y = 1,
(32)
and let u(x, y) = sin'rrx sinh~-y. The effectivity index (for a regular three-direction mesh as shown in Fig. 3) is given in Table 7. Comparing Tables 2 and 7 we see a similar behaviour influenced only mildly by the structure of the solution. The limiting value for the effectivity index seems to be different for the mesh 1/1. For the other meshes this effect is less visible. For these particular meshes, this might be explained by a superconvergence phenomenon.
Adaptively constructed meshes So far, we have restricted our experiments to uniform meshes. We will now consider nonuniform meshes either based on an adaptive procedure or randomly refined. We use problem (32) again. We have constructed three sequences ( A - C ) of meshes. Each sequence is obtained by refining a part of the elements of the previous mesh according to a refinement indicator and adjusting neighbouring elements. In sequences A and B the refinement indicator was an error estimator (with different parameters), in sequence C the indicator was a random number. Typical meshes of each sequence are shown in Fig. 6. The effectivity indices are given in Table 8. Comparing Tables 8 and 7 we see that different refinement strategies lead to different effectivity indices, although the difference is not big. Once more we see a range of the effectivity indices. We note that for these general meshes superconvergence results cannot be applied.
I. Babu,~ka et al. / Quality assessment of the a-posteriori error estimation
298
/l//2/
,G,'/\
Y// y//yy
///
/
/ / /
..-v z
/
(b) ---
//
•
/
///
/
"
k
/,"
\ /
,
~ "
\\\
/
Fig. 6. Nonuniform meshes for problem (32).
Effectiuity indices f o r problems with unsmooth solutions So far, we have considered problems with a s m o o t h solution. L e t u s now e x a m i n e the b e h a v i o u r of the error estimator for solutions with singularities. T h e o r e m s 1 and 2 indicate a
L Babu~ka et a L / Quality assessment of the a-posteriori error estimation
299
Table 8 Effectivity indices for problem (32) and nonuniform meshes Sequence A
Sequence B
Sequence C
Elements
Effectivity index
Elements
Effectivity index
Elements
Effectivity index
8 19 47 94 181 323 593 1020 1773
0.841 0.944 0.970 1.008 1.018 1.041 1.043 1.060 1.051
8 19 38 67 107 201 355 607 991
0.841 0.944 0.941 0.980 0.949 1.013 1.039 1.032 1.056
8 23 45 96 200 435 964 2129 4588
0.841 0.906 0.964 1.008 0.979 0.990 0.972 1.008 0.992
2 a
a(x,y) =
3
1
4
a(x,y) = 1 a ( x , y ) = Fig. 7. Coefficient a for problem (33).
decrease of the effectivity index, w h e n the singularity of the exact solution grows, and also that this effect diminishes when adaptive meshes are used. T o this purpose we consider the domain O = ( - 1, 1) × ( - 1, 1), partitioned into 4 square subdomains as shown in Fig. 7. Let u be the solution of
-V(a(x, y)Vu) =0
in O,
0u
a(x, Y)-~n=g(x,
where a(x, y) is constant in each quadrant chosen so that U(X, y ) = r a ( c i
COS(Ot~)÷s i
y)
on F,
(33)
(a(x, y ) = 1 or a(x, y ) = 5, see Fig. 7). 5 is
sin(a~o)),
(34)
where (r, ~0) are polar coordinates with center in the origin. For a = 0.5 and a = 0.1 we get 5 = 3 + v~- = 5 . 8 2 8 . . . and 5 --- 166.447 . . . . respectively.
Table 9 Effectivity indices for problem (33) (u = r~fQo)) and 3-direction meshes Number of elements
a = 0.5
a = 0.1
32 128 512 2048 8192
0.658 0.694 0.711 0.719 0.722
0.520 0.512 0.524 0.536
300
I. Babu~ka et al. / Qualiw assessment o f the a-posteriori error estimation
//LK_
~V
\\// ,/ \ \ / / // /. \\ \\\\/////~;'i_ \\///// /J" / V /
/\\\\
/
\\X
~ / / F S/ / ~/ ~
J~ /// \\ \\ \\ \\ \\
,. ;:7 ~ ~ , \ )
F"
(al
I
]
]
"
"\
[
Fig. 8. Nonuniform meshes for problem (33).
Table 9 shows the effectivity index for uniform meshes. (In each quadrant the meshes are uniform three-direction meshes of mesh ratio 1/1 as in Fig. 3. The orientation of the meshes is changed in subdomains 2 and 4 to preserve the symmetry of the solution.) Although the estimator in Theorem 1 includes the ratio of a, if influences only elements on the interfaces and hence this factor can be neglected. The decrease of the effectivity index is clearly related to the necessity to use a higher degree p in (16) for the estimate. Hence the unsmoothness of the solution decreases the effectivity index. This has also been observed in Table 6. We have seen in Table 8 that nonuniformity of the mesh does not influence the effectivity index too much. On the other hand an adaptive procedure refines the mesh in the places where the solution in unsmooth and hence a lower p can be used. We have created two sequences of adaptively refined meshes (A and B). The meshes B1 and B2 are less refined than the meshes A1 and A2. The refinement in meshes B is more concentrated around the origin than in meshes A. Examples of each sequences are displayed in Fig. 8. If the mesh is properly refined, a lower p can be used and the effectivity index improves. Table 10 confirms this conclusion (for c~ = 0.5).
L Babu~ka et al. / Quality assessment of the a-posteriori error estimation
301
Table 10 Effectivity indices for problem (33) (u = r~f(~)) and adapted meshes Sequence A
Sequence B
Elements
Effectivity index
Elements
Effectivity index
32 87 194 370 725 1297 2346
0.658 0.712 0.738 0.761 0.787 0.807 0.828
32 64 96 168 244 363 485 604 768 965
0.658 0.719 0.787 0.830 0.880 0.918 0.960 0.959 0.975 0.995
The effect of the unsmoothness of the solution can also be observed if the singular part o f the solution is removed. So we consider problem (33) in the domain l ) = O - [ - 0.5, 0.5] x [ - 0 . 5 , 0.5] and prescribe a(au/an)=g o n aO. T a b l e 11 shows the effectivity indices for this problem using uniform meshes as for Table 9. Again, a = 0.5. W e note that the effectivity index is now very similar to problems (29), (30)." The oscillation
effect
It is often assumed, that the error of the interpolant is close to the error fo the finite element method. To demonstrate that this assumption is in general incorrect, we show that the error of the finite element method could oscillate. To this end we show in Fig. 9 the typical error behaviour for the problem discussed in the previous section with uniform meshes and a = 0.5. Displayed are two patches of regular three-direction meshes at the same location - - f a r from the singularity--in O. Figure 9 also shows the error indicators (in parentheses). In Fig. 10 we illustrate graphically the error behaviour for the entire mesh where the mesh has 363 elements. (An element is black if the energy norm error in this element exceeds a certain threshold. Here, the threshold is 4.08 x 10-4.) We see that the error indicator shows approximately the average error for adjacent elements and hence does not exhibit the oscillation of the error of the finite element solution. Note that the interpolant does not show the oscillation. This oscillation effect is caused by the pollution of the error stemming from the singularity of strength r ~/2. If we enforce a smooth solution by removing the center part of the domain, as in problem (33), the results given in Fig. 11 show that the oscillation disappears. (We use elements from the same location as in Fig. 9; the numbers are therefore directly comparable.)
Table 11 Effectivity indices, for problem (33) (u = r~f(~) on O
a7) and 3-direction meshes
Number of Elements
a = 0.5
24 96 384 1536 6144
0.964 1.027 1.059 1.070 1.073
L Babu~ka et al. / Quality assessment o f the a-posteriori error estimation
302
(1.2~/
2.2(-3) /
2.5(-3) / (1.4~)/
/ 1.o(-3)
/
1.1(-3)
/ (1.2(-3))/ / (1.4(-3)) 2.4(-3) / 2.2(-3) / (1.0~)/
/
(1.2U
/9.7(-4) (1.0(-3)) / /
5.6(-4) /
5.9(-4) /
5.6(-4) /
5.9(-4)
(3.2~ (3.°U / 2.7(-4) /// 2.6(-4) (3.o(-4)) / (3.2(-4))
(2.7U
/ 2.6(-4)
/ i.i(-3) (1.2(-3))
small element side: 1/8
///
(2.7(-4))
(3.o~ / 2.6(-4)
/ (3.0(-4))
small element side: 1/16
values in (..): error estimator Fig. 9. Oscillation of the finite e l e m e n t error.
The oscillation effect is weaker, if a 4:0.5 or the mesh is constructed adaptively. We show in Fig. 12 the local results for a = 0.8 and in Fig. 13 the results for a = 0.5 and adaptively refined meshes. In Fig. 14 we show the error behaviour for an adaptively refined mesh and the threshold 2.0 × 10 -3. In the next section we will examine this effect. Analysis of the oscillation phenomenon
The oscillation phenomenon we have observed in the previous section occurs also in the case of the following problem on O = [ - 1, 1] × [0, 1] - A u = 0 inO, u = 0 onFD={(x,y)lO<~x<~l,y=O }, Ou --=g On
on F N = F - F ~,
Fig. 10. Finite e l e m e n t e r r o r ( a = 0.5, u n i f o r m meshes).
I. Babu~ka et at / Quality assessment o f the a-posteriori error estimation
(3.2~/ (3.o~// 2.8(-4)
1.3(-3) / 1.1(-3) / (1.4~)/ (1.2~/
/ //
/ 1.3(-3)
1.1(-3) (1.2(-3)) /
(1.4(-3))
1.0(-3)/ 1.2(-3)/ (1.1~)/ (1.3~)/ / //
1.0(-3) / 1.2(-3) (1.1(-3))/ (1.3(-3))
303
2.4(-4)
/
//
2.6(-4) / (3.0(-4)) / /
2.5(-4) /
(2.T~/
/ /"
small element side: 1/8
2.9(-4) (3.2(-4))
2.8(-4) / (3.0~/
2.5(-4) / (2.7(-4)) /
2.8(-4) (3.2(-4))
small element side: 1/16
Fig. 11. Finite element error and estimators (in parentheses) for a problem without a singularity (a = 0.5, adaptively refined mesh).
with g selected so that u =r 1/2 sin(½q~)-- Im z 1/2, with z = x + iy. We consider now a uniform 3-direction mesh of the type shown in Fig. 3. Then the finite element solution has essentially the form (35) where IIR(.,.)II H' < o(h) and G(x, y) is a piecewise linear function, which coincides in the nodal points of the mesh with the function G*(x, y) =Ar -1/2 sin ½Oh =Ah Im Z-1/2 (36) Us(X, y) =u'(x, y) +G(x, y) +R(x, y),
The function G*(x, y) is the adjoint singularity function. It describes the pollution error, which is of the same order (in II" IIE) as the error of the interpolant. The error--measured in
6.0(-4) /
7.0(-4) /
(5-1U (5.9~ / / 5.1(-4) (6.2(-4)),l / / 6.1(-4) (?.0(-4)) 5.9(-4) / (5.6~ 6.7(-4) (6.4~// / /
4.6(-4) (5.4(-4))
/
/
5.3(-4)
(6.2(-4))
small element side: 1/8
1.4(-4) / (l'3U
1.6(-4) / (1.3~/
/
/
1.3(-4) / 1.4(-4) (1.5(-4)) i / (1.6(-4)) 1.4(-4) / 1.5(-4) /
(1.4~)/ (1.5~/ / (1.4(-4)) 1.2(-4) / / (1.5(-4)) 1.3(-4) / small element side: 1/16
Fig. 12. Finite element error and estimators (in parentheses) for problem with singularity ( a = 0.8, uniform mesh).
1. Babugka et al. / Quality assessment o f the a-posteriori error estimation
304
1.4(-3) / / ,//
1.8(-3) /
1.0(.-3) / (1.3(-3)) /
1.4(-3)
/
1.5(-3) (1.7(-3))
1.4(-3)
/
(1"5 U
/ 1.o(-3) J
(1.3(-3))
/
1.4(-3)
// (1.6(-3))
s m a l l e l e m e n t side: 1 / 8 Fig. 13. Finite element error and estimators (in parentheses) for problem with singularity (a = 0.5, adaptively refined mesh).
Fig.
14. Finite
element error (a =0.5).
for
adapted
mesh
the energy norm--is locally (i.e. on single elements) of order h e and is governed by the second derivatives of u. The contributions of G to the energy error are of the same order, but are governed by the first derivative of G. Now observe that in our case both errors are governed by the second derivative of Im z ~/2. This leads to a cancellation and consequently to the oscillation phenomenon. To demonstrate this, we consider a patch of two triangles as in Fig. 15 with local coordinates ~:, rt with the origin in the lower left vertex of the patch. Using (35), (36) we can assume that on the patch: u = a(sc2 - .q2) + b~Zr/+ linear functions + higher-order terms, G = ( a a ~ + ab~7)h + constant term + higher-order terms,
where a is related to A in (36). Neglecting R in (35), we get by simple computations Ilell~l
(a2+bZ)(l+Za+3a2)+a
Ile1122
( a2 + b 2 ) ( 1 - 2a + 3a2) + a 2
2
(II e II Ei is the energy norm of the error in element i) and we see that we have in the case without pollution (a = 0) He 1121/[1 e ]]~2 = 1,
7/
ii
Fig. 15. Mesh for local analysis of errors.
I. Babu~ka et al. / Quality assessment of the a-posteriori error estimation
305
while t h e r a t i o is ~ 1 in t h e p r e s e n c e o f t h e p o l l u t i o n error. This explains o u r o b s e r v a t i o n in t h e p r e v i o u s section. L e t us n o t e t h a t t h e said effect occurs w h e n t h e p o l l u t i o n e r r o r a n d t h e e r r o r o f the i n t e r p o l a n t have t h e s a m e strength. This is in a g r e e m e n t with o u r results in the previous section, w h e r e the oscillation p h e n o m e n o n was only visible for t h e singularity a = 0.5.
Conclusions S u m m a r i z i n g t h e results t h a t we have p r e s e n t e d in t h e p r e v i o u s sections we d r a w t h e following conclusions. (1) T h e e s t i m a t o r p e r f o r m s a e x p e c t e d f r o m its p r o p e r t i e s given in T h e o r e m s 1 a n d 2. (2) T h e effectivity i n d e x d e p e n d s on t h e mesh, especially on the a n g l e o f t h e triangles. H e n c e m e s h e s w i t h o u t small angles a r e p r e f e r a b l e . T h e n o t i o n o f t h e angle d e p e n d s on t h e differential operator. (3) T h e effectivity index is n o t to sensitive to t h e t o p o l o g y of the m e s h ( e x c e p t for t h e m i n i m a l angle). (4) T h e effectivity index can b e l a r g e r o r s m a l l e r t h a n 1 d e p e n d i n g on t h e c h a r a c t e r o f t h e s o l u t i o n as p r e d i c t e d in T h e o r e m s 1 a n d 2. (5) If t h e s o l u t i o n has s i n g u l a r b e h a v i o u r , the quality o f t h e e r r o r e s t i m a t o r d e t e r i o r a t e s . T h e d e t e r i o r a t i o n c a n b e a v o i d e d by a d a p t i v e r e f i n e m e n t , as follows f r o m T h e o r e m 1.
Acknowledgements T h e w o r k o f t h e first a u t h o r (I.B.) was p a r t i a l l y s u p p o r t e d by t h e O N R u n d e r g r a n t N00014-90-J-1030. T h e w o r k o f t h e o t h e r a u t h o r s (L.P. a n d R . R . ) was p a r t i a l l y s u p p o r t e d by t h e N S F u n d e r g r a n t CCRo88-20279.
References [1] I. BABU~KA,R. DUR.~N and R. RODRfGUEZ, "Analysis of the efficiency of an a-posteriori error estimator for linear triangular elements", SIAM J. Numer. Anal., to be published. [2] I. BABU~KAand W.C. RHEINBOLDT,"Error estimates for adaptive finite element computations", SIAMJ. Numer. Anal. 15, pp. 736-755, 1978. [3] 1. BABU.~KAand W.C. RHEINBOLDT,"A posteriori error estimates for the finite element method", Int. J. Numer. Methods Eng. 12, pp. 1597-1615, 1978. [4] I. BAaU~V:Aand W.C. RHEINBOLDT,"Reliable error estimation and mesh adaptation for the finite element method", in: Computational Methods in Nonlinear Mechanics, edited by J.T. Oden, North-Holland, Amsterdam, 1978. [5] I. BASU.~KAand R. RODRiGUEZ, "The problem of the selection of an a-posteriori error indicator based on smoothing techniques", Tech. Rep. BN 1126, University of Maryland at College Park, Inst. for Phys. Sci. and Technol., 1991. [6] I. BABU.~KAand M. SURI, "The h - p version of the finite element method with quasiuniform meshes", Model Math. Anal. Numer. (RAIRO) 21, pp. 199-238, 1987. [7] I. BAnu.~r,A and D. Yu, "Asymptotically exact a-posteriori error estimator for biquadratic elements", Finite Elements in Analysis and Design 3, pp. 341-354. [8] I. BAaU~KA,O.C. ZIENrdEWlCZ,J. GAGO and E.R. OE A. OLIVEIRA(eds.), Accuracy Estimates and Adaptiue Refinements in Finite Element Computations, J. Wiley, Chichester, UK, 1986. [9] P.L. BAEHMANN,M.S. SHEPHARD and J.E. FLAHERTY, "A-posteriori error estimation for triangular and tetrahedral quadratic elements using interior residuals", Tech. Rep. SCOREC Rep 14-1990, Rensselaer Polytechnic Institute, Troy, NY, 1990.
306
1. Babu,~ka et al. / Quality assessment of the a-posteriori error estimation
[10] R.E. BANK and B.D. WELFERT,"A-posteriori error estimates for the Stokes equation: A comparison", Comput. MethodsAppl. Mech. Eng. 82, pp. 323-340, 1990. [11] R. DURAN, M.A. MUSCHmTT~ and R. RODRI¢;UEZ, "'On the asymptotic exactness of error estimators for triangular finite elements", Numer. Math. 59, pp. 107-127, 1991. [12] R.E. EWING, "A-posteriori error estimation", Comput. Methods" Appl. Mech. Eng. 82, pp. 59-72, 1990. [13] J.T. ODEN, "Adaptive finite element methods for problems in solid and fluid mechanics", in: Finite Elements, Theory and Applications, edited by D.L. DwovH<, M.Y. HUSSAINI and R.G. VOIG'I', Springer Verlag, New York, 1988. [14] J.T. ODEN, L. DEMKOWITZ, L. RA('HOWITZ and T.A. WESTERMANN, "Towards a universal h - p adaptive finite element strategy. Part 2: A-posteriori error estimation", Comput. Methods AppL Mech. Eng. 77, pp. 113-180. 1989. [15] J.T. ODEN, L. DEMKOWITZ, L. R A c n o w n z and T.A. WESTERMANN, "A-posteriori error analysis in finite elements: The residual method for symmetrizable problems with applications to compressible Euler and Navier-Stokes equations", Comput. Methods Appl. Mech. Eng. 82, pp. 183-204, 1990. [16] L. PLANK, E. STE1N and D. BISCHOFF, "Accuracy and adaptivity in the numerical analysis of thin-walled structures", Comput. Methods Appl. Mech. Eng. 82, pp. 223-256, 1990. [17] M.S. SHEPHARD, Q. NIu and P.L. BAEHMANN, "Some results using projectors for error indication and estimation", in: Adaptice Methods in Partial Differential Equations, edited by J.E. FLAHERT',', P.J. PASt,OW, M.S. SHEPrtARD and J.D. VASILAKIS,Philadelphia, PA, pp. 83-99, 1989. [18] T. STROUBOUUS and K.A. HAOUE, "Recent experiences with error estimation and adaptivity. Part I. Review o1 error estimators for scalar elliptic equations", Comput. Methods Appl. Mech. Eng., to be published, 1992. [19] T. STROUBOULIS and K.A. HAOUE, "Recent experiences with error estimation and adaptivity. Part II. Error estimation for h-adaptive approximations on grids of triangles and quadrilaterals", Comput. Methods Appl. Mech. Eng., to be published, 1992. [20] R. VERFORTrt, "A-posteriori error estimators for the Stokes equation", Numer. Math. 55, pp. 309-325, 1989. [21] L.B. WArtLBIN, "Local behaviour in finite element methods", in: Handbook of Numerical Ana{vsis, edited by P.G, CIARt.ET and J.L. LIONS, North-Holland, Amsterdam, Vol. 2, pp. 353 522, 1991. [22] J.Z. ZHU and O.C. ZIENKIEWlCZ, "Adaptive techniques in the finite element method", Commun. AppL Numer. Methods 4, pp. 197-204, 1988. [23] O.C. ZXENK1EWlC'Z and J.Z. ZHtJ, " A simple error estimator and adaptive procedure for practical engineering analysis", Int. J. Numer. Methods" Eng. 24, pp. 337-357, 1987. [24] O.C. ZIENKIEWlCZ and J.Z. ZHU, "The three R's of engineering analysis and error estimation and adaptivity", Comput. Methods Appl. Mech. Eng. 82, pp. 95-114, 1990.
285
Elsevier FINEL 250
Quality assessment of the a-posteriori error estimation in finite elements Ivo Babu~ka, Lothar Plank a and Rodolfo Rodriquez a,b a Institute for Physical Science and Technology and Department of Mathematics, University of Maryland, College Park, MD 20740, USA b Departamento de Matemdtica, Facuhad de Ciencias Exactas, Universidad Nacional de La Plata, C.C. 172, 1900 La Plata, Argentina Received October 1991 Revised March 1992 Abstract. The paper addresses results of numerical experimentation with an a-posteriori error estimator. The
paper concentrates on the problem of the selection of benchmarks which is directed to a numerical verification of basic theoretical features of the estimator.
Introduction
During recent years much interest has been focused on the design of a-posteriori error estimates in the finite element method, their experimental verification and their use for adaptive procedures. We refer for example to [1-24] and further citations there. The a-posteriori error estimators are often derived by purely heuristic reasoning and without mathematical analysis. If a mathematical analysis is made, it adresses most often only asymptotic properties and special meshes. Usually the estimators are numerically analyzed (verified) on a set of examples (benchmarks) which are selected at best by heuristic reasoning. Often these benchmark computations are used for the design of a correction factor which is used to improve the quality of the estimators for this set of benchmarks. This pragmatic approach to test and compare the estimators by benchmark computations has serious shortcomings. The quality of an estimator is usually very sensitive to the structure of the solution and the meshes used. The performance of an estimator can be very different. It may depend on whether the meshes are nearly translation invariant or general, or whether the meshes are general or such that the error indicators are nearly equal on all elements (i.e. the meshes are equilibrated), or whether the solution is smooth, has point singularities or is generally unsmooth, etc. Some of these aspects are adressed in [5]. Hence the benchmark computations could motivate misleading conclusions. Each benchmark computation should be a representative of a more or less precisely defined class of computational problems and the conclusions used for this class of problems. In this paper we will address various aspects of the relation between theoretical understanding of an estimator and the benchmark selection. We will adress a certain estimator which is used in practice, and discuss its properties. We do not intend to compare it with any other error estimator. Correspondence to: Ivo Babu~ka, Institute for Physical Science and Technology and Department of Mathematics, University of Maryland, College Park, MD 20740, USA. 0168-874X/92/$05.00 © 1992 - Elsevier Science Publishers B.V. All rights reserved
286
1. Babu~ka et al. / Quality assessment of the a-posteriori error estimation
Model problem ket /2 ~ E2 be a b o u n d e d polygonal domain and let 1" be its boundary. A s s u m e further that .0 = U i~_ t~i, where gl i is a polygonal domain (with b o u n d a r y a/2 i) and that /2i c~ S'lj = 0 for i 4:j. We will be interested in solving the following problem:
-div(aVu)=f 0u a--=g an
in/2,
on I~:,
u=0
on I)~,
l't~Ul~.=l',
(1)
where n is the o u t w a r d normal to /2. We assume that f ~ H ~(/2), g ~ H~(1)~), a >/a 0 >~ 0 is constant on /2i, i = 1. . . . . m. We will assume that either m e a s ( F o) ~ 0 or, if F D = 9i. f a r d x d y + frg ds = 0 and u ( A ) = 0, A ~ . (Because of the assumption about f and g. u ~ H~(/2), a > 1, and so u ( A ) = 0 makes sense.) We will consider the weak solution of (1). To this end we d e n o t e by Hk(/2) the usual Sobolev space and let H,l{u
= 0 on
F u r t h e r we let
B ( u , /')
=-~f
SU 0t! 0U ~)l') ff - - - - + - - - dx dy ax ax
ay a y ,
(2)
be the bilinear form defined on H(}(/2) × H~I(~) and
X(t!) = f~ft' d x d y +
frgL,
ds
(3)
be a linear functional on H(}(/2). T h e n the weak solution u of (1) in H(I exists, is unique and
B(u,c)=F(t')
Vc~H,I.
(4)
If 1"0 = ~, we consider H I ( / 2 ) instead of H01(/2) and assume u ( A ) = 0, which is well-defined, as we have in fact higher regularity for u. Let us consider a family ~ of triangular meshes M on /2 satisfying a uniform minimal angle condition. We d e n o t e by 3 - the closed triangles of the mesh M, U 3 - = ~ . W e will assume that if J - ~ M then i n t e r i o r ( Y ) c / 2 t for some l, i.e. 3 - lies in a domain where a is a constant. To every mesh M we associate a p a r a m e t e r t(M)>1 O, for example t ( M ) = N where N is the n u m b e r of the nodes that do not belong to 1"D or t ( M ) = h where h = m a x ( d ( J ) ) and d ( 3 - ) denotes the diameter of ~ . W e will always assume that h = m a x ( d ( 3 j ) ) ~ 0 as t ( M ) ~ O. F u r t h e r let S ( M ) be the set of all continuous functions o n / 2 whose restriction on a g ~ M is a linear function. Obviously S( M ) c Hl(g2). Let So(M) = H~(12) C~S ( M ) and S A ( M ) = {u S ( M ) , u ( A ) = 0}. T h e n the finite element solution u s ~ So(M) satisfies
B(us, c)=B(u,u)
Vc~S0(M),
(5)
w h e r e u is the exact solution of (1). If 1"D = ~, we replace So(M) by SA(M). We will consider a family of finite element solutions u s associated to ~t" and assume that II u - US(M)tl E = B(u US(M), U -- US(M))1/2 ~ 0 as t ( M ) ~ O. By (5) we associate to any solution u its finite element solution u s = us(M, u) and an error estimator ~ ( u , M ) that will be defined in the next section. In fact ~ d e p e n d s on u s and M but because of (5) we write ~ ( u , M). This estimator approximates the energy n o r m of the error, i.e. ~ ( u , M ) = II e [[ E = (B(e, e)) 1/2, where e = u - u s . -
287
I. Babugka et al. / Quality assessment o f the a-posteriori error estimation
Error estimator
Let M ~atv, ~ / , ~ t ~ M be two triangles with common side 3"l,i as shown in Fig. 1. By nt, i we denote the outer normal from ~ t to 5~; we will also write ~ t = ° ~ i n and ~ / = 3-out, as ~in is the triangle for which the normal nt,/ is the outward one. Let us now define for every side 3"l,i
[[a~. ](x) = (aVus) l~root .nl, i - (aVus) l:~, .nl,i,
(6)
the jump of a(aUs/Ont, i) across the side 3"t,i. This value is independent of the choice of nt, ~, i.e. [[a(OUs/anl,i)](x) = ~a(bUs/OUi,t)]l(x) and so we will write ~a(OUs/On)](x). For each triangle 3- ~ M let
R : = - [div(aVu) I ~-+ div(aVus) l : ] + f I : .
(7)
We may write for any v ~ H d the residual equation
B(e,v)=
~ f e:vdxdy+ 5r~M
g-a
vds+ ~ J ~ a T l l v d s ,
"~
yEG
3'
(a)
Otl
where G stands for the set of all interior edges of M and n for the outward normal to 12. (If
Fo = 0 eqn. (8) is the same, but v ~ Hi(12).) For any 3 - ~ M we define 1
i----(fR:
11:=
dx d y
(9)
where t3-1 is the area of 3". (1I: is the L 2 ( 3 -) projection of the interior residual onto the set of constants). Further for each 3' ~ FN, where 3' is a side of a triangle ,Y, we define
1
aUs I
l-l, = T i f,( g - -Tnn l dS,
(10)
where 13" I is the length of 3'. ( H , is the L2(y) projection of the boundary residual onto the set of constants.) For each edge y of the mesh let
{
[[2(~Us/On)]]
j~=
H~,
if y ~ G,
if 3'c_rN,
(11)
if y c F o ,
and for each triangle 3 - ~ M we define the error indicator r / : : n : :=
13-1
"'S~dxdy+7
~
13/I a--i--~ds
,
(12)
Fig. 1. Scheme of neighbouring elements.
L Babu~ka et aL / Quality assessment of the a-posteriori error estimation
288
which approximates the energy error on ~5r. Further we define the error estimator ~ ( u , M ) :=
r/.~-
•,
(13)
" 17 ~ M
where K is a suitable constant which will be defined later. It is convenient to write g*(u,
M ) :=
Y[ r/.~
(14)
.'f ~ M
for the choice K = 1. Let us note that the error indicator is meaningful for more general f and g than we have assumed. Nevertheless our generality is sufficient for practical purposes. Obviously the set of problems we are discussing is precisely described.
Quality of the error estimator
In the previous section, we defined the error estimator g~(u, M ) which approximates IIe IIE. In the following we will discuss the relation between g'(u, M ) and I1e II E, especially the effectivity index of the estimator ~: = ~ ( u , M ) / l l e l l t . : . Before addressing this problem we introduce some notation which we will need in the following. Let S P ( M ) be the space of all continuous function on J2 whose restriction on 3~M is a polynomial of degree p. Let S P ( M ) = S p h i l l ( s O ) and u~ ~ S g ( M ) be the finite element solution of problem (1), i.e. B(u}',v)=B(u,v)
Vv~Sg(M).
(If F D = fl, we use S P ( M ) instead of S g ( M ) . ) Obviously Uls' is the p-version solution of problem (1). For more about the h-p-version we refer to [6]. We have for any k > 1:
II U~s'
-
u
II ~ ~ C p - ¢ ~ - j) II u It Hkh min{p'k- I}.
(15)
with C dependent on k, but not on p and h. Remark. Although in (15) the maximal element side is used, an analogous element by element estimate holds. Now we are able to state the basic theorems about N * ( u , M). (See [1] for the detailed proof and the next section for an outline of the main ideas.) Theorem 1. For any integer p there is a constant K f depending only on the minimal angle ce.i o f the triangles used and on p, such that
Il e HF. ~ K f ~ * ( u, M ) + Hu - u~ ll E + O( h 3/2)
(16)
and the constant K f grows logarithmically with p. Theorem 2. that
There is a constant K l > 0, depending only on %7 and on the jumps of a, such
~ * ( u , M ) < ~ K l llell E + O( h3/2).
(17)
L Babugka et al. / Quality assessment o f the a-posteriori error estimation
289
Table 1 Constants C~ and C~- for various angles a :
a
C:
7.5 15.0 22.5 30.0 37.5 45.0 52.5 60.0
0.051 0.072 0.087 0.099 0.108 0.115 0.119 0.121
C.~
t
p=2
p=4
p=6
p=8
2.390 1.682 1.363 1.169 1.035 0.939 0.875 0.850
3.306 2.341 1.918 1.670 1.508 1.400 1.334 1.309
3.660 2.609 2.156 1.895 1.727 1.615 1.547 1.522
3.988 2.839 2.343 2.058 1.876 1.757 1.684 1.657
The constants K~ and K~' are computable. We have K~ = max C~-,
(18)
3~M
K 1 = ( max C~r)Q,,
(19)
Qa =
(20)
where
max
:ng-'*O
(al:/al~,).
For 2 ~
(21)
3.45 sin-l/2( ot:/2) <~C~.~ 5.85 sin-1/2( a : / 2 ) ,
(22)
and
where t~9. is the minimal angle of 3 . In Table 1 we give the values of C~- and C~-. Let us note that for h small only very few elements have a side on the interface 092; we can in fact use Qa = 1. Remark. We have dealt only with the Laplace equation. In this case we use the minimal angle % r of the triangle. For a general differential equation we have to replace it by ot~which is the minimal angle of the triangle after a local transformation of the general equation to Laplace's equation (this is always possible). The t e r m O ( h 3/2) in (16) and (17) depends only on II f tl nl(n) and II g II nkr), but not on u. For any u the term II u - uff II E can be made arbitrarily small by (15) using sufficiently large p. We note that the smoothness of u is governed by the smoothness of f and g, transition of boundary conditions and smoothness of F. It is unessential that the estimates (16) and (17) have an asymptotic character. For example if f = const on every 3 - and g = const on every triangle side of FN, the term O(h 3/2) disappears. Hence as t ( M ) ~ O, the term O(h 3/2) is negligible in general. Theorems 1 and 2 indicate that the error estimator ~ * ( u , M ) can be either smaller or bigger than the true error. Further the error estimator loses accuracy when a 3 ~ 0 and for unsmooth u the reliability of the estimator decreases because a higher p is needed in the estimator. The coefficient x, which scales the error estimator (~'(u, M ) = K ~ ' * ( u , M)), can be calculated in different ways. We can use the geometrical average of the bounds in (18), (19).
290
1. Babugka et al. / Quality assessment of the a-posteriori error estimation
Another possibility is to choose K such that the error estimator is asymptotically exact for uniform meshes of equilateral triangles and f = 0. In this case we get x = ( 8 ~ ) - 1 / 2 = 0.26864 which will be used for all examples. Theorems 1 and 2 indicate the expected theoretical properties of the estimator well. The structure of Theorems 1 and 2 clearly suggests the set of benchmark experiments that should be performed.
Outline of the proof of Theorems 1 and 2 For simplicity assume that f and g are constant on each g resp l ~ F N. (The general case is discussed in [1].) With ep := u~ - u S E S(p, we have I l e l l E < ][ep]lE+ ]]U--Uf][E
(23)
II epll2: = B ( ep, ep) = B(uJs' - u, % ) + B( u - u s, % ) = B ( e, % ) .
(24)
and
For any continuous function v let v I denote the linear interpolant of v on 3 . Using (8), we have:
y"fJ'(e"-4) dsl
l,ep[,~=B(e,ep)=.,ic. ~ [fY :x(ep-elp) ~,2M
rlJ
~-T~(ep-elp)
1 +~,c~ 3
<
Y'. J~M
dxdy
.......
Va l j ~-jl
/ca..
r" ( e t , - e I ) ds
2) )/ 2
~JC~-II ep II t.:,J,
where
1 -~lfT(v-
CP~-=
sup
Vl) d x d y I 2 + ' ~
'
"~e,,\P,,
1
I~-
f/ (U - UI) dsl 2
) 1/2
Ica~
f lV v l
2 dx dy
Here, Pp and P0 denote the space of all polynomials of degree p or 0 on 5r, respectively. Hence ]l epll E < ~ K ~ g ' * ( u , M ) , with K f = sup C~-.
J~M
Using (23) we get (16) without the term O(h3/2). This is related to the simplification assumption about f and g made in this section. Let us now address inequality (17). For any triangle Y as shown in Fig. 2, let ¢0 be the cubic bubble function vanishing on a3- with ¢0(Q0) = 1, where Q0 is the barycenter of 3-. Let
L Babu~ka et al. / Quality assessment o f the a-posteriori error estimation
291
Fig. 2. Notation for triangular elements.
further ~ ° ( 3 - ) be the space of quadratic functions vanishing at the vertices of .Y- and let {~°t}3=1 be the canonical basis of this space (i.e. q~i(Qj)= ~ij, i = 1, 2, 3, j = 1. . . . . 6) and let V: be the space spanned by {~ot}~=0. Set w ~ = E3=oWiq~i~ V: with coefficients w0, Wl, w2, (-03 chosen such that /i2
dx d y = [.Y-I J[ y "': dx dy, a ftJtoo: ds = •
where
j2 li 11i [ f 7 - ' T d s ,
(25)
i = 1, 2, 3,
(26)
li[aJt i
[a]t, is the average of a over the triangles sharing the edge li, i.e. / =
[a]ti
| /
x
I 2(al~'i"+al:°u') ~a[:
if li = a~in ~ &Y-°ut'
(27)
if lic&Y'NF.
The coefficients w 0. . . . . o)3 satisfying (25) and (26) are OJi
t°O =
2[a]l i '
a-----~--:
V
2 i= 1 [ a ] l i
li
).
Since for any interior edge l the value wv) are the same on both triangles 3- sharing the edge, the function w composed from w : on every triangle is continuous and hence oo~ H~I'2). Using (23), (24) we get: ( g ~ * ( u , M ) ) 2=
]~ r/2 3-~M
17z " ' ~ d x d y + ~ 1 E I11 ,: a leo:
j2
]
=E
13-If
=E
[a]t /j~o f J l f ~ o : dx dy+ ~1 ,~.~a~aljl: , : ds ]
,~~M
ds
~-~M
-<
where [a]t Co =
max
1~&fr,3-~M
a I ~-
Since oJ ~ H i we get from (8)
(g'*(u, M))2<~QB(e, w) < C~ Ilell Elleo lIE.
(28)
I. Babu~ka et al. / Quality assessment o f the a-posteriori error estimation
292
3 Now let for r = {r i}i=o
C9=
i l Vt.,r l I12
sup
"
[Irll.~O
dx dy
3
r 2
,
1
o + 5 E
ri-
i=1
where v~ E E~ is the function satisfying f~v r d x dy = I J l r o , T h e r e f o r e w z is the function v r corresponding to
I lilJz, [a]li
IJ-l/Tj
r° -
- - - ,
r i --
alj
,
fliUr ds = I li[r i, i = 1, 2, 3.
i=1,2,3
and so
I1~o~ IIu2 . , z - a [ . ~ f lIVVrtl 2 d x d y
2 ,
~ c;~,l~( =c.~( IJ
a]: 2
211'2j
+ Z121 3
1{ 2
i=l
[a]6
( a~ll~, )
l
01:7
a
f.,2p,2 ~< t _ . j ~ a
2
7"/aw,
where C~ =
alj
max ,c~J,~M [a],"
Hence II o)I[ e ~< { max " ,Y-~M
C~-}C'~'*(u, M)
and so, using (28), we obtain ~*(u,
M ) ~< Qa( max C ) ) II etl ~, J~M
where
Qa = C.C~
a[:7 = max - .gem a [ ~,
and we get (17).
Experimental results In this section we will experimentally verify the reliability of various conclusions that can be m a d e from T h e o r e m s 1 and 2 for the given estimator.
Influence of mesh topology T h e o r e m 1 was proven for general meshes without any restrictions on the mesh topology. H e n c e , it can be expected that particular m e s h e s with different topologies could influence the
I. Babu~ka et aL / Quality assessment of the a-posteriori error estimation
// //
W
293
//
mesh ratio 1/1
//// //// mesh ratio 1/2
Fig. 3. Regular three-direction mesh.
effectivity index in the same order of magnitude as indicated by the theoretical upper and lower bounds. In addition the influence of the angle has to be visible (and could also depend on the topology). To check the validity of this conclusion for a smooth solution we consider the following problem: For 12 = (0, 1) x (0, 1), let u be the solution of -Au =f
in 12,
u =0
on ~e.
(29)
We choose f such that u = sin wx sin -try.
Three-direction meshes Let us consider the three-directional meshes which are obtained by uniform refinement of a basic mesh of the type 1/n, n = 1, 2, 4, 8 as shown in Fig. 3. In Table 2 we report the effectivity index for these meshes. We also report the u p p e r and lower bounds for the effectivity index from Theorems 1, 2, when p = 2 was used.
Criss-cross meshes Let us consider now the "criss-cross" meshes as shown in Fig. 4. The effectivity indices for these meshes are given in Table 3.
Diamond-shaped meshes Finally, we consider the "diamond-shaped" meshes as shown in Fig. 5. The effectivity indices are given in Table 4.
Table 2 Effectivity index for problem (29) (u = sinrrx sin~ry) and 3-direction meshes Number of elements
Mesh
8 16 32 64 128 256 512 1024 2048 4096
0.852
Upper bound Lower bound
2.34 0.29
1/1
1/2
1/4
1/8
0.993 1.070
1.289 1.183
1.147
1.804 1.503
1.252 1.180
2.068 1.571
1.277 1.192
2.145 1.594
1.287 2.87 0.21
2.169 3.87 0.15
5.38 0.12
1. Babu~ka et al. / Quality assessment of the a-posteriori error estimation
294
/ \
,
\/\
\/\/ /\/\
mesh ratio 1/1
, mesh ratio 1/2
Fig. 4. Regular criss-cross mesh.
Table 3 Effectivity index for problem (29) (u = sin'rrx sin'rry) and "criss-cross'" meshes Number of elements
Mesh
8 16 32 64 128 256 512 1024 2048 4096
1.35(I
Upper bound Lower bound
1/1
1/2
1/4
I/8
{L939 1.1136
1.1}49 1.053
1.077
1.349 1.179
1.{185 1.088
1.518 1.214
1.093 1.I19(I
1.562 1.223
1.I195 2.34 I).29
2.87 {}.21
1.574 3.87 0.15
5.38 I). 12
Comparing the tables we can make the following conclusions, which are in agreement with Theorems 1, 2: (1) For h ~ 0 the effectivity index stays in the theoretical bounds. In addition we see that in our concrete case the effectivity index converges to a limiting value. This is not surprising, because the mesh is translation invariant and the solution can be periodically extended. In this case various results about superconvergence can be used [21]. (2) The effectivity index depends on the topology. The error estimator is not asymptotically exact. (3) The dependence of the effectivity index on the angle is clearly visible.
mesh ratio 1/1
mesh ratio 1/2 Fig. 5. Regular "diamond-shaped" mesh.
L Babu~ka et al. / Quality assessment of the a-posteriori error estimation
295
Table 4 Effectivity index for problem (29) (u = sinrrx sinlry) and "diamond-shaped" meshes Number of elements
Mesh
6 12 20 24 40 48 72 80 144 160 272 288 544 576 1056 1088 2112 2176 4160 4224 8448
1.465
Upper bound Lower bound
2.34 0.29
1/1
1/2
1/4
1/8
0.750 1.017 0.523 1.026 0.385 1.107 1.111 1.188 0.959 1.131 1.558 1.234 1.708 1.138 1.740 1.246 2.274 1.139 1.793 2.515 2.87 0.21
3.87 0.15
5.38 0.12
(4) The theoretical bounds for the effectivity index are relatively pessimistic. This is because these bounds are independent of the topology and the solution. Better bounds in our case could be obtained by exploiting the superconvergence effect [5].
Influence of the solution structure Theorems 1 and 2 have been shown for a general solution. We expect that the effectivity index may depend on the solution structure. To check this conclusion consider the solution of the following problem on g2 = (0, 1) × (0, 1): -Au=f au --=0
by
ing2, for y = 0 a n d
u=0 y=l
for x = 0 and x = l, (30)
with f such that u = sin ~ x . We will use meshes as in Fig. 3 with the orientation shown there (mesh 1 / 2 and 1 / 4 ) and with the opposite orientation ( 2 / 1 and 4 / 1 ) . From Table 5 we see in fact that the effectivity index depends on the relation between the mesh and the structure of the solution (especially the relation between the second derivatives and the mesh orientation). Hence in general we have to deal with uncertainties in the effectivity index. (No correction factor could improve the situation.) We note that our factor K = 0.2684... could be changed, but in general (for all funotions) the range in the effectivity index could not be improved. Again, we see that Theorems 1 and 2 reliably predict the outcome of the numerical experiments.
296
I. Babugka et al. / Quality assessment 0["the a-posteriori error estimation
Table 5 Effectivity index for problem (30) (u = sin'n'x) and 3-direction meshes Number of elements
Mesh
8 16 32 64 128 256 512 1024 2048 4096
1.013
Upper bound Lower bound
2.34 0.29
1/1
l/2
1/4
4/1
2/1
1.564
[).717
1.1117
2.258
0.5117
1.598
/).783
1.131
2.273
11.553
1.608
il.8011
1.137
2.278
().566
1.611
11.804
1.139
2.279
11.569
1,612
0.806
2.87 0.21
3.87 0.t5
2.87 0.21
3.~7 1/.15
Influence o f the smoothness o f the solution T h e o r e m s 1 a n d 2 c o n t a i n h i g h e r - o r d e r t e r m s w h i c h we n e g l e c t e d in t h e r e p o r t e d t h e o r e t i c a l b o u n d s in T a b l e s 2 - 5 w h e n c o m p u t i n g t h e b o u n d s . T h e s e t e r m s i n f l u e n c e t h e c o n v e r g e n c e o f t h e e f f e c t i v i t y i n d e x as h --* 0. T o c h e c k t h e e f f e c t o f t h e s e t e r m s w e will c o n s i d e r o n c e m o r e p r o b l e m (29), u s i n g n o w t h e s o l u t i o n
u ( x , y ) = sin a'rrx sin a w y .
(31)
( W e use r e g u l a r t h r e e - d i r e c t i o n m e s h e s , as in Fig. 3 w i t h m e s h r a t i o 1 / 1 . ) We see from Table 6 that the insufficient smoothness influences the effectivity index and that the effectivity index can behave erratically. Nevertheless the bounds--where higher order terms are neglected--are still valid. W e s e e in g e n e r a l t h a t t h e e f f e c t i v i t y i n d e x decreases with decreasing smoothness of the solution. The lower bound goes down with i n c r e a s i n g p , b u t t h e u p p e r b o u n d is n o t i n f l u e n c e d . H e n c e w e c a n e x p e c t t h a t t h e e f f e c t i v i t y i n d e x will go d o w n ( a s s u m i n g h e u r i s t i c a l l y t h a t t h e e f f e c t i v i t y i n d e x m o v e s as t h e m i d d l e o f t h e interval). T h e s a m e e f f e c t will b e s e e n later. O n c e m o r e w e s e e g o o d a g r e e m e n t w i t h t h e c o n c l u s i o n s s t e m m i n g f r o m T h e o r e m 1 a n d t h a t w e c a n n o t h o p e for t h e e x i s t e n c e o f a universal correction factor.
Table 6 Effectivity index for problem (31) (u = sinawx sina'rry) and 3-direction meshes Number of elements
a= 1
a=2
a=4
8 32 128 512 2048 8192
0.852 1.071/ 1.147 1.180 1.192
1/.833 0.876 1.080 1.161 1.187
0.438 0.932 0.887 1.090 1.168
Upper bound Lower bound
2.34 0.29
2.34 0.29
2.34 0.29
a = 1l
0.379 0.5311 0.7811 1.022 1.148 2.34 0.29
L Babu~ka et al. / Quality assessment o f the a-posteriori error estimation
297
Table 7 Effectively index for problem (32) (u(x, y) = sinarx sinh'rry) and 3-direction meshes Number of elements
Mesh 1/1
1/2
1/4
1/8
0.841
8
16 32 64 128 256 512 1024 2048 4096
0.958 0.978
1.269 1.102
1.040
1.795 1.457
1.166 1.065
2.042 1.534
1.190 1.071
2.136 1.560
1.199
Upper bound Lower bound
2.34 0.29
2.87 0.21
2.166 3.87 0.15
5.38 0.12
Harmonic solutions The error estimator contains terms related to the residual and the jumps, respectively. To isolate the influence of the jump term, we consider the following problem (again i n / 2 = (0, 1) x (0, 1)): -Au=0
inO,
u=0
forx=O,x=landy=O,
0u
-- =g On
for y = 1,
(32)
and let u(x, y) = sin'rrx sinh~-y. The effectivity index (for a regular three-direction mesh as shown in Fig. 3) is given in Table 7. Comparing Tables 2 and 7 we see a similar behaviour influenced only mildly by the structure of the solution. The limiting value for the effectivity index seems to be different for the mesh 1/1. For the other meshes this effect is less visible. For these particular meshes, this might be explained by a superconvergence phenomenon.
Adaptively constructed meshes So far, we have restricted our experiments to uniform meshes. We will now consider nonuniform meshes either based on an adaptive procedure or randomly refined. We use problem (32) again. We have constructed three sequences ( A - C ) of meshes. Each sequence is obtained by refining a part of the elements of the previous mesh according to a refinement indicator and adjusting neighbouring elements. In sequences A and B the refinement indicator was an error estimator (with different parameters), in sequence C the indicator was a random number. Typical meshes of each sequence are shown in Fig. 6. The effectivity indices are given in Table 8. Comparing Tables 8 and 7 we see that different refinement strategies lead to different effectivity indices, although the difference is not big. Once more we see a range of the effectivity indices. We note that for these general meshes superconvergence results cannot be applied.
I. Babu,~ka et al. / Quality assessment of the a-posteriori error estimation
298
/l//2/
,G,'/\
Y// y//yy
///
/
/ / /
..-v z
/
(b) ---
//
•
/
///
/
"
k
/,"
\ /
,
~ "
\\\
/
Fig. 6. Nonuniform meshes for problem (32).
Effectiuity indices f o r problems with unsmooth solutions So far, we have considered problems with a s m o o t h solution. L e t u s now e x a m i n e the b e h a v i o u r of the error estimator for solutions with singularities. T h e o r e m s 1 and 2 indicate a
L Babu~ka et a L / Quality assessment of the a-posteriori error estimation
299
Table 8 Effectivity indices for problem (32) and nonuniform meshes Sequence A
Sequence B
Sequence C
Elements
Effectivity index
Elements
Effectivity index
Elements
Effectivity index
8 19 47 94 181 323 593 1020 1773
0.841 0.944 0.970 1.008 1.018 1.041 1.043 1.060 1.051
8 19 38 67 107 201 355 607 991
0.841 0.944 0.941 0.980 0.949 1.013 1.039 1.032 1.056
8 23 45 96 200 435 964 2129 4588
0.841 0.906 0.964 1.008 0.979 0.990 0.972 1.008 0.992
2 a
a(x,y) =
3
1
4
a(x,y) = 1 a ( x , y ) = Fig. 7. Coefficient a for problem (33).
decrease of the effectivity index, w h e n the singularity of the exact solution grows, and also that this effect diminishes when adaptive meshes are used. T o this purpose we consider the domain O = ( - 1, 1) × ( - 1, 1), partitioned into 4 square subdomains as shown in Fig. 7. Let u be the solution of
-V(a(x, y)Vu) =0
in O,
0u
a(x, Y)-~n=g(x,
where a(x, y) is constant in each quadrant chosen so that U(X, y ) = r a ( c i
COS(Ot~)÷s i
y)
on F,
(33)
(a(x, y ) = 1 or a(x, y ) = 5, see Fig. 7). 5 is
sin(a~o)),
(34)
where (r, ~0) are polar coordinates with center in the origin. For a = 0.5 and a = 0.1 we get 5 = 3 + v~- = 5 . 8 2 8 . . . and 5 --- 166.447 . . . . respectively.
Table 9 Effectivity indices for problem (33) (u = r~fQo)) and 3-direction meshes Number of elements
a = 0.5
a = 0.1
32 128 512 2048 8192
0.658 0.694 0.711 0.719 0.722
0.520 0.512 0.524 0.536
300
I. Babu~ka et al. / Qualiw assessment o f the a-posteriori error estimation
//LK_
~V
\\// ,/ \ \ / / // /. \\ \\\\/////~;'i_ \\///// /J" / V /
/\\\\
/
\\X
~ / / F S/ / ~/ ~
J~ /// \\ \\ \\ \\ \\
,. ;:7 ~ ~ , \ )
F"
(al
I
]
]
"
"\
[
Fig. 8. Nonuniform meshes for problem (33).
Table 9 shows the effectivity index for uniform meshes. (In each quadrant the meshes are uniform three-direction meshes of mesh ratio 1/1 as in Fig. 3. The orientation of the meshes is changed in subdomains 2 and 4 to preserve the symmetry of the solution.) Although the estimator in Theorem 1 includes the ratio of a, if influences only elements on the interfaces and hence this factor can be neglected. The decrease of the effectivity index is clearly related to the necessity to use a higher degree p in (16) for the estimate. Hence the unsmoothness of the solution decreases the effectivity index. This has also been observed in Table 6. We have seen in Table 8 that nonuniformity of the mesh does not influence the effectivity index too much. On the other hand an adaptive procedure refines the mesh in the places where the solution in unsmooth and hence a lower p can be used. We have created two sequences of adaptively refined meshes (A and B). The meshes B1 and B2 are less refined than the meshes A1 and A2. The refinement in meshes B is more concentrated around the origin than in meshes A. Examples of each sequences are displayed in Fig. 8. If the mesh is properly refined, a lower p can be used and the effectivity index improves. Table 10 confirms this conclusion (for c~ = 0.5).
L Babu~ka et al. / Quality assessment of the a-posteriori error estimation
301
Table 10 Effectivity indices for problem (33) (u = r~f(~)) and adapted meshes Sequence A
Sequence B
Elements
Effectivity index
Elements
Effectivity index
32 87 194 370 725 1297 2346
0.658 0.712 0.738 0.761 0.787 0.807 0.828
32 64 96 168 244 363 485 604 768 965
0.658 0.719 0.787 0.830 0.880 0.918 0.960 0.959 0.975 0.995
The effect of the unsmoothness of the solution can also be observed if the singular part o f the solution is removed. So we consider problem (33) in the domain l ) = O - [ - 0.5, 0.5] x [ - 0 . 5 , 0.5] and prescribe a(au/an)=g o n aO. T a b l e 11 shows the effectivity indices for this problem using uniform meshes as for Table 9. Again, a = 0.5. W e note that the effectivity index is now very similar to problems (29), (30)." The oscillation
effect
It is often assumed, that the error of the interpolant is close to the error fo the finite element method. To demonstrate that this assumption is in general incorrect, we show that the error of the finite element method could oscillate. To this end we show in Fig. 9 the typical error behaviour for the problem discussed in the previous section with uniform meshes and a = 0.5. Displayed are two patches of regular three-direction meshes at the same location - - f a r from the singularity--in O. Figure 9 also shows the error indicators (in parentheses). In Fig. 10 we illustrate graphically the error behaviour for the entire mesh where the mesh has 363 elements. (An element is black if the energy norm error in this element exceeds a certain threshold. Here, the threshold is 4.08 x 10-4.) We see that the error indicator shows approximately the average error for adjacent elements and hence does not exhibit the oscillation of the error of the finite element solution. Note that the interpolant does not show the oscillation. This oscillation effect is caused by the pollution of the error stemming from the singularity of strength r ~/2. If we enforce a smooth solution by removing the center part of the domain, as in problem (33), the results given in Fig. 11 show that the oscillation disappears. (We use elements from the same location as in Fig. 9; the numbers are therefore directly comparable.)
Table 11 Effectivity indices, for problem (33) (u = r~f(~) on O
a7) and 3-direction meshes
Number of Elements
a = 0.5
24 96 384 1536 6144
0.964 1.027 1.059 1.070 1.073
L Babu~ka et al. / Quality assessment o f the a-posteriori error estimation
302
(1.2~/
2.2(-3) /
2.5(-3) / (1.4~)/
/ 1.o(-3)
/
1.1(-3)
/ (1.2(-3))/ / (1.4(-3)) 2.4(-3) / 2.2(-3) / (1.0~)/
/
(1.2U
/9.7(-4) (1.0(-3)) / /
5.6(-4) /
5.9(-4) /
5.6(-4) /
5.9(-4)
(3.2~ (3.°U / 2.7(-4) /// 2.6(-4) (3.o(-4)) / (3.2(-4))
(2.7U
/ 2.6(-4)
/ i.i(-3) (1.2(-3))
small element side: 1/8
///
(2.7(-4))
(3.o~ / 2.6(-4)
/ (3.0(-4))
small element side: 1/16
values in (..): error estimator Fig. 9. Oscillation of the finite e l e m e n t error.
The oscillation effect is weaker, if a 4:0.5 or the mesh is constructed adaptively. We show in Fig. 12 the local results for a = 0.8 and in Fig. 13 the results for a = 0.5 and adaptively refined meshes. In Fig. 14 we show the error behaviour for an adaptively refined mesh and the threshold 2.0 × 10 -3. In the next section we will examine this effect. Analysis of the oscillation phenomenon
The oscillation phenomenon we have observed in the previous section occurs also in the case of the following problem on O = [ - 1, 1] × [0, 1] - A u = 0 inO, u = 0 onFD={(x,y)lO<~x<~l,y=O }, Ou --=g On
on F N = F - F ~,
Fig. 10. Finite e l e m e n t e r r o r ( a = 0.5, u n i f o r m meshes).
I. Babu~ka et at / Quality assessment o f the a-posteriori error estimation
(3.2~/ (3.o~// 2.8(-4)
1.3(-3) / 1.1(-3) / (1.4~)/ (1.2~/
/ //
/ 1.3(-3)
1.1(-3) (1.2(-3)) /
(1.4(-3))
1.0(-3)/ 1.2(-3)/ (1.1~)/ (1.3~)/ / //
1.0(-3) / 1.2(-3) (1.1(-3))/ (1.3(-3))
303
2.4(-4)
/
//
2.6(-4) / (3.0(-4)) / /
2.5(-4) /
(2.T~/
/ /"
small element side: 1/8
2.9(-4) (3.2(-4))
2.8(-4) / (3.0~/
2.5(-4) / (2.7(-4)) /
2.8(-4) (3.2(-4))
small element side: 1/16
Fig. 11. Finite element error and estimators (in parentheses) for a problem without a singularity (a = 0.5, adaptively refined mesh).
with g selected so that u =r 1/2 sin(½q~)-- Im z 1/2, with z = x + iy. We consider now a uniform 3-direction mesh of the type shown in Fig. 3. Then the finite element solution has essentially the form (35) where IIR(.,.)II H' < o(h) and G(x, y) is a piecewise linear function, which coincides in the nodal points of the mesh with the function G*(x, y) =Ar -1/2 sin ½Oh =Ah Im Z-1/2 (36) Us(X, y) =u'(x, y) +G(x, y) +R(x, y),
The function G*(x, y) is the adjoint singularity function. It describes the pollution error, which is of the same order (in II" IIE) as the error of the interpolant. The error--measured in
6.0(-4) /
7.0(-4) /
(5-1U (5.9~ / / 5.1(-4) (6.2(-4)),l / / 6.1(-4) (?.0(-4)) 5.9(-4) / (5.6~ 6.7(-4) (6.4~// / /
4.6(-4) (5.4(-4))
/
/
5.3(-4)
(6.2(-4))
small element side: 1/8
1.4(-4) / (l'3U
1.6(-4) / (1.3~/
/
/
1.3(-4) / 1.4(-4) (1.5(-4)) i / (1.6(-4)) 1.4(-4) / 1.5(-4) /
(1.4~)/ (1.5~/ / (1.4(-4)) 1.2(-4) / / (1.5(-4)) 1.3(-4) / small element side: 1/16
Fig. 12. Finite element error and estimators (in parentheses) for problem with singularity ( a = 0.8, uniform mesh).
1. Babugka et al. / Quality assessment o f the a-posteriori error estimation
304
1.4(-3) / / ,//
1.8(-3) /
1.0(.-3) / (1.3(-3)) /
1.4(-3)
/
1.5(-3) (1.7(-3))
1.4(-3)
/
(1"5 U
/ 1.o(-3) J
(1.3(-3))
/
1.4(-3)
// (1.6(-3))
s m a l l e l e m e n t side: 1 / 8 Fig. 13. Finite element error and estimators (in parentheses) for problem with singularity (a = 0.5, adaptively refined mesh).
Fig.
14. Finite
element error (a =0.5).
for
adapted
mesh
the energy norm--is locally (i.e. on single elements) of order h e and is governed by the second derivatives of u. The contributions of G to the energy error are of the same order, but are governed by the first derivative of G. Now observe that in our case both errors are governed by the second derivative of Im z ~/2. This leads to a cancellation and consequently to the oscillation phenomenon. To demonstrate this, we consider a patch of two triangles as in Fig. 15 with local coordinates ~:, rt with the origin in the lower left vertex of the patch. Using (35), (36) we can assume that on the patch: u = a(sc2 - .q2) + b~Zr/+ linear functions + higher-order terms, G = ( a a ~ + ab~7)h + constant term + higher-order terms,
where a is related to A in (36). Neglecting R in (35), we get by simple computations Ilell~l
(a2+bZ)(l+Za+3a2)+a
Ile1122
( a2 + b 2 ) ( 1 - 2a + 3a2) + a 2
2
(II e II Ei is the energy norm of the error in element i) and we see that we have in the case without pollution (a = 0) He 1121/[1 e ]]~2 = 1,
7/
ii
Fig. 15. Mesh for local analysis of errors.
I. Babu~ka et al. / Quality assessment of the a-posteriori error estimation
305
while t h e r a t i o is ~ 1 in t h e p r e s e n c e o f t h e p o l l u t i o n error. This explains o u r o b s e r v a t i o n in t h e p r e v i o u s section. L e t us n o t e t h a t t h e said effect occurs w h e n t h e p o l l u t i o n e r r o r a n d t h e e r r o r o f the i n t e r p o l a n t have t h e s a m e strength. This is in a g r e e m e n t with o u r results in the previous section, w h e r e the oscillation p h e n o m e n o n was only visible for t h e singularity a = 0.5.
Conclusions S u m m a r i z i n g t h e results t h a t we have p r e s e n t e d in t h e p r e v i o u s sections we d r a w t h e following conclusions. (1) T h e e s t i m a t o r p e r f o r m s a e x p e c t e d f r o m its p r o p e r t i e s given in T h e o r e m s 1 a n d 2. (2) T h e effectivity i n d e x d e p e n d s on t h e mesh, especially on the a n g l e o f t h e triangles. H e n c e m e s h e s w i t h o u t small angles a r e p r e f e r a b l e . T h e n o t i o n o f t h e angle d e p e n d s on t h e differential operator. (3) T h e effectivity index is n o t to sensitive to t h e t o p o l o g y of the m e s h ( e x c e p t for t h e m i n i m a l angle). (4) T h e effectivity index can b e l a r g e r o r s m a l l e r t h a n 1 d e p e n d i n g on t h e c h a r a c t e r o f t h e s o l u t i o n as p r e d i c t e d in T h e o r e m s 1 a n d 2. (5) If t h e s o l u t i o n has s i n g u l a r b e h a v i o u r , the quality o f t h e e r r o r e s t i m a t o r d e t e r i o r a t e s . T h e d e t e r i o r a t i o n c a n b e a v o i d e d by a d a p t i v e r e f i n e m e n t , as follows f r o m T h e o r e m 1.
Acknowledgements T h e w o r k o f t h e first a u t h o r (I.B.) was p a r t i a l l y s u p p o r t e d by t h e O N R u n d e r g r a n t N00014-90-J-1030. T h e w o r k o f t h e o t h e r a u t h o r s (L.P. a n d R . R . ) was p a r t i a l l y s u p p o r t e d by t h e N S F u n d e r g r a n t CCRo88-20279.
References [1] I. BABU~KA,R. DUR.~N and R. RODRfGUEZ, "Analysis of the efficiency of an a-posteriori error estimator for linear triangular elements", SIAM J. Numer. Anal., to be published. [2] I. BABU~KAand W.C. RHEINBOLDT,"Error estimates for adaptive finite element computations", SIAMJ. Numer. Anal. 15, pp. 736-755, 1978. [3] 1. BABU.~KAand W.C. RHEINBOLDT,"A posteriori error estimates for the finite element method", Int. J. Numer. Methods Eng. 12, pp. 1597-1615, 1978. [4] I. BAaU~V:Aand W.C. RHEINBOLDT,"Reliable error estimation and mesh adaptation for the finite element method", in: Computational Methods in Nonlinear Mechanics, edited by J.T. Oden, North-Holland, Amsterdam, 1978. [5] I. BASU.~KAand R. RODRiGUEZ, "The problem of the selection of an a-posteriori error indicator based on smoothing techniques", Tech. Rep. BN 1126, University of Maryland at College Park, Inst. for Phys. Sci. and Technol., 1991. [6] I. BABU.~KAand M. SURI, "The h - p version of the finite element method with quasiuniform meshes", Model Math. Anal. Numer. (RAIRO) 21, pp. 199-238, 1987. [7] I. BAnu.~r,A and D. Yu, "Asymptotically exact a-posteriori error estimator for biquadratic elements", Finite Elements in Analysis and Design 3, pp. 341-354. [8] I. BAaU~KA,O.C. ZIENrdEWlCZ,J. GAGO and E.R. OE A. OLIVEIRA(eds.), Accuracy Estimates and Adaptiue Refinements in Finite Element Computations, J. Wiley, Chichester, UK, 1986. [9] P.L. BAEHMANN,M.S. SHEPHARD and J.E. FLAHERTY, "A-posteriori error estimation for triangular and tetrahedral quadratic elements using interior residuals", Tech. Rep. SCOREC Rep 14-1990, Rensselaer Polytechnic Institute, Troy, NY, 1990.
306
1. Babu,~ka et al. / Quality assessment of the a-posteriori error estimation
[10] R.E. BANK and B.D. WELFERT,"A-posteriori error estimates for the Stokes equation: A comparison", Comput. MethodsAppl. Mech. Eng. 82, pp. 323-340, 1990. [11] R. DURAN, M.A. MUSCHmTT~ and R. RODRI¢;UEZ, "'On the asymptotic exactness of error estimators for triangular finite elements", Numer. Math. 59, pp. 107-127, 1991. [12] R.E. EWING, "A-posteriori error estimation", Comput. Methods" Appl. Mech. Eng. 82, pp. 59-72, 1990. [13] J.T. ODEN, "Adaptive finite element methods for problems in solid and fluid mechanics", in: Finite Elements, Theory and Applications, edited by D.L. DwovH<, M.Y. HUSSAINI and R.G. VOIG'I', Springer Verlag, New York, 1988. [14] J.T. ODEN, L. DEMKOWITZ, L. RA('HOWITZ and T.A. WESTERMANN, "Towards a universal h - p adaptive finite element strategy. Part 2: A-posteriori error estimation", Comput. Methods AppL Mech. Eng. 77, pp. 113-180. 1989. [15] J.T. ODEN, L. DEMKOWITZ, L. R A c n o w n z and T.A. WESTERMANN, "A-posteriori error analysis in finite elements: The residual method for symmetrizable problems with applications to compressible Euler and Navier-Stokes equations", Comput. Methods Appl. Mech. Eng. 82, pp. 183-204, 1990. [16] L. PLANK, E. STE1N and D. BISCHOFF, "Accuracy and adaptivity in the numerical analysis of thin-walled structures", Comput. Methods Appl. Mech. Eng. 82, pp. 223-256, 1990. [17] M.S. SHEPHARD, Q. NIu and P.L. BAEHMANN, "Some results using projectors for error indication and estimation", in: Adaptice Methods in Partial Differential Equations, edited by J.E. FLAHERT',', P.J. PASt,OW, M.S. SHEPrtARD and J.D. VASILAKIS,Philadelphia, PA, pp. 83-99, 1989. [18] T. STROUBOUUS and K.A. HAOUE, "Recent experiences with error estimation and adaptivity. Part I. Review o1 error estimators for scalar elliptic equations", Comput. Methods Appl. Mech. Eng., to be published, 1992. [19] T. STROUBOULIS and K.A. HAOUE, "Recent experiences with error estimation and adaptivity. Part II. Error estimation for h-adaptive approximations on grids of triangles and quadrilaterals", Comput. Methods Appl. Mech. Eng., to be published, 1992. [20] R. VERFORTrt, "A-posteriori error estimators for the Stokes equation", Numer. Math. 55, pp. 309-325, 1989. [21] L.B. WArtLBIN, "Local behaviour in finite element methods", in: Handbook of Numerical Ana{vsis, edited by P.G, CIARt.ET and J.L. LIONS, North-Holland, Amsterdam, Vol. 2, pp. 353 522, 1991. [22] J.Z. ZHU and O.C. ZIENKIEWlCZ, "Adaptive techniques in the finite element method", Commun. AppL Numer. Methods 4, pp. 197-204, 1988. [23] O.C. ZXENK1EWlC'Z and J.Z. ZHtJ, " A simple error estimator and adaptive procedure for practical engineering analysis", Int. J. Numer. Methods" Eng. 24, pp. 337-357, 1987. [24] O.C. ZIENKIEWlCZ and J.Z. ZHU, "The three R's of engineering analysis and error estimation and adaptivity", Comput. Methods Appl. Mech. Eng. 82, pp. 95-114, 1990.