Stable multiscale bases and local error estimation for elliptic problems

MATHEMATICS ELSEVIER

Applied Numerical Mathematics 23 (1997) 2147

Stable multiscale bases and local error estimation for elliptic problems Stephan Dahlke

a,1

Wolfgang Dahmen

a,*,2,

Reinhard Hochmuth

b,3,

Reinhold Schneider c

Institut fiir Geometrie und Praktische Mathematik, RWTH Aachen, Templergraben 55, 52056 Aachen, Germany b Fachbereich Mathematik und Informatik, Institutfiir Mathematik I, Freie Universitiit Berlin, Arnimallee 2-6, 14195 Berlin, Germany c Fachbereich Mathematik, Technische Hochschule Darmstadt, Schloflgartenstrafle 7, 64289 Darmstadt, Germany

Received 22 December 1995; revised 4 September 1996

Abstract This paper is concerned with the analysis of adaptive multiscale techniques for the solution of a wide class of elliptic operator equations covering, in principle, singular integral as well as partial differential operators. The central objective is to derive reliable and efficient a-posteriori error estimators for Galerkin schemes which are based on stable multiscale bases. It is shown that the locality of corresponding multiresolution processes combined with certain norm equivalences involving weighted sequence norms of wavelet coefficients leads to adaptive space refinement strategies which are guaranteed to converge in a wide range of cases, again including operators of negative order. Keywords: Stable multiscale bases; Norm equivalences; Elliptic operator equations; Galerkin schemes;

A-posteriori error estimators; Convergence of adaptive schemes

1. Introduction The increasing importance of adaptive techniques in large scale computation is reflected by a vast amount of recent literature on this topic primarily in connection with finite element schemes (see, e.g., [1-4,7,26,32,35]). How to fit such adaptive techniques into the context of stable splittings for multilevel Schwarz type preconditioners for elliptic problems has been briefly indicated in [30]. On * Corresponding author. E-mail: [email protected]. I The work of this author has been supported by Deutsche Forschungsgemeinschaft (Da 360/1-1). 2The work of this author has been supported in part by Deutsche Forschungsgemeinschaft (Da 117/8-2). 3The work of this author has been supported by the Graduiertenkolleg 'Analyse und Konstruktion in der Mathematik' funded by Deutsche Forschungsgemeinschaft. 0168-9274/97/$17.00 Copyright © 1997 Elsevier Science B.V. All rights reserved PII S0168-9274(96)00060-8

22

S. Dahlke et aL / Applied Numerical Mathematics 23 (1997) 21-47

the other hand, there have been several attempts to apply wavelet concepts to the solution of differential and more generally pseudo-differential equations. As the above mentioned multilevel schemes these concepts hinge upon making successive corrections of current solutions when progressing to finer scales of discretization. However, the wavelet methodology differs from conventional finite element or finite difference discretizations and associated solvers in that direct use of bases is made which span appropriate complements between successive trial spaces. Moreover, loosely speaking, the representation of elliptic operators and their inverses relative to wavelet bases is nearly diagonal. The benefit gained from such concepts has, however, usually be paid for by considerable efforts in constructing suitable wavelet bases. On the other hand, several cases have been studied where such bases are available and have proved to offer significant advantages. For instance, pre-wavelets have been shown to yield robust preconditioners in combination with sparse grid discretizations for two- and three-dimensional anisotropic problems [25]. Divergence-free wavelets with small support (realizing any degree of accuracy and whose construction is essentially the same for any spatial dimension) have been applied to the Stokes problem [34]. More generally, there are examples of wavelet bases with built-in Ladyzhenskaya-Babu~ka-Brezzi condition for various types of saddle point problems [15]. Finally, such bases give rise to matrix compression techniques when dealing with discretizations of pseudo-differential or singular integral operators as they arise, for instance, in connection with boundary element methods. Wavelet bases defined on two-dimensional manifolds in It~3 satisfying all the requirements which guarantee optimal compression and convergence rates have now become available [19]. This is an important aspect, since for such equations on two-dimensional manifolds the fact that corresponding system matrices are not sparse is the major computational bottleneck. So far the compression techniques apply to Galerkin or collocation schemes based on essentially uniformly refined trial spaces. In view of the above remarks, it is natural to complement the sparse (or nearly sparse) wavelet representation of operators by sparse representations of the solutions. In fact, it is well known that singular behavior of a function is directly reflected by large wavelet coefficients near the singularities while in regions of high regularity the wavelet coefficients have very small modulus. As a consequence often only relatively few significant wavelet coefficients suffice to represent a function with good accuracy. Detecting these coefficients during a solution process and confining all calculations essentially to these coefficients suggests itself as a natural adaptive solution concept. Apparently this differs in many respects from adaptive techniques in connection with multigrid or other multilevel methods based on finite difference or finite element discretizations. The above mentioned instances of promising stable multiscale bases for PDEs as well as integral equations have stirred our interest in substantiating the still somewhat vague and heuristically founded adaptive potential of multiscale bases. Therefore the objective of this paper is to show that once suitable multiscale bases are available • a unified adaptive strategy converges for a wide class of elliptic operator equations covering differential as well as integral operators, • while neither requiring any a priori knowledge about possible types of singularities nor assuming the validity of the saturation property. While in conventional finite element or boundary element schemes the typical a-posteriori error indicators depend very much on the specific underlying discretization unified means here that essentially the same type of criterion applies throughout the above mentioned range of cases independent of the specific wavelet. In fact, local error estimators arise essentially as coefficients of corresponding multiscale expansions of residuals. It is also worth stressing that adaptive refinements in the present

S. Dahlke et al./ Applied Numerical Mathematics 23 (1997) 21--47

23

setting do not require handling any geometric problems of constructing suitable meshes. The refinement process concerns only the selection of those basis functions which capture the significant part of the searched object and thus involves essentially the efficient administration of index sets. Moreover, whereas in most studies of adaptive schemes the so-called saturation property is assumed (see, e.g., [3,7,29]) or its validity can only be established under relatively strong assumptions on the data we will show that under rather weak assumptions on the right hand sides this property is indeed enforced by the adaptive scheme. So assumptions on the unknown solution can be reduced to a relatively very weak assumption on accessible data which is actually the most one can hope for. A particular further motivation for us is that for integral equations the understanding of local error estimators is comparatively less developed than for partial differential equations. For instance in [9] lower and upper a-posteriori error estimates for the energy norm are derived for first-kind integral equations but relative to quasi-uniform meshes. In the study of Haar-wavelet-based multilevel schemes for single layer potentials [29] the saturation property is assumed. To our knowledge the results of this paper offer for the first time a provably convergent adaptive scheme based on reliable and efficient a-posteriori error estimators also for integral operators. By this we mean that without any restrictions on the mesh the current error is bounded from above and below by expressions involving computable local quantities. The layout of the paper is as follows. In Section 2 we describe the class of problems to which the subsequent analysis applies. As mentioned before operators of positive as well as negative order are covered. Ellipticity in this context means that the operator under consideration maps a certain Hilbert space in a boundedly invertible manner into another Hilbert space. We also assume that the Schwarz kernel of the operator has certain asymptotic properties matched in many cases of interest and that Galerkin schemes are stable. We deliberately choose such an abstract setting in order to bring out clearly the role and relevance of various functional analytic ingredients. The analysis is essentially based on two pillars, namely (i) the validity of equivalences of Sobolev norms and weighted sequence norms for wavelet coefficients in a range depending on approximation properties and regularity of the underlying basis; (ii) compression of matrices in wavelet representation. Both effects distinguish wavelet-based discretizations in an essential way from conventional finite element or finite difference discretizations. Specifically the fact that the range of norm equivalences cover Sobolev spaces of positive and negative order seems to be crucial for advancing the analysis so that, for instance, the saturation assumption can be dispensed with. Combining (i) with the ellipticity yields asymptotically sharp lower and upper bounds for residuals and thereby for the current error of a Galerkin approximation in an energy-like norm. The fact that these bounds can be reduced to quantities which within a given tolerance still provide lower and upper bounds but are computable at reasonable cost is facilitated by (ii). At this point we benefit from earlier investigations of wavelet methods for pseudo-differential equations (see, e.g., [17,18]) by which we can sufficiently quantify the fact that the wavelet representation of a wide class of operators is almost diagonal. We find it particularly promising that on one hand wavelet discretizations lead to sparse representations of operators which usually give rise to densely populated matrices while on the other hand, due to their built in analysing capabilities they suggest also a very natural adaptive selection of trial spaces. Therefore we briefly review in Section 3 some results on multiscale decompositions, some related stability concepts and resulting norm equivalences. In principle, these facts are known. Since, as

S. Dahlke et al. /Applied Numerical Mathematics 23 (1997) 21-47

24

mentioned above, they play a crucial role in the present context and since the results in the form needed here are only very recent we decided to include this material. So those readers who are familiar with these facts should simply skip over this section. Section 4 is devoted to the derivation of a-posteriori error estimators in terms of wavelet coefficients of the residual with the aid of (i) while as mentioned before the reduction of these bounds to computable quantities is facilitated by (ii). Special cases of such estimates for univariate second-order two-point boundary value problems already appear in [5]. Furthermore we propose an adaptive procedure and prove its convergence under very mild assumptions on the right hand side data. We emphasize that no a-priori assumptions on the solution such as the occurrence of certain types of singularities or the saturation property are needed. Although technically completely different the procedure is similar in spirit to that in [23] for piecewise linear finite elements for the Laplace equation. In this sense the present investigation has been to some extent inspired by the work in [5,23]. While the a-posteriori error estimates are not confined to symmetric problems symmetry is used though in the convergence proof. We conclude Section 4 with briefly indicating that also in this context the assumption of symmetry can be relaxed. Incidentally the discussion of a typical example concerning dominating convection sheds some light on the limitations of such arguments. While all assertions remain valid in an asymptotic sense certain constants are seen increase with growing convection which eventually excludes the applicability of the estimates in any discretization regime of practical relevance. Of course, one has to question the suitability of classical Galerkin schemes for problems of this type. This paper focusses mainly on the adaptive selection of a discretization in a Galerkin context not on the solution of the corresponding systems of equations. On the other hand, as in the multigrid context, both issues are closely intertwined. Therefore we briefly address this issue in Section 5, in particular, with regard to preconditioning.

2. The scope of problems and basic assumptions We will be concerned with linear operator equations

Au = f

(2.1)

where A will be assumed to be a boundedly invertible operator from some Hilbert space H1 into another Hilbert space He, i.e., IIAUIIH2 ~,~ ]IU[IHI,

zt

E H1,

(2.2)

where "a ~ b" means that both quantities can be uniformly bounded by some constant multiple of each other. Likewise "~<" indicates inequalities up to constant factors. We will write out such constants explicitly only when their value matters. To give an idea of the scope of problems we have in mind, we will briefly indicate a few representative examples. An important case arises when

a(u,v) = (Au, v),

u,v E Hi,

(2.3)

defines a symmetric bilinear form on H1. Here ellipticity means that a(u,u)

N

Ilull ,,

u e H1,

(2.4)

S. Dahlke et al. / Applied Numerical Mathematics 23 (1997) 21--47

25

which is easily seen to imply (2.2). Of course, the simplest examples of this type arise when Y2 C ]~d d a2/Ox~ is the Laplacian is some bounded domain, Au = - A u or Au = - A u + cu where A = ~--~j=l and c > 0. Here H 1 ~ H01 (~'-~) o r H 1 z H1 (J'~) while//2 = H-I(Y2), (Hi(y2)) ' (which denotes the dual of H 1(Y2)), respectively. As usual for m E No, Hm(Y2) denotes the Sobolev space of square integrable functions on Y2 whose ruth-order weak derivatives belong to L2(Y2) and H~(Y2) is the closure of C~(Y2) with respect to the 1[-I[Hm-norm. For s E ~ the spaces HS(Y2) are then as usual obtained by interpolation or duality. In order to solve an exterior boundary value problem AU=0

in]~3\$2,

U=f

onSY2,

it is tempting to transform this problem into a boundary integral equation. For instance, the indirect method leads to the equation

A u = f,

A = I + 2K,

where

Ku(x) : : [½ -

+

1 Jfny.(x-y) :

It(y)

dy

at2

and O~(x) denotes the interior angle at x E 0Y2 located on an edge of 0Y2. One can then recover the solution U of the exterior boundary value problem from the solution u of the integral equation on the domain OY2. In this case it is known that Hi = / / 2 = H ° = L2(F), F := aY2. Alternatively, the direct method yields an integral equation of the first kind

Au = V u = (½ + K ) f,

Vu(x) = f 4 7 uO) f i x - Yl dy. OY2

It is known that this fits into the above framework with H1 = H-1/2(F), H2 = H1/2(F). Due to the fact that V has order different from zero the problem of preconditioning corresponding stiffness matrices arises as in the case of differential operators [18] while on the other hand this formulation is better suited for handling different kinds of boundary conditions or transmission problems. At this point we will not discuss any of these concrete cases further but rather collect next those assumptions and ingredients which will matter in our analysis. According to the above examples, let Y2 be some domain or more generally a manifold and let L2(Y2) be the space of square integrable functions on Y2. We will always assume that either Hi C L2(~2) c/-/2 or/-/2 C L2(Y2) C H1 where inclusion is to be understood as continuous imbedding and Ha = H~ is the dual space of H1 relative to the dual form (., .) induced by the inner product in L2(Y2). We emphasize that it is not essential in what follows that (2.1) is a scalar problem. Correspondingly Hi could be products of Sobolev spaces. To simplify the exposition we will assume though throughout the remainder of this paper that H1 = H t, H2 = H -t for some t E R where H t denotes a Sobolev space and H - t will always stand for its dual space (which does agree with H-t(Y2) when Y2 has no boundary or when H t = H~(O)). For simplicity we will also write I1" lit = I1" lira. Relation (2.2) then reads

clllAull_t <~ Ilullt ~< c211Aull-t,

u c n t,

(2.5)

S. Dahlke et al. / Applied Numerical Mathematics 23 (1997) 21-47

26

and we will refer to the positive constants c1, c2 later again. Moreover, we will assume that A is still continuous in a Sobolev scale near t, i.e., there exists some 7- > 0 such that

IIA ll-t+

v

c H t+s, 0 <<.Isl

~

(2.6)

As pointed out above we will admit operators with Schwarz kernels of global support

(Av)(x) = I K(x,y)v(y)dy, ~2

where we require that whenever d + p + lal + I/3[ > 0

10 0 K(x,y)l

dist(x,Y) -(a+p+l~l+l~l)

(2.7)

holds with constants depending only on the multi-indices a,/3 E Za+ where Z d := {a = ( a l , . . . , ad) E zd: ai ~> 0, i = 1 , . . . , d}. Estimates of the type (2.7) are known to hold for a wide range of cases including classical pseudo-differential operators and Calder6n-Zygmund operators (see, e.g., [18,33]). Thus the single and double layer potential operator above as well as classical differential operators fall into this category. Finally, we will assume that Galerkin schemes for (2.1) are stable. When A is positive definite and selfadjoint this is the case for any trial space. More generally, in the framework of pseudo-differential operators injectivity of A, coercivity of the real part of the principal part of the symbol and certain approximation properties of the trial spaces also imply stability of corresponding Galerkin schemes even when A is not symmetric [18].

3. Muitiresolution Our goal is to develop adaptive Galerkin methods for the approximate solution of (2.1). However, in contrast to conventional finite element discretizations we will work with trial spaces that do not only exhibit the usual approximation properties and good localization but in addition lead to expansions of any element in the underlying Hilbert spaces in terms of multiscale or wavelet bases with certain stability properties. We will show that these stability properties ultimately stated in terms of norm equivalences for Sobolev spaces will indeed allow us to improve on previous theoretical investigations for the above scope of problems. To correspond to the above range of applications we formulate the relevant facts for the following general framework. These results are essentially known (cf. [14,18]) but for the convenience of the reader we include a brief summary of the relevant facts. Suppose H is a Hilbert space (of functions defined on 12, say) with inner product (.,-). Throughout this section orthogonality will always be understood relative to this inner product. Again typical examples are H = L2(f2), H = Hs((2) or products of such spaces. Let S = {Sj}j= o be a sequence of closed nested subspaces of H whose union is dense in H. In all cases of practical relevance the spaces Sj are spanned by single scale bases ~j = {¢j,k: k EIj} which are uniformly stable, i.e.,

Ilcllt=(ij)

~

Ckej,k

(3.1)

H

uniformly in j C No. Here we denote as usual I1" l[2 = (, ") and [[cll2=(ij) =

Ick[ 2-

S. Dahlkeet al. / AppliedNumericalMathematics23 (1997)21-47 Successively updating a current approximation in stable bases

27

Sj-l to a better one in Sj can be facilitated if

Jj} for some complement Wj of Sj_I in Sj are available. Defining for convenience ~0 := ~0, W0 := So, any vn = 2kEIn CkOn,k E Sn has then an alternative multiscale representation k~j = {¢j,k: k E

n

5=0 kEJj which corresponds to the direct sum decomposition n

sn =Gw . j=0

Of course, there is a continuum of possible choices of such complements. Orthogonal decompositions would correspond to wavelets. However, orthogonality often interferes with locality and the actual computation of orthonormal bases might be too expensive. Moreover, in certain applications orthogonal decompositions are actually not best possible [18]. The essential constraint on the choice of Wj is that

jEN0

forms a

Riesz-basis of H, i.e., every v E H has a unique expansion o~

(3.2) j=0 kEJj such that llvlIu ~

~

~

j=0

kEJj

where ~ = {~j,k: k E

I(v,¢j,k)l:

,

v E H,

(3.3)

Jj, j E No} forms a biorthogonal system

(¢j,k,~j',k') = 5j,j,Sk,k,,

j,j' E No,

k E

gj,

k' E Jj,,

(3.4)

and is in fact also a Riesz-basis for H (cf. [14]). To explain one aspect why this is important let Tn denote the transformation that takes the coefficients dj,k in the multiscale representation of vn into the coefficients Ck of the single scale representation. It corresponds to the synthesis part of the fast wavelet transform. In fact, it is known that the Riesz basis property of gt is equivalent to Tn being well conditioned, i.e.,

IITntl, llT:l[I = O ( 1 ) ,

n -+ ec,

(3.5)

where I1" II denotes the spectral norm [13,141. With such a pair of biorthogonal bases ~ and ~ one can associate canonical truncation projectors n

Qnv := Z j=0

n

Z (v'~j,k)ff)j,k' kEJj

Otnv := Z j=0

Z (v,~yj,k)~)j,k kEgj

S. Dahlke et al. / Applied Numerical Mathematics 23 (1997) 21~17

28

which are obviously adjoints of each other. Of course, when ~/' is a Riesz-basis then the Qn and hence their adjoints Q~n are uniformly bounded in H. Denoting by Sn the range of Q~n we have therefore two sequences S and ,5 of nested closed subspaces Sj and Sj, respectively, whose union is easily seen to be dense in H [13]. While the Riesz-basis property of k~ implies the existence of a biorthogonal Riesz-basis ~ as well as the uniform boundedness of the projectors Qn and Q~, the converse is known not to be true in general [14]. Additional conditions that do ensure the Riesz-basis property for a general Hilbert space setting have been established in [14]. Here we are only interested in their specialization to the particular case H = L2(~2) (where as above Y2 is either a closed surface or a domain in R d) and denote by H 8 a corresponding scale of Sobolev spaces. What turns out to matter is that both S and ,,~ should have some approximation and regularity properties which can be stated in terms of the following pair of estimates. There exists some "7 > 0 such that the inverse estimate

Ilvnllgs(m < 2nSIIvnllr=(m,

vn E & ,

(3.6)

holds for s < "7. Moreover, there exists some m ~> "7 such that the direct estimate inf I I v - v llLz(m < 2-SnllvllHs(

vn ESn

),

v ~ H*(Y2),

(3.7)

holds for s ~< m. Such estimates are known to hold for every finite element or spline space. For instance, for piecewise linear finite elements one has "7 = 3/2, m = 2. It will be convenient to introduce the following notation. Let OO

J := {)~ = (j,k): k E J j , j E No} = U ({J} x Jj), j=O and define I,Xl := j

if )~ E Jj.

The following result from [14] will play a central role in the subsequent analysis. Theorem 3.1. Suppose that ko = {Ca: ,X E J } and ~ = {~;~: ,k E J } are biorthogonal collections in L2( f2) and that the associated sequence of projectors { QJ}j=0 is uniformly bounded. If both S and satisfy (3.6) and (3.7) relative to some "7, "7' > O, "7 <~m, "7' <~m', then 1/2

IIvlIH*(a2)<'(Z22PqSI(V,~bA)I 2) \ XEJ

,

s e (-'7','7), (3.8)

1/2

,

s E (-'7,'7'),

v E H*(,Q).

" )~EJ

Moreover, the projectors Qj, Q~ are uniformly bounded in Hs(g2), s E (-7','7) and s E (-'7,'7'), respectively. R e m a r k 3.1. Instead of powers of 2 in (3.6)-(3.8) we could have used powers of a for some a > 1 which reflects the subdivision rate of successive refinement levels. Since this entails no essential differences we will stick in the following with halving mesh sizes.

S. Dahlke et al. / Applied NumericalMathematics23 (1997) 21-47

29

For more information about the construction of multiscale bases ~, ~ with the above properties the reader is referred to [8,16,19]. Throughout the remainder of this paper we will assume that ~ and satisfy the assumptions of Theorem 3.1 where for t as in (2.5) d 7 > t + ~,

d 7 t > - t + ~,

(3.9)

so that there exist positive constants 0 < c3, c4 < oc for which

c3(Z2-etlal(v~' )~)lz)/'e)~cJ

<<"I[v'l-t <<"c4(z2-etlXl(v~' )~)le)l/e-acJ

Finally, for our applications it will be important to work with assume that diam(supp qSn,k), diam(supp ~n,k) ~ 2 - n ,

(3.10,

local bases, i.e., we will always

n E N.

(3.11)

Furthermore, it is desirable that the ~ , k , ~n,k have the same property diam(supp~,k),diam(suppz~,k)

~ 2 -~,

n E N.

(3.12)

4. Multiscale error estimates

4.1. Some preliminary remarks In the symmetric case (2.3), (2.4) it is natural to estimate the accuracy of the Galerkin solution with respect to the energy norm

Iq" II := a(., 01/2 Since by (2.4), this norm is equivalent to II-IIH, = II" lit we will employ this latter norm most of the time in the general (possibly non-symmetric) situation. To motivate the subsequent development we briefly recall the most common starting point of adaptive strategies for selfadjoint elliptic problems. The basic idea is very simple (see, e.g., [3,7]). Suppose that u t, u" are Galerkin approximations to the solution u of (2.1) from spaces S ~, S" respectively, where u" is a more accurate solution, i.e., S t C S". R e m a r k 4.1. When A is selfadjoint positive definite one has for u t, u", as above

IW' - u'11 ~< Ilu - u'll,

(4.1)

Moreover, one has for some/3 < 1

I1~ - u'll ~< (1 -/32)-1/211u" - u'll

(4.2)

if and only if

Ilu - u"ll ~
(4.3)

The assertion is a trivial consequence of the orthogonality of the error u - u" to u t - u" so that Ilu - u'll 2 = biu - u"ll 2 + Ilu" - u'll 2

30

S. Dahlke et al. / Applied Numerical Mathematics 23 (1997) 21-47

The relation (4.3) is often called saturation property requiring that the new approximation u" is strictly better than the previous one u'. Thus if the saturation property holds the quantity Ilu' - u"ll provides a computable lower and upper bound for the tree error Ilu - u'll. Such a-posteriori bounds are called efficient and reliable. In the finite element context S ~ could be a trial space corresponding to a refined mesh or a trial space with the same mesh as S ~ but containing higher order trial functions. In either case the expression u ~ - u" can be evaluated efficiently typically by solving local problems (see, e.g., [7]). The resulting bounds are then usually comprised of sums of local terms whose precise form, however, depends strongly on the particular discretization at hand. In the following we will also study energy estimates or more generally estimates relative to the norm I1" lit but for the above general multiscale basis setting and for the scope of problems described before. Our goal is to develop computable efficient and reliable error bounds without assuming the saturation property. We will show that under relatively mild assumptions on the given data these a-posteriori error estimates lead to an adaptive strategy which can be proved to converge. Making crucial use of the norm equivalences based on wavelet expansions these estimates will be formulated in terms of wavelet coefficients and thus take a rather unified form for a wide range of cases.

4.2. Residual estimates The spaces Sj satisfying estimates of the form (3.7) typically correspond to uniform refinements of underlying meshes. It is well known that such uniform refinements do not provide economical approximations when the approximated object exhibits singularities. Therefore our goal is to determine possibly small subspaces of the form

SA := span {¢~: A E A}, which are adapted to the particular solution u of (2.1). More precisely, given some tolerance eps, we seek a possibly small set A C J such that the Galerkin approximation UA E SA, defined by ( A U A , V) = (f, V),

V E SA,

(4.4)

satisfies

Ilu - UAllt <<.eps. This will be done by successively updating A based on appropriate a-posteriori estimates of a current Galerkin approximation UA. To this end, we follow first standard lines and exploit the fact that, in principle, one can evaluate the residual

rA = AUA -- f = A(UA -- u)

(4.5)

and that, on account of (2.5), we have

clllrzll-t ~< I l u -

UAIIt <<.cmllrAIl-t.

(4.6)

Recall also from Section 2 that the Galerkin scheme is assumed to be stable for the spaces SA (which is for instance the case when A is selfadjoint and positive definite). In terms of the projectors

Q z v := ~ (v, ~A)¢A AEA

S. Dahlkeet al. /Applied NumericalMathematics23 (1997) 21-47

31

this can be conveniently reformulated as

IIQ~AAVAII_t "~ IlvZlit, vz E SA.

(4.7)

Next we will make essential use of the norm equivalences in Theorem 3.1, recalling that under the assumption (3.9), IIrAlI-t can be estimated by weighted sequence norms of the wavelet coefficients of rA. In fact, by (3.10), we have

C3( ~ 2-2t[Al[(ra,¢A)12) 1/2 <~ [[rzll-t <~C4( ~ 2-2tlAll(rz'¢A)[ 2) 1/2 AEJ\A AEJ\A

(4.8)

where we have used that, since UA is a Galerkin solution,

rA = E(rA,~3A)/~A = E (rA'~A)~X" AEJ XEJ\A Thus, combining (4.6) and (4.8), we obtain, in principle, an efficient and reliable error bound. To understand these expressions a little better let us abbreviate the residual weights

and note that inserting the expansion

UA : E ~A'~bA' X'E A of the Galerkin solution UA yields the representation = 2-tl;~l f x -

E

(a¢:e,¢~)ux,,

(4.9)

XtEA where f:~ := (f, ¢>,) are the wavelet coefficients of the right hand side f relative to the dual wavelets ~ . In view of (4.6) and (4.8), the quantities

EJ\A ~

are asymptotically sharp lower and upper

bounds for the error ]]u - UAllt. Obviously, the/2-norm of all the residual weights is of little practical use yet since this involves infinitely many terms. The objective of our following considerations is to replace them by computable expressions which still form within certain controllable tolerances lower and upper bounds for the error ][u - Udllt. To this end, note that the ~x are influenced by two sources, namely (i) by the wavelet coefficients of the right hand side f which reflect its smoothness and approximability, and (ii) by the behavior of the current approximate solution UA and properties of the operator A. It will help to estimate both effects separately. We will begin with the second one.

4.3. A basic estimate We will first show that for almost all A E J the sums ~--~,eA(A¢;~,, ~b~)u:~, can actually be neglected. To this end, it is well known that for a large class of operators A the (A¢~,, CA) decay when either the levels I,Xl, I: '1 or the supports $2~, I2A, of the wavelets ¢~, ~b~,, respectively, are far apart. This

S. Dahlke et al. / Applied Numerical Mathematics 23 (1997) 21-47

32

phenomenon has been carefully studied and quantified in the context of matrix compression for the fast solution of linear systems arising from boundary integral discretizations [6,17,18,33,36]. In the following we shall make crucial use of these results which provide the key to reducing the estimates for the residual to finitely many computable terms. In fact, the decay of the entries (A~:~,, ~:~) is typically governed by the following basic estimate on the entries 2-1lAI-tA'll~ 2-(IA'I+IAI)tI(ACA" CA)] < (1 + 2min(IAl,lA'l)dist(/2A, Y2A,))d+2m'+2t' (4.10) whose validity we shall assume in the following (see Lemma 4.1 below). Here again 2t = p is the order of A, cr > 0 is some fixed real number depending on the regularity of the wavelets and m t is a positive integer which typically represents the order of exactness of S or the order of vanishing moments of the ~A (which have to be suitably interpreted depending on the type of the underlying domain Y2, see, e.g., [19]). In fact, if ~ is a stable multiscale basis for L2(Y2) where Y2 is a domain in R d, ~ is said to have vanishing moments of order m / if

P(x)¢),(x)dx=O,

,~CJ\J0,

(4.11)

holds for all polynomials P of degree less than m I. Since estimates of the form (3.7) usually imply that the approximation spaces contain polynomials of the order corresponding to the highest approximation order and since the ¢~ are orthogonal to SI~I-1 the moment conditions are indeed closely related to the order of exactness of the dual multiresolution sequence ,S. The fact that polynomials appear in (4.1 l) is actually not essential. What matters is that (4.1 l) holds for a finitely dimensional family of functions which locally approximate any smooth function well. Accordingly, this notion can be modified. For instance, when the domain I2 is a manifold represented by a family of smooth parametric mappings ~;i : [] --+ Oi C ~2, where [] is a fixed parameter domain, one can work with the condition

f P(x)¢;~(ni(x)) dx :- O, A e J,

(4.12)

[] for any polynomial P of degree less than m / (see [19,20]). Let us now briefly outline some circumstances under which estimates of the form (4.10) hold. To begin with a simple but instructive case let A -- - A , t = l, H 1 = H01(12), ~ c I~d, and suppose that the V¢)~ have vanishing moments of order m I, i.e.,

P(x).

V¢;~(x) dx = 0,

~2

for any vector valued function whose components are polynomials of degree less than m ~. Then one has for any such polynomial and [)~l > IA~I (A~,;v, ¢;~)=

=[

~< max I V ¢ , v ( x ) - P ( x ) l / I V ¢ ; ~ ( x ) l xE~2), J S2

P¢ ll dx.

S. Dahlke et al. / Applied Numerical Mathematics 23 (1997) 21-47

33

Thus when V~b;~, is still H61der continuous with exponent or' > 0 elementary calculations show that the first factor can be estimated by a constant times 2-1~1°'21~'1(1+~')2dl;¢1/2 while the second factor is bounded by a constant times 2-dlA12dl&l/22 I~1 so that overall one obtains in this case ',

{ 21;~I+IX'I2(d/2÷~')(I;¢I-tXl), if ~ , N ~;~, ~ 1~, < 0, if f2;~ N f2~, = 0.

Obviously this is a special case of (4.10). The above estimate is still crude. In fact, a much stronger decay occurs when I; 1 is much larger than I)~'1 and ~2x does not intersect the singular support of ~ , [33]. Recall from (2.7) that we admit more generally operators with Schwarz kernels of global support. The following result has been established in [17,18] with the aid of relations of the type (4.11) or (4.12) combined with (2.7). It confirms the validity of (4.10) for the range of problems under consideration. L e m m a 4.1. Assume that (2.7) holds and that k~ has vanishing moments of order m' in the sense of (4.11), (4.12), respectively. Moreover, suppose that for r from (2.6) cr > 0 satisfies (see (3.9)) cr <~ r,

t + cr < 7,

t - tr > - 7 ' .

(4.13)

Then for d + 2t + 2m I > O, d/2 + m' + t >~ tr and any )~, ~' C J one has

2-11~l-lX'llo 2-(l~'l+l~l)tl(A~bx, , ~b~)[ ~< (1 + 2min(IXhlX'l)dist(f2x, f2a,)) a+2m'+2t"

(4.14)

Proof. It has been shown in [18,33,36] that under the above assumptions on k~ 2-(l~l+l;¢l)(d/2+m ') I(AWx',~b~,) I < (dist(/2x, ax,))a+2m,+2t

holds when f2x ~ f2a,. Straightforward calculations show that this implies 2-11N-IX'll(cl/2+m' +t) 2-(I;¢l+l;q)t[(A~x" ~bx)l < (1 + 2min(IXl,lX'l)dist(~x, f2;v)) d+2m'+2t'

provided that dist(~2)~, ~2~,) ~> 2 -min(IN'l'VI) and )~, ~' E g \ do. Since we have assumed that d/2 + m' + t ~> ~r we have confirmed (4.14) for ~, )~' E a \ a0 and dist(f2)~, Dx,) > 2 - min(IN,p''[). Now let dist(f2~, ~?x,) < 2 -min(l~l&'l) and suppose without loss of generality that > I '1, i.e., ~ a \ a0, )~' ~ g. Using Schwarz' inequality and the continuity of A (2.6) yields I(A~;~,, ~bx) I ~< IIA

, II-t+ ll ll -

<

I1,+ II llt- .

(4.15)

On account of (4.13), we can apply now the norm equivalence (3.8) to each factor on the right hand side of (4.15) which, upon using biorthogonality, yields

[(AWa,, a)l

2t(1~1+1~'1)2~(I;¢1-1~1)"

This confirms the remaining part of the assertion and completes the proof.

[]

Note that we have exploited above the flexibility of choosing the order m' of vanishing moments sufficiently large which is possible in the biorthogonal setting.

34

S. Dahlke et al. / Applied Numerical Mathematics 23 (1997) 21-47

4.4. Computable bounds Here estimates of the form (4.10) or (4.14) should be viewed as representatives of a certain decay property which suffices to prove the principal fact that the bounds in (4.8) can be reduced to finitely many local quantities. The reasoning is closely related to matrix compression techniques based on estimates of the form (4.14) [18]. Lemma 4.2. Assume that (4.10) holds with cr > d/2. Let 0 < ¢$ < a - d/2 be fixed and choose for E > 0 positive numbers El, E2 > 0 such that E21(m'+t) + 2 -Us~ ~< E.

(4.16)

For ~ E J and e define the influence set

:--

e J: J i l l - I 'll

-
Then there exists a constant c5 depending only on ~, rn ~, t, the constant c3 and the stability constants in (4.7) of the Galerkin scheme such that the quantities ea :=

~ E J\A,

(4.17)

2-2t/'kl[e'kl2)1/2~c5EllQSlfll-t"

(4.18)

Z

(A~b)e,¢a)ua,,

X'EA\J;~,~

satisfy

(Z AEJ\A

Proof. It has been shown in [18,33], that the infinite symmetric matrix R e = (a~,,v)~,,;v~j defined by a~,A, =

2-11;q-I~'ll~ (1 + 2min(lXt'lX'l)dist(/2A,J'~A'))d+2m'+2t'

A E J \ A, At E J \ Jx,e,

O,

otherwise,

satisfies

IIR II

E,

(4.19)

where II • II denotes here the spectral norm and the constant in (4.19) depends only on the constants Cl, c2, c3, c4 and t, a. By definition (4.17) we obtain 2-2tl)'l[e;~l2-- y ~ 2 -2tlAI AEJ\A

~EJ\A

<~ Z XEJ\A

~

(A¢)e, CA)UA' 2

~'EA\Jx,e

(

Z

2-(l'Xl+l'Vl)t](A¢'v,¢'x)l 2tl'vllu.vI)

AtEA\J;~,e

Therefore we infer from (4.14) and (4.19) that

35

S. Dahlke et al. /Applied Numerical Mathematics 23 (1997) 21-47

2-2tlXl levi 2 Z )~EJ\A

Y~

~

)~EJ\A MEA\J~, e

2-11Al-lA'll~ ..... (1 + 2min(I)q'l~'l)dist(f2,X, J'2M))d+2m'+2t

)~'EA\J.x,~

i2

.X'EA\J~,,~

< 2 2 t 2 ,,o g IlUallt ~ S211QtAAUAll2_t = s 2 1 l Q A f [ [ _ t ,

where we have used (4.7) and the fact that (4.20)

QIAAUA = Q'Af in the last step.

[]

R e m a r k 4.2. It is obvious from the above proof that the quantity I[Q'Afll-t in (4.18) can be replaced by either IluAl[t or Ilfll-t with a modified constant c5. The latter choice has the advantage of being independent of A but the disadvantage that in a strict sense it is not computationally accessible while the first two choices can be estimated via the coefficients f),, u;~, A E A and the corresponding norm equivalences. Lemma 4.2 shows that up to a controllable tolerance the portion e;~ in the residual weight (iA can be neglected which suggests replacing for a given e > 0 and a given A C J the quantities ~i;~ in (4.8) by

d;~(A,c) = d)~ := 2 -tIll f;~ -

Z

(A¢;~,,¢~)u~, ,

)~ E J \ A.

(4.21)

,k' EANJx,~

It is clear that also the sums (~-])~EJ\A d2) 1/2 are not finite yet. However, consider the set

NA,e := {)~ E J \ A: A @ J~,e ¢ 0}

(4.22)

which consists of all indices in the complement of A whose influence set intersects A. In this sense Na,e describes the significant margin around A which reflects the smearing by the operator A of the contributions from A relative to the tolerance e. Alternatively one has

NA,~ = U{J)e,~ \ A: .V E A } so that

#NA,~ < c~.

(4.23)

Since by definition J~,~ M A = 0 whenever/~ E J \ (A U NA,e) we infer from (4.21) that d~(A,e) = 2-tl~llf:~ I for )~ E J \ (AUNA,e).

(4.24)

Hence all but finitely many of the terms d),(A, ~) depend only on the right hand side f. This dependence, in turn, can readily be interpreted as an approximation error. In fact, since f E H -t the series E)~EJ 2-2t1~1 Ira[ 2 converges so that }--:~;~EJ\ (NA,eUA)2-2tl)~t If)~[2 can be made arbitrary small by choosing A appropriately. In fact, the contribution of f to (~-~;~EJ\(NA,~UA)d2) 1/2 is just

36

S. Dahlke et al. / Applied Numerical Mathematics 23 (1997) 21-47

IIf -

2-2tI~IIA[ 2

Q~AUNA,,fII2-t "

inf

AEJ\(NA,eUA)

Ilf - vll2-t

vC'SAuNA,e

~< in..f IIf-vJJ2-t •

(4.25)

vESA

In the following we will assume that quantitative estimates of the form

(

2-2elO~llfml2

csIIQSdll_, ~< c~llfll-t

~< c7 in f I I f - v i i - ,

AEJ\A

(4.26)

vCSA

are accessible for any A C J (see Remark 4.2). Such contribution can be subsumed under the perturbation of order e provided that SA admits a sufficiently good approximation to f. This is precisely the problem of adaptively approximating an explicitly given function or distribution. Thus any singularities of f will be reflected by the initial choice of a starting set A. The above estimates of the quantities e~ show how much this information is smeared due to the pseudo-locality of the elliptic operator which becomes accessible through the multiscale representation. If f is very smooth the contribution of (4.25) will be negligible and the adaptive choice of larger sets/1 will be dominated by the behavior of the current approximation UA and the action of A. Lemma 4.2 and (4.25) suggest to define the quantities

a;~(A, e) = a>, := 2 -tIll

(Ag,;~,, ~:~)u:~, ,

)~ E J \ A,

(4.27)

A' CAMJ>,,e

noting that, in view of (4.22), a:~(A,e)=O

.k~NA,~.

for,~EJ\A,

(4.28)

We will show next that these quantities still provide sufficiently accurate error bounds.

4.5. A-posteriori error bounds We will show next that the quantities a;~(A, e) still give rise to a new a-posteriori estimate which is up to any chosen tolerance efficient and reliable. For the special case of second-order two-point boundary value problems a similar result was obtained by S. Bertoluzza who was as far as we know the first to establish an a-posteriori estimate of the following type [5]. T h e o r e m 4.1. Under the assumptions in Lemma 4.2 one has

~

Ilu--uAIIt <<.C2C4

a:4A,

)=

+c'sEllfll-t+c7 in f I I f - v l l - t

)~CNA,E

,

(4.29)

VCSA

as well as ( ~-~ ax(A,e)e) 1/2 ,~EN A,e

1 C1C3

Ilu - UAIIt + c

llfll-t + c7 inf IIf - v[l-t. vE,-S~A

(4.30)

S. Dahlkeet al. / AppliedNumericalMathematics23 (1997)21--47

37

Proof. Both inequalities follow from (4.8), the triangle inequality, and Lemma 4.2. In fact, we apply (4.17), (4.21) and (4.27) to obtain 3;~ <~ d;~(A, e) +

2-tl~'llem I ~< a:~(A,e) + 2-tl:'llem I +

2-tl;qlf;~ I.

Inserting this into (4.8) gives an upper bound for the residual lira I1-~ which, in view of (4.6), provides

I,U--UAllt ~C2C4{( Z a)~(A.c)2)l/2+ ( Z 2-2tlAIleAI2)I/2 ACJ\A AEJ\A

Since by (4.27), Y'~J\A a)~(A, e) 2 = ~-]~NA.~ ax(A, e) 2, (4.26) and (4.18) in Lemma 4.2 yield the estimate (4.29). Likewise, because again by triangle inequality, also a), (.4, e) <<.5),+2-tl)'lle),l+2-tl)qlf), I, employing the lower bounds in (4.6) and (4.8), a completely analogous reasoning provides (4.30). [] Recall that, given A, our basic strategy is to find a possibly small set A C J containing A so that the error of the Galerkin approximation is reduced by some fixed factor. To this end the following observation is useful. Lemma 4.3. Suppose that A c A C J and let u S be the Galerkin solution with respect to S S = span{,;~: A E A}.

Then we have ( Z

aA(A'~)2) 1/2

1 CI C3

A~AnNA,~

IluS

--

UAllt q- c~llfil-t + c7 inf IIf - vll-c.

vESA

Proof. For A E A we have

(AO;~,, ¢:~)u:~,--- (Auz, ~)~) = (A(uA -- uS), z~)~) + f~, MEA

and therefore, by (4.17) and (4.21), d~(A,z) = 2 -tl;~l f:~ - Z

(A¢~,,~,)u~, + e~ < 2-tI~'II(A(uA - uff),¢:0 1+ 2-*l~lle~ I.

MEA

Because of (4.6) and (4.8) we obtain

y ~ 2-2tl;~ll(A(uA- us),¢),) I2 ~< c3211A(uz

-

us/ll2_

)~ES\A

Thus (4.18) provides

d~(A,¢) 2

A\A

<~ - - I l u S - uAIl~ + 5¢IIQA/II_,.

C3C1

1

[[UA--us]l 2.

(4.31)

38

S. Dahlke et al. /Applied Numerical Mathematics 23 (1997) 21-47

Since by definition (4.27) of a~(A, e) one has

lax(a,e)l <<.Id~(a,e)l + 2-tl~llf;~ ] the assertion follows from (4.26).

[]

Deriving sharp lower and upper a-posteriori bounds for operators with negative order (cf. [24,31,37]) is considered to be a difficult task. To our knowledge no efficient and reliable estimators for arbitrarily refined meshes have so far been known. For example, Carstensen and Stephan proved in [10] aposteriori estimates without deriving lower bounds. The results in [9] require quasi-uniform meshes. It is remarkable, that within our approach operators with negative order seem to have to some extent even advantages over operators with positive order, because of their smoothing property, i.e., (A~2~,, ~ ) decays faster for I.XI--+ o~ if A~b;~, is smoother.

4.6. A convergent adaptive strategy The next step is to use the a-posteriori error estimates for an adaptive refinement strategy. Although totally different in a technical sense and with regard to the whole setting the results are similar in spirit to those by D0rfler [23] who considers adaptive refinement of piecewise linear finite elements for Poisson's equation in two dimensions. In the present setting refinement simply means adding properly selected basis functions ~bA to the current solution space. We will describe assumptions, that guarantee an improvement for the approximate solution after the refinement step without assuming the saturation property. Our next goal is to use (4.31) for selecting a set A containing A such that the saturation property is automatically guaranteed to hold. So far all our previous estimates did not require any symmetry assumptions on the operator A. For the next step, however, it seems that more information about A is needed. The simplest setting would again be the symmetric case, i.e., a(.,-) := (A., .)

(4.32)

defines a symmetric bilinear form such that II" ]I := a(., .)1/2 satisfies

c8[I. It

II-[It

c911. II.

(4.33) We will formulate and prove the next result for this simple setting and will indicate later how to weaken the assumptions somewhat. Theorem 4.2. Suppose that (4.32) and (4.33) hold and let eps > 0 be a given tolerance. Fix any 8* E (0, 1) and define l-O* • ClC3 2C2C4~] Choose any #* > 0 such that 1 O* Iz*Ce <~ 2(2 0")¢2c 4 and set #*eps e .--

Ce :=

(

--

1

+

(4.34)

(4.35)

-

2dsIIf[[-t"

(4.36)

39

S. Dahlke et al. /Applied Numerical Mathematics 23 (1997) 21-47

Suppose that A c J is chosen so that 1

,

c7 inf [ I f - v[[-t < ~# eps.

(4.37)

VCSA

Then whenever A C J, A C A is chosen so that

( z AEAANA,s

A

,c

there exists a constant ~ E (0, 1), depending only on the constant 0", the constants in (4.33) and the constants ci, i = 1,..., 4 such that either

(4.39)

I l u - ~,xll < ~ l l u - uatl or

(

)1/2(

Z

a)~(A'e)2

\1/2

=

AENA,e

E


a.x(A,e) 2)

AEJ\A

Proof. We first assume, that

Ilu- UAllt ~> ep_._ss c~

(4.40)

where the constant C~ > 0 is defined by (4.34). When A satisfies (4.38) we infer from (4.31) and (4.29) that Z aA(A'e)2 -c'5~llfll-~-c7inf Ilf-vll-~ ~nN~,~ VC~A inf IIf-vll-t} ClC 3 ( (1 - 0") {1 C2C--' ~ - Ilu--umllt--C'sellfll-t--C7 VCSA

[[uX--UAl[t)clc3

- c ~ l l f l l - t - c7 in f ClC 3

1o, k, C2C4

vESA

Ilf vll-t~/ -

II~ -- ZtAllt - (2 -

O*)c'se[lfll-t-

Ilf - vll-t vCSA

(2 - 0")c7 inf

)

•

Employing (4.36) and (4.37), yields I l u ~ - uz[It

~ clc3(~Hu-

uzNt - (2 -O*)p* eps).

Moreover, taking (4.40) into account, we conclude that [[U~ -- UAIlt ~ ClC3

(1-0"

---C2C4

(2 -

o*)~*c~

>/clc3(1 - 0") i[u _ uzllt, 2c2c4

)

Ilu -

UAllt

(4.41)

40

S. Dahlke et al. / Applied Numerical Mathematics 23 (1997) 21--47

where we have used that by (4.35), C2%7

\ C2C4

"2-C22-C44

-2~22-C44"

By (4.33) and (4.41), we obtain

I1~*~- uzll >1

CLC3C8(1-- 0") Ilu - uall.

(4.42)

2c2c4c 9

At this point we exploit symmetry by applying Remark 4.1 which confirms the assertion (4.39) with t~ =

W~

- ~,

where c6

On the other hand, (c~

Ilu \1/2

a~(A,~.)2)

CLC3C8(1 - 0") :=

2c2c4c9

UAl[t< eps/Ce yields, in view of (4.30), (4.36) and (4.37) ~ _ _1 Ilu - UAllt + c~ellfll-t + C7 in f Ill -- vii-, Cl C3

vESA

~ < - -eps +/,t*eps = cle3C~

(1/(CLC3)+>*Ce)eps C~

Taking (4.34) and (4.35) into account, we see that

(~/(clc3) + t**c~)
(Z

a)~(A'c)2) 1/2 < eps

)~ENA.e

which completes the proof. Note that (~-~ENA, e

[]

a~(A,g)2) 1/2 < eps yields, by (4.29),

(4.36) and (4.37) (4.43)

Ilu - uallt <~c2c4( 1 q- #*)eps. R e m a r k 4.3. By Remark 4.2, the term

Ilfll-t

in (4.36) can be replaced by

IluAllt

or

( ~EA22'AItIuA'2)1/2 SO that (4.36) can be replaced by >*eps

2zlaltluale)'/2"

C, (~-'~AEA

The constants change then in an obvious way. The above result may be formulated in terms of the following algorithm.

(4.44)

S. Dahlke et al. / Applied Numerical Mathematics 23 (1997) 21--47

41

Assumptions: We assume, that the constants Cl, ce, c3, ca, c~, c7, c8, co or estimates for these constants are known.

Initialization: Fix 0* E (0, 1) and the desired accuracy eps. Compute Ce, #* according to (4.34), (4.35).

Algorithm Step 1: Compute e according to (4.36). Step 2: Determine an index set A c J such that c7 inf IIf -

vN-t <

½t

*eps-

vESA

Step 3: Compute the Galerkin solution UA with respect to A. Step 4: Compute

( )1/2~

a),(A,e) 2

:=

(4.45)

._-. .kENA,e

If r]A,e < eps Stop, accept ua as solution which satisfies (4.43). Otherwise, go to Step 5. Step 5: Determine an index set A, A c A c J such that Z

a;~(A'e)2

~> (1--O*)?']A,e.

Set A - + A and go to Step 3. R e m a r k 4.4. Again we could have used the computable quantities IIQ'afll-t or IlUAlit in the definition of e instead of Nfll-t which would require an additional evaluation of these terms in each step as well as a possible change of e. Shooting directly for the final desired accuracy might require starting in the above algorithm with a rather refined set A. To better balance the influence of f and A it may therefore be preferable to view the above algorithm as one loop in a scheme of the following type where we assume again the respective initializations of the above algorithm. (I) Choose the final desired accuracy e p s / a n d some initial eps 0 > 0 (which could be much larger than epsf). Put eps = eps 0. (II) Apply the algorithm with eps. (III) If eps ~< epsf Stop, accept UA as the approximate solution. Otherwise set eps/2 --+ eps and go to (II).

4.7. Some comments on the role of symmetry Let us point out finally that symmetry is a convenient but not quite necessary assumption in the above context. We will briefly indicate one way of weakening this hypotheses. In many cases there is a symmetric Ht-elliptic bilinear form s(-, .) such that the difference

k(u, v) := a(u, v) - s(u, v),

u, v E H t,

(4.46)

42

S. Dahlke et al. / Applied Numerical Mathematics 23 (1997) 21-47

is a compact bilinear form, which can be estimated by a weaker norm than 11-lit. Under this assumption the bilinear form s(., .) defines a norm

Ilull := ~u),

u C n t,

which is equivalent to

I1 lit,

i.e.,

Ilull ~ Ilullt.

(4.47)

Since for A C J, A c

Ilu- uall 2 -Ilu-

usII = -[[ua

--

usII 2

= 8 ( U -- UA, U -- U A ) -- 8 ( U -- U S , U -- ZtS) -- 8 ( U A -- U S , UA -- U S ) -~ 28 (U -- U S , U S -- U A ) ,

the Galerkin orthogonality a ( u -- u S , UA -- u S )

= O,

provides - u ~ , UA -- uS),

s ( u - u s , UA -- u ~ ) = - k ( u

(4.48)

so that

Ilu - usII 2 = Ilu -

UAll 2 -

IlUA - us]]2 + 2 k ( u - u S , u A - u s ) .

(4.49)

Now we assume further that ]k(zt-

~(llu- usII 2 + Ilu-

Iz~, UA -- U S ) I ~

uzll 2)

(4.50)

holds for some ~7 > 0. We emphasize that under the given assumptions on the bilinear form k(., .) this is not an assumption on the particular solution u but rather on the approximation properties of the spaces S A which should imply that the error u - UA measured in a weaker norm induced by k eventually is dominated by the error in the energy norm. Thus in addition to the right hand side f this may impose an additional condition on the choice of the initial index sets A (see (4.37)) namely that A contains sufficiently many low frequencies. This point will be illustrated by the example discussed below. As in the proof of Theorem 4.2 we conclude from (4.41) and (4.47) that

Ilu~- UAll 2 ~

~'211u -

uall 2

for some n¢ E (0, 1), which, on account of (4.49) and (4.50), yields /

Ilu - usII ~< ~l

/~t2 + 2r/ 1 - 27

Ilu - u a l l .

(4.51)

If r/is small enough one obtains IIu - usII ~< ~"llu - um[I for some ~;~t C (0, 1), i.e., we can again confirm the validity of the saturation property with respect to A and .~.

S. Dahlke et a l . / Applied Numerical Mathematics 23 (1997) 21-47

43

Let us consider next a typical example, namely a simple boundary value problem on a bounded Lipschitz domain ~ C R a with boundary F = ~ ~2 involving a symmetric elliptic second-order partial differential operator and a first-order term destroying symmetry. Let x = ( x l , . . . , Xd) E ~d and A the partial differential operator defined by A=-A+~x~,

~>0.

The weak formulation of the homogeneous Dirichlet problem Au(x) = f(x),

x e D,

ul0n --0, for a given function f E H-l(S2) reads a ( u , v ) := (Vu, V v ) + / 3 ( S x , u , v ) = ( f , v ) ,

v • HI(O ).

The bilinear form s ( u , v ) := (Vu, Vv),

u , v • g~(f2),

is symmetric and Ho1(~2)-elliptic. Furthermore we have

k(u,v) := Z(Ox,U,v). Since Ik( u - uX, uA -- us)l ~
II(A')-'vll,÷s

< c(~s)llvll_,+s,

v • H-'+s(~),

holds for some s • (0, 1]. Employing the Aubin-Nitsche trick we infer from the estimate

[l(I-Qz)vlll

<~ ~'llvll~÷~,

v • H1+~(~2),

(4.53)

for some ~?t > 0, that Ilu - uall~-~ <~ ~'c(~, s)llu - uzlll as well as

I1~ - ~ l l ~ - s ~< ~'c(~, sDIlu - usIl~.

(4.54) (4.55)

Thus, by (4.52), (4.54) and (4.55), we obtain with ~7 := (3/2)0'3c(3, s)

usIl~

[k(u - u~,/z z -/z~') ] ~ ~(llu + Ilu - uallff). Of course, when {¢:~: I~1 < j} c_ A standard estimates (see, e.g., [14,18]) show that r/' in (4.53) behaves like 2 -sj. Thus, in principle, how ever large fl may be, including uniformly refined spaces up to a sufficiently high level of discretization, one can always ensure the validity of the estimate (4.50). In fact, on a sufficiently high scale the discrete operator is elliptic. However, for large ~ this may exclude the analysis of adaptive refinements for that range of discretization levels which are of practical relevance. Similar arguments apply for boundary integral equations and even for systems of boundary integral equations with inherent asymmetries. For more details about such examples we refer to [1 t,28].

44

S. Dahlke et al. / Applied Numerical Mathematics 23 (1997) 21-47

5. Solving the linear systems So far we have only discussed the adaptive selection of subspaces which permit approximate solutions of sufficient accuracy at the expense of possibly few degrees of freedom. On the other hand, the adaptive choice of subspaces is closely related to the efficient solution of corresponding linear systems. To describe this, let

a a := ((A~bA',~bA))~y~A denote the stiffness matrix relative to the wavelet basis. By the same arguments as used in [18] for the full spaces it follows from (4.7) that also cond2 (DAtAADA t) ~ 1,

(5.1)

where (D~)~,~, := 2sl~l~,~,. In fact, defining 2~l;q (v, ¢:0 ¢;~, ACJ

one clearly has, in view of (3.4),

SsQA = QAZs,

~v,~-I = ~ ' - s .

(5.2)

Thus Theorem 3.1 says that

IISsvll Ilvlls+ , (--y',-y). Setting WA := ZtVA for VA E SA, (5.3), (5.2), and (4.7) imply IIWAII0 ,'~ lira[it ~

IlQ'Aavall_t ~

!

(5.3)

!

IlZ_tQaZQzZ_twzllo,

which means that the operator Z~_tQ~AAQAZ_t and their inverses are uniformly bounded on L2($2). t t m QAS-t which Moreover, it is easy to see that DAtAADA t is the matrix representation of S_tQ A confirms (5.1). The above reasoning also shows that the operators

A~\ z := QA\AAQX\A are for A C A C J well-conditioned as long as the difference in levels IAI occurring in A \ A remains bounded. Now suppose that we have determined an approximate solution UA E SA and A is the next updated index set generated by the above algorithm. Note that this requires determining the entries of AX\ A. The above comments suggest now to solve the problem

(Aux, v ) = (f,v),

y e s S,

iteratively using the starting value vX

:=

VA~)A~

~eX where VA : =

{

u;~, i f A E A , w;~, ifA C A \ A ,

S. Dahlke et al. /Applied Numerical Mathematics 23 (1997) 21-47

45

and

WA\ A -- ~ WA~.k, AEA\A

UA = ~ UA~A, AEA

are the solutions to

QA\AAWS\A : QA\A f,

QtAAUA = Q~Af,

respectively. As pointed out above w s \ A can be approximated well by a few iterations on a (relatively small) system essentially without spending any effort on preconditioning. Finally, the matrix compression techniques studied in [18,33] readily apply to the matrix representations of the operators AA when A has a global Schwarz kernel.

6. Conclusions Exploiting properties of stable multiscale bases, especially norm equivalences for certain ranges of Sobolev norms, we have shown that certain a-posteriori error estimators in terms of wavelet coefficients of residuals are efficient and reliable. Moreover, they give rise to adaptive schemes which are guaranteed to converge for a wide range of elliptic operator equations including those of negative order without assuming a priori the validity of the saturation property. To our knowledge these are the first results of this type for the latter class of problems.

References [1] I. Babu~ka and A. Miller, A feedback finite element method with a posteriori error estimation: Part I. The finite element method and some basic properties of the a posteriori error estimator, Comput. Methods Appl. Mech. Engrg. 61 (1987) 1-40. [2] I. Babugka and W.C. Rheinboldt, Error estimates for adaptive finite element computations, SIAM J. Numer. Anal. 15 (1978) 736-754. [3] R.E. Bank and A. Weiser, Some a posteriori error estimates for elliptic partial differential equations, Math. Comp. 44 (1985) 283-301. [4] E Bastian and G. Wittum, Adaptive multigrid methods: The UG concept, in: Adaptive Methods--Algorithms, Theory and Applications, Notes on Numerical Fluid Mechanics 46 (Vieweg, Braunschweig, 1994) 17-37. [5] S. Bertoluzza, A posteriori error estimates for the wavelet Galerkin method, AppL Math. Lett. 8 (1995). [6] G. Beylkin, R. Coifman and V. Rokhlin, The fast wavelet transform and numerical algorithms, Comm. Pure AppL Math. 44 (1991) 141-183. [7] E Bomemann, B. Erdmann and R. Kornhuber, A posteriori error estimates for elliptic problems in two and three space dimensions, SIAM J. Numer. Anal. 33 (1996) 1188-1204. [8] J.M. Carnicer, W. Dahmen and J.M. Pefia, Local decomposition of refinable spaces and wavelets, Appl. Comput. Harmonic Anal. 3 (1996) 127-153. [9] C. Carstensen, Efficiency of a posteriori BEM-error estimates for first-kind integral equations on quasiuniform meshes, Math. Comp. 65 (1996) 69-84. [ 10] C. Carstensen and E.E Stephan, Adaptive boundary element methods, Preprint, Universitat Hannover (1993). [l l] M. Costabel, Boundary integral operators on Lipschitz domains: elementary results, SIAM J. Math. Anal. 19 (1988) 613-626.

46

S. Dahlke et al. / Applied Numerical Mathematics 23 (1997) 21-47

[12] M. Costabel and E.E Stephan, A direct boundary integral equation method for transmission problems, J. Math. Appl. (1985) 367--413. [13] W. Dahmen, Some remarks on multiscale transformations, stability and biorthogonality, in: P.J. Laurent, A. Le Mthaut6 and L.L. Schumaker, eds., Curves and Surfaces II (AKPeters Wellesley, Boston, 1994) 157-188. [14] W. Dahmen, Stability of multiscale transformations, J. Fourier Anal. Appl. 2 (1996) 341-361. [15] W. Dahmen, A. Kunoth and K. Urban, A wavelet-Galerkin method for the Stokes problem, Computing 56 (1996) 259-302. [16] W. Dahmen, A. Kunoth and K. Urban, Biorthogonal spline wavelets on the interval---Stability and moment conditions, Preprint, RWTH Aachen (1996). [ 17] W. Dahmen, S. PrOssdorf and R. Schneider, Wavelet approximation methods for pseudodifferential equations II: Matrix compression and fast solution, Adv. Comput. Math. 1 (1993) 259-335. [ 18] W. Dahmen, S. Prtissdorf and R. Schneider, Multiscale methods for pseudo-differential equations on smooth manifolds, in: C.K. Chui, L. Montefusco and L. Puccio, eds., Proceedings of the International Conference on Wavelets: Theory, Algorithms, and Applications (Academic Press, 1994) 385-424. [19] W. Dahmen and R. Schneider, Wavelets on manifolds, in preparation. [20] W. Dahmen and R. Schneider, Multiscale methods for boundary integral equations: Adaptive quadrature, in preparation. [21] R.A. DeVore, Degree of nonlinear approximation, in: C.K. Chui, L.L. Schumaker and J.D. Ward, eds., Approximation Theory VI (Academic Press, 1989) 175-201. [22] R.A. DeVore, B. Jawerth and V. Popov, Compression of wavelet decompositions, Amer. J. Math. 114 (1992) 737-785. [23] W. Dtrfler, A convergent adaptive algorithm for Poisson's equation, SlAM J. Numer. Anal. 33 (1996) 1106-1124. [24] B. Faermann, Lokale a-posteriori-Fehlersch~itzer bei der Diskretisierung von Randintegralgleichungen, Dissertation, Christian-Albrechts-Universi~'t Kiel, Germany (1993). [25] M. Griebel and P. Oswald, Tensor product type subspace splittings and multilevel iterative methods for anisotropic problems, Adv. Comput. Math. 4 (1995) 171-206. [26] P. Hansbo and C. Johnson, Adaptive finite element methods in computational mechanics, Comput. Methods Appl. Mech. Engrg. 101 (1992) 143-181. [27] S. Hildebrandt and E. Wienholtz, Constructive proofs of representation theorems in separable Hilbert spaces, Comm. Pure Appl. Math. 17 (1964) 369-373. [28] R. Hochmuth, A posteriori estimates and adaptive schemes for transmission problems, Preprint, RWTH Aachen (1996). [29] P. Mund and E.P. Stephan, Adaptive multilevel boundary element methods for the single layer potential in II~3, Preprint, Universit~t Hannover (1995). [30] P. Oswald, Multilevel Finite Element Approximation, Teubner Skripten zur Numerik (Teubner, Stuttgart, 1994). [31] E. Rank, Adaptive h-, p- and hp-versions for boundary element methods, Internat. J. Numer. Methods Engrg. 28 (1989) 1335-1349. [32] U. RUde, On the V-cycle of the fully adaptive multigrid method, in: Adaptive Methods--Algorithms, Theory and Applications, Notes on Numerical Fluid Mechanics 46 (Vieweg, Braunschweig, 1994) 251-260. [33] R. Schneider, Multiskalen- und Wavelet-Matrixkompression: Analysisbasierte Methoden zur Lt~sung groBer vollbesetzter Gleichungssysteme, Habilitationsschrift, TH Darmstadt, Germany (1995). [34] K. Urban, On divergence-free wavelets, Adv. Comput. Math. 4 (1995) 51-82. [35] R. Verfiihrt, A posteriori error estimation and adaptive mesh refinement techniques, J. Comput. Appl. Math. 50 (1994) 67-83.

S. Dahlke et al. / Applied Numerical Mathematics 23 (1997) 21-47

47

[36] T. von Petersdorff and C. Schwab, Fully discrete multiscale Galerkin BEM, Research Rept. No. 95-08, Seminar Angew. Mathematik, ETH ZUrich (1995). [37] W.L. Wendland and D. Yu, Adaptive BEM for strongly elliptic integral equations, Numer. Math. 53 (1988) 539-558. [38] J. Wloka, Partial Differential Equations (Cambridge University Press, Cambridge, 1987).