ICM 2002 · Vol. I · 607–620
arXiv:math/0212414v1 [math.NA] 1 Dec 2002
Adaptive Methods for PDE’s Wavelets or Mesh Re?nement?
Albert Cohen?
Abstract Adaptive mesh re?nement techniques are nowadays an established and powerful tool for the numerical discretization of PDE’s. In recent years, wavelet bases have been proposed as an alternative to these techniques. The main motivation for the use of such bases in this context is their good performances in data compression and the approximation theoretic foundations which allow to analyze and optimize these performances. We shall discuss these theoretical foundations, as well as one of the approaches which has been followed in developing e?cient adaptive wavelet solvers. We shall also discuss the similarities and di?erences between wavelet methods and adaptive mesh re?nement. 2000 Mathematics Subject Classi?cation: 65N50, 41A25, 41A46, 42C40. Keywords and Phrases: Adaptivity, Mesh re?nement, Wavelets, Multiscale, Methods, Nonlinear approxmimation.
1. Introduction
Among those relevant phenomenons which are modelled by partial di?erential or integral equations, countless are the instances where the mathematical solutions exhibit singularities. Perhaps the most classical examples are elliptic equations on domains with re-entrant corners, or nonlinear hyperbolic systems of conservation laws. While such singularities are sources of obvious theoretical di?culties— classical solutions should be abandonned to the pro?t of weak solutions—they are also an obstacle to the convergence of numerical approximation methods, in the sense that they deteriorate the rate of decay of the error with respect to the size of the discrete problem : achieving a prescribed accuracy will typically require ?ner resolution and therefore heavier computational cost and memory storage, in comparison to the approximation of smooth solutions. Let us remark that singularities
? Laboratoire Jacques-Louis Lions, Universit? Pierre et Marie Curie, Paris, France. E-mail: e cohen@ann.jussieu.fr
608
Albert Cohen
often have a physical relevance : they represent the concentration of stress in elasticity, boundary layers in viscous ?uid ?ows, shock waves in gas dynamics... It is therefore a legitimous requirement that they should be accurately resolved by the numerical method. In this context, the use of adaptive methods, appears as a natural solution to improve the approximation at a reasonable computational cost. Here, the word adaptivity has a twofold meaning : (i) the discretization is allowed to be re?ned only locally, in particular near the singularities of the solution, and (ii) the resolution algorithm uses information gained during a given stage of the computation in order to derive a new re?ned discretization for the next stage. The most typical example is adaptive mesh re?nement based on a-posteriori error estimates in the ?nite element context. While these methods have proved to be computationally successful, the theory describing their advantages over their non-adaptive counterpart is far from being complete. In particular, the rate of convergence of the adaptive algorithm, which describes the trade-o? between the accuracy and complexity of the approximation, is not clearly understood. In recent years, wavelet bases have been proposed as an alternative to adaptive mesh re?nement, motivated by their good performances in data (more speci?cally image) compression. In wavelet-based adaptive schemes, the set of basis functions which describe the approximate solution is updated at each stage of the computation. Intuitively, the selection of the appropriate basis functions plays a similar role as the selection of the mesh points in adaptive ?nite element methods, and one could therefore expect similar performances from both approaches. On a more rigorous level, a speci?c feature of the wavelet approach is the emergence of a sound theoretical setting which allows to tackle fundational questions such as the rate of convergence of the adaptive method. The goal of this paper is to give some elements of comparison between adaptive wavelet and mesh re?nement methods from this perspective. We shall ?rst describe in §2 a general setting which leads us in §3 to a ?rst comparison between wavelets and adaptive ?nite elements from the point of view of approximation theory. We discuss in §4 the relation between these results and adaptive algorithms for PDE’s. After recalling in §5 the classical approach in the ?nite element context, we present in §6 an adaptive wavelet strategy which has been applied to various problems, and discuss its fundational speci?cities. Finally, we shall conclude in §7 by pointing out some intrinsic shortcoming of wavelet-based adaptive methods.
2. A general framework
Approximation theory is the branch of mathematics which studies the process of approximating general functions by simple functions such as polynomials, ?nite elements or Fourier series. It plays therefore a central role in the accuracy analysis of numerical methods. Numerous problems of approximation theory have in common the following general setting : we are given a family of subspaces (SN )N ≥0 of a
Adaptive Methods for PDE’s Wavelets or Mesh Re?nement? normed space X, and for f ∈ X, we consider the best approximation error σN (f ) := inf
g∈SN
609
f ?g
X.
(1)
Typically, N represents the number of parameters which are needed to describe an element in SN , and in most cases of interest, σN (f ) goes to zero as this number tends to in?nity. If in addition σN (f ) ≤ CN ?s for some s > 0, we say that f is approximated at rate s. Given such a setting, the central problem of approximation theory is to characterize by some analytic (typically smoothness) condition those functions f which are approximated at some prescribed rate s > 0. Another important problem is how to design simple approximation procedures f → fN ∈ ΣN which avoid solving the minimization problem (1), while remaining near optimal in the sense that f ? fN
X
≤ CσN (f ),
(2)
for some constant C independent of N and f . As an example, consider approximation by ?nite element spaces Vh de?ned from regular conforming partitions Th of a domain ? ? I d into simplices with R uniform mesh size h. The approximation theory for such spaces is quite classical, see e.g. [12], and can be summarized in the following way. If W t,p denotes the classical Sobolev space, consisting of those functions f ∈ Lp such that Dα f ∈ Lp for |α| ≤ t, we typically have f ∈ W t+r,p ? inf
g∈Vh
f ?g
W t,p
≤ Chr
(3)
provided that Vh is contained in W t,p and that Vh has approximation order larger than t + r, i.e. contains all polynomials of degree strictly less than t + r. Such classical results also hold for fractional smoothness. We can express them in terms of the number of parameters, remarking that N := dim(Vh ) ? h?d , so that if we set X = W t,p and SN := Vh with h := N ?1/d , we have obtained f ∈ W t+r,p ? σN (f ) ≤ CN ?r/d . (4)
We have thus identi?ed an analytic condition which ensures the rate s = r/d. Note that this is not a characterization (we only have an implication), yet a deeper analysis shows that an “if and only if” result holds if we slightly modify the notion of Sobolev smoothness (using Besov classes, see [13]). In summary, the rate of approximation in W s,p is governed by the approximation order of the Vh spaces, the dimension d and the level of smoothness of f measured in Lp . Let us ?nally remark that near-optimal approximation procedures can be obtained if we can ?nd a sequence of ?nite element projectors PN : X → SN such that PN X→X ≤ K with K independent of N : in this case, we simply take fN = PN f and remark that f ? fN X ≤ (1 + K)σN (f ). In the following we shall address the same questions in the cases of adaptive ?nite element and wavelet approximation. As we shall see, a speci?c feature to such cases is that the spaces ΣN are not linear vector spaces.
610
Albert Cohen
3. Adaptive ?nite elements and wavelets
In the adaptive ?nite element setting, the number of parameters N is proportional to the number of triangles, but for a given budget N the partition T and the ?nite element space VT are allowed to be locally re?ned in a way which depends on the function f to be approximated. It is therefore natural to de?ne the approximation spaces SN as SN := ∪#(T )≤N VT . (5)
It should be well understood that the SN are not linear vector spaces (the sum of two elements does not in general fall in SN when their triangulation do not match) but any g ∈ SN is still described by O(N ) parameters, which encode both its triangulation T and its coordinates in VT . The requirement of adaptivity has thus led us to the concept of nonlinear approximation. Wavelet bases o?er another track toward nonlinear adaptive approximation. The simplest prototype of a wavelet basis is the Haar system. Let us describe this system in the case of expanding a function f de?ned on [0, 1] : the ?rst component in this expansion is simply the average of f , i.e. the orthogonal projection f, e0 e0 onto the function e0 = χ[0,1] . The approximation is then re?ned into the average of f on the two half intervals of equal size. This re?nement amounts in adding the orthogonal projection f, e1 e1 onto the function e1 = χ[0,1/2] ? χ[1/2,1] . Iterating this re?nement process, we see that the next components have the same form as e1 up to a change of scale : at re?nement level j, we are adding the projection onto the functions ψj,k (x) = 2j/2 ψ(2j x ? k), k = 0, · · · , 2j ? 1, (6) where ψ = e1 . Since all these functions are orthogonal to the previous ones, letting j go to +∞, we obtain the expansion of f into an orthonormal basis of L2 ([0, 1]) f=
λ∈?
fλ ψλ ,
(7)
with fλ := f, ψλ . In the above notation λ concatenates the scale and space parameters j and k, and ? is the set of all indices (including also the ?rst function e0 ). In order to keep track of the scale j corresponding to an index λ = (j, k) we shall use the notation |λ| = j. More general wavelet systems in one or several space dimension are built from similar nested approximation processes, involving e.g. spline functions or ?nite elements in place of piecewise constant functions (see [23] or [13] for a general presentation). This brief description suggests that a natural construction of adaptive wavelet approximations is obtained by using only a limited set of indices λ as the scale |λ| grows, which depends on the function to be approximated and typically corresponds to those wavelets whose supports are close to its singularities. It is therefore natural to de?ne the approximation spaces SN as the set of all N terms combinations SN := {
λ∈Λ
dλ ψλ ; #(Λ) ≤ N }.
(8)
Adaptive Methods for PDE’s Wavelets or Mesh Re?nement?
611
Again this is obviously not a linear space, since we allow to approximate a function by choosing the best N terms which di?er from one function to another. Note that we still do have SN + SN = S2N . Both adaptive ?nite element and wavelet framework have obvious similiarities. However, the answer to the two basic questions raised in the previous section—what are the properties of f which govern the decay of σN (f ) and how to compute in a simple way a near optimal approximation of f in ΣN —is only fully understood in the wavelet framework. Concerning the ?rst question, a striking result by DeVore and his collaborators [24] is the following : with X := W t,p , best N -term wavelet approximation satis?es f ∈ W t+r,q ? σN (f ) ≤ CN ?r/d , (9)
with q and r connected by the relation 1/q = 1/p + r/d, assuming that the multiresolution approximation spaces associated to the wavelet basis are in W t,p and contain the polynomials of degree strictly less than t + r. Such an estimate should be compared with the linear estimate (4) : the same convergence rate is governed by a much weaker smoothness assumption on f since q < p (as in the linear case, an “i? and only if” result can be obtained up to slight technical modi?cations in the statement of (9)). This result gives a precise mathematical meaning to the spatial adaptation properties of best N -term wavelet approximation: a function f having isolated discontinuity, has usually a smaller amount of smoothness s + t when measured in Lp than when measured in Lq with 1/q = 1/p+t/d, and therefore σN (f ) might decrease signi?cantly faster than εN (f ). The answer to the second question is given by a result due to Temlyakov : if f = λ∈? dλ ψλ , and if we measure the approximation error in X = W t,p , a near optimal strategy when 1 < p < ∞ consists in the thresholding procedure which retains the N largest contributions dλ ψλ X : if ΛN is the corresponding set of indices, one can prove that there exists C > 0 independent of N and f such that f?
λ∈ΛN
dλ ψλ
X
≤ CσN (f ).
(10)
This fact is obvious when X = L2 using the orthonormal basis property. It is a remarkable property of wavelet bases that it also holds for more general function spaces. In summary, thresholding plays for best N -term wavelet approximation an analogous role as projection for linear ?nite element approximation. In the adaptive ?nite element framework, a similar theory is far from being complete. Partial answers to the basic questions are available if one chooses to consider adaptive partitions with shape constraints in terms of a uniform bound on the aspect ratio of the elements
K∈T
max [Diam(K)]d /vol(K) ≤ C.
(11)
Such a restriction means that the local re?nement is isotropic, in a similar way to wavelets. In such a case, we therefore expect a rate of approximation similar to (9). Such a result is not available, yet the following can be proved [13] for Lagrange
612
Albert Cohen
?nite elements of degree m : if for any given tolerance ε > 0, one is able to build a partition T = T (ε) of cardinality N = N (ε) such that on each K ∈ T the local error of approximation by polynomials satis?es ε ≤ inf f ? p 2 p∈Πm
W t,p (K)
≤ ε,
(12)
then we can build global approximants fN ∈ VT ? SN such that f ∈ W t+r,q ? f ? fN
X
≤ CN ?r/d ,
(13)
with q and r connected by the relation 1/q = 1/p + r/d and assuming s + t < m. The e?ective construction of T (ε) is not always feasible, in particular due to the conformity constraints on the partition which does not allow to connect very coarse and very ?ne elements without intermediate grading. However, this result shows that from an intuitive point of view, the adaptive ?nite element counterpart to wavelet thresholding amounts in equilibrating the local error over the partition. One can actually use these ideas in order to obtain the estimate (9) for adaptive ?nite elements under the more restrictive assumption that 1/q < 1/p + r/d. Let us ?nally mention that the approximation theory for adaptive ?nite elements without shape constraints is an open problem.
4. Nonlinear approximation and PDE’s
Nonlinear approximation theory has opened new lines of research on the theory of PDE’s and their numerical discretization. On the one hand, it is worth revisiting the regularity theory of certain PDE’s for which the solutions develop singularities but might possess signi?cantly higher smoothness in the scale of function spaces which govern the rate of nonlinear approximation in a given norm than in the scale which govern the rate of linear approximation in the same norm. Results of this type have been proved in particular for elliptic problems on nonsmooth domains [20] and for scalar 1D conservation laws [25]. These results show that if u is the solution of such equations, the rate of decay of σN (u) is signi?cantly higher for best N -term approximation than for the projection on uniform ?nite element spaces, therefore advocating for the use of adaptive discretizations of such PDE’s. On the other hand these results also provide with an ideal benchmark for adaptive discretizations of the equation, since σN (u) represents the best accuracy which can be achieved by N parameters. In the wavelet case these parameters are typically the N largest coe?cients of the exact solution u. However, in the practice of solving a PDE, these coe?cients are not known, and neither is the set Λ corresponding to the indices of the N largest contributions dλ ψλ . It is therefore needed to develop appropriate adaptive resolution strategies as a substitute to the thresholding procedure. Such strategies aim at detecting the indices of the largest coe?cients of the solutions and to compute them accurately, in a similar way that adaptive mesh re?nement strategies aim at contructing the optimal mesh for ?nite element approximation. In both contexts, we could hope for an algorithm which builds approximations uN ∈ ΣN such that u ? uN X is bounded up to a ?xed
Adaptive Methods for PDE’s Wavelets or Mesh Re?nement?
613
multiplicative constant by σN (u) for a given norm of interest, but this requirement is so far out of reach. A more reasonable goal is that the adaptive strategy exhibits the optimal rate of approximation : if σN (u) ≤ CN ?s for some s > 0, then u?uN X ≤ CN ?s up to a change in the constant. Another requirement is that the adaptive algorithm should be scalable, i.e. the number of elementary operations in order to compute uN remains proportional to N . Let us ?nally remark that the norm · X for which error estimates can be obtained is often dictated by the nature of the equation (for example X = H 1 in the case of a second order elliptic problem) and that additional di?culties can be expected if one searches for estimates in a di?erent norm.
5. The classical approach
The classical approach to numerically solving linear and nonlinear partial differential or integral equations F (u) = 0 by the ?nite element method is typically concerned with the following issues : (c1) Well-posedness of the equation, i.e. existence, uniqueness and stability of the solution. (c2) Discretization into a ?nite element problem FT (uT ) = 0 by the Galerkin method with uT ∈ VT , analysis of well-posedness and of the approximation error u ? uT X . (c3) Numerical resolution of the ?nite dimensional system. (c4) Mesh re?nement based on a-posteriori error estimators in the case of adaptive ?nite element methods. Several di?culties are associated to each of these steps. First of all, note that the well-posedness of the ?nite element problem is in general not a consequence of the well-posedness of the continuous problem. Typical examples even in the linear case are saddle point problems. For such problems, it is well known that, for Galerkin discretizations to be stable, the ?nite element spaces for the di?erent solution components have to satisfy certain compatibility conditions (LBB or InfSup condition), which are also crucial in the derivation of optimal error estimates. Thus the discrete problem does not necessarily inherit the “nice properties” of the original in?nite dimensional problem. Concerning the numerical resolution of the discrete system, a typical source of trouble is its possible ill-conditioning, which interferes with the typical need to resort on iterative solvers in high dimension. An additional di?culty occuring in the case of integral equations is the manipulation of matrices which are densely populated. Finally, let us elaborate more on the adaptivity step. Since more than two decades, the understanding and practical realization of adaptive re?nement schemes in a ?nite element context has been documented in numerous publications [1, 2, 3, 27, 33]. Key ingredients in most adaptive algorithms are a-posteriori error estimators which are typically derived from the current residual F (uT ) : in the case where the Frechet derivative DF (u) is an isomorphism between Banach function spaces X to Y , one can hope to estimate the error u ? uT X by the evaluation of
614
Albert Cohen
F (uT ) Y . The rule of thumb is then to decompose F (uT ) Y into computable local error indicators ηK which aim to describe as accurately as possible the local error on each element K ∈ T . In the case of elliptic problems, these indicators typically consist of local residuals and other quantities such as jumps of derivatives across the interface between adjacent elements. A typical re?nement algorithm will subdivide those elements K for which the error indicator ηK is larger than a prescribed ? tolerance ε resulting in a new mesh T . Note that this strategy is theoretically in accordance with our remarks in §3 on adaptive ?nite element approximation, since it tends to equilibrate the local error. Two other frequently used strategies consist in re?ning a ?xed proportion of the elements corresponding to the largest ηK , or the smallest number of elements K for which the ηK contribute to the global error up to a ?xed proportion. It is therefore hoped that the iteration of this process from an initial mesh T0 will produce optimal meshes (Tn )n≥0 in the sense that the associated solutions un := uTn ∈ VTn converge to u at the optimal rate : σN (u) ≤ CN ?r ? u ? un
X
≤ C[#(Tn )]?r ,
(14)
up to a change in the constant C. Unfortunately, severe obstructions appear when trying to prove (14) even in the simplest model situations. One of them is that ηK is in general not an estimate by above of the local error, reducing the chances to derive the optimal rate. For most adaptive re?nement algorithms, the theoretical situation is actually even worse in the sense that it cannot even be proved that the re?nement step actually results in a reduction of the error by a ?xed amount and that un converges to u as n grows. Only recently [26, 30] have proof of convergence appeared for certain type of adaptive ?nite element methods, yet without convergence rate and therefore no guaranteed advantage over their non-adaptive counterparts.
6. A new paradigm
Wavelet methods vary from ?nite element method in that they can be viewed as solving systems that are ?nite sections of one ?xed in?nite dimensional system corresponding to the discretization of the equation in the full basis. This observation has led to a new paradigm which has been explored in [15] for linear variational problems. It aims at closely intertwining the analysis—discretization—solution process. The basic steps there read as follows : (n1) Well-posedness of the variational problem. (n2) Discretization into an equivalent in?nite dimensional problem which is well posed in ?2 . (n3) Devise an iterative scheme for the ?2 -problem that exhibits a ?xed error reduction per iteration step. (n4) Numerical realization of the iterative scheme by means of an adaptive application of the involved in?nite dimensional operators within some dynamically updated accuracy tolerances. Thus the starting point (n1) is the same. The main di?erence is that one aims at staying as long as possible with the in?nite dimensional problem. Only at the very
Adaptive Methods for PDE’s Wavelets or Mesh Re?nement?
615
end, when it comes to applying the operators in the ideal iteration scheme (n4), one enters the ?nite dimensional realm. However, the ?nite number of degrees of freedom is determined at each stage by the adaptive application of the operator, so that at no stage any speci?c trial space is ?xed. The simplest example is provided by the Poisson equation ??u = f on a domain ? with homogeneous boundary conditions, for which the variational formu1 lation in X = H0 reads : ?nd u ∈ X such that a(u, v) = L(v), for all v ∈ X, (15)
with a(u, v) := ? ?u?v and L(v) := ? f v. The well-posedness for a data f ∈ X ′ = H ?1 is ensured by the Lax-Milgram lemma. In the analysis of the wavelet discretization of this problem, we shall invoke the fact that wavelet bases provide norm equivalence for Sobolev spaces in terms of weighted ?2 norms of the coe?cients : if u = λ uλ ψλ , one has u
2 Hs
?
λ
uλ ψλ
2 Hs
?
λ
22s|λ| |uλ |2 .
(16)
We refer to [13] and [21] for the general mechanism allowing to derive these equivalences, in particular for Sobolev spaces on domains with boundary conditions such 1 as H0 . Therefore, if we renormalize our system in such a way that ψλ X = 1, we obtain the norm equivalence u 2 ? U 2, (17) X where U := (uλ )λ∈? and obtains · denotes the ?2 norm. By duality, one also easily f
2 X′
? F
2
,
(18)
with F := ( f, ψλ )λ∈? . The equivalent ?2 system is thus given by AU = F, (19)
where A(λ, ?) = a(ψλ , ψ? ) is a symmetric positive de?nite matrix which is continuous and coercive in ?2 . In this case, a converging in?nite dimensional algorithm can simply be obtained by the Richardson iteration U n := U n?1 + τ (F ? AU n?1 ) (20)
with 0 < τ < 2[λmax (A)]?1 and U 0 = 0, which guarantees the reduction rate U ? U n ≤ ρ U ? U n?1 with ρ = max{1 ? τ λmin (A), τ λmax (A) ? 1}. Note that renormalizing the wavelet system plays the role of a multiscale preconditioning, similar to multigrid yet operated at the in?nite dimensional level. At this stage, one enters ?nite dimensional adaptive computation by modifying the Richardson iteration up to a prescribed tolerance according to U n := U n?1 + τ (COARSE(F, ε) ? APPROX(AU n?1 , ε)) (21)
where F ? COARSE(F, ε) ≤ ε and AU ? APPROX(AU n?1 , ε) ≤ ε, and the U n are now ?nite dimensional vector supported by adaptive sets of indices Λn . The
616
Albert Cohen
procedure COARSE, which simply corresponds to thresholding the data vector F at a level corresponding to accuracy ε, can be practically achieved without the full knowledge of F by using some a-priori bounds on the size of the coe?cients f, ψλ , exploiting the local smoothness of f and the oscillation properties of the wavelets. The procedure APPROX deserves more attention. In order to limitate the spreading e?ect of the matrix A, one invokes its compressibility properties, namely the possibily to truncate it into a matrix AN with N non-zero entries per rows and columns in such a way that A ? AN
?2 →?2
≤ CN ?s .
(22)
The rate of compressibility s depends on the available a-priori estimates on the o?-diagonal entries A(λ, ?) := ? ?ψλ ?ψ? which are consequences of the smoothness and vanishing moment properties of the wavelet system, see [14]. Once these properties are established, a ?rst possibility is thus to choose APPROX(AU n?1 , ε) = AN U n?1 (23)
with N large enough so that accuracy ε is ensured. Clearly the modi?ed iteration (21) satis?es U ?U n ≤ ρ U ?U n?1 +2τ ε, and therefore ensures a ?xed reduction 2τ rate until the error is of the order 1?ρ ε, or until the residual F ? AU n is of order ε. A natural idea is therefore to update dynamically the tolerance ε, which is ?rst set to 1 and divided by 2 each time the approximate residual COARSE(F, ε)? APPROX(AU n?1 , ε) is below [ 2τ A +3]ε (which is ensured to happen after a ?xed 1?ρ number of steps). We therefore obtain a converging adaptive strategy, so far without information about the convergence rate. It turns out that the optimal convergence rate can also be proved, with a more careful tuning of the adaptive algorithm. Two additional ingredients are involved in this tuning. Firstly, the adaptive matrix vector multiplication APPROX has to be designed in a more elaborate way than (23) which could have the e?ect of in?ating too much the sets Λn . Instead, one de?nes for a ?nite length vector V
j 2τ A 1?ρ
APPROX(AV, ε) =
l=0
A2j?l [V2l ? V2l?1 ]
(24)
where VN denotes the restriction of V to its N largest components (with the notation j V1/2 = 0), and j is the smallest positive integer such that the residual l=0 A ? A2j?l V2l ?V2l?1 + A V ?V2j is less than ε. In this procedure, the spreading of the operator is more important on the largest coe?cients which are less in number, resulting in a signi?cant gain in the complexity of the outcome. Secondly, additional coarsening steps are needed in order to further limitate the spreading of the sets Λn and preserve the optimal rate of convergence. More precisely, the procedure COARSE is applied to U n with a tolerance proportional to ε, for those n such that ε will be updated at the next iteration.
Adaptive Methods for PDE’s Wavelets or Mesh Re?nement?
617
With such additional ingredients, it was proved in [15] that the error has the optimal rate of decay in the sense that σN (u) ≤ CN ?s ? u ? un
X
? U ? U n ≤ C[#(Λn )]?s ,
(25)
and that moreover, the computational cost of producing un remains proportional to #(Λn ). It is interesting to note that this strategy extends to non-elliptic problems such as saddle-points problems, without the need for compatibility conditions, since one inherits the well-posedness of the continuous problem which allows to obtain a converging in?nite dimensional iteration, such the Uzawa algorithm or a gradient descent applied to the least-square system (see also [15, 19]). The extension to nonlinear variational problems, based on in?nite dimensional relaxation or Newton iterations, has also been considered in [16]. It requires a speci?c procedure for the application of the nonlinear operator in the wavelet coe?cients domain which generalizes (23). It should also be mentioned that matrix compressibility also applies in the case of integral operators which have quasi-sparse wavelet discretizations. Therefore several of the obstructions from the classical approach—conditioning, compatibility, dense matrices—have disappeared in the wavelet approach. Let us ?nally mention that the coarsening steps are not really needed in the practical implementations of the adaptive wavelet method (for those problems which have been considered so far) which still does exhibit optimal convergence rate. However, we do not know how to prove (25) without these coarsening steps. There seems to be a similar situation in the ?nite element context : it has recently been proved in [10] that (14) can be achieved by an adaptive mesh re?nement algorithm which incorporates coarsening steps, while these steps are not needed in practice.
7. Conclusions and shortcomings
There exist other approaches for the development of e?cient wavelet-based adaptive schemes. In particular, an substantial research activity has recently been devoted to multiresolution adaptive processing techniques, following the line of idea introduced in [28, 29]. In this approach, one starts from a classical and reliable scheme on a uniform grid (?nite element, ?nite di?erence or ?nite volume) and applies a discrete multiresolution decomposition to the numerical data in order to compress the computational time and memory space while preserving the accuracy of the initial scheme. Here the adaptive sets Λn are therefore limited within the resolution level of the uniform grid where the classical scheme operates. This approach seems more appropriate for hyperbolic initial value problems [17, 22], in which a straightforward wavelet discretization might fail to converge properly. It should again be compared to its adaptive mesh re?nement counterpart such as in [6, 5]. Let us conclude by saying that despite its theoretical success, in the sense of achieving for certain classes of problems the optimal convergence rate with respect to the number of degrees of freedom, the wavelet-based approach to adaptive numerical simulation su?ers from three major curses.
618
Albert Cohen
The curse of geometry : while the construction of wavelet bases on a rectangular domains is fairly simple—one can use tensor product techniques and inherit the simplicity of the univariate construction—it is by far less trivial for domains with complicated geometries. Several approaches have been proposed to deal with this situation, in particular domain decomposition into rectangular patches or hierarchical ?nite element spaces, see [13, 21], and concrete implementations are nowaday available, but they result in an unavoidable loss of structural simplicity in comparison to the basic Haar system of §3. The curse of data structure : encoding and manipulating the adaptive wavelet approximations U n to the solution means that we both store the coe?cients and the indices of the adaptive set Λn which should be dynamically updated. The same goes for the indices of the matrix A which are used in the matrix-vector algorithm (24) at each step of the algorithm. This dynamical adaptation, which requires appropriate data structure, results in major overheads in the computational cost which are observed in practice : the numerical results in [4] reveal that while the wavelet adaptive algorithm indeed exhibits the optimal rate of convergence and slightly outperforms adaptive ?nite element algorithms from this perspective, the latter remains signi?cantly more e?cient from the point of view of computational time. The curse of anisotropy : adaptive wavelet approximation has roughly speaking the same properties as isotropic re?nement. However, many instances of singularities such as boundary layers and shock waves, have anisotropic features which suggests that the re?nement should be more pronounced in one particular direction. From a theoretical point of view, the following example illustrate the weakness of wavelet bases in this situation : if f = χ? with ? ? I d a smooth domain, then R the rate of best N -term approximation in X = L2 is limited to r = 1/(2d ? 2) and therefore deteriorates as the dimension grows. Wavelet bases should therefore be reconsidered if one wants to obtain better rates which take some advantage of the geometric smoothness of the curves of dicsontinuities. On the adaptive ?nite element side, anisotropic re?nement has been considered and practically implemented, yet without a clean theory available for the design of an optimal mesh. The signi?cance of wavelets in numerical analysis remains therefore tied to these curses and future breakthrough are to be expected once simple and appropriate solutions are proposed in order to deal with them.
References
[1] Babushka, I. and W. Reinhbolt (1978) A-posteriori analysis for adaptive ?nite element computations, SIAM J. Numer. Anal. 15, 736–754. [2] Babuˇka, I. and A. Miller (1987), A feedback ?nite element method with as posteriori error estimation: Part I. The ?nite element method and some basic properties of the a-posteriori error estimator, Comput. Methods Appl. Mech.
Adaptive Methods for PDE’s Wavelets or Mesh Re?nement?
619
Engrg. 61, 1–40. [3] Bank, R.E. and A. Weiser (1985), Some a posteriori error estimates for elliptic partial di?erential equations, Math. Comp., 44, 283–301. [4] Barinka, A., T. Barsch, P. Charton, A. Cohen, S. Dahlke, W. Dahmen, K. Urban (1999), Adaptive wavelet schemes for elliptic problems—Implementation and numerical experiments, IGPM Report # 173 RWTH Aachen, SIAM J. Sci. Comp. 23, 910–939. [5] Berger, M. and P. Collela (1989) Local adaptive mesh re?nement for shock hydrodynamics, J. Comp. Phys. 82, 64–84. [6] Berger, M. and J. Oliger (1984) Adaptive mesh re?nement for hyperbolic partial di?erential equations, J. Comp. Phys. 53, 482–512. [7] Bertoluzza, S. (1995) A posteriori error estimates for wavelet Galerkin methods, Appl. Math. Lett. 8, 1–6. [8] Bertoluzza, S. (1997) An adaptive collocation method based on interpolating wavelets, in Multiscale Wavelet Methods for PDEs, W. Dahmen, A. J. Kurdila, P. Oswald (eds.), Academic Press, 109–135. [9] Bihari, B. and A. Harten (1997) Multiresolution schemes for the numerical solution of 2–D conservation laws, SIAM J. Sci. Comput. 18, 315–354. [10] Binev, P., W. Dahmen and R. DeVore (2002), Adaptive ?nite element methods with convergence rates, preprint IGPM-RWTH Aachenm, to appear in Numerische Mathematik. [11] Canuto, C. and I. Cravero (1997), Wavelet-based adaptive methods for advection-di?usion problems, Math. Mod. Meth. Appl. Sci. 7, 265–289. [12] Ciarlet, P.G. (1991), Basic error estimates for the ?nite element method, Handbook of Numerical Analysis, vol II, P. Ciarlet et J.-L. Lions eds., Elsevier, Amsterdam. [13] Cohen, A. (2000), Wavelets in numerical analysis, Handbook of Numerical Analysis, vol. VII, P.G. Ciarlet and J.L. Lions, eds., to appear in 2003 as a book“Numerical analysis of wavelet methods”, Elsevier, Amsterdam. [14] Cohen, A., W. Dahmen and R. DeVore (2000), Adaptive wavelet methods for elliptic operator equations—convergence rate, Math. Comp. 70, 27–75. [15] Cohen, A., W. Dahmen and R. DeVore (2002), Adaptive wavelet methods for operator equations—beyond the elliptic case, Found. of Comp. Math. 2, 203– 245. [16] Cohen, A., W. Dahmen and R. DeVore (2002), Adaptive wavelet methods for nonlinear variational problems, preprint IGPM-RWTH Aachen, submitted to SIAM J. Num. Anal. [17] Cohen, A., S.M. Kaber, S. Mueller and M. Postel (2002), Fully adaptive multiresolution ?nite volume schemes for conservation laws, Math. of Comp. 72, 183–225. [18] Cohen, A. and R. Masson (1999), Wavelet adaptive methods for elliptic problems—preconditionning and adaptivity, SIAM J. Sci. Comp. 21, 1006–1026. [19] Dahlke, S., W. Dahmen and K. Urban (2001), Adaptive wavelet methods for saddle point problems—optimal convergence rates, preprint IGPM-Aachen [20] Dahlke, S. and R. DeVore (1997), Besov regularity for elliptic boundary value
620
Albert Cohen
problems, Communications in PDEs 22, 1–16. [21] Dahmen, W. (1997), Wavelet and multiscale methods for operator equations, Acta Numerica 6, 55–228. [22] Dahmen, W., B. Gottschlich-M¨ ller and S. M¨ ller (2001), Multiresolution u u Schemes for Conservation Laws, Numerische Mathematik 88, 399–443. [23] Daubechies, I. (1992), Ten lectures on wavelets, SIAM, Philadelphia. [24] DeVore, R. (1997), Nonlinear Approximation, Acta Numerica 51–150. [25] DeVore, R. and B. Lucier (1990), High order regularity for conservation laws, Indiana J. Math. 39, 413–430. [26] D¨r?er, W. (1996), A convergent adaptive algorithm for Poisson’s equation, o SIAM J. Num. Anal. 33, 1106–1124 [27] Eriksson, K., D. Estep, P. Hansbo, and C. Johnson (1995), Introduction to adaptive methods for di?erential equations, Acta Numerica 4, Cambridge University Press, 105–158. [28] Harten, A. (1994), Adaptive multiresolution schemes for shock computations, J. Comp. Phys. 115, 319–338. [29] Harten, A. (1995), Multiresolution algorithms for the numerical solution of hyperbolic conservation laws, Comm. Pure and Appl. Math. 48, 1305–1342. [30] Morin, P., R. Nocetto and K. Siebert (2000), Data oscillation and convergence of adaptive FEM, SIAM J. Num. Anal. 38, 466–488. [31] Petrushev, P. (1988), Direct and converse theorems for spline and rational approximation and Besov spaces, in Function spaces and applications, M. Cwikel, J. Peetre, Y. Sagher and H. Wallin, eds., Lecture Notes in Math. 1302, Springer Verlag, Berlin, 363–377. [32] Temlyakov, V. (1998), Best N -term approximation and greedy algorithms, Adv. Comp. Math. 8, 249–265. [33] Verf¨ rth, R. (1994), A-posteriori error estimation and adaptive mesh re?neu ment techniques, Jour. Comp. Appl. Math. 50, 67–83.