A PDE approach to fractional diffusion: a posteriori error analysis
Long Chen, Ricardo H. Nochetto, Enrique Otárola, Abner J. Salgado
AA PDE APPROACH TO FRACTIONAL DIFFUSION:A POSTERIORI ERROR ANALYSIS
LONG CHEN, RICARDO H. NOCHETTO, ENRIQUE OT ´AROLA,AND ABNER J. SALGADO
Abstract.
We derive a computable a posteriori error estimator for the α -harmonic extension problem, which localizes the fractional powers of ellipticoperators supplemented with Dirichlet boundary conditions. Our a posteriorierror estimator relies on the solution of small discrete problems on anisotropic cylindrical stars . It exhibits built-in flux equilibration and is equivalent to theenergy error up to data oscillation, under suitable assumptions. We design asimple adaptive algorithm and present numerical experiments which reveal acompetitive performance. Introduction
The objective of this work is the derivation and analysis of a computable, efficientand, under certain assumptions, reliable a posteriori error estimator for problemsinvolving fractional powers of the Dirichlet Laplace operator ( − ∆) s with s ∈ (0 , R n ( n ≥
1) with boundary ∂ Ω, s ∈ (0 , f : Ω → R be given. We shall be concerned with the following problem:find u such that(1.1) ( − ∆) s u = f, in Ω , u = 0 , on ∂ Ω . One of the main difficulties in the study of problem (1.1) is that the fractionalLaplacian is a nonlocal operator; see [5, 7, 24]. To localize it, Caffarelli and Silvestreshowed in [7] that any power of the fractional Laplacian in R n can be realized as anoperator that maps a Dirichlet boundary condition to a Neumann-type conditionvia an extension problem on the upper half-space R n +1+ . For a bounded domain Ω,the result by Caffarelli and Silvestre has been adapted in [4, 8, 35], thus obtaining anextension problem which is now posed on the semi-infinite cylinder C = Ω × (0 , ∞ ). Date : Version of November 15, 2018.2000
Mathematics Subject Classification.
Key words and phrases.
Fractional diffusion; finite elements; a posteriori error estimates; non-local operators; nonuniformly elliptic equations; anisotropic elements; adaptive algorithm.LC has been supported by NSF grants DMS-1115961, DMS-1418732, and DOE prime award a r X i v : . [ m a t h . NA ] N ov L. CHEN, R.H. NOCHETTO, E. OT´AROLA, AND A.J. SALGADO
This extension is the following mixed boundary value problem:(1.2) div ( y α ∇ U ) = 0 , in C , U = 0 , on ∂ L C , ∂ U ∂ν α = d s f, on Ω × { } , where ∂ L C = ∂ Ω × [0 , ∞ ) is the lateral boundary of C , and d s is a positive normal-ization constant that depends only on s ; see [5, 7] for details. The parameter α isdefined as(1.3) α = 1 − s ∈ ( − , , and the so-called conormal exterior derivative of U at Ω × { } is(1.4) ∂ U ∂ν α = − lim y → + y α U y . We will call y the extended variable and the dimension n + 1 in R n +1+ the ex-tended dimension of problem (1.2). The limit in (1.4) must be understood in thedistributional sense; see [5, 7, 8] for more details. As noted in [4, 7, 8, 35], thefractional Laplacian and the Dirichlet to Neumann operator of problem (1.2) arerelated by d s ( − ∆) s u = ∂ U ∂ν α in Ω . Based on the ideas presented above, the following simple strategy to find thesolution of (1.1) has been proposed and analyzed by the last three authors in [29]:given a sufficiently smooth function f we solve (1.2), thus obtaining a function U = U ( x (cid:48) , y ); setting u : x (cid:48) ∈ Ω (cid:55)→ u ( x (cid:48) ) = U ( x (cid:48) , ∈ R , we obtain the solutionof (1.1). The results of [29] provide an a priori error analysis which combinesasymptotic properties of Bessel functions with polynomial interpolation theory onweighted Sobolev spaces. The latter is valid for tensor product elements which maybe graded in Ω and exhibit a large aspect ratio in y (anisotropy) which is necessaryto fit the behavior of U ( x (cid:48) , y ) with x (cid:48) ∈ Ω and y >
0. The resulting a priori errorestimate is quasi-optimal in both order and regularity for the extended problem(1.2). These results are summarized in Section 3.The main advantage of the algorithm described above, is that we are solvingthe local problem (1.2) instead of dealing with the nonlocal operator ( − ∆) s ofproblem (1.1). However, this comes at the expense of incorporating one moredimension to the problem, thus raising the question of computationally efficiency.A quest for the answer has been the main drive in our recent research programand motivates the study of a posteriori error estimators and adaptivity. The latteris also motivated by the fact that the a priori theory developed in [29] requires f ∈ H − s (Ω) and Ω convex. If one of these conditions is violated the solution U ( x (cid:48) , y ) may have singularities in the direction of the x (cid:48) -variables and thus exhibitfractional regularity. As a consequence, quasi-uniform refinement of Ω would notresult in an efficient solution technique and then an adaptive loop driven by an aposteriori error estimator is essential to recover optimal rates of convergence.In this work we derive a computable, efficient and, under certain assumptions,reliable a posteriori error estimator and design an adaptive procedure to solve prob-lem (1.2). As the results of [29] show, meshes must be highly anisotropic in theextended dimension y if one intends for the method to be optimal. For this reason,it is imperative to design an a posteriori error estimator which is able to deal with POSTERIORI ERROR ANALYSIS FOR FRACTIONAL DIFFUSION 3 such anisotropic behavior. Before proceeding with our analysis, it is instructive tocomment about the anisotropic a posteriori error estimators and analysis advocatedin the literature.A posteriori error estimators are computable quantities, i.e., they may depend onthe computed solution, mesh and data, but not on the exact solution. They provideinformation about the quality of approximation of the numerical solution. They areproblem-dependent and may be used to make a judicious mesh refinement in orderto obtain the best possible approximation with the least amount of computationalresources. For isotropic discretizations, i.e., meshes where the aspect ratio of allcells is bounded independently of the refinement level, the theory of a posteriorierror estimation is well understood. Starting with the pioneering work of Babuˇskaand Rheinbolt [3], a great deal of work has been devoted to its study. We refer to[1, 37] for an overview of the state of the art. However, despite of what might beclaimed in the literature, the theory of a posteriori error estimation on anisotropicdiscretizations, i.e., meshes where the cells have disparate sizes in each direction, isstill in its infancy.To the best of our knowledge the first work that attempts to deal with anisotropica posteriori error estimation is [34]. In this work, a residual a posteriori errorestimator is introduced and allegedly analyzed on anisotropic meshes. However,such analysis relies on assumptions on the exact and discrete solutions and on themesh, which are neither proved nor there is a way to explicitly enforce them in thecourse of computations; see [34, §
6, Remark 3]. Subsequently, in [21] the conceptof matching function is introduced in order to derive anisotropic a posteriori errorindicators. The correct alignment of the grid with the exact solution is crucialto derive an upper bound for the error. Indeed, this upper bound involves thematching function, which depends on the error itself and then it does not providea real computable quantity; see [21, Theorem 2]. For similar works in this directionsee [23, 22, 28]. In [33], the anisotropic interpolation estimates derived in [16] areused to derive a Zienkiewicz–Zhu type of a posteriori error estimator. However, asproperly pointed out in [33, Proposition 2.3], the ensuing upper bound for the errordepends on the error itself, and thus, it is not computable .In our case, since the coefficient y α in (1.2) either degenerates ( s < /
2) or blowsup ( s > / y and the nonuniform coefficient y α , upon consideringlocal problems on cylindrical stars . The solutions of these local problems allow usto define a computable and anisotropic a posteriori error estimator which, undercertain assumptions, is equivalent to the error up to data oscillations terms. Inorder to derive such a result, a computationally implementable geometric conditionneeds to be imposed on the mesh, which does not depend on the exact solutionof problem (1.2). This approach is of value not only for (1.2), but in general foranisotropic problems since rigorous anisotropic a posteriori error estimators are notavailable in the literature.The outline of this paper is as follows. Section 2 sets the framework in whichwe will operate. Notation and terminology are introduced in § § § L. CHEN, R.H. NOCHETTO, E. OT´AROLA, AND A.J. SALGADO [29]. The need for a new approach in a posteriori error estimation is examinedin Section 4, where we show that the standard approaches either do not workor produce suboptimal results. This justifies the introduction of our new errorestimator on cylindrical stars . Section 5 is the core of this work and is dedicated tothe development and analysis of our new error estimator. After some preliminarysetup carried out in §§ § § Notation and preliminaries
Notation.
Throughout this work Ω is an open, bounded and connected do-main of R n , n ≥
1, with polyhedral boundary ∂ Ω. We define the semi-infinitecylinder with base Ω and its lateral boundary, respectively, by C := Ω × (0 , ∞ ) , ∂ L C := ∂ Ω × [0 , ∞ ) . Given Y > C Y := Ω × (0 , Y ) . Thelateral boundary ∂ L C Y is defined accordingly.Throughout our discussion we will be dealing with objects defined in R n +1 andit will be convenient to distinguish the extended dimension. A vector x ∈ R n +1 ,will be denoted by x = ( x , . . . , x n , x n +1 ) = ( x (cid:48) , x n +1 ) = ( x (cid:48) , y ) , with x i ∈ R for i = 1 , . . . , n + 1, x (cid:48) ∈ R n and y ∈ R .If X and Y are normed vector spaces, we write X (cid:44) → Y to denote that X iscontinuously embedded in Y . We denote by X (cid:48) the dual of X and by (cid:107) · (cid:107) X thenorm of X . The relation a (cid:46) b indicates that a ≤ Cb , with a constant C thatdoes not depend on a or b nor the discretization parameters. The value of C mightchange at each occurrence.2.2. The fractional Laplace operator.
Our definition is based on spectral the-ory. For any f ∈ L (Ω), the Lax-Milgram Lemma provides the existence anduniqueness of w ∈ H (Ω) that solves − ∆ w = f in Ω , w = 0 on ∂ Ω . The operator ( − ∆) − : L (Ω) → L (Ω) is compact, symmetric and positive, whenceits spectrum { λ − k } k ∈ N is discrete, real, positive and accumulates at zero. Moreover,there exists { ϕ k } k ∈ N ⊂ H (Ω), which is an orthonormal basis of L (Ω) and satisfies(2.1) − ∆ ϕ k = λ k ϕ k in Ω , ϕ k = 0 on ∂ Ω . Fractional powers of the Dirichlet Laplace operator can then be defined for w ∈ C ∞ (Ω) by(2.2) ( − ∆) s w = ∞ (cid:88) k =1 λ sk w k ϕ k , POSTERIORI ERROR ANALYSIS FOR FRACTIONAL DIFFUSION 5 where w k = ´ Ω wϕ k . By density ( − ∆) s can be extended to the space(2.3) H s (Ω) = (cid:40) w = ∞ (cid:88) k =1 w k ϕ k : ∞ (cid:88) k =1 λ sk w k < ∞ (cid:41) = H s (Ω) , s ∈ (0 , ) ,H / (Ω) , s = ,H s (Ω) , s ∈ ( , . The characterization given by the second equality is shown in [25, Chapter 1]. For s ∈ (0 ,
1) we denote by H − s (Ω) the dual of H s (Ω).2.3. The Caffarelli-Silvestre extension problem.
To exploit the Caffarelli-Silvestre result [7], or its variants [4, 6, 8], we need to deal with a nonuniformlyelliptic equation. To this end, we consider weighted Sobolev spaces with the weight | y | α , α ∈ ( − , D ⊂ R n +1 , we then define L ( | y | α , D ) to be the space of allmeasurable functions defined on D such that (cid:107) w (cid:107) L ( | y | α ,D ) = ˆ D | y | α w < ∞ . Similarly we define the weighted Sobolev space H ( | y | α , D ) = (cid:8) w ∈ L ( | y | α , D ) : |∇ w | ∈ L ( | y | α , D ) (cid:9) , where ∇ w is the distributional gradient of w . We equip H ( | y | α , D ) with the norm(2.4) (cid:107) w (cid:107) H ( | y | α ,D ) = (cid:16) (cid:107) w (cid:107) L ( | y | α ,D ) + (cid:107)∇ w (cid:107) L ( | y | α ,D ) (cid:17) / . Since α ∈ ( − ,
1) we have that | y | α belongs to the so-called Muckenhoupt class A ( R n +1 ); see [15, 17, 27, 36]. This, in particular, implies that H ( | y | α , D ) equippedwith the norm (2.4), is a Hilbert space and the set C ∞ ( D ) ∩ H ( | y | α , D ) is dense in H ( | y | α , D ) (cf. [36, Proposition 2.1.2, Corollary 2.1.6], [20] and [17, Theorem 1]).We recall now the definition of Muckenhoupt classes; see [27, 36]. Definition 2.5 (Muckenhoupt class A ) . Let ω be a weight and N ≥
1. We say ω ∈ A ( R N ) if(2.6) C ,ω = sup B (cid:18) B ω (cid:19) (cid:18) B ω − (cid:19) < ∞ , where the supremum is taken over all balls B in R N .If ω belongs to the Muckenhoupt class A ( R N ), we say that ω is an A -weight,and we call the constant C ,ω in (2.6) the A -constant of ω .To study problem (1.2) we define the weighted Sobolev space(2.7) ◦ H L ( y α , C ) = (cid:8) w ∈ H ( y α , C ) : w = 0 on ∂ L C (cid:9) . As [29, (2.21)] shows, the following weighted Poincar´e inequality holds:(2.8) ˆ C y α w (cid:46) ˆ C y α |∇ w | , ∀ w ∈ ◦ H L ( y α , C ) . Then, the seminorm on ◦ H L ( y α , C ) is equivalent to the norm (2.4). For w ∈ H ( y α , C ), we denote by tr Ω w its trace onto Ω × { } , and we recall that the traceoperator tr Ω satisfies, (see [29, Proposition 2.5], [8, Proposition 2.1])(2.9) tr Ω ◦ H L ( y α , C ) = H s (Ω) , (cid:107) tr Ω w (cid:107) H s (Ω) ≤ C tr Ω (cid:107) w (cid:107) ◦ H L ( y α , C ) . L. CHEN, R.H. NOCHETTO, E. OT´AROLA, AND A.J. SALGADO
Let us now describe the Caffarelli-Silvestre result and its extension to boundeddomains; see [7, 35]. Given f ∈ H − s (Ω), let u ∈ H s (Ω) be the solution of ( − ∆) s u = f in Ω. We define the α -harmonic extension of u to the cylinder C , as the function U ∈ ◦ H L ( y α , C ), solution of problem (1.2), namely( − ∆) s u = d s ∂ U ∂ν α in Ω , where d s = 2 − s Γ(1 − s )Γ( s ) . Finally, we must mention that(2.10) C tr Ω ≤ d − / s . Indeed, given ψ ∈ ◦ H L ( y α , C ) we define Ψ ∈ ◦ H L ( y α , C ) as the solution of − div( y α ∇ Ψ) = 0 , in C , Ψ = 0 , on ∂ L C , Ψ = tr Ω ψ on Ω × { } . It is standard to show that Ψ is the minimal norm extension of tr Ω ψ . Moreover,separation of variables gives d s (cid:107) tr Ω ψ (cid:107) H s (Ω) = (cid:107)∇ Ψ (cid:107) L ( y α , C ) , [8, Proposition 2.1].Therefore (cid:107) tr Ω ψ (cid:107) H s (Ω) = 1 d s (cid:107) Ψ (cid:107) ◦ H L ( y α , C ) ≤ d s (cid:107) ψ (cid:107) ◦ H L ( y α , C ) . Estimate (2.10) will be useful to obtain an upper bound of the error by the estima-tor. 3.
A priori error estimates
In an effort to make this contribution self-contained here we review the mainresults of [29], which deal with the a priori error analysis of discretizations of prob-lem (1.1). This will also serve to make clear the limitations of this theory, therebyjustifying the quest for an a posteriori error analysis. To do so in this section, andthis section only, we will assume the following regularity result, which is valid if,for instance, the domain Ω is convex [18](3.1) (cid:107) w (cid:107) H (Ω) (cid:46) (cid:107) ∆ x (cid:48) w (cid:107) L (Ω) , ∀ w ∈ H (Ω) ∩ H (Ω) . Since C is unbounded, problem (1.2) cannot be directly approximated with finite-element-like techniques. However, as [29, Proposition 3.1] shows, the solution U ofproblem (1.2) decays exponentially in the extended variable y so that, by truncat-ing the cylinder C to C Y and setting a vanishing Dirichlet condition on the upperboundary y = Y , we only incur in an exponentially small error in terms of Y [29,Theorem 3.5].Define ◦ H L ( y α , C Y ) = (cid:8) v ∈ H ( y α , C Y ) : v = 0 on ∂ L C Y ∪ Ω × { Y } (cid:9) . Then, the aforementioned problem reads: find v ∈ ◦ H L ( y α , C Y ) such that(3.2) ˆ C Y y α ∇ v ∇ φ = d s (cid:104) f, tr Ω φ (cid:105) H − s (Ω) × H s (Ω) , for all v ∈ ◦ H L ( y α , C Y ), where (cid:104)· , ·(cid:105) H − s (Ω) × H s (Ω) denotes the duality pairing between H − s (Ω) and H s (Ω), which is well defined as a consequence of (2.9).If U and v denote the solution of (1.2) and (3.2), respectively, then [29, Theorem3.5] provides the following exponential estimate (cid:107)∇ ( U − v ) (cid:107) L ( y α , C ) (cid:46) e −√ λ Y / (cid:107) f (cid:107) H − s (Ω) , POSTERIORI ERROR ANALYSIS FOR FRACTIONAL DIFFUSION 7 where λ denotes the first eigenvalue of the Dirichlet Laplace operator and Y is thetruncation parameter.In order to study the finite element discretization of problem (3.2) we must firstunderstand the regularity of the solution U , since an error estimate for v , solutionof (3.2), depends on the regularity of U as well [29, § U is much worse in theextended direction, namely (cid:107) ∆ x (cid:48) U (cid:107) L ( y α , C ) + (cid:107) ∂ y ∇ x (cid:48) U (cid:107) L ( y α , C ) (cid:46) (cid:107) f (cid:107) H − s (Ω) , (3.3) (cid:107) U yy (cid:107) L ( y β , C ) (cid:46) (cid:107) f (cid:107) L (Ω) , (3.4)where β > α + 1. This suggests that graded meshes in the extended variable y play a fundamental role. In fact, estimates (3.3)–(3.4) motivate the construction ofa mesh over C Y as follows. We first consider a graded partition I Y of the interval[0 , Y ] with mesh points(3.5) y k = (cid:18) kM (cid:19) γ Y , k = 0 , . . . , M, where γ > / (1 − α ) = 3 / (2 s ). We also consider T Ω = { K } to be a conformingand shape regular mesh of Ω, where K ⊂ R n is an element that is isoparametricallyequivalent either to the unit cube [0 , n or the unit simplex in R n . The collectionof these triangulations T Ω is denoted by T Ω . We construct the mesh T Y as thetensor product triangulation of T Ω and I Y . In order to obtain a global regularityassumption for T Y , we assume that there is a constant σ Y such that if T = K × I and T = K × I ∈ T Y have nonempty intersection, then(3.6) h I h I ≤ σ Y , where h I = | I | . It is well known that this weak regularity condition on the meshallows for anisotropy in the extended variable (cf. [14, 29]). The set of all triangu-lations of C Y that are obtained with this procedure and satisfy these conditions isdenoted by T . Figure 1 shows an example of this type of meshes in three dimensions.For T Y ∈ T , we define the finite element space(3.7) V ( T Y ) = (cid:8) W ∈ C ( C Y ) : W | T ∈ P ( K ) ⊗ P ( I ) ∀ T ∈ T Y , W | Γ D = 0 (cid:9) . where Γ D = ∂ L C Y ∪ Ω × { Y } is called the Dirichlet boundary; the space P ( K ) is P ( K ) — the space of polynomials of total degree at most 1, when the base K ofan element T = K × I is simplicial. If K is an n -rectangle P ( K ) stands for Q ( K )— the space of polynomials of degree not larger than 1 in each variable. We alsodefine U ( T Ω ) = tr Ω V ( T Y ), i.e., a P ( K ) finite element space over the mesh T Ω .The Galerkin approximation of (3.2) is given by the unique function V T Y ∈ V ( T Y ) such that(3.8) ˆ C Y y α ∇ V T Y ∇ W = d s (cid:104) f, tr Ω W (cid:105) H − s (Ω) × H s (Ω) , ∀ W ∈ V ( T Y ) . Existence and uniqueness of V T Y immediately follows from V ( T Y ) ⊂ ◦ H L ( y α , C Y )and the Lax-Milgram Lemma. It is trivial also to obtain a best approximationresult `a la Cea. This best approximation result reduces the numerical analysis ofproblem (3.8) to a question in approximation theory which in turn can be answeredwith the study of piecewise polynomial interpolation in Muckenhoupt weighted
L. CHEN, R.H. NOCHETTO, E. OT´AROLA, AND A.J. SALGADO
Figure 1.
A three dimensional graded mesh of the cylinder(0 , × (0 , Y ) with 392 degrees of freedom. The mesh is con-structed as a tensor product of a quasi-uniform mesh of (0 , with cardinality 49 and the image of the quasi-uniform partitionof the interval (0 ,
1) with cardinality 8 under the mapping (3.5).Sobolev spaces; see [29, 30]. Exploiting the Cartesian structure of the mesh ispossible to handle anisotropy in the extended variable, construct a quasi interpolantΠ T Y : L ( C Y ) → V ( T Y ), and obtain (cid:107) v − Π T Y v (cid:107) L ( y α ,T ) (cid:46) h K (cid:107)∇ x (cid:48) v (cid:107) L ( y α ,S T ) + h I (cid:107) ∂ y v (cid:107) L ( y α ,S T ) , (cid:107) ∂ x j ( v − Π T Y v ) (cid:107) L ( y α ,T ) (cid:46) h K (cid:107)∇ x (cid:48) ∂ x j v (cid:107) L ( y α ,S T ) + h I (cid:107) ∂ y ∂ x j v (cid:107) L ( y α ,S T ) , with j = 1 , . . . , n +1; see [29, Theorems 4.6–4.8] and [30] for details. However, since U yy ≈ y − α − as y ≈
0, we realize that U / ∈ H ( y α , C Y ) and the second estimate isnot meaningful for j = n + 1. In view of estimate (3.4) it is necessary to measurethe regularity of U yy with a stronger weight and thus compensate with a gradedmesh in the extended dimension. This makes anisotropic estimates essential.Notice that T Y = M T Ω , and that T Ω ≈ M n implies T Y ≈ M n +1 .Finally, if T Ω is shape regular and quasi-uniform, we have h T Ω ≈ ( T Ω ) − /n . Allthese considerations allow us to obtain the following result; see [29, Theorem 5.4]and [29, Corollary 7.11]. Theorem 3.9 (a priori error estimate) . Let T Y ∈ T be a tensor product grid, whichis quasi-uniform in Ω and graded in the extended variable so that (3.5) holds. If V ( T Y ) is defined by (3.7) and V T Y ∈ V ( T Y ) is the Galerkin approximation definedby (3.8) , then we have (cid:107) U − V T Y (cid:107) ◦ H L ( y α , C ) (cid:46) | log( T Y ) | s ( T Y ) − / ( n +1) (cid:107) f (cid:107) H − s (Ω) , where Y ≈ log( T Y ) . Alternatively, if u denotes the solution of (1.1) , then (cid:107) u − V T Y ( · , (cid:107) H s (Ω) (cid:46) | log( T Y ) | s ( T Y ) − / ( n +1) (cid:107) f (cid:107) H − s (Ω) . Remark 3.10 (domain and data regularity) . The results of Theorem 3.9 hold trueonly if f ∈ H − s (Ω) and the domain Ω is such that (3.1) holds. POSTERIORI ERROR ANALYSIS FOR FRACTIONAL DIFFUSION 9 A posteriori error estimators: the search for a new approach
The function U , solution of the α -harmonic extension problem (1.2), has asingular behavior on the extended variable y , which is compensated by consideringanisotropic meshes in this direction as dictated by (3.5). However, the solution U , may also have singularities in the direction of the x (cid:48) -variables and thus exhibitfractional regularity, which would not allow us to attain the almost optimal rateof convergence given by Theorem 3.9. In fact, as Remark 3.10 indicates, it isnecessary to require that f ∈ H − s (Ω) and that the domain has the property(3.1) in order to have an almost optimal rate of convergence. If any of thesetwo conditions fail singularities may develop in the direction of the x (cid:48) -variables,whose characterization is as yet an open problem; see [29, § Residual estimators.
Simply put, the so-called residual error estimators usethe strong form of the local residual as an indicator of the error. To obtain thestrong form of the equation, integration by parts is necessary. Let us consider anelement T ∈ T Y and integrate by parts the term ˆ T y α ∇ V T Y · ∇ W = ˆ ∂T W y α ∇ V T Y · ν − ˆ T div( y α ∇ V T Y ) W, where ν denotes the unit outer normal to T . Since α ∈ ( − ,
1) the boundaryintegral is meaningless for y = 0. As we see, even the very first step (integrationby parts) in the derivation of a residual a posteriori error estimator fails! At thispoint there is nothing left to do but to consider a different type of estimator.4.2. Local problems on stars over isotropic refinements.
Inspired by [2,26] we can construct, over shape regular meshes, a computable error estimatorbased on the solution of small discrete problems on stars. Its construction andanalysis is similar to the developments of Section 5 so we shall not dwell on thisany further. Since we consider shape regular meshes, such estimator is equivalentto the error up to data oscillation without any additional conditions on the mesh,but under some suitable assumptions; see § ,
1) and s ∈ (0 , f ( x (cid:48) ) = π s sin( πx (cid:48) ), so that u ( x (cid:48) ) = sin( πx (cid:48) ), and the solution U to (1.2) is U ( x (cid:48) , y ) = 2 − s π s Γ( s ) sin( πx (cid:48) ) K s ( πy ) , where K s denotes the modified Bessel function of the second kind; see [29, § α -harmonic extension we are solving a twodimensional problem so the optimal rate of convergence in the H ( y α , C )-seminorm that we expect is O ( T − . Y ). Figure 2 shows the experimental rate of convergenceof this algorithm for the cases s = 0 . s = 0 . O (cid:16) T − s/ Y (cid:17) and coincides with the suboptimal one obtained with quasi-uniform refinement; see[29, § cylindrical stars togetherwith a new anisotropic error estimator, which will treat the x (cid:48) -coordinates and theextended direction, y , separately. −2 −1 Degrees of Freedom (DOFs) E rr o r DOFs − . s = 0 . − . s = 0 . Figure 2.
Computational rate of convergence T Y ) − s/ for anisotropic adaptive algorithm for n = 1, s = 0 . s = 0 . A posteriori error estimators: cylindrical stars
It has proven rather challenging to derive and analyze a posteriori error esti-mators over a fairly general anisotropic mesh. For this reason, we introduce an implementable geometric condition which will allow us to consider graded meshesin Ω in order to compensate for possible singularities in the x (cid:48) -variables, whilepreserving the anisotropy in the extended direction, necessary to retain optimalorders of approximation. We thus assume the following condition over the familyof meshes T : there exists a positive constant C T such that for every mesh T Y ∈ T (5.1) h Y ≤ C T h z (cid:48) , for all the interior nodes z (cid:48) of T Ω , where h Y denotes the largest mesh size in the y direction, and h z (cid:48) ≈ | S z (cid:48) | /n ; see § h z (cid:48) and S z (cid:48) .We remark that this condition is satisfied in the case of quasi-uniform refinementin the variable x (cid:48) , which is a consequence of the convexity of the function involvedin (3.5). In fact, a simple computation shows(5.2) h Y = y M − y M − = Y M γ (cid:0) ( M ) γ − ( M − γ (cid:1) ≤ γ Y M , where γ > / (1 − α ) = 3 / (2 s ). We must reiterate that this mesh restriction is fullyimplementable. We refer the reader to Section 6 for more details on this. POSTERIORI ERROR ANALYSIS FOR FRACTIONAL DIFFUSION 11
Remark 5.3 ( s -independent mesh grading) . We point out that the term γ = γ ( s )in (5.2) deteriorates as s becomes small because γ > / (2 s ). However, a modifiedmesh grading in the y -direction has been proposed in [12, § s ∈ (0 , h Y ≤ C Y /M where C does not depends on s .5.1. Preliminaries.
Let us begin the discussion on a posteriori error estimationwith some terminology and notation. Given a node z on the mesh T Y , we exploitthe tensor product structure of T Y , and we write z = ( z (cid:48) , z (cid:48)(cid:48) ) where z (cid:48) and z (cid:48)(cid:48) arenodes on the meshes T Ω and I Y respectively.Given a cell K ∈ T Ω , we denote by N ( K ) and ◦ N ( K ) the set of nodes and interiornodes of K , respectively. We set N ( T Ω ) = (cid:91) K ∈ T Ω N ( K ) , ◦ N ( T Ω ) = (cid:91) K ∈ T Ω ◦ N ( K ) . Given T ∈ T Y , we define N ( T ) and ◦ N ( T ) accordingly, i.e., as the set of nodes andinterior and Neumann nodes of T , respectively. Similarly, we define ◦ N ( T Y ) and N ( T Y ). Any discrete function W ∈ V ( T Y ) is uniquely characterized by its nodalvalues on the set ◦ N ( T Y ). Moreover, the functions φ z ∈ V ( T Y ), z ∈ ◦ N ( T Y ), suchthat φ z ( w ) = δ z w for all w ∈ N ( T Y ) are the canonical basis of V ( T Y ), i.e., W = (cid:88) z ∈ ◦ N ( T Y ) W ( z ) φ z . The functions { φ z : z ∈ ◦ N ( T Y ) } are the so-called shape functions of V ( T Y ).Analogously, given a node z (cid:48) ∈ ◦ N ( T Ω ), we also consider the discrete functions ϕ z (cid:48) ∈ U ( T Ω ) = tr Ω V ( T Y ) defined by ϕ z (cid:48) ( w (cid:48) ) = δ z (cid:48) w (cid:48) for all w (cid:48) ∈ N ( T Ω ). The set { ϕ z (cid:48) : z (cid:48) ∈ ◦ N ( T Ω ) } is the canonical basis of U ( T Ω ).The shape functions { φ z : z ∈ N ( T Y ) } satisfy two properties which will proveuseful in the sequel. First, we have the so-called partition of unity property , i.e.,(5.4) (cid:88) z ∈ N ( T Y ) φ z = 1 in ¯ C Y . Second, for any z ∈ ◦ N ( T Y ), the corresponding shape function φ z belongs to V ( T Y )whence we have the so-called Galerkin orthogonality , i.e.,(5.5) ˆ C Y y α ∇ ( v − V T Y ) ∇ φ z = 0 . The partition of unity property also holds for the shape functions { ϕ z (cid:48) : z (cid:48) ∈ N ( T Ω ) } : (cid:88) z (cid:48) ∈ N ( T Ω ) ϕ z (cid:48) = 1 in ¯Ω . Given z (cid:48) ∈ N ( T Ω ) and the associated shape function ϕ z (cid:48) , we define the extended shape function ˜ ϕ z (cid:48) by ˜ ϕ z (cid:48) ( x (cid:48) , y ) = ϕ z (cid:48) ( x (cid:48) ) (0 , Y ) ( y ). These functions satisfy thefollowing partition of unity property:(5.6) (cid:88) z (cid:48) ∈ N ( T Ω ) ˜ ϕ z (cid:48) = 1 in ¯ C Y . Given z (cid:48) ∈ N ( T Ω ), we define the star around z (cid:48) as S z (cid:48) = (cid:91) K (cid:51) z (cid:48) K ⊂ Ω , and the cylindrical star around z (cid:48) as C z (cid:48) := (cid:91) { T ∈ T Y : T = K × I, K (cid:51) z (cid:48) } = S z (cid:48) × (0 , Y ) ⊂ C Y . Given an element K ∈ T Ω we define its patch as S K := (cid:83) z (cid:48) ∈ K S z (cid:48) . For T ∈ T Y itspatch S T is defined similarly. Given z (cid:48) ∈ N ( T Ω ) we define its cylindrical patch as D z (cid:48) := (cid:91) {C w (cid:48) : w (cid:48) ∈ S z (cid:48) } ⊂ C Y . For each z (cid:48) ∈ N ( T Ω ) we set h z (cid:48) := min { h K : K (cid:51) z (cid:48) } .5.2. Local weighted Sobolev spaces.
In order to define the local a posteriorierror estimators we first need to define some local weighted Sobolev spaces.
Definition 5.7 (local spaces) . Given z (cid:48) ∈ N ( T Ω ) and its associated cylindricalstar C z (cid:48) , we define W ( C z (cid:48) ) = (cid:8) w ∈ H ( y α , C z (cid:48) ) : w = 0 on ∂ C z (cid:48) \ Ω × { } (cid:9) . The space W ( C z (cid:48) ) defined above is Hilbert due to the fact that the weight | y | α belongs to the class A ( R n +1 ); see Definition 2.5. Moreover, as the following resultshows, a weighted Poincar´e-type inequality holds and, consequently, the semi-norm (cid:57) w (cid:57) C z (cid:48) = (cid:107)∇ w (cid:107) L ( y α , C z (cid:48) ) defines a norm on W ( C z (cid:48) ); see also [29, § Proposition 5.8 (weighted Poincar´e inequality) . Let z (cid:48) ∈ N ( T Ω ) . If the function w ∈ W ( C z (cid:48) ) , then we have (5.9) (cid:107) w (cid:107) L ( y α , C z (cid:48) ) (cid:46) Y (cid:57) w (cid:57) C z (cid:48) . Proof.
By density [36, Corollary 2.1.6], it suffices to reduce the considerations to asmooth function w . Given x (cid:48) ∈ S z (cid:48) , we have that w ( x (cid:48) , Y ) = 0 so that w ( x (cid:48) , y ) = − ˆ Y y ∂ y w ( x (cid:48) , ξ ) d ξ. Multiplying the expression above by | y | α , integrating over C z (cid:48) , and using the CauchySchwarz inequality, we arrive at ˆ C z (cid:48) | y | α | w ( x (cid:48) , y ) | d x (cid:48) d y ≤ ˆ C z (cid:48) | y | α (cid:18) ˆ Y | ξ | α | ∂ y w ( x (cid:48) , ξ ) | d ξ ˆ Y | ξ | − α d ξ (cid:19) d x (cid:48) d y = ˆ Y | y | α d y ˆ Y | ξ | − α d ξ ˆ C z (cid:48) | ξ | α | ∂ y w ( x (cid:48) , ξ ) | d x (cid:48) d ξ ≤ C , | y | α Y ˆ C z (cid:48) | y | α | ∂ y w ( x (cid:48) , y ) | d x (cid:48) d y, where in the third inequality we used that y α ∈ A ( R n +1 ). In conclusion, (cid:107) w (cid:107) L ( y α , C z (cid:48) ) (cid:46) Y (cid:107) ∂ y w (cid:107) L ( y α , C z (cid:48) ) , which is (5.9). (cid:3) POSTERIORI ERROR ANALYSIS FOR FRACTIONAL DIFFUSION 13
Remark 5.10 (anisotropic weighted Poincar´e inequality) . Let z (cid:48) ∈ N ( T Ω ). If w ∈ W ( C z (cid:48) ), then by extending the one-dimensional argument in the proof ofProposition 5.8 to a n -dimensional setting, we can also derive (cid:107) w (cid:107) L ( y α , C z (cid:48) ) (cid:46) h z (cid:48) (cid:107)∇ x (cid:48) w (cid:107) L ( y α , C z (cid:48) ) . An ideal a posteriori error estimator.
Here we define an ideal a poste-riori error estimator on anisotropic meshes which is not computable . However, itprovides the intuition required to define a discrete and computable error indicator,as explained in § ζ z (cid:48) ∈ W ( C z (cid:48) ) to be the solution of(5.11) ˆ C z (cid:48) y α ∇ ζ z (cid:48) ∇ ψ = d s (cid:104) f, tr Ω ψ (cid:105) H − s (Ω) × H s (Ω) − ˆ C z (cid:48) y α ∇ V T Y ∇ ψ, for all ψ ∈ W ( C z (cid:48) ). The existence and uniqueness of ζ z (cid:48) ∈ W ( C z (cid:48) ) is guaranteed bythe Lax–Milgram Lemma and the weighted Poincar´e inequality of Proposition 5.8.The continuity of the right hand side of (5.11), as a linear functional in W ( C z (cid:48) ),follows from (2.9) and the Cauchy-Schwarz inequality. We then define the globalerror estimator(5.12) ˜ E T Y = (cid:88) z (cid:48) ∈ N ( T Ω ) ˜ E z (cid:48) / , in terms of the local error indicators(5.13) ˜ E z (cid:48) = (cid:57) ζ z (cid:48) (cid:57) C z (cid:48) . The properties of this ideal estimator are as follows.
Proposition 5.14 (ideal estimator) . Let v ∈ ◦ H L ( y α , C Y ) and V T Y ∈ V ( T Y ) solve (3.2) and (3.8) , respectively. Then, the ideal estimator ˜ E T Y , defined in (5.12) – (5.13) , satisfies (5.15) (cid:107)∇ ( v − V T Y ) (cid:107) L ( y α , C Y ) (cid:46) ˜ E T Y , and for all z (cid:48) ∈ N ( T Ω )(5.16) ˜ E z (cid:48) ≤ (cid:107)∇ ( v − V T Y ) (cid:107) L ( y α , C z (cid:48) ) . Proof. If e T Y := v − V T Y denotes the error, then for any w ∈ ◦ H L ( y α , C Y ) we have ˆ C Y y α ∇ e T Y ∇ w = d s (cid:104) f, tr Ω w (cid:105) H − s (Ω) × H s (Ω) − ˆ C Y y α ∇ V T Y ∇ w = d s (cid:104) f, tr Ω ( w − W ) (cid:105) H − s (Ω) × H s (Ω) − ˆ C Y y α ∇ V T Y ∇ ( w − W )= (cid:88) z (cid:48) ∈ N ( T Ω ) d s (cid:104) f, tr Ω [( w − W ) ˜ ϕ z (cid:48) ] (cid:105) H − s (Ω) × H s (Ω) − ˆ C z (cid:48) y α ∇ V T Y ∇ [( w − W ) ˜ ϕ z (cid:48) ]for any W ∈ V ( T Y ), where to derive the expression above we have used Galerkinorthogonality (5.5) and the partition of unity property (5.6).Notice that for each z (cid:48) ∈ N ( T Ω ) the function ( w − W ) ˜ ϕ z (cid:48) ∈ W ( C z (cid:48) ). Indeed, if z (cid:48) is an interior node,(5.17) ( w − W ) ˜ ϕ z (cid:48) | ∂ C z (cid:48) \ Ω ×{ } = 0 because of the vanishing property of the shape function ϕ z (cid:48) on ∂S z (cid:48) together withthe fact that w = W = 0 on Ω × { Y } ; otherwise, if z (cid:48) is a Dirichlet node then w = W = 0 on ∂ C z (cid:48) , and we get (5.17).Set W = Π T Y w , where Π T Y denotes the quasi-interpolation operator introducedin [29, § (cid:57) ( w − W ) ˜ ϕ z (cid:48) (cid:57) C z (cid:48) as follows: (cid:57) ( w − Π T Y w ) ˜ ϕ z (cid:48) (cid:57) C z (cid:48) (cid:46) ˆ C z (cid:48) y α |∇ ( w − Π T Y w ) | ˜ ϕ z (cid:48) + ˆ C z (cid:48) y α | w − Π T Y w | |∇ x (cid:48) ˜ ϕ z (cid:48) | (cid:46) (cid:57) w (cid:57) D z (cid:48) . To bound the first term above we use the local stability of Π T Y [29, Theorems 4.7and 4.8] together with the fact that 0 ≤ ˜ ϕ z (cid:48) ≤ x ∈ C Y . For the second termwe resort to the local approximation properties of Π T Y [29, Theorems 4.7 and 4.8](5.18) ˆ C z (cid:48) y α | w − Π T Y w | |∇ x (cid:48) ˜ ϕ z (cid:48) | (cid:46) h z (cid:48) (cid:16) h z (cid:48) (cid:107)∇ x (cid:48) w (cid:107) L ( D z (cid:48) ,y α ) + h Y (cid:107) ∂ y w (cid:107) L ( D z (cid:48) ,y α ) (cid:17) (cid:46) (cid:57) w (cid:57) D z (cid:48) , where we used that |∇ x (cid:48) ˜ ϕ z (cid:48) | = |∇ x (cid:48) ϕ z (cid:48) | (cid:46) h − z (cid:48) together with (5.1).Set ψ z (cid:48) = ( w − Π T Y w ) ˜ ϕ z (cid:48) ∈ W ( C z (cid:48) ) as test function in (5.11) to obtain ˆ C Y y α ∇ e T Y ∇ w = (cid:88) z (cid:48) ∈ N ( T Ω ) ˆ C z (cid:48) y α ∇ ζ z (cid:48) ∇ ψ z (cid:48) (cid:46) (cid:88) z (cid:48) ∈ N ( T Ω ) (cid:57) ζ z (cid:48) (cid:57) C z (cid:48) (cid:57) w (cid:57) D z (cid:48) (cid:46) (cid:88) z (cid:48) ∈ N ( T Ω ) (cid:57) ζ z (cid:48) (cid:57) C z (cid:48) / (cid:107)∇ w (cid:107) L ( y α , C Y ) = ˜ E T Y (cid:107)∇ w (cid:107) L ( y α , C Y ) , where we used that (cid:57) ψ z (cid:48) (cid:57) C z (cid:48) (cid:46) (cid:57) w (cid:57) D z (cid:48) and the finite overlapping property of thestars S z (cid:48) . Since w is arbitrary, we set w = e T Y and obtain (5.15).Finally, inequality (5.16) is immediate:˜ E z (cid:48) = (cid:57) ζ z (cid:48) (cid:57) C z (cid:48) = ˆ C z (cid:48) y α ∇ ζ z (cid:48) ∇ ζ z (cid:48) = ˆ C z (cid:48) y α ∇ e T Y ∇ ζ z (cid:48) ≤ (cid:107)∇ e T Y (cid:107) L ( y α , C z (cid:48) ) (cid:57) ζ z (cid:48) (cid:57) C z (cid:48) . This concludes the proof. (cid:3)
Remark 5.19 (anisotropic meshes) . Examining the proof of Proposition 5.14, werealize that a critical part of (5.18) consists of the application of inequality (5.1),which we recall reads: h Y ≤ C T h z (cid:48) for all z (cid:48) ∈ N ( T Ω ). Therefore, Proposition 5.14shows how the resolution of local problems on cylindrical stars allows for anisotropicmeshes on the extended variable y and graded meshes in Ω. The latter enables usto compensate possible singularities in the x (cid:48) -variables. Remark 5.20 (relaxing the mesh condition) . Owing to the anisotropy of the mesh,condition (5.1) can be violated only near the top of the cylinder. Near the bottomof the cylinder, the size of the elements will be much smaller than h z (cid:48) . If one couldprove that the error v − V T Y decays exponentially (as the exact solution v does)an examination of the proof of Proposition 5.14 reveals that condition (5.1) canbe removed. Proving this decay, however, requires local pointwise error estimateswhich are not availabe and are currently under investigation. POSTERIORI ERROR ANALYSIS FOR FRACTIONAL DIFFUSION 15
A computable a posteriori error estimator.
Although Proposition 5.14shows that the error estimator ˜ E T Y is ideal, it has an insurmountable drawback: foreach node z (cid:48) ∈ N ( T Ω ), it requires knowledge of the exact solution ζ z (cid:48) to the localproblem (5.11) which lies in the infinite dimensional space W ( C z (cid:48) ). This makes thisestimator not computable. However, it provides intuition and establishes the basisto define a discrete and computable error estimator. To achieve this, let us nowdefine local discrete spaces and local computable error indicators, on the basis ofwhich we will construct our global error estimator. Definition 5.21 (discrete local spaces) . For z (cid:48) ∈ N ( T Ω ), define the discrete space W ( C z (cid:48) ) = (cid:8) W ∈ C ( ¯ C Y ) : W | T ∈ P ( K ) ⊗ P ( I ) ∀ T = K × I ∈ C z (cid:48) ,W | ∂ C z (cid:48) \ Ω ×{ } = 0 (cid:9) . where, if K is a quadrilateral, P ( K ) stands for Q ( K ) — the space of polynomialsof degree not larger than 2 in each variable. If K is a simplex, P ( K ) correspondsto P ( K ) ⊕ B ( K ) where where P ( K ) stands for the space of polynomials of totaldegree at most 2, and B ( K ) is the space spanned by a local cubic bubble function.We then define the discrete local problems: for each cylindrical star C z (cid:48) we define η z (cid:48) ∈ W ( C z (cid:48) ) to be the solution of(5.22) ˆ C z (cid:48) y α ∇ η z (cid:48) ∇ W = d s (cid:104) f, tr Ω W (cid:105) H − s (Ω) × H s (Ω) − ˆ C z (cid:48) y α ∇ V T Y ∇ W, for all W ∈ W ( C z (cid:48) ). We also define the global error estimator(5.23) E T Y = (cid:88) z (cid:48) ∈ N ( T Ω ) E z (cid:48) , in terms of the local error indicators(5.24) E z (cid:48) = (cid:57) η z (cid:48) (cid:57) C z (cid:48) . We next explore the connection between the estimator (5.23) and the error. Wefirst prove a local lower bound for the error without any oscillation term and freeof any constant.
Theorem 5.25 (localized lower bound) . Let v ∈ ◦ H L ( y α , C Y ) and V T Y ∈ V ( T Y ) solve (3.2) and (3.8) respectively. Then, for any z (cid:48) ∈ N ( T Ω ) , we have (5.26) E z (cid:48) ≤ (cid:107)∇ ( v − V T Y ) (cid:107) L ( y α , C z (cid:48) ) . Proof.
The proof repeats the arguments employed to obtain inequality (5.16). Let z (cid:48) ∈ N ( T Ω ), and let η z (cid:48) and E z (cid:48) be as in (5.22) and (5.24). Then, E z (cid:48) = (cid:57) η z (cid:48) (cid:57) C z (cid:48) = ˆ C z (cid:48) y α ∇ η z (cid:48) ∇ η z (cid:48) = ˆ C z (cid:48) y α ∇ e T Y ∇ η z (cid:48) ≤ (cid:57) e T Y (cid:57) C z (cid:48) (cid:57) η z (cid:48) (cid:57) C z (cid:48) , which concludes the proof. (cid:3) Remark 5.27 (strong efficiency) . The oscillation and constant free lower bound(5.26) implies a strong concept of efficiency: the relative size of E z (cid:48) dictates meshrefinement regardless of fine structure of the data f , and thus works even in thepre-asymptotic regime.To obtain an upper bound we must assume the following. Conjecture 5.28 (operator M z (cid:48) ) . For every z (cid:48) ∈ N ( T Ω ) there is a linear operator M z (cid:48) : W ( C z (cid:48) ) → W ( C z (cid:48) ) such that, for all w ∈ W ( C z (cid:48) ) , satisfies: • For every cell K ∈ T Ω such that K ⊂ S z (cid:48) (5.29) ˆ K ×{ } tr Ω ( w − M z (cid:48) w ) = 0 , • For every cell T ⊂ C z (cid:48) and every W ∈ V ( T Y )(5.30) ˆ T y α ∇ ( w − M z (cid:48) w ) ∇ W = 0 . • Stability (5.31) (cid:57) M z (cid:48) w (cid:57) C z (cid:48) (cid:46) (cid:57) w (cid:57) C z (cid:48) , where the hidden constant is independent of the discretization parameters but maydepend on α . We next introduce the so-called data oscillation . For every z (cid:48) ∈ N ( T Ω ), wedefine the local data oscillation as(5.32) osc z (cid:48) ( f ) := d s h sz (cid:48) (cid:107) f − f z (cid:48) (cid:107) L ( S z (cid:48) ) where f z (cid:48) | K ∈ R is the average of f over K , i.e.,(5.33) f z (cid:48) | K := K f. The global data oscillation is defined as(5.34) osc T Ω ( f ) := (cid:88) z (cid:48) ∈ N ( T Ω ) osc z (cid:48) ( f ) . We also define the total error indicator (5.35) τ T Ω ( V T Y , S z (cid:48) ) := (cid:0) E z (cid:48) + osc z (cid:48) ( f ) (cid:1) / ∀ z (cid:48) ∈ N ( T Ω ) , which will be used to mark elements for refinement in the adaptive finite elementmethod proposed in Section 6. Let K T Ω = { S z (cid:48) : z (cid:48) ∈ N ( T Ω ) } and, for any M ⊂ K T Ω , we set(5.36) τ T Ω ( V T Y , M ) := (cid:88) S z (cid:48) ∈ M τ T Ω ( V T Y , S z (cid:48) ) / . With the aid of the operators M z (cid:48) and under the assumption that Conjec-ture 5.28 holds, we can bound the error by the estimator, up to oscillation terms. Theorem 5.37 (global upper bound) . Let v ∈ ◦ H L ( C Y , y α ) and V T Y ∈ V ( T Y ) solve (3.2) and (3.8) , respectively. If f ∈ L (Ω) and Conjecture 5.28 holds, then the totalerror estimator τ T Ω ( V T Y , K T Ω ) , defined in (5.36) satisfies (5.38) (cid:107)∇ ( v − V T Y ) (cid:107) L ( y α , C Y ) (cid:46) τ T Ω ( V T Y , K T Ω ) . Proof.
Let e T Y = v − V T Y denote the error and let ψ z (cid:48) = ( w − Π T Y w ) ˜ ϕ z (cid:48) ∈ W ( C z (cid:48) ),for any w ∈ ◦ H L ( y α , C Y ), where Π T Y is the quasi-interpolation operator introduced POSTERIORI ERROR ANALYSIS FOR FRACTIONAL DIFFUSION 17 in [29, §
4] and [30, § (cid:57) ψ z (cid:48) (cid:57) C z (cid:48) (cid:46) (cid:57) w (cid:57) D z (cid:48) obtained as partof the proof of Proposition 5.14. Following the proof of Proposition 5.14, we have ˆ C Y y α ∇ e T Y ∇ w = (cid:88) z (cid:48) ∈ N ( T Ω ) ˆ C z (cid:48) y α ∇ e T Y ∇ ψ z (cid:48) = (cid:88) z (cid:48) ∈ N ( T Ω ) ˆ C z (cid:48) y α ∇ e T Y ∇M z (cid:48) ψ z (cid:48) − (cid:88) z (cid:48) ∈ N ( T Ω ) ˆ C z (cid:48) y α ∇ e T Y ∇ ( ψ z (cid:48) − M z (cid:48) ψ z (cid:48) ) . We now examine each term separately. First, for every z (cid:48) ∈ N ( T Ω ) we have M z (cid:48) ψ z (cid:48) ∈ W ( C z (cid:48) ), whence the definition of the discrete local problem (5.22) yields (cid:88) z (cid:48) ∈ N ( T Ω ) ˆ C z (cid:48) y α ∇ e T Y ∇M z (cid:48) ψ z (cid:48) = (cid:88) z (cid:48) ∈ N ( T Ω ) ˆ C z (cid:48) y α ∇ η z (cid:48) ∇M z (cid:48) ψ z (cid:48) ≤ (cid:88) z (cid:48) ∈ N ( T Ω ) E z (cid:48) / (cid:88) z (cid:48) ∈ N ( T Ω ) (cid:57) ψ z (cid:48) (cid:57) C z (cid:48) / (cid:46) E T Y (cid:107)∇ w (cid:107) L ( y α , C Y ) , where in the last inequality we used the stability assumption (5.31) of the operator M z (cid:48) , the inequality (cid:57) ψ z (cid:48) (cid:57) C z (cid:48) (cid:46) (cid:57) w (cid:57) D z (cid:48) , and the finite overlapping property of thestars C z (cid:48) .Second, for any z (cid:48) ∈ N ( T Ω ), we use the conditions (5.29) and (5.30) imposed onthe operator M z (cid:48) to derive ˆ C z (cid:48) y α ∇ e T Y ∇ ( ψ z (cid:48) − M z (cid:48) ψ z (cid:48) ) = d s ˆ S z (cid:48) ( f − f z (cid:48) ) tr Ω ( ψ z (cid:48) − M z (cid:48) ψ z (cid:48) ) , where we used that ∇ V T Y is constant on every T . Moreover, since f z (cid:48) is the L ( S z (cid:48) )projection onto piecewise constants of f we have that, for any (cid:37) such that (cid:37) | K ∈ R ˆ S z (cid:48) ( f − f z (cid:48) ) tr Ω ( ψ z (cid:48) − M z (cid:48) ψ z (cid:48) ) = ˆ S z (cid:48) ( f − f z (cid:48) ) (tr Ω ( ψ z (cid:48) − M z (cid:48) ψ z (cid:48) ) − (cid:37) ) ≤ (cid:107) f − f z (cid:48) (cid:107) L ( S z (cid:48) ) (cid:107) tr Ω ( ψ z (cid:48) − M z (cid:48) ψ z (cid:48) ) − (cid:37) (cid:107) L ( S z (cid:48) ) . After suitably choosing (cid:37) , and using a standard interpolation-type estimate, we get ˆ S z (cid:48) ( f − f z (cid:48) ) tr Ω ( ψ z (cid:48) − M z (cid:48) ψ z (cid:48) ) (cid:46) h sz (cid:48) (cid:107) f − f z (cid:48) (cid:107) L ( S z (cid:48) ) (cid:107) tr Ω ( ψ z (cid:48) − M z (cid:48) ψ z (cid:48) ) (cid:107) H s (Ω) . Consequently, (cid:88) z (cid:48) ∈ N ( T Ω ) ˆ C z (cid:48) y α ∇ e T Y ∇ ( ψ z (cid:48) − M z (cid:48) ψ z (cid:48) ) ≤ (cid:88) z (cid:48) ∈ N ( T Ω ) C tr Ω d s h sz (cid:48) (cid:107) f − f z (cid:48) (cid:107) L ( S z (cid:48) ) (cid:57) ψ z (cid:48) − M z (cid:48) ψ z (cid:48) (cid:57) C z (cid:48) (cid:46) osc T Ω ( f ) (cid:107)∇ w (cid:107) L ( y α , C Y ) , where we applied the trace inequality (2.9), the estimate (2.10) on C tr Ω , the sta-bility assumption (5.31) of M z (cid:48) , the bound (cid:57) ψ z (cid:48) (cid:57) C z (cid:48) (cid:46) (cid:57) w (cid:57) D z (cid:48) , and the finiteoverlapping property of the stars C z (cid:48) .Collecting the above estimates, we obtain the asserted bound (5.38). (cid:3) Remark 5.39 (role of oscillation) . The Definition 5.21 of W ( C z (cid:48) ) is meant toprovide enough degrees of freedom for existence of the operator M z (cid:48) satisfying(5.29)–(5.31). This leads to a solution free oscillation term (5.32). Otherwise, if we were not able to impose (5.30), then the oscillation (5.32) should be supplementedby the term osc z (cid:48) ( V T Y ) := (cid:107) y α ∇ V T Y − σ z (cid:48) (cid:107) L ( y − α , C z (cid:48) ) , with σ z (cid:48) being the local average of y α ∇ V T Y , for (5.38) to be valid. This termcannot be guaranteed to be of higher order due to the presence of the weight y − α .In fact, computations show that, for s > (or α < osc z (cid:48) ( V T Y )is of lower order than the actual error estimator E z (cid:48) unless a stronger mesh gradingin the extended direction is employed to control this term. This can be achieved,for instance, by taking the largest grading parameter γ in (3.5) for ± α , namely γ > / (1 − | α | ). Numerical experiments show no degradation of the convergencerate, expressed in terms of degrees of freedom, for the ensuing meshes T Y withstronger anisotropic refinement. Remark 5.40 (construction of M z (cid:48) ) . The construction of the operator M z (cid:48) isan open problem. We design the local space W ( C z (cid:48) ) in order to provide enoughdegrees of freedom for the existence of the operator M z (cid:48) satisfying (5.29)–(5.31).The numerical experiments of Section 6 provide consistent computational evidencethat the upper bound (5.38) is valid without osc z (cid:48) ( V T Y ), and thus indirect evidenceof the existence of M z (cid:48) with the requisite properties (5.29)–(5.31).6. Numerical experiments
Here we explore computationally the performance and limitations of the a poste-riori error estimator introduced in § SOLVE → ESTIMATE → MARK → REFINE . Design of AFEM.
The modules in (6.1) are as follows: • SOLVE : Given a mesh T Y we compute the Galerkin solution of (3.2): V T Y = SOLVE ( T Y ) . • ESTIMATE : Given V T Y we calculate the local error indicators (5.24) and the localoscillations (5.32) to construct the total error indicator of (5.35): (cid:8) τ T Ω ( V T Y , S z (cid:48) ) (cid:9) S z (cid:48) ∈ K T Ω = ESTIMATE ( V T Y , T Y ) . • MARK : Using the so-called D¨orfler marking [13] (bulk chasing strategy) withparameter 0 < θ ≤
1, we select a set M = MARK ( (cid:8) τ T Ω ( V T Y , S z (cid:48) ) (cid:9) S z (cid:48) ∈ K Ω , V T Y ) ⊂ K T Ω of minimal cardinality that satisfies τ T Ω ( V T Y , M ) ≥ θτ T Ω ( V T Y , K T Ω ) . • REFINE : We generate a new mesh T (cid:48) Ω by bisecting all the elements K ∈ T Ω contained in M based on the newest-vertex bisection method; see [32, 31]. Wechoose the truncation parameter as Y = 1 + log( T (cid:48) Ω ) to balance the approx-imation and truncation errors; see [29, Remark 5.5]. We construct the mesh I (cid:48) Y by the rule (3.5), with a number of degrees of freedom M sufficiently large sothat condition (5.1) holds. This is attained by first creating a partition I (cid:48) Y with POSTERIORI ERROR ANALYSIS FOR FRACTIONAL DIFFUSION 19 M ≈ ( T (cid:48) Ω ) /n and checking (5.1). If this condition is violated, we increase thenumber of points until we get the desired result. The new mesh T (cid:48) Y = REFINE ( M ) , is obtained as the tensor product of T (cid:48) Ω and I (cid:48) Y .6.2. Implementation.
The AFEM (6.1) is implemented within the MATLAB © software library i FEM [11]. The stiffness matrix of the discrete system (3.8) isassembled exactly, and the forcing boundary term is computed by a quadratureformula which is exact for polynomials of degree 4. The resulting linear system issolved by using the multigrid method with line smoother introduced and analyzedin [12].To compute the solution η z (cid:48) to the discrete local problem (5.22), we loop aroundeach node z (cid:48) ∈ N ( T Ω ), collect data about the cylindrical star C z (cid:48) and assem-ble the small linear system (5.22) which is solved by the built-in direct solver ofMATLAB © . All integrals involving only the weight and discrete functions are com-puted exactly, whereas those also involving data functions are computed element-wise by a quadrature formula which is exact for polynomials of degree 7.For convenience, in the MARK step we change the estimator from star-wise toelement-wise as follows: We first scale the nodal-wise estimator as E z (cid:48) / ( S z (cid:48) ) andthen, for each element K ∈ T Ω , we compute E K := (cid:88) z (cid:48) ∈ K E z (cid:48) . The scaling is introduced so that (cid:80) K ∈ T Ω E K = (cid:80) z (cid:48) ∈ N ( T Ω ) E z (cid:48) . The cell-wise dataoscillation is now defined as osc K ( f ) := d s h sK (cid:107) f − ¯ f K (cid:107) L ( K ) , where ¯ f K denotes the average of f over the element K . This quantity is computedusing a quadrature formula which is exact for polynomials of degree 7.Unless specifically mentioned, all computations are done without explicitly en-forcing the mesh restriction (5.1), which shows that this is nothing but an artifactin our proofs; see Remark 5.20. How to remove this assumption is currently underinvestigation. Nevertheless, computations (not shown here for brevity) show thatoptimality is still retained if one imposes (5.1).For the cases where the exact solution to problem (3.8) is available, the error iscomputed by using the identity (cid:107)∇ ( v − V T Y ) (cid:107) L ( y α , C Y ) = d s ˆ Ω f tr Ω ( v − V T Y ) , which follows from Galerkin orthogonality and integration by parts. Thus, we avoidevaluating the singular/degenerate weight y α and reduce the computational cost.The right hand side of the equation above is computed by a quadrature formulawhich is exact for polynomials of degree 7. On the other hand, if the exact solutionis not available, we introduce the energy E ( w ) = 12 ˆ C Y y α |∇ w | − d s (cid:104) f, tr Ω w (cid:105) H − s (Ω) × H s (Ω) . where w ∈ ◦ H L ( y α , C Y ) and f ∈ H − s (Ω). Consequently, Galerkin orthogonalityimplies E ( V T Y ) − E ( v ) = 12 (cid:107)∇ ( v − V T Y ) (cid:107) L ( y α , C Y ) , where v ∈ ◦ H L ( y α , C Y ) and V T Y ∈ V ( T Y ) denote the solution to problems (3.2) and(3.8), respectively. We remark that for a discrete function W the energy E ( W ) canbe computed simply using a matrix-vector product. Then, an approximation of theerror in the energy norm can be obtained by computing (cid:16) E ( V T Y ) − E ( V T ∗ Y )) (cid:17) , where V T ∗ Y is the Galerkin approximation to v on a very fine grid T ∗ Y , which isobtained by a uniform refinement of the last adaptive mesh.We also compute the effectivity index of our a posterior error estimator definedas the ratio of the estimator and the true error, i.e., τ T Ω ( V T Y , K T Ω ) (cid:107)∇ ( v − V T Y ) (cid:107) L ( y α , C Y ) , as well as the aspect ratio of elements Th K h I ∀ T = K × I ∈ T Y . Smooth and compatible data.
The purpose of this numerical example is toshow how the error estimator (5.23)–(5.24) based on the solution of local problemson anisotropic cylindrical stars allows us to recover optimality of our AFEM. Werecall that adaptive isotropic refinement cannot be optimal; see § , and f ( x , x ) = sin(2 πx ) sin(2 πx ). Then, the exact solution to (1.1)is given by u ( x , x ) = 8 − s π s sin(2 πx ) sin(2 πx ).The asymptotic relation (cid:107)∇ ( U − V T Y ) (cid:107) L ( y α , C Y ) ≈ T Y ) − / is shown in Fig-ure 3 which illustrates the quasi-optimal decay rate of our AFEM driven by theerror estimator (5.23)–(5.24) for all choices of the parameter s considered. Theseexamples show that anisotropy in the extended dimension is essential to recoveroptimality of our AFEM.6.4. Smooth but incompatible data.
The numerical example presented in [29, § f ∈ H − s (Ω) is indeed necessary to obtain an optimal rate of convergence with aquasiuniform mesh in the x (cid:48) -direction. The heuristic explanation for this is thata certain compatibility between the data and the boundary condition is necessary.The results of [29, § , and f = 1. Then, for s < , we have that f / ∈ H − s (Ω)whence we cannot invoke the results of Theorem 3.9. Nevertheless, as the resultsof Figure 4 show, we recover the optimal rate of convergence. POSTERIORI ERROR ANALYSIS FOR FRACTIONAL DIFFUSION 21 −2 −1 Degrees of Freedom (DOFs) E rr o r DOFs − s = 0 . s = 0 . s = 0 . s = 0 . −1 Degrees of Freedom (DOFs) E s t i m a t o r DOFs − s = 0 . s = 0 . s = 0 . s = 0 . Figure 3.
Computational rate of convergence for our anisotropicAFEM on the smooth and compatible right hand side of § n = 2 and s = 0 .
2, 0 .
4, 0 . s = 0 .
8. The left panel showsthe decrease of the error with respect to the number of degrees offreedom, whereas the right one that for the total error estimator. Inall cases we recover the optimal rate T Y ) − / . The aspect ratios,averaged over x (cid:48) , of the cells on the bottom layer [0 , y ] in the finestmesh are: 1 . × , 6 . × , 396 and 36 .
2, respectively. Theaverage effectivity indices are 1 . , . , . , .
62, respectively. −1 Degrees of Freedom (DOFs) E rr o r DOFs − s = 0 . s = 0 . s = 0 . s = 0 . −1 Degrees of Freedom (DOFs) E s t i m a t o r DOFs − s = 0 . s = 0 . s = 0 . s = 0 . Figure 4.
Computational rate of convergence for our anisotropicAFEM on the smooth but incompatible right hand side of § n = 2 and s = 0 .
2, 0 .
4, 0 . s = 0 .
8. The left panel showsthe decrease of the error with respect to the number of degrees offreedom, whereas the right one that for the total error indicator.In all cases we recover the optimal rate T Y ) − / . Notice that,for s < , the right hand side f = 1 / ∈ H − s (Ω) and so a quasiu-niform mesh in Ω does not deliver the optimal rate of convergence[29, § x (cid:48) , of the cells on thebottom layer [0 , y ] in the finest mesh are: 2 . × , 6 . × ,387 and 33 .
4, respectively. The average effectivity indices are 1 . .
41, 1 .
51, 1 .
67, respectively.6.5.
L-shaped domain with compatible data.
The result of Theorem 3.9 re-lies on the assumption that the Laplace operator − ∆ x (cid:48) on Ω supplemented with Dirichlet boundary conditions possesses optimal smoothing properties, i.e., (3.1).As it is well known [18], (3.1) holds when Ω is a convex polyhedron. Let us thenconsider the case when this condition is not met.We consider the classical L-shaped domain Ω = ( − , \ (0 , × ( − ,
0) and f ( x , x ) = sin(2 πx ) sin( πx ), which is a smooth and compatible right hand side,i.e., f = 0 on ∂ Ω, for all values of s ∈ (0 , −1 Degrees of Freedom (DOFs) E rr o r DOFs − s = 0 . s = 0 . s = 0 . s = 0 . −1 Degrees of Freedom (DOFs) E s t i m a t o r DOFs − s = 0 . s = 0 . s = 0 . s = 0 . Figure 5.
Computational rate of convergence for our anisotropicAFEM on the smooth and compatible right hand side over anL-shaped domain of § n = 2 and s = 0 .
2, 0 .
4, 0 . s = 0 .
8. The left panel shows the decrease of the error with respectto the number of degrees of freedom, whereas the right one that forthe total error indicator. In all cases we recover the optimal rate T Y ) − / . As the exact solution of this problem is not known,we compute the rate of convergence by comparing the computedsolution with one obtained on a very fine mesh. The aspect ratios,averaged over x (cid:48) , of the cells on the bottom layer [0 , y ] in the finestmesh are: 2 . × , 1 . × , 667 and 56 .
7, respectively. Theaverage effectivity indices are 1 .
73, 1 .
77, 1 .
73, 1 .
74, respectively.6.6.
L-shaped domain with incompatible data.
Let us combine the singularityintroduced by the data incompatibility of § § − , \ (0 , × ( − , f = 1. As the results of Figure 6 show, we again recover the optimal rate ofconvergence for all possible cases of s . As the exact solution of this problem is notknown, we compute the rate of convergence by comparing the computed solutionwith one obtained on a very fine mesh.We display some meshes in Figure 7. As expected, when s < the data f = 1 isincompatible with the equation and this causes a boundary layer on the solution.To capture it, the AFEM refines near the boundary. In contrast, when s > the refinement is only near the reentrant corner, since there is no boundary layeranymore. POSTERIORI ERROR ANALYSIS FOR FRACTIONAL DIFFUSION 23 −1 Degrees of Freedom (DOFs) E rr o r DOFs − s = 0 . s = 0 . s = 0 . s = 0 . −1 Degrees of Freedom (DOFs) E s t i m a t o r DOFs − s = 0 . s = 0 . s = 0 . s = 0 . Figure 6.
Computational rate of convergence for our anisotropicAFEM on the smooth but incompatible right hand side over anL-shaped domain of § n = 2 and s = 0 .
2, 0 .
4, 0 . s = 0 .
8. The left panel shows the decrease of the error with respectto the number of degrees of freedom, whereas the right one that forthe total error indicator. In all cases we recover the optimal rate T Y ) − / . As the exact solution of this problem is not known,we compute the rate of convergence by comparing the computedsolution with one obtained on a very fine mesh. The aspect ratios,averaged over x (cid:48) , of the cells on the bottom layer [0 , y ] in the finestmesh are: 2 . × , 1 . × , 632 and 55 .
6, respectively. Theaverage effectivity indices are 1 .
36, 1 .
40, 1 .
44, 1 .
49, respectively.
Figure 7.
Adaptively graded mesh for an L-shaped domain withincompatible data (see § s = 0 . s = 0 . s < the data f = 1 is incompatible withthe equation and this causes a boundary layer on the solution. Tocapture it, our AFEM refines near the boundary. In contrast, when s > the refinement is only near the reentrant angle, since thereis no boundary layer anymore.6.7. Discontinuous coefficients.
In this example we compute the fractional pow-ers of a general elliptic operator with piecewise discontinuous coefficients. We in-voke the formulas derived by Kellogg [19] (see also [32, § − , , f ( x , x ) = (( x ) − x ) −
1) and L w = − div x (cid:48) ( A ∇ x (cid:48) w ) , with A = (cid:40) (cid:37) I , x x > , I , x x ≤ , where I denotes the identity tensor and (cid:37) = 161 . barely in H (Ω).On the other hand, as the results of [29, §
7] show, the problem: find u such that L s u = f, in Ω u | ∂ Ω = 0also admits an extension, which can be discretized with the techniques and ideaspresented in [29, § T Y ) − / for all the considered cases. −1 Degrees of Freedom (DOFs) E rr o r DOFs − s = 0 . s = 0 . s = 0 . s = 0 . −1 Degrees of Freedom (DOFs) E s t i m a t o r DOFs − s = 0 . s = 0 . s = 0 . s = 0 . Figure 8.
Computational rate of convergence for our anisotropicAFEM on the problem with discontinuous coefficients describedin § n = 2 and s = 0 .
2, 0 .
4, 0 . s = 0 .
8. The leftpanel shows the decrease of the error with respect to the numberof degrees of freedom, whereas the right one that for the total errorindicator. In all cases we recover the optimal rate T Y ) − / . Asthe exact solution of this problem is not known, we compute therate of convergence by comparing the computed solution with oneobtained on a very fine mesh. The aspect ratios, averaged over x (cid:48) , of the cells on the bottom layer [0 , y ] in the finest mesh are:3 . × , 1 . × , 719 and 69 .
1, respectively. The averageeffectivity indices are 1 .
55, 1 .
64, 1 .
69, 1 .
71, respectively.6.8.
The role of oscillation.
Let us explore the role of oscillation as discussed inRemark 5.39. We consider the same example as in § P ( K ) = P ( K ). In this case the discrete local space W ( C z (cid:48) ) does not have enoughdegrees of freedom to impose (5.30). Consequently, in order to have (5.38) theoscillation (5.32) must be supplemented by the term(6.2) osc z (cid:48) ( V T Y ) = (cid:107) y α ∇ V T Y − σ z (cid:48) (cid:107) L ( C z (cid:48) ,y − α ) , where σ z (cid:48) is the local average of y α ∇ V T Y .Figure 9 shows the experimental rates of convergence obtained by our AFEMfor the total error where the oscillation term (5.32) is supplemented with (6.2). Aswe can see, especially for s = 0 .
8, the results are not optimal. To remedy this we
POSTERIORI ERROR ANALYSIS FOR FRACTIONAL DIFFUSION 25 −1 Degrees of Freedom (DOFs) E s t i m a t o r DOFs − s = 0 . s = 0 . s = 0 . s = 0 . −1 Degrees of Freedom (DOFs) E s t i m a t o r DOFs − s = 0 . s = 0 . s = 0 . s = 0 . Figure 9.
Experimental rate of convergence for the example of § P ( K ) = P ( K ) in Definition 5.21. In this case,the oscillation (5.32) must be supplemented by (6.2) to attain anupper bound. The expression (6.2) possesses a singularity andweight which cannot be captured with the mesh grading necessaryfor optimal approximation orders; this yields suboptimal decayrates (left figure). However, a graded mesh with γ > / (1 − | α | ) in(3.5) gives optimal decay rates (right figure).notice that a graded mesh is able to capture the singular behavior of (6.2). Indeed,the underlying weight in this expression is y − α , so that a grading parameter in (3.5)of γ > / (1 − | α | ) would yield an optimal decay for (6.2). For this reason, setting γ > / (1 − | α | ), we are able to obtain an optimal decay rates for both the errorand the oscillation (6.2) and all values of s . The results shown in Figure 9 confirmthis. References [1] M. Ainsworth and J.T. Oden.
A posteriori error estimation in finite element analysis . Pureand Applied Mathematics (New York). Wiley-Interscience, New York, 2000.[2] I. Babuˇska and A. Miller. A feedback finite element method with a posteriori error estimation.I. The finite element method and some basic properties of the a posteriori error estimator.
Comput. Methods Appl. Mech. Engrg. , 61(1):1–40, 1987.[3] I. Babuˇska and W.C. Rheinboldt. Error estimates for adaptive finite element computations.
SIAM J. Numer. Anal. , 15(4):736–754, 1978.[4] C. Br¨andle, E. Colorado, A. de Pablo, and U. S´anchez. A concave–convex elliptic probleminvolving the fractional laplacian.
Proceedings of the Royal Society of Edinburgh, Section: AMathematics , 143:39–71, 2013.[5] X. Cabr´e and Y. Sire. Nonlinear equations for fractional Laplacians ii: Existence, uniquenessand qualitative properties of solutions. arXiv:1111.0796v1, 2011.[6] X. Cabr´e and J. Tan. Positive solutions of nonlinear problems involving the square root ofthe Laplacian.
Adv. Math. , 224(5):2052–2093, 2010.[7] L. Caffarelli and L. Silvestre. An extension problem related to the fractional Laplacian.
Comm. Part. Diff. Eqs. , 32(7-9):1245–1260, 2007.[8] A. Capella, J. D´avila, L. Dupaigne, and Y. Sire. Regularity of radial extremal solutions forsome non-local semilinear equations.
Comm. Part. Diff. Eqs. , 36(8):1353–1384, 2011.[9] C. Carstensen and S.A. Funken. Fully reliable localized error control in the FEM.
SIAM J.Sci. Comput. , 21(4):1465–1484 (electronic), 1999/00.[10] J.M. Casc´on and R.H. Nochetto. Quasioptimal cardinality of AFEM driven by nonresidualestimators.
IMA J. Numer. Anal. , 32(1):1–29, 2012. [11] L. Chen. i FEM: An integrated finite element methods package in matlab. Technical report,University of California at Irvine, 2009.[12] L. Chen, R.H. Nochetto, E. Ot´arola, and A.J. Salgado. Multilevel methods for nonuniformlyelliptic operators. arXiv:1403.4278, 2014.[13] W. D¨orfler. A convergent adaptive algorithm for Poisson’s equation.
SIAM J. Numer. Anal. ,33(3):1106–1124, 1996.[14] R.G. Dur´an and A.L. Lombardi. Error estimates on anisotropic Q elements for functions inweighted Sobolev spaces. Math. Comp. , 74(252):1679–1706 (electronic), 2005.[15] E. B. Fabes, C.E. Kenig, and R.P. Serapioni. The local regularity of solutions of degenerateelliptic equations.
Comm. Part. Diff. Eqs. , 7(1):77–116, 1982.[16] L. Formaggia and S. Perotto. Anisotropic error estimates for elliptic problems.
Numer. Math. ,94(1):67–92, 2003.[17] V. Gol (cid:48) dshtein and A. Ukhlov. Weighted Sobolev spaces and embedding theorems.
Trans.Amer. Math. Soc. , 361(7):3829–3850, 2009.[18] P. Grisvard.
Elliptic problems in nonsmooth domains , volume 24 of
Monographs and Studiesin Mathematics . Pitman (Advanced Publishing Program), Boston, MA, 1985.[19] R.B. Kellog. On the Poisson equation with intersecting interfaces.
Applicable Anal. , 4:101–129, 1974/75.[20] A. Kufner and B. Opic. How to define reasonably weighted Sobolev spaces.
Comment. Math.Univ. Carolin. , 25(3):537–554, 1984.[21] G. Kunert. An a posteriori residual error estimator for the finite element method onanisotropic tetrahedral meshes.
Numer. Math. , 86(3):471–490, 2000.[22] G. Kunert. A local problem error estimator for anisotropic tetrahedral finite element meshes.
SIAM J. Numer. Anal. , 39(2):668–689 (electronic), 2001.[23] G. Kunert and R. Verf¨urth. Edge residuals dominate a posteriori error estimates for linearfinite element methods on anisotropic triangular and tetrahedral meshes.
Numer. Math. ,86(2):283–303, 2000.[24] N.S. Landkof.
Foundations of modern potential theory . Springer-Verlag, New York, 1972.[25] J.-L. Lions and E. Magenes.
Non-homogeneous boundary value problems and applications.Vol. I . Springer-Verlag, New York, 1972.[26] P. Morin, R.H. Nochetto, and K.G. Siebert. Local problems on stars: a posteriori errorestimators, convergence, and performance.
Math. Comp. , 72(243):1067–1097, 2003.[27] B. Muckenhoupt. Weighted norm inequalities for the Hardy maximal function.
Trans. Amer.Math. Soc. , 165:207–226, 1972.[28] S. Nicaise. A posteriori error estimations of some cell centered finite volume methods fordiffusion-convection-reaction problems.
SIAM J. Numer. Anal. , 44(3):949–978, 2006.[29] R.H. Nochetto, E. Ot´arola, and A.J. Salgado. A pde approach to fractional diffusionin general domains: A priori error analysis.
Found. Comput. Math. , pages 1–59, 2014.DOI:10.1007/s10208-014-9208-x.[30] R.H. Nochetto, E. Ot´arola, and A.J. Salgado. Piecewise polynomial interpolation in Muck-enhoupt weighted Sobolev spaces and applications. arXiv:1402.1916, 2014.[31] R.H. Nochetto, K.G. Siebert, and A. Veeser. Theory of adaptive finite element methods: anintroduction. In
Multiscale, nonlinear and adaptive approximation , pages 409–542. Springer,Berlin, 2009.[32] R.H. Nochetto and A. Veeser. Primer of adaptive finite element methods. In
Multiscale andadaptivity: modeling, numerics and applications , volume 2040 of
Lecture Notes in Math. ,pages 125–225. Springer, Heidelberg, 2012.[33] M. Picasso. An anisotropic error indicator based on Zienkiewicz-Zhu error estimator: appli-cation to elliptic and parabolic problems.
SIAM J. Sci. Comput. , 24(4):1328–1355, 2003.[34] K.G. Siebert. An a posteriori error estimator for anisotropic refinement.
Numer. Math. ,73(3):373–398, 1996.[35] P.R. Stinga and J.L. Torrea. Extension problem and Harnack’s inequality for some fractionaloperators.
Comm. Part. Diff. Eqs. , 35(11):2092–2122, 2010.[36] B.O. Turesson.
Nonlinear potential theory and weighted Sobolev spaces , volume 1736 of
Lec-ture Notes in Mathematics . Springer-Verlag, Berlin, 2000.[37] R. Verf¨urth.
A Review of A Posteriori Error Estimation and Adaptive Mesh-RefinementTechniques . John Wiley, 1996.
POSTERIORI ERROR ANALYSIS FOR FRACTIONAL DIFFUSION 27 (L. Chen)
Department of Mathematics, University of California at Irvine, Irvine, CA92697, USA
E-mail address : [email protected] (R.H. Nochetto) Department of Mathematics and Institute for Physical Science andTechnology, University of Maryland, College Park, MD 20742, USA
E-mail address : [email protected] (E. Ot´arola) Department of Mathematics, University of Maryland, College Park,MD 20742, USA and Department of Mathematical Sciences, George Mason University,Fairfax, VA 22030, USA.
E-mail address : [email protected] (A.J. Salgado) Department of Mathematics, University of Tennessee, Knoxville, TN37996, USA
E-mail address ::