A posteriori estimates for finite element approximations

Chapter 4

A posteriori estimates for finite element approximations

Finite element method is one of the main numerical tools in modern numerical analysis. From the mathematical point of view, this method has its origin in the works of W. Ritz [187] and R. Courant [57]. A detailed exposition of its mathematical basis and applications can be found in numerous books (see, e.g., P. Ciarlet [53], G. Streng and G. Fix [204]). In what follows, we assume that readers are familiar with foundations of the finite element method and focus our attention on a posteriori error estimates. This subject has been intensively investigated in the last decades and generated a vast amount of publications. The goal of this chapter is to give a concise overview of several methods developed for controlling errors in terms of energy norms and other quantities. We restrict ourselves with consideration of only one class of linear elliptic problems. On this paradigm, we try to emphasize principal mathematical ideas that form the basis of each method. References to publications that contain a detailed exposition of various approaches to a posteriori error estimation for FEM are given at the end of the chapter. 39

40

4.1

C H A P T E R 4. A POSTERIORI ESTIMATES FOR FEM

Methods ative

4.1.1

based

norm

upon

of the

Approximation

estimates

equation

errors

and

of the neg-

residual

residuals

We begin by recalling basic relations between residuals and errors that hold for systems of linear simultaneous equations. Let Jt E M ~• det Jt ~ 0, and f be a given n-vector. Then, there exists a vector u E I~n satisfying the system Au + f = O. (4.1.1) Assume that v is an approximation of u. We have A ( v - u) - A v + /

and, therefore,

IIv - ull ~ iI,A-x IIIL,Av+ fll,

(4.1.2)

This is the simplest residual type estimate that provides an upper bound of the error e : - v - u in terms of the norm of the residual

r'-Av+ Upper and lower estimates of obtained if the quantities ,~min

--

min

f.

Ilell evaluated

II,Ayll

u~ R'~ IlYll yr

and

in the norm of r can easily be

)kmax

--

max

IIAyll

u~R" Ilyl[ ur

are known. Since jte = r, we see that

IIAell /~min ~

Ilrll

llell -- Ilell ~

)~max

(4.1.3)

and, therefore, --1

-1

Amaxllr[I ~ Ilell <__Aminl[rl["

(4.1.4)

Since u is a solution, we have

ilA~ll iluJl

il/ll ~ ll~JJ

/~max

(4.1.5)

,Xmaxllfll ~ Ilull ~ )~minllfll"

(4.1.6)

,~min ~_

--

and -1

--1

4.1. RESIDUAL BASED ESTIMATES

41

The estimates (4.1.4) and (4.1.6) imply the basic error relation: ,~min /~max

Ilrll < llv- ull < ,~max Ilrll [lYll Ilull - )~min Ilfll"

(4.1.7)

The quantity ~)kmin is called the condition number of the matrix ,,4 and is denoted by Cond.A. Now, we represent (4.1.7) in the form (Cond ,4)- 1 ~

___ IIv - ull < Cond A !!rl[.

Ilull

-

(4.~.s)

II/ll

The estimate (4.1.8) reads as follows: the relative value of the error is bounded from above and below by the value of the residual normalized with respect to the norm of f. The factors depend on the condition number of the matrix ,,4. Some generalizations of these estimates can be found in [10]. In principle, the above consideration can also be applied to a wider set of linear problems, where ,4 E Z:(X, Y) is a coercive linear operator acting from a Banach space X to another space Y and f is a given element of Y. However, if ,4 is related to a boundary-value problem, then one should properly define the spaces X and Y and find a practically meaningful analog of the estimate (4.1.8).

4.1.2

R e s i d u a l t y p e error e s t i m a t e s for a p p r o x i m a t e l u t i o n s of b o u n d a r y - v a l u e p r o b l e m s

so-

Let .,4: X ~ Y be a linear elliptic operator. Consider the boundary-value problem ~tu+f=0 inf,, u=u0 on0f~. (4.1.9) Assume that v E X is an approximation of u. Then, we should measure the error in X and the residual in Y, so that the principal form of the a estimate is IIv - ullx < CII.Av + / l l Y , (4.1.10) where the constant C is independent of v. The key question is as follows:

Which spaces X and Y should we choose for a particular boundary-value problem ? Consider the classical problem -Au = f u = 0

in fl, on 0f~,

f e L2(f~),

(4.1.11) (4.1.12)

42

CHAPTER 4. A POSTERIORI ESTIMATES FOR FEM

whose generalized solution satisfies the relation

f V u . Vw dx - fn f w dx

Vw E Vo := ~1 (12)

(4.1.13)

and is subject to the energy estimate (4.1.14)

I]Vu[12,~ ~_ Cn]lfll2,~,

where Cf~ is a constant in the Friederichs inequality. Assume that the approximation v E V0 possesses extra regularity properties, so that Av E L2(12). Then, from (4.1.13) it follows that

/

V ( u - v) . Vw dx - / f ~ ( f + Av)w dx,

Vw E Vo.

(4.1.15)

Setting here w = u - v, we obtain the estimate

IIV(u- v)ll2,~ ~ C~][f + AvI[2,~.

(4.~.16)

The right-hand side of (4.1.16) is formed by the L2-norm of the residual of (4.1.11). However, this estimate may be useless for approximations constructed by conforming approximations of the energy space Vo. Indeed, the general convergence theory can guarantee that a sequence {Vk } of such approximations converges to an exact solution in H 1, but no convergence of higher derivatives could be observed. Therefore, the quantity IIAvk + fl] does not necessarily converges to zero, unlike the left-hand side of (4.1.16), which decreases and tends to zero. This means that the consistency (the key property of any practically meaningful estimate) is lost. Now, the following question arises: Which norm of the residual leads to a consistent estimate of the error in the energy norm? To answer it, we replace (4.1.13) by the integral identity n V u . Vw dx = ( f , w),

V w E I/o,

(4.1.17)

where (f, w) denotes the value of f E H -I on w E V0. If f E L2(12), then the right-hand side of (4.1.17) is expressed by the integral and we arrive at (4.1.13). However, as we have seen, the set of linear continuous functiona]s defined on V0 is wider than the set of the "regular" functionals presented as Lebesgue integrals of summable functions. It also includes the functionals from H -I. Since < I f IIIVwll2,n, (4.1.18)

4.1.

43

RESIDUAL BASED ESTIMATES

we set w = u in (4.1.17) and find that

IlVull2,~ ~ Ifl.

(4.1.19)

This estimate gives an upper bound of the energy norm of u in terms of the norm of f in H -x and suggests an answer to the question stated. Assume that v E V0 is an approximation of u. By (4.1.13), we have ff~ V(u - v) . V w dx - s ( f w - V v " V w ) d x .

(4.1.20)

The right-hand side of (4.1.20) is a linear continuous functional defined on o

the space H 1(f~), so that it belongs to H-l(f~). Thus, (4.1.20) has the form a V ( u - v). Vw dx - (Av + f, w},

(4.1.21)

where f + A v e H -1 (f~). Setting w = v - u, we obtain the estimate

LIV(u - v)112,n ~ I f + my I.

(4.1.22)

Note that If+Avl=

sup

I(f+mv'~}[

sup

Ira v ( ~ - , ) .

o

V~l

IIv~il~,a

and we arrive at the relation ]IV(u - v)]]2,f~ -- I f + Av l,

(4.1.23)

which shows that the quantity I f + A v k l tends to zero for any sequence {Vk} that tends to u in H 1. Thus, the estimate (4.1.23) is consistent. From (4.1.23) it follows the estimate

IlY + ~ll(-~),a < IIV(~- ~)ll~,a < (c~ + 1)~/~11/+ ~vllc-~),a, which yields error bounds in terms of another negative norm. Similar estimates can be derived for any uniformly elliptic equation of the divergent type. Indeed, let

d o( i,j=l

ou)

C H A P T E R 4. A P O S T E R I O R I E S T I M A T E S FOR FEM

44

where aij(x) -- aji(x) are measurable bounded functions satisfying the condition Aminlr/I 2 < aij(x)rli~Tj < Am~xlr/I2,

V~ E I~n, x E ~2

(4.1.24)

with positive constants )~min and )~max independent of x and ~. Henceforth, we will also use a the compact notation .Au = div AVu,. Assume t h a t f E L2(fl) and consider the boundary-value problem

Au + f - O,

in ~2,

u - 0

on 0~,

whose generalized solution u E V0 is defined by the integral identity

/

A V u . V w dx - / ~

f w dx,

Vwe Vo.

(4.1.25)

Let v E V0 be an approximation of u. Then,

/

A V ( u - v) . V w d x - f ~ ( f w - A V v . V w ) d x ,

Vw e Yo.

The right-hand side of this relation is a bounded linear functional on Vo. Therefore, it determines a distribution, which is f + div(AVv) e H -1. Hence, we have the relation

~ A V ( u - v) . V w dx - (f + div(AVv), w),

Vw e Vo.

Setting w - u - v, we derive the estimate

-1 llV( It -- V)ll2,gt --< )~min I f "+- div(AVv)].

(4.1.26)

Next,

,~l- +

j~., ~ : ~v ., ~ u ~ '

=

sup

-

sup

I (f + d i v ( A V u ) , ~) i

< )~maxllV(u - v)l12,~.

(4.~.27)

45

4.2. EVALUATION OF NEGATIVE NORMS Combining (4.1.26) and (4.1.27), we obtain

--1 --1 '~max I n(v) l ~ [[V(~ - v)ll2,a < "~min IR(v)l,

(4.1.28)

where

R(v) = f + div(AVv) e H-l(f~) is a residual understood in the sense of distributions. The estimate (4.1.28) shows that upper and lower bounds of the error can be evaluated in terms of the negative norm of R(v). It presents a certain analog of the estimate (4.1.8) for approximate solutions of elliptic type boundary-value problems. There is a simple case in which the estimate of IR(v) l is easy to find. Assume that R(v) c LP(f~) with p >_ 1. Then

In(v) l _< sup If the space Lq(12) with q - ~ the estimate

[]R(v)Ilp,a II~olIq,a

(4.1.29)

is embedded in H1 (f~), then (4.1.29) implies

IR(v)l < Cq(a)llR(v)llp,a, where Cq(f~) is the constant in the inequality o

II~llq,a <_ c~(a)llV~ll~,a,

v~ e H i ( a ) .

However, estimates of such a type have only a conventional meaning, because in the majority of cases approximate solutions obtained by finite element approximations cannot guarantee such a regularity of the residual.

4.2

E v a l u a t i o n of n e g a t i v e n o r m s of residuals

To use the estimate (4.1.28) one needs either to find IR(v) l or to obtain a sharp estimate of this quantity. A wealth of papers devoted to the residual method is, in fact, focused on the way of overcoming the difficulties arising in the process. For continuous and piecewise smooth approximations, there are two main approaches to finding such estimates. One group encompasses the methods estimating I R(v) I directly by using Galerkin orthogonality property and interpolation estimates for finite element patches. They provide explicit error bounds represented by sums of residuals on the elements and

C H A P T E R 4. A P O S T E R I O R I E S T I M A T E S F O R F E M

46

jumps over interelement boundaries. These sums have weight coefficients that depend not only on the entire domain, but also on the shape of subdomains. In the other group of methods, the quantity JR(v) [ is estimated by solving auxiliary boundary-value problems whose right-hand side is formed by the residual. Below, we demonstrate ideas of both approaches with a paradigm of the problem (4.1.25). 4.2.1

Estimation

of the

residual

in a one-dimensional

problem To clarify the idea of the first method, we consider a simple one-dimensional problem: find u such that (au')' + f = O,

(4.2.1)

u(O) = u(X) = O,

(4.2.2)

where a - a(x) > 0 and f = f ( x ) are given functions. Let I denote the interval (0, 1), f E L2(I), and a(x) E C(I). A weak statement of the problem (4.2.1)-(4.2.1) is as follows: find u E Vo(I) such that au'w' dx -

f w ' dx,

Vw e Vo(I) := H 1(I).

(4.2.3)

Divide I into a number of subintervals Ii = (xi, xi+l), where xo - O, xg+l = 1, and [xi+l - x i [ = hi. Assume that we have an aproximate solution

O

v E H 1 (I) that is smooth on any interval Ii. In this case, In(v)] =

-

sup weVo(Z),

sup

f~ ( - a v ' w ' + f w ) d x =

IIw'll2,z EN=0 fli ri(v)w dx + EN~=Ia(x,)w(xi)j(v'(x~))

(4.2.4)

IIw'll ,z

where j(r

:= r

+ 0) - r

is the "jump-function" (note that j(r ous at x0) and =

- 0) - 0 if the function r is continu-

+ f

is the residual on Ii. The supremum in (4.2.4) is taken over all functions from Vo(I) other than zero function. To simplify further notation, we henseforth

4.2. E V A L U A T I O N OF N E G A T I V E N O R M S

47

will not write this explicitely, but always assume that zero functions are eliminated in quotient type relations. To find an upper bound of the right-hand side of (4.2.4), we restrict the class of admissible approximations and hereafter assume that v is the Galerkin approximation obtained on a finite-dimensional subspace V0h formed by piecewise affine continuous functions. Let v is the Galerkin approximation Uh E Voh, which obeys the relation O~UhW h

dx =

f wh dx

VWh E VOh.

(4.2.5)

Then,

I n(uh)l =

sup ~vo(z)

fd (--au~h( w -- 7rhW)' + f ( w -- 7rhW)) dx Iiw'112,I

where 7rh : Yo --4 Yoh is the interpolation operator defined by the conditions rhV E Voh, 7rhV(O) = rrhV(1) = 0 and

rrhV(Xi) = v(xi),

Vxi, i = 1, 2,..., N.

(4.2.6)

Integrating by parts, we obtain

IR(uh) l =

sup ~Vo(I)

{ EN=~fx, r~(uh)(w - ~rhw)dx + IIw'll2,z 4- EN=I a(xi)(w(xi) - rhW(Xi))j(u~h(Xi)) flW'lr J

By (4.2.6), w(xi) - rhW(Xi) = 0 SO that the second sum vanishes. The first one admits the estimate N

i=O

N

i

i=0

It can be represented in a more convenient form, provided that the estimates

[lw- ~h~ll2,z, ~ c~llw'll2,z,

(4.2.7)

hold with constants ci dependent on I~. Then, we obtain the estimate N

i=O

N

i

i=0

<-

E i=0

ci llri(uh)]]2,1`

)lw'll2,z,

(4.2.8)

48

C H A P T E R 4. A P O S T E R I O R I E S T I M A T E S F O R F E M

which implies the desired upper bound N

)~/e 2

[R(uh)[ <

2

(4.2.9)

Ecillri(uh)[12,ii i=0

This bound is the sum of local residuals with weights given by the interpolation constants ci. In this simple problem, the constants ci can easily be found. Indeed, let 7 / b e a constant that satisfies the condition

I1~'i1~,,, IIw-~hwll 22,Ii

inf ~-~

_

(4.2 10)

> 7~.

Then, - t / 2 w' IIw- ~hwll~,z, <_ ri II 112,I, and one can set ci

Xi+I

L

i

=

")'I

-x/2 Note that 9

xi+l

Iw'l 2 dx

-

I(w - 7rhW)' + (Trhw)'l 2 dx,

,/xi

where (Trhw)' is constant on (xi, xi+l). In view of (4.2.6),

L

zi+' (w - 7rhW)'(lrhw)' dx - 0 i

and

Lxi+l

Xi+l

Iw'l 2 dx >

[(w

cXi

--

7rhW)'[2 dx.

i

Thus, we have

inf weHl(i,)

L~i{ 1 ]wlI2 d X fx,+l

JT, i

[W -- 7rhW[ 2

>

inf

>

dx

-- w E H I ( I i )

-

fx,+l i( w _ rhW),12 dx ~x~ > [ x i + l [W 7rhW[ 2 dx Jxi

f~'+~ [r/'[2 dx > -

so t h a t 7 i -

7r2/h2 and c i - hi/Tr.

inf

,,z/ f~'+~

l~TI2 dx

71"2

h~'

4.2. E V A L U A T I O N OF N E G A T I V E N O R M S

4.2.2

Estimation lems

49

o f t h e r e s i d u a l for l i n e a r e l l i p t i c p r o b -

Consider the problem (4.1.25). Let f~ be represented as a union of elements Ti. For the sake of simplicity, we assume that f~ is a polygonal domain and f~ = t_Ji=1Ti. Assume that Yoh consists of piecewise affine continuous functions. The respective Galerkin approximation Uh satisfies the relation - -

N

- -

~ A V u h . Vwh dx - f f f w h dx,

VWh E Voh,

(4.2.11)

where

Voh -- {Wh C Vo l wh C p l ( T i ) , Ti E ~h}. We have

I R(uh)l -

sup ~Vo

f" (fw -

A V u h . Vw) dx = IlVwll2,~

= sup

f ~ ( f ( w - rhW) -- A V u h . V(w - rhW)) dx

~Vo

IlVwll2,~

where 7rh is the Clement's interpolation operator (cf. (2.3.18)-(2.3.19)). Let Oh -- A V u h . Then, for any simplex Ti

f T ah . V ( w - 7rhw) dx -

-- fOTi (ah'v)(W -- rhw)ds -- fTi divah(w -- Trhw) dx' where v is the unit outward normal to OTi. Since on the boundary W--rhW -0, we obtain

I R(uh) I -- sup ~Vo

{ E~=~ fri (div ah + f)(w - ~rhw)dz ILVwlI2,~

+

+ E~=xEj~, fE,~llVwll2,~J(~h'V'-J)(w~hW)ds} . In the above relation, j(ah'Vij) denotes the jumps of ah'Vij on Eij and vii is the unit normal vector to Eij outward to Ti if i < j. Note that

T (div orh + f ) ( w - 7rhw)dx _< [1div (7h +

fll2,T,IIw-

<_ clill div ah + fll2,T~ diam(Ti)llwlll,2,~N(T~),

7rhWil2,Ti

CHAPTER 4. A POSTERIORI ESTIMATES FOR FEM

50

where I.gN(Ti) is the patch of elements having common nodes with 7'/ (see (2.3.20)-(2.3.21)). Then, we estimate the first sum as follows (cf. (2.1.2): N

~I /T~(divah + f)(w-Trhw)dx < dl

c~i diam(T~)2ll div Oh + Y112,T,

Z

Ilwil~,2,. ,

(4.2.12)

i--1

where the constant dl depends on the maximal number of elements in the set wN(Ti). For the second one, we have N

i=1

N

j>i N

ij N

< ~ ~ IIJ(ah'V~IlI2,E,~C2'~IE'jlX/211wlIx,2,~(T,)< i=1

j>i

< d2 Z~C~jlE~IIIJ((rh'V~j)II~'E'J i=1

IIwlIx,2,a"

(4.2.13)

j>i

In the above relations, the constant d~. depends on the maximal number of elements in the set wE(T~). The relations (4.2.12) and (4.2.13) imply the estimate

I R(uh) I < Co

c2i diam(Ti)211 div ah + II[2,T,

~

+

i=1

(4.2.14) i=1

j>i

Here Co - max(dl,d2)v/1 + C~, where C~ comes from (2.2.1). The righthand side of (4.2.14) is the sum of locally defined quantities r/(T~) multiplied by constants depending on properties of the chosen splitting ~'h (e.g., on the maximal number of elements in a neighborhood of a node, the minimal angle of the triangulation, etc.). For quasi-uniform meshes all generic constants cli have approximately the same value and can be replaced by a single constant cl. If the constants c2~j are also estimated by a single constant c2, then (4.2.14) implies an estimate 1/2

'R(uh)' <-C~ ( ~

'

(4.2.15)

4.2. E V A L U A T I O N OF N E G A T I V E N O R M S

51

where r/2(T~) = c2 diam(Ti)2l[ div ah + fli~,T~ + -2c~ ~

iEijlllj(ah, vij)ll 22 , E ~ y

"

EijET~

(4.2.16) Two parts of (4.2.16) are the residual norm on Ti and integrals related to the jumps on the edges. The multiplier 1/2 arises, because any interior edge is common for two elements. It is worth noting that for a strongly nonhomogeneous mesh such an estimate may lead to a considerable overestimation o f [ R ( u h ) [ . The sharper estimate (4.2.14) involves the constants cli and c2ij whose evaluation requires the solution of the following variational problems:

inf IIwlI ,2, (T,)diam(Ti)

(4.2.17)

inf

(4.2.18)

and

eWo IIw- hwll2, ,, 4.2.3

Estimation

of the residual

norm

by implicit

meth-

ods For a conforming approximation Uh, we have the basic relation

a A V ( u - Uh) " V w dx - (R(uh), w),

Vw e Vo.

(4.2.19)

Hence, the error e - u - Uh is a solution of the problem div A V e + R(uh) = 0 e= 0

in fl,

(4.2.20)

on 0f~,

(4.2.21)

where R(uh) = f + d i v A V u h E H -1 (fl) and the first equation is understood in the sense of distributions. Formally, this idea can be extended to a wide class of linear problems. Indeed, if ,4 : X --+ Y is a linear operator for which we consider the problem

.Au + f = 0 u = 0

in f~,

(4.2.22)

on 0[2,

(4.2.23)

in ~,

(4.2.24)

on OFt,

(4.2.25)

then, for any v E X,

Ae + R(v) - 0 e= 0

CHAPTER 4. A POSTERIORI ESTIMATES FOR FEM

52

where R(v) = f + .Av. However, the evaluation of the error by directly using (4.2.24)-(4.2.25) may be a difficult task, because, in general, R(v) is a distribution. Moreover, the accuracy of an approximate solution of (4.2.24)-(4.2.25) obviously affects the accuracy of the error estimation, so that a new error estimation problem arises. In the literature focused on practical applications of this method, it is often suggested to split the functional (R(uh), w) into a number of functionals defined as solutions of local subproblems. A version of such an approach is described below. Let ~ be a union of nonoverlapping domains ~i, i = 1, 2, ...N. Assume that for each Qi we have found a function ui such that N


(4.2.26)

i--1 f~i

Since A is a positive definite symmetric matrix, we supply L2(12i, II~'~) with a new norm

Ay(x). y(x) dx

[[[ Y II!~,'-

,

i

which is equivalent to the standard one. inequality, we derive the estimate

Then, in view of the Cauchy

N

I t- Z

ILlv ~, ili~, ni v ~ iti,, <

i--1

<111v ~ --

Ilia

Ili vu~ III2f ~ i

<- - c Ilvwll2

i:1

tll Vu~ III2

,f~

i:1

(4.2.27) where c -

V/Amax. Hence,

I R(~h)l - sup I
~/

I< ~

IIIW~ II1~,

(4.2.2S)

Now, the main task is to find the functions ui. Let us define them as solutions of local subproblems with Neumann type boundary conditions on patches of elements 12i with a common node i. Denote the common edge of

4.2. EVALUATION OF NEGATIVE NORMS

53

fti and f~j by Fij and F/o := Ofti N Oft. For each fti we solve the problem: find ui E Hl(fti), such that ui - 0 on Fi0 and

AVui . Vwdx = f fwdx+ N

AVuh. Vwdx,

Vw E Vo(fti).

(4.2.29)

In this relation, g is a function defined on all Fij and it is assumed that the integral over Fij is equal to zero if Fij - 0. The space Vo(f~i) is defined as follows. If F 0 / ~ 0 then

Vo(f~i) - {w e Hl(f~i) I w - 0 on Fi0, }. If F0i = 0, then Vo(fti) is the space H 1(f~i). The integer-valued function aji is defined by the rule / - / +1' i f / > j , aj - 1 , if j > / ,

( O,

ifi-j.

Usually, the function g is constructed on the basis of the information encompassed in Uh. For example, if the domains Fti coincide with the elements of a certain finite element sampling, then g can be determined by computing fluxes of a respective finite element solution. Each internal boundary F ij generates two integrals with equal absolute values and opposite signs. Therefore, the sum of all integrals does not contain such terms, and we obtain N

f AVui" Vwdx i = l l"tl

: / f w d x - / AVuh" VWdx = (R(uh), W), f~

Vw E V0(f~). (4.2.30)

f~

Hence, (4.2.28) provides an upper bound of the error in terms of the norms of ui. Consider a function ~ : ft -+ 1I~ that coincides with ui(x) if x E f/i. Assume that the functions ui preserve continuity on the boundaries ['ij and the function ~(x) belongs to Hl(ft). Then (4.2.30) reads

f

f~

AY('5 + Uh) " Vwdx = f fwdx f~

Vw E Vo(f~).

(4.2.31)

54

C H A P T E R 4. A P O S T E R I O R I E S T I M A T E S FOR FEM

The relation (4.2.31) means that u = ui + Uh on f~i. Therefore, ui = u - Uh and by these functions we also obtain local errors. However, this formally simple procedure contain certain technical difficulties. One of them is that for internal domains the function g (which defines the Neumann type boundary conditions of the local subproblems) cannot be taken arbitrarily. This follows from the fact that the Neumann problem may be unsolvable if the external data do not satisfy an additional condition. For the problem (4.2.29) this condition is as follows:

N

/ .fdx + J~Xr~igds-O.aj f~i

=

(4.2.32)

""

Therefore, it is required a special procedure (often called the "equilibration procedure") that transforms g in order to satisfy (4.2.32) on each element. After that, exact solutions of local problems must be found and then used to compute the upper bound of the error. In practice, the accuracy of such an estimate depends on the choice of g and on the accuracy of the computed approximations of ui.

4.3

Error estimation methods using adjoint problems

Global error estimates give a general idea on the quality of an approximate solution and stopping criteria. However, from the viewpoint of engineering purposes, this information is not sufficient. In many cases, it is important to have special (problem-oriented) measures of the error associated with certain subdomains, lines, and points. Characteristics of such a type can be constructed by means of specially selected linear functionals ~s, s = 1, 2, ...M. Having computable estimates of the quantities (es, u - Uh), one could judge the quality of a finite element approximation Uh in reference to the desired property. Moreover, such estimates also lead to integral type estimates of the error. A sensible example of such a quantity is given by the functional

(e,

-

- f

(,

u) dx,

where the weight function ~0 is positive in w C C ~ and vanishes in f~ \ w. Such a quantity gives information on the behavior of v - u in w. One can also construct other quantities characterizing v - u along certain lines

4.3.

METHODS

USING

ADJOINT

PROBLEMS

55

and also at some points of 12. In other words, by choosing ~Oo, one can create a wide spectrum of "quantities of interest" that complement the information obtained by global error estimates. This direction is often called g o a l - o r i e n t e d error estimation methods (see, e.g., [3, 152]), because g is chosen to fit a special goal. Certainly, an estimate of the desired quantity is easy to obtain with the help of the obvious relation I{e,

-

h)l < Ilellll

-

hllv,

provided that the global error I I v - u l l v and the norm I[gll are estimated. However, in many cases, such an estimate will strongly overestimate the quantity of interest. A posteriori estimates of the errors evaluated in terms of linear functionals are derived by attracting the a d j o i n t boundary-value problem whose right-hand side is formed by the functional g. Let us represent this idea in the simplest form. Consider a system of linear simultaneous equations A?~ -- f~

where A is a positive definite matrix and f is a given vector. Let v be an approximate solution. We need to estimate the quantity g. (u - v), where g is a given vector. Define u t by the relation A*ut

= g

in which A* is the matrix adjoint to A. Then, g . (u - v) = A*ut

. u-

g . v = f . ut - g . v = (f - Av) . ut

Hence, an estimate of g. (u - v) can be obtained fairly easily provided that ut is defined. Certainly, this quantity cannot completely characterize the error (e.g., g. ( u - v) vanishes if the residual R ( v ) = f - A v is orthogonal to u t ) . Therefore, it is desirable to obtain estimates for various functionals gs, which amounts to solving several adjoint problems. In this section, we briefly describe a version of the approach to the a posteriori error estimation for finite element approximations based on the consideration of an adjoint problem. As before, we show this with a the paradigm of the problem (4.1.25). The reader interested in a deep investigation in this subject can find a thorough exposition of these methods in the book [26] (see also [27] and other publications given in the literature comments at the end of this chapter).

CHAPTER 4. A POSTERIORI ESTIMATES FOR FEM

56

Consider again the problem (4.1.25). Now, we do not assume that A is a symmetric matrix and denote its adjoint matrix by A*. Also, it is convenient to denote the solution of (4.1.25) by uf, i.e

/AVus.Vwdz-ffwdx, The

adjoint

Vwe Vo(f).

problem is to find ue E V0(f~) such that

f A*Vul.Vwdx- f twdx, Vwe Vo(n). Let f~ be divided into a number of elements Ti, i = 1,2,...N. Given approximations on the elements, we define a finite-dimensional subspace Voh C Vo(f~) and the Galerkin approximations Ufh and Uth as functions satisfying the relations

f AVufh.VWhdX-/fwhdx,

VWhE Voh,

/ A*Vulh.Vwh dx-- /f ewhdx,

~CwhE Voh.

Since

f e(uf - Ufh)dx -- fn A*Vut.

V(uf -

Ufh)dx

and

f A*Vueh"

V(uf -

Ufh)dz -- O,

we arrive at the relation (4.3.1) whose right-hand side is expressed in the form N

f AV(uf - Ufh)" V(ut

-

Ulh)dx -

i--1 Ti

div (AV(uf i--1

Ti

1/

+~

OTi

j(vi"

Ufh)) (ul -- Ulh) dx+

}

AV(uf - Ufh)) (ue -- Ulh) ds .

4.3. M E T H O D S USING A D J O I N T PROBLEMS

57

This relation implies the estimate N

I,.,l(u$ - UIh)dx-~{lldiv

AV(ul - UIh)II2.T, IlU,- UZhlI2,T, +

i=1

+89 IIJ(v," AV(uS - ~'Zh))II~,oT, Iluz -- UZ~II=,oT, } -N

-- ~

{ IIf + div AVuzhlI=,T, I1~'~ -- uzhlI=,T, +

i=1

+~ IIJ(vi "AV(ulh))II2,0T~ llul -- U~hll2,0T~

}.

(4.3.2/

We see that the principal terms of this estimate are the same as in (4.2.16), but the weights are given by the norms of ul - Uih. Their values depend on the properties of the adjoint problem. Assume that the latter has a sufficiently regular solution, and Ulh is constructed by piecewise affine continuous approximations. Then the norms llul- UlhllT~ and l l u l - U~hll2,0T~ are estimated by the quantities halull2,2,T~ with c~ - 1 and 1/2 and multipliers ci and cij, respectively. In this case, we obtain an estimate with constants defined by the standard H 2 -+ Voh interpolation operator whose evaluation is easier than that of the constants arising in the H I -+ Voh interpolation procedures. Moreover, if for any e E L2(~) the solution ul lies in H 2, then it is easy to derive an estimate of llu$ - UIhll2,~. Indeed,

II~s - ~'shll~ n '

-

sup

~L~-

II~ll~,~

By (4.3.2), we estimate the numerator of this fraction as a sum of element residuals and interelement jumps with multipliers ci and ~ij and the global multiplier Ilu~112,2,~. The latter norm is subject to the regularity estimate

Ilu~ll~,~,,~ < ~ I1~11~,~, which can also be viewed as a stability relation in the adjoint problem. Then, Ilgl12,~ occurs in the numerator and the fraction is reduced to an estimate of the form similar to (4.2.14), which, however, has other interpolation constants and includes a new "stability" constant. Estimates of such a type for various problems are analyzed in [112, 113, 114, 68] and other papers cited therein. In [3, 152] and some other publications, it was suggested another way of evaluating errors in terms of linear functionals. It is based on the relation (4.3.1) and the algebraic identity

58

CHAPTER

4.

A POSTERIORI

ESTIMATES

FOR FEM

which is valid for any symmetric bilinear form b and ~ E IR. By this identity, the right-hand side of (4.3.1) is expressed in the form 14f oA*V(r et + ~ ) " V(r

+

where et = u t - Uth and e I = u I - U f h . Further, these two integrals are estimated by the equilibrated residual method based upon solving local problems with equilibrated Neumann type boundary conditions (cf. 2.3.2). In Section 5, we present a new estimation method for goal-oriented functionals that uses the post-processing of approximate solutions ut and u I. Further, in Chapter 6, these functionals will also be estimated by functional type a posteriori error bounds.

4.4

M e t h o d s based upon post-processing of finite element approximations

Finite element solutions are usually smooth at the interior points of elements, but in the entire domain they have a restricted regularity. For example, the solutions obtained by piecewise affine continuous approximations have piecewise constant gradients, i.e., they belong to L~176 In the vast majority of cases, this class of functions is much wider than the one that contains gradients of the exact solution. This fact is easy to demonstrate by considering the problem (4.1.25). Indeed, if / E L2(~), then the computed vector-valued function ah -- A V u h lies in L~(~), but a = A V u lies in H(~, div). Moreover, if u e H k ( ~ ) , k _ 2, and the coefficients of A are smooth, then a E H k-1. In more general terms, the situation that often arises in finite element methods is as follows. For a conforming approximation Uh E Vh, the function T u h (where T is a certain linear operator) lies in the space U. However, a priori estimates of the solution guarantee that TuEUCU.

Assume that we have a computationally inexpensive continuous mapping G such that G(Tvh) E U,

VVh E Vh.

Then, we may hope that the function G(Tvh) that possesses a priori regularity properties would be much closer to Tu than Wvh. These arguments form

4.4. POST-PROCESSING

59

TUh~

U

Wu

J

GTUh

Figure 4.4.1: Post-processing operator and error estimation. the basis of various post-processing algorithms that change a computed solution in accordance with some a priori knowledge of properties of the exact solution. Such algorithms are widely used in advanced numerical methods. If the error caused by violations of a priori regularity properties is dominant and the post-processing operator G is properly constructed, then IIGTuh - Tul] << [ITuh - T u l l .

Schematically, such a situation is illustrated in Fig. 4.4.1. In this case, the explicitly computable norm [ I G T u h - Tuh[I can be used to evaluate upper and lower bounds of the error. Indeed, assume that there is a positive number a < 1 such that for the mapping T the estimate [ I G T u h -- T u l l ~ c~ []Tuh - Tu[]

(4.4.1)

holds. Then, we have (1 -

c~)IlTuh - Tull _< I I T u h

--

Tu]l-

IlGTuh - Tu]l

< [IGTuh -- Tuhll ~ IlGTuh - Tull + IlTuh - Tull <_ < (1 + ~ ) I I T ~ h - Tull.

(4.4.2)

Thus, we observe that if the image of Tuh computed by a certain postprocessing operator G is much closer to the exact solution than Tuh (i.e., if a << 1), then IlTuh - Tull ~ IlGTuh -- Tuhll.

(4.4.3)

In this case, the right-hand side of (4.4.3) can be used as an error indicator. The quality of this indicator is the better the smaller the value of a.

CHAPTER 4. A POSTERIORI ESTIMATES FOR FEM

60

In practice, various types of post-processing operators are used. Below, we consider some of them for the problem (4.1.25), where T is either the gradient operator or the operator A V.

4.4.1

Post-processing

by averaging

Post-processing operators are often constructed by averaging Wuh on finite element patches or on the entire domain f~. I n t e g r a l a v e r a g i n g on p a t c h e s If Tuh is a square summable function, then post-processing operators are obtained by various averaging procedures. Let f~i be a patch of finite elements, i.e.,

-fli - U Tij,

j - 1, 2, ...Mi.

Let pk (fli, ~ ) be a subspace of U that consists of vector-valued polynomial functions of degrees less than or equal to k. Define gi C Pk(f~i,I~'*) as a solution of the following minimization problem: inf

f]g-

gEP~(Di,Rn) J fli

Tuh] 2 dx.

(4.4.4)

The minimizer gi can be used to define the values of an averaged function at some points (nodes). Further, these values are utilized by a prolongation procedure that defines an averaged function

GTuh " fl --+ ~. Consider the simplest case. Let T be the operator V and Uh be a piecewise affine continuous function. Then,

VUh E pO (fli, i~n)

on Tij.

We denote the respective values of Vuh by (Vuh)ij. Set k - 0 and find

gi E po such that / ]gi - VUh]2 dx - gEP~

f l g - Vuhl2 dx -

f~i

fli

Mi =

inf

Mi

]gl2lfh]- 2g . Z ( V u h ) i j l T i j l + ~

j=l

j=l

](Vuh)q[2lTq]

} .

(4.4.5/

4.4. POST-PROCESSING

61

It is easy to see that gi is given by a weighted sum of ( V U h ) i j , namely, Mi Iq'l l

gi - E

"~(Vuh)ij.

(4.4.6)

j--1

Now, we define the value of G(Vuh)(Xi) as gi. Repeat this procedure for all nodes and define the vector-valued function GV(uh) by the piecewise affine prolongation of these values. If the mesh is regular and all the quantities ITij] are equal, then (4.4.6) reads Mi

1

(4.4.7)

=

j--1

Various averaging formulas of this type are represented in the form Mi

Mi

g i - Z Aij(Vuh)ij, j=l

Z

Aij - 1,

(4.4.8)

j=l

where the quantities Aij are the weight factors. For internal nodes, they may be taken in accordance with (4.4.6) or defined by the rule

2r ' where [~ijI is the radian measure of the angle of Tij associated with the node i. However, if a node belongs to the boundary, then it is better to choose special weights. Their values depend on the mesh and on the type of the boundary. The reader can find a detailed consideration of this question in [101]. D i s c r e t e a v e r a g i n g on p a t c h e s Consider the problem mi

inf

Z

[g(xs) -

Tuh(xs)[ e

gEP~(~i) s = l

where the points x s are specially selected in f~i. Usually, the points x s are the so-called superconvergent points (see below). Let gi E ~k(f~i) be the minimizer of this problem. Then, we define the values of GTuh at some points and further use an inexpensive prolongation procedure to define GTuh at any point.

62

CHAPTER

4. A P O S T E R I O R I E S T I M A T E S F O R F E M

If k - 0 , then (cf. (4.4.6)) 1

mi

(4.4.9)

gi - - - E T u h ( x s ) . mi s=l

Global averaging In many cases, an efficient averaging operator is obtained if local minimization problems on patches are replaced by a global problem. Assume that T u b is a square summable function and solve the problem: find ~h E Vh (12) C U such that I]gh -- Tuhll2~ --

inf Ilgh -- Tuh[12~ 9 gh ~Y~ (n)

(4.4.10)

The function gh can be viewed as G T u h . Very often gh is a better image of Tu than the functions obtained by local procedures. Moreover, on can await that mathematical justifications of the methods based on advanced averaging procedures can be performed under weaker assumptions what makes them applicable to a wider class of problems (see, e.g., [43]). 4.4.2

Superconvergence

Let Uh be a Galerkin approximation of u computed on Vh. For piecewise affine approximations of the problem (4.1.25), we have the estimate IIV(u - Uh)ll2,n <__clh.

(4.4.11)

Here, cl depends on the H2-norm of the exact solution and on properties of the sampling 7~. Under apropriate regularity assumptions similar estimates also hold for the maximum-norm. Such estimates show the maximal convergence rates that might be generally observed the type of approximations under consideration. However, it was discovered (see, e.g., [155, 216, 229]) that in many cases this rate may be higher. For example, the estimate [u(xs)--Uh(Xs)l <_ C h 2+~

a > 0

(4.4.12)

may hold at a certain point xs, which is called a superconvergent point for the sequence {Uh}. Certainly, the location of superconvergent points (and their existence) strongly depends on the structure of 7~. Analogously, an operator G possesses a superconvergence property in w C ~2 if I I V u - GVuh[12,~ <_ c2h TM,

(4.4.13)

63

4.4. P O S T - P R O C E S S I N G

where the constant c2 depends on higher norms of u and the structure of 7~. For the problem (4.1.25) estimates of such a type can be found, e.g., in [101, 126, 127]. For example, let ~ E n~n, n = 2, be a domain with piecewise smooth boundary, A E C2(~,Mn• u E H3(nl), where w CC f~l CC ~ and w with its neighborhood is covered by congruent triangles. If the operator G is constructed in accordance with (4.4.6), then (4.4.13) holds with a = 1 ([127]). Nowadays it is a matter of common knowledge that convergent sequences of Galerkin approximations often possess certain superconvergence properties (see, e.g., [18, 213, 216, 219]). By exploiting these properties, one can usually construct a simple post-processing operator G satisfying (4.4.1). In this estimate, the value of a depends on h, which improves the quality of the error indication as h decreases. For this reason, error indicators of such a type are widely used in engineering computations (see, e.g., [17, 29, 30, 224, 226, 227, 228]). 4.4.3

Post-processing

by equilibration

For a solution of the problem (4.1.25), it is a priori known that div a + f = 0,

(4.4.14)

where a = A V u . This suggests an idea to construct an operator G such that d i v ( G ( A V u h ) ) + f = O.

(4.4.15)

If G possesses additional properties (linearity, boundedness), then we may hope that the function G A V u h is closer to a than A V u h and use the quantity ]]AVuh - GAVuhl] as an error indicator. This idea can be applied to an important class of problems A*Tu + f = 0,

Tu = AAu,

(4.4.16)

where r is a positive definite operator, h is a linear continuous operator, and A* is the adjoint operator. Note that in continuum mechanics, equations of the type (4.4.14) are referred to as the equilibrium equations. Therefore, it is natural to call an operator G (whose range consists of functions satisfying these equations) an equilibration operator. One way for constructing such an operator is as follows. Assume that we know the functions a~, i = 0, 1, 2, ...M, that satisfy the conditions div a0 + f = 0,

divai = 0,

i = 1, 2, .., M.

CHAPTER 4. A POSTERIORI ESTIMATES FOR FEM

64

In this case, G may be defined as the projection to the subspace

Span(ao,al,...,aM). It should be outlined that by G we can also obtain a guaranteed upper bound of the error. This fact is very easy to establish in the simplest case of A = ]I. Indeed, V ( u - Uh) "V ( u - uh)dx = / ( a fl

- VUh) . V ( u - uh)dx,

fl

where a is a function in H(div, ft) satisfying (4.4.14). Then, we arrive at the estimate IIV(u - Uh)ll _

inf

a: d i v a + f - - 0

[la - V u h l l ,

(4.4.17)

which shows that

[IV(u--uh)ll <_ IIG(AVuh)- Vuhll.

(4.4.18)

This type a posteriori error estimates have been proposed and investigated in [128, 129, 130] and some other papers. They lead to estimators that measure approximation errors by computing errors in constitutive relations. The major difficulty of this approach consists of finding a suitable projection a the set of functions exactly satisfying (4.4.14). In general, this may be an uneasy task, which makes a straightforward employment of (4.4.18) too expensive. In practice, various procedures that generate functions satisfying (4.4.15) or functions very close to equilibrated fields are used. Equilibrium finite element approximations in linear elasticity can be constructed with the help of the method exposed in [111].

4.5 4.5.1

Error estimation of goal-oriented quantities by gradient averaging techniques General

scheme

Post-processing procedures creates effective a posteriori estimates not only for integral error estimates, but also for errors evaluated by goal-oriented functionals (cf. Section 4.3). In the present section, we describe such a method suggested in [123, 124] for a class of linear elliptic type problems. Below, we consider it for the problem whose generalized solution is defined as a function u E V0 that meets the relation ( ~ u , Aw) = (f, w),

Vw E V0.

(4.5.1)

4.5. GOAL-ORIENTED ERROR ESTIMATES

65

In (4.5.1), .,4 E s U), A e s U), V is a Banach space, Vo is a subspace of V, f E Vo*, U is a Hilbert space equipped with scalar product (., .), and (-, .) denotes the duality pairing of V0 and Vo*. We assume that the operators ,4 and A satisfy the conditions

c~llyll~ ~ (~4y,y) ~ c211yll~, Vy ~ u, c311wlIy ~ IIAwllv, Vw ~ Vo,

(4.5.2) (4.5.3)

Henceforth, the problem (4.5.1) is called the primal problem (or Problem

P). Let e be an element of Vo*. Our aim is to derive estimates of the quantity (~,u - ~),

where ~ is an arbitrary element from V0. For this purpose, we use the adjoint problem (or Problem :Pa): find v E Vo such that (A*Av, Aw) = (~, w)

Vw E Vo,

(4.5.4)

where A* is the operator adjoint to ,4, i.e., (.Ay, z) = (y, .,4*z) for all y, z E U. If Jt is a self-adjoint operator, then both problems are associated with a functional of the type 1

J(w) = -~(AAw, Aw) + (#, w). Under the conditions (4.5.2) and (4.5.3), this functional has a unique minimizer on Vo for any # C Vo*. P r o p o s i t i o n 4.5.1. Let u and v be solutions of problems P and Pa, re-

spectively. Then, (~, u - 5> := E(5, ~) = Eo(5, ~) + E1 (5, ~),

(4.5.5)

Eo(~, V) = (y, ~ - (AA~, A~

(4.5.6)

E~ (~, ~) - (AA(u - ~), A(v - ~)).

(4.5.7)

where

and

C H A P T E R 4. A P O S T E R I O R I E S T I M A T E S F O R F E M

66

Proof. By (4.5.4), we have

(e, u - ~) - (A*Av, A(u - ~)). Since (A*Av, A(u - ~)) - (AA(u - ~), Av), where u satisfies (4.5.1), we find that (e, u - ~) - (AA(u - ~), h(v - v-)) + (f, U) - (AA~, Av-). El In (4.5.5), the term E0(~, ~') is explicitly computable, whereas the term E1 (~, v') contains unknown solutions of the problems (P) and (:Pa). Let Vh and Vr be two finite-dimensional subspaces of Vo, and let ~ - Uh and ~ " - v,, where Uh and v, are solutions of the problems (AAuh, Awh) = ( f , Wh) (A*Av,, Aw~) - (~, w,)

vwh ~ yh, Vwr E v,.

(4.5.8) (4.5.9)

In the particular case Vh -- VT, the relation (4.5.8) implies that Eo(Uh, vr) = 0 so that there remains the term containing the product of the (unknown) energy errors. On the contrary, for noncoinciding spaces both terms Eo and E1 are present. Moreover, it is easy to show that the term Eo dominates if is close to v. Indeed, if ~" -~ v in V, then A(v -v-') -+ 0 in U, so that E1 (~, v~ -+ O. However, Eo(~, v--')+ E1 (~,~') = ( g , u - ~) # 0. Hence, if ~ is sufficiently close to v, then the directly computable term Eo(~, v-') contains the major part of the estimated quantity. The case of nonhomogeneous boundary conditions in Problem :P is reduced to the case under consideration. Indeed, assume that we find a function u ~ Yo + u 0 = {w ~ V ] w = w o + u o ,

uo e V, wo ~ Yo}

that satisfies (4.5.1). Set u - u - uo and represent Problem P as follows: O

find u E Vo such that O

(AA'~,Aw) - (.f, w)

Vw E Vo,

4.5. GOAL-ORIENTED ERROR ESTIMATES where

67

O

(I,

= (I,

- (Ahuo, hw).

Let g E go + uo be an approximation of u. Then, w = g - uo E go is an approximation of ~. By (4.5.5) we obtain 0

0

0

(e,u--~,) = (~,u- w)= (f,v-'~- (,AA w, Av-")+ El(W, v-"). Certainly, finding a sharp approximation of the adjoint problem may require high computational costs. A more economical way is to use an approximate solution of the adjoint problem having approximately the same quality as the approximate solution of the primal one and to recover the unknown functions Au and Av by some post-processing techniques. Let Gh and Gr be averaging operators defined on Vh and Vr, respectively. Replace E(Uh, Vr) by the directly computable functional

E(uh, v~) := Eo(Uh, vr) + Ez (Uh, Vr),

(4.5.10)

where

E1 (Uh, v~) = (A(Gh(AUh) -- Auh), G~(Av~) - Av~). If the operators Gh and G~ perform a proper recovery of Auh andAv~, then it is natural to expect that the difference between Ez (Uh, v~) and E(uh, vr) is given by higher order terms and, thus, the latter quantity can successfully be used instead of Ez. For a class of boundary-value problems, this observation was theoretically and numerically justified in [123]. Below, we present a concise summary of these results. It is worth noting that the efficiency of the above approach strongly improves if one is interested not in a single solution of the primal problem but analyzes a set of approximate solutions obtained for various boundary conditions and right-hand sides (such a situation often arises in the engineering design). In this case, the adjoint problem must be solved only once, and its solution U can further be used in the estimate (4.5.5) for various primal problems.

4.5.2

Example

Let ~ be a bounded domain in n~2 whose boundary consists of two measurable nonintersecting parts 01 ~ and 02~. Consider the problem div(AVu) + f = 0

U----0 AVu. v = g

in

~,

Oil 01~'~, on 02~,

68

C H A P T E R 4. A P O S T E R I O R I E S T I M A T E S FOR F E M

where v is the outward normal to 0f~ and A = {aij (x) }i,j=l 2 is a symmetric positive definite matrix with smooth entries. Assum that f E L2(~),

g E L2(02~),

Now, the problem (4.5.1) has the form (4.5.11) n

n

o2n

In this case, U = L 2 (~, II~2 ), V = H 1(~), V o = { v E H 1(~) I v = 0 o n 0 1 ~ } , and Problem :Pa is to find v E Vo such that f A V v . V w d x = (g.,w) n

Vw E Vo.

(4.5.12)

By(4.5.5), we obtain (4.5.13)

~(u - ~) = Eo(~, v-) + E1 (~, v-), where

f

f~

Oat2

12

and El(u, v~ = / AV(v - ~) " V (u - ~)dx. n Let Vh and Vr be two finite-dimensional subspaces in 17o, constructed by piecewise affine finite element approximations on the triangulations 7~ and Tr, respectively, and Uh and Vr be solutions of the following two problems: /AVuh.Vwhdx-/fwhdx+/gwhds, n

12 A V v r " V w r d x = (~,Wr),

~WhEVh, (4.5.14)

o2~ VWr ~. Vr.

(4.5.15)

4.5. GOAL-ORIENTED ERROR ESTIMATES

69

Define Gh" L ~ (~h, lilt2) ~ H 1(fib, ~2) and G~- 9 L ~ (fL-, ~2) -+ H~ (~2r, ~ 2 ) as gradient averaging operators on patches satisfying (4.4.6). functional (4.5.10) has the form

Now, the

E(wh, w~.) -- Eo(Wh, wT) + E1 (Wh, wT),

(4.5.16)

where

/ flh

Let us qualitatively estimate the asymptotic order of the terms E0, El, and El. We have

IEo(uh,,~)i

A ( V u - Vuh). Vv~ dx <_

-f~

< iiAlio~,~iiVv~i12,~llV(u

- uh)iL2,~ <

< l i A l ] ~ , , ( l l V ( , - v~)li2,~ + liVvii~,~)i]V(u - uh)ll2,,.

Since Vr --+ V in V, the term in brackets is bounded and we see that

IEo(uh, Vr)l <_-chiul2,2,f~,

(4.5.17)

provided that u E H2(f~). If the solutions of Problems ~' and ~'a are sufficiently regular, then the terms E1 (Uh, v~.) and E1 (Uh, vr) are asymptotically smaller ones. Indeed,

tEl (Uh, vr)l <---ch'riul2,2,~ Iv62,2,~ and

HE1(Uh, vr)t --

A(Vuh - Gh(VUh))" (Vvr - GT(Vvr))dx <_ _< IlAlicc,e liVuh - Gh(VUh)l[2,e liVv~ - G~(Vv~)l[2,e.

CHAPTER 4. A POSTERIORI ESTIMATES FOR FEM

7O Since

IlWh - ah(wh)ll~,a

< I I W - Gh(Wh)lt~,n + I l W h - W l I ~ , , ,

IlVur - Gr(Vur)ll~,a < IIW' - C~(W'.)II~,a + IlVu~ - W'll~,a,

the term E1 has the same order hT if both gradient averaging operators possess superconvergence properties. In this case, Eo(Uh, vr) contains the major part of the quantity (e, u - Uh) (except the case Vh =- V~., in which this term vanishes). It remains to compare E1 (Uh, vr) with El(Uh, vr). This question was considered in [123] under additional assumptions necessary to guarantee superconvergence properties of the gradient averaging operators. There it was shown that for sufficiently small h and T, . v

IE1 (Uh, vr) -- E1 (Uh, vr)l -<_ c(hT~ + Th~ + hT(h + T) 89

+ T2) + #(h, T),

where m is any positive integer greater than 2 and #(h, r) contains higher order terms. The asymptotic behavior of the error indicator E suggests the following numerical strategy. 1. Define T~ that is sufficiently fine and accounts particular features of the functional ~. For example, if

(e, u - ~> = f ~o(u - ~) dx, where the weight ~o is a locally supported function, i.e., supp~ C w CC fl, then the mesh should have more degrees of freedom in a neighborhood of w. Solve Problem 7~a on Sir and find v~; 2. Define Vh, solve the primal problem, and find Uh; 3. Compute Eo(Uh, vr) and use post-processed values of VUh and Vvr to evaluate E1 (Uh, vr). It is worth noting that if we have a series of approximate solutions ukh, k 1, 2, ...m, computed for different right-hand sides and boundary conditions, then the operation of in part 1 must be performed only once.

4.5. GOAL-ORIENTED ERROR ESTIMATES 4.5.3

71

E r r o r e s t i m a t e s in t e r m s o f s e m i n o r m s

E s t i m a t e s of u - ~ computed on a set of linear functionals can be used for an evaluation of this difference in some seminorms. Let w be a chosen collection of simplexes and {~1, ~o2,..., ~Od} be a set of functions t h a t vanish in fl \ ca. By the m e t h o d described above, we can e s t i m a t e the quantities

I~,. "= J ~os(x)(u(x) - ~(x)) dx. o,)

However, we are also interested in estimating a certain norm of u - ~, for example, the norm

Ilu - ~ll~,~ = f I~(~) - ~(~)1 ~ d~. Note t h a t

f ~?(u ~)dx -

l i ~ - ~11~,~ -

~up

~

and the s u p r e m u m is attained if ~ = u - ~. Thus, if the difference u - ~ is known to belong to a certain set S(w) C L2(w), then the problem is reduced to the evaluation of the seminorm

f v(u - ~)dx lu-~lz'=

sup

~

In general, a seminorm does not provide complete information on the error, because the relation [u - ulz = 0 does not mean t h a t u = ~. Nevertheless, seminorms m a y give a useful unformation if E contains sufficiently large a m o u n t of linearly independent functions. Let E = Span{~ol, ~02, ..., ~d}. Then, d

E ai f ~oi(u - ~)dx Iu - u l ~ =

sup

i=1

~

1.

aiaj f ~oi~ojdx i,j=l

w

(4.5.18)

72

C H A P T E R 4. A P O S T E R I O R I E S T I M A T E S F O R F E M

This problem is equivalent to the following one" find the quantity /~2(W, E)

--

inf B a . a ,

(4.5.19)

aEQ

where a = (al, a2, ..., ad), I -- (Iv, , Iv2,... , Ivy), Q-{aE~d

l a'l-1},

{bij}i,j=l,d,

bij - / ~ i ~ j

and

B-

dx.

o)

The problem (4.5.19) has a minimax form a2 (w,.=.) -

inf sup{B a. a + A[(a. I) - 1]} >_ XER > sup inf ,rB a a + A[(a I)

aERa

AER aERaL

"

.

1]} _

(4.5.20) .

Assume that the functions ~i are chosen in such a way that B is a positive definite matrix. Then (4.5.20) holds as equality. The minimization problem that stands on the right-hand side in (4.5.20) has a solution a - - ~)~B - 1 i . Therefore,

A2

}

a2(w,E) -- supper - - ~ - B - 1 I 9 I and

1

Iu - ~ l ~ = ~ ( ~ , s )

A

1 - B-1I. I'

= ( B - 1 I 9 I)89

We see that the value of ]u - ~]~ is not difficult to find, provided that the functions {~s} are properly chosen and the respective quantities Iv. are defined. Let us show that if ~ E E(w) and u is a sufficiently smooth function, then [u - ul-: gives an estimate of Iiu - u[[2,o+. Let E(w) be a subspace of L2(w) and ~ be an arbitrary function from E(w). Then

f ((u - ~)~ -4- ( ~ - ~)~) dx I 1 ~ - ~11~,~ -

sup ~ ,~i=(~)

< -

Ilvl12,~

f (~ - ~)~ dx

<

sup

oJ

+ I]~ - ~II~,=.

(4.11)

73

4.5. G O A L - O R I E N T E D E R R O R E S T I M A T E S Since ~ - ~ E E(w), we have

$(~Ii~

~11:,~,

-

< sup - ,~(~)

"

~),7 d~

11,711~,~

+

I1~

-

~11~,~,.

By rearranging again the numerator of the fraction on the right-hand side, we arrive at the estimate

f ((~ - u)rl + (u - u)rl) dx Ilu

- ~11~.,~ -

sup w ,~z(~)

q - I l u - ~11~,~ <

11,7112,~,

f (u ~)o dx -

_< sup

"

,~--(~)

I1~11~,~

+

211u- ~ll2,,,, - l u -

~lz +

2llu- ~ll2,.,,

in which ~ is an arbitrary function from E(w). Hence, we see that

(4.5.21)

25(u, .=.(w)),

I1~- ~11~,~ __ lu - ~1= + where 5(u,E(w)) =

inf

I1~- ~il2,~

is the distance between u and 2(w). If E(f~) is a set of polynomials and u is a smooth function then 5(u, .=.(w)) can be estimated with the help of well-known results of approximation theory (see, e.g., [53]). By similar arguments, we obtain

f ((u - ~)~7 + ( ~ - ~)~7) dx

I1~- ~11~.,~

-

sup ,~(~)

~'

>

Iluli~,~

_>

f (~ - "5)rTdx sup

~

~112,~.

-I1~-

Since ~ - ~ belongs to E(w), we find that

f ( ( ~ - u)~7 + (u - "5)~7) dx

I1~-

oJ

~112,~, > sup f(u

_

sup ,~--(~)

-

f (~ - u)rl dx

"5)~7 dx

o)

Ilull~,~

sup ,~L~(~)

to

I1~11~.,~

-II~-

= I~ - ~lz - 211~ - ~11~,~,

~ll~,~ = v~

c

s(w).

74

CHAPTER

4. A P O S T E R I O R I

ESTIMATES

FOR FEM

From here, it follows that Ilu - ~]12,,. _> ]u - ~]~ - 2(~(u, E(w)).

(4.5.22)

By (4.5.21) and (4.5.22), we conclude that the error arising if JJu - UhJJ2,~ is replaced by lu - Uhl=_ depends on the regularity of u and approximation properties of S(w).

4.6

N o t e s for the Chapter

A posteriori error estimation for finite element approximations is the main subject of the books by M. Aintworth and T. Oden [3], I. Babu~ka and T. Stroboulis [18], W. Bangerth and R. Rannacher [26], K. Eriksson, D. Estep, P. Hansbo and C. Johnson [6S], and R. Verfiirth [213]. They present various approaches, results of numerical experiments and a wide list of references. Below, we present a concise overview of some other publications related to the material considered in Chapter 4. Explicit residual method was originated by I. Babu~ka and W. C. Rheinboldt [15, 16]. Later it was investigated and extended by many authors. Among other publications associated with residual based methods we refer to M. Ainsworth and J. T. Oden [1], C. Carstensen [42], C. Carstensen and S.A. Funken [44, 45] (in these papers the authors consider ways of computing bounds of the interpolation constants in the residual based error estimator), C. Carstensen and R. Verfiirth [46], K. Eriksson and C. Johnson [69, 70], C. Johnson and P. nansbo [110]. Estimates for problems with biharmonic operator are analysed in A. Charbonneau, K. Dossou and R. Pierre [47]. A posteriori estimates taking into account the influence of of the non-discretized part of the domain on the approximation error are considered in W. Dhrfler and M. Rumpf [60]. Estimates in the Lee-norm can be found, e.g., in S. W. Brady and A. R. Elcrat [31]. Adaptive methods for convection-diffusion problems are considered in C. Johnson [109]. A posteriori error estimators based on solving local Neumann type problems on each element were proposed in R.E. Bank and A. Weiser [24]. Its asymptotic properties were investigated in R. Duran and R. Rodriguez [63]. This error estimation method was in the focus of many researches. We present several publications, in which readers will find other references: I. Babu~ka, T. Strouboulis, C. S. Upadhyay and S. K. Gangaraj [19], I. Babu~ka, F. Ihlenburg, T. Strouboulis and S. K. Gangaraj [14], E. Stein and S. Ohnimus [202], E. Stein, F.J. Bartold, S. Ohnimus and M. Schmidt [203]. Probably earliest results on superconvergence were established in the papers by L. A. Oganesjan and L. A. Ruchovets [155] and M. Zl~mal [229].

4.6. N O T E S FOR THE C H A P T E R

75

Also, we refer to I. Hlavs and M. Kfi~ek [101], M. Ki~ek and P. Neittaanm~iki [127], P. Neittaanm~iki and M. Kii~ek [148]. In the books L. Wahlbin [216] and M. Ki~ek, P. Neittaanm~iki and R Stenberg [125] readers can find a detailed exposition of the subject and other references. Surveys on superconvergence are presented in M. Ki~ek and P. Neittaanm~iki [126, 127] and J. R. Whiteman and G. Goodsell [219]. Simple and efficient error indicators based upon gradient averaging were suggested in O. C. Zienkiewicz and J. Z. Zhu [226]. This and other indicators were analyzed in many papers. Here, we refer to, e.g., M. Ainsworth, J. Z. Zhu, A. W. Craig and O. C. Zienkiewicz [5], I. Babu~ka and R. Rodriguez [17], R. D. Lazarov [136], R. Rodrigues [189], O. C. Zienkiewicz and J. Z. Zhu [227, 228], O. C. Zienkiewicz, B. Boroomand and J.Z. Zhu [224]. In C. Carstensen and S. Bartels [43], it is shown that each averaging technique leads to a certain a posteriori error estimate. Estimates based upon global averaging of gradients were considered in B.-O. Heimsund, X.-C. Tai and J. Wang [115]. In J. Wang [217], a generalization of the patch recovery technique of O. C. Zienkiewicz and J. Z. Zhu is considered. It consists of approximating an approximate solution data by a certain smooth function that is computed by the least-squares method. A posteriori estimates that express approximation errors by errors in the equilibrium relations were suggested in P. Ladev eze and D. Leguillon [128] (see also P. Ladev eze, J.-P. Pelle and Ph. Rougeot [129], P. Ladev eze and Ph. Rougeot [130]). Post-processing procedures based upon equilibrating of recovered stresses were considered in B. Boroomand and O. C. Zienkiewicz [29, 30]. They result in stress fields that satisfy an equilibrium condition in a weak form. A method of equilibration is presented in P. Destuynder and B. Mdtivet [59]. This paper also contain numerical tests, in which a posteriori error estimates obtained by equilibration are compared with other estimates. A posteriori error estimates are intended not only for the verification of the accuracy of numerical solutions. They must also provide a useful information for mesh adaptation algorithms. These questions are analyzed, e.g., in R. E. Bank and R. K. Smith [23], G. F. Carey and D. L. Humphrey [41], R. LShner, K. Morgan and O. C. Zienkiewicz [139], J. T. Oden, L. Demkowicz, T. Strouboulis and P. Devloo [150], J. Peraire and A. T. Patera [159], R. Verfiirth [213]. Methods of mesh adaptation are considered in the book by W. C. Rheinboldt [167]. The reader can find a detailed exposition of the a posteriori error estimation methods based upon using adjoint boundary-value problems in, e.g., W. Bangerth and R. Rannacher [26], R. Becker and R. Rannacher [27]. Applications of this techniques to elastic and elasto-plastic problems

76

CHAPTER 4. A POSTERIORI ESTIMATES FOR FEM

were considered in R. Rannacher and F. T. Suttmeier [165, 166], and to transport problems in P. Houston, R. Rannacher, E. Siili [103]. Applications to problems in the theory of optimal control were considered in R. Becker, H. Kapp and R. Rannacher [28]. Error estimation methods for goal-oriented quantities were considered, e.g., in M. Ainsworth and J .T .Oden [3] and J. T. Oden and S. Prudhomme [152]. Finally, it is necessary to say a few words about extensions of the residual based approach to nonlinear problems. In this case, the correlation between errors and residuals is not so obvious. This fact is easy to discover even in elementary examples. Let, for example a real-valued function r have an isolated root Xo. If r is a smooth function, then 1 ,, -

+

- x0) + 5r

-

+

. . .

and we see that [ I x - Xo[[ < [r162

(4.6.1)

provided that r ~ 0 and ] I x - x01] is so small that the third and other higher terms of the expansion are negligible with respect to the linear one. However, for a "distant" point 5 in which r is small the estimate (4.6.1) may be wrong because of the evident violation of the above locality condition. This simple example shows that, in general, error estimation methods based upon calculation of residuals require a deep investigation of the corresponding operator (mapping) at a neighborhood of the exact solution. Roughly speaking, a suitable estimate can be obtained if the operator considered admits a proper linearization at x0. Such an estimate is valid in a ball centered at xo whose radius may be a priori unknown and, therefore, should be defined by an additional investigation. For these reasons, a posteriori error estimates for the abstract problem

F(x) = 0 , where F : X -+ Z is a continuous mapping of a Banach space X to a Banach space Z are usually derived under additional assumptions on F (see, e.g., J. Pousin and J. Rappaz [161, 162] and R. Verfiirth [212, 213]). In particular, F must be a C 1 mapping, which possesses a unique solution u in a certain ball B(u, 5). Moreover, the differential DF(u) must be an isomorphism in a neighborhood of u. An extension of this approach to parabolic problems is exposed in R. Verffirth [213]). Residual type a posteriori error estimates for finite element approximations of elliptic obstacle problems were obtained in Z. Chen and R. H.

4.6. N O T E S FOR THE C H A P T E R

77

Nochetto [48]. Their estimator contains terms of four types. Two of them are the same as in residual type estimates for linear elliptic problems and two others ("obstacle oscillations" and "obstacle jump residuals") reflect the behavior of an approximate solutions near the coincidence set. Various a posteriori error estimates for finite element approximations of variational inequalities were suggested by M. Ainsworth, J. T. Oden and C. Y. Lee [4], R. Hoppe and R. Kornhuber [102], R. Kornhuber [122], A. Veeser [210]. In J. Medina, M. Picasso and J. Rappaz [141], a priori and a posteriori error estimates are presented for nonlinear diffusion-convection problems. An aproaches to a posteriori error estimation for elasto-plastic problems is exposed in M. Schulz and W. L. Wendland [199, 200]. A posteriori error estimates for finite element approximations of problems in the theory of fluids are considered in, e.g., M. Ainsworth and J. T. Oden [2], M. Amara, M. Ben Younes and C. Bernardi [6], R. E. Bank and B. D. Welfert [25], C. Johnson and R. Rannacher [112], C. Johnson, R. Rannacher and M. Boman [113], J. T. Oden, W. Wu and M. Ainsworth [154], C. Padra [158], W. Strouboulis and J. W. Oden [205], R. Verfiirth [211].