A PDE approach to fractional diffusion in general domains: a priori error analysis
AA PDE APPROACH TO FRACTIONAL DIFFUSION IN GENERALDOMAINS: A PRIORI ERROR ANALYSIS ∗ RICARDO H. NOCHETTO † , ENRIQUE OT ´AROLA ‡ , AND
ABNER J. SALGADO § Abstract.
The purpose of this work is the study of solution techniques for problems involv-ing fractional powers of symmetric coercive elliptic operators in a bounded domain with Dirichletboundary conditions. These operators can be realized as the Dirichlet to Neumann map for a degen-erate/singular elliptic problem posed on a semi-infinite cylinder, which we analyze in the frameworkof weighted Sobolev spaces. Motivated by the rapid decay of the solution of this problem, we pro-pose a truncation that is suitable for numerical approximation. We discretize this truncation usingfirst degree tensor product finite elements. We derive a priori error estimates in weighted Sobolevspaces. The estimates exhibit optimal regularity but suboptimal order for quasi-uniform meshes.For anisotropic meshes, instead, they are quasi-optimal in both order and regularity. We presentnumerical experiments to illustrate the method’s performance.
Key words.
Fractional diffusion; finite elements; nonlocal operators; degenerate and singularequations; second order elliptic operators; anisotropic elements.
AMS subject classifications.
1. Introduction.
Singular integrals and nonlocal operators have been an activearea of research in different branches of mathematics such as operator theory andharmonic analysis (see [56]). In addition, they have received significant attentionbecause of their strong connection with real-world problems, since they constitute afundamental part of the modeling and simulation of complex phenomena that spanvastly different length scales.Nonlocal operators arise in a number of applications such as: boundary controlproblems [31], finance [21], electromagnetic fluids [48], image processing [36], materialsscience [8], optimization [31], porous media flow [25], turbulence [5], peridynamics [55],nonlocal continuum field theories [32] and others. Therefore the domain of definitionΩ could be rather general.To make matters precise, in this work we shall be interested in fractional powersof the Dirichlet Laplace operator ( − ∆) s , with s ∈ (0 , R n ( n ≥ ∂ Ω. Given s ∈ (0 ,
1) and a smooth enough function f , find u such that (cid:40) ( − ∆) s u = f, in Ω ,u = 0 , on ∂ Ω . (1.1)Our approach, however, is by no means particular to the fractional Laplacian. Insection 7 we will discuss how, with little modification, our developments can be appliedto a general second order, symmetric and uniformly elliptic operator. ∗ This work is supported by NSF grants: DMS-1109325 and DMS-0807811. AJS is also supportedby NSF grant DMS-1008058 and an AMS-Simons Grant. EO is supported by the Conicyt-FulbrightFellowship Beca Igualdad de Oportunidades. † Department of Mathematics and Institute for Physical Science and Technology, University ofMaryland, College Park, MD 20742, USA. [email protected] , ‡ Department of Mathematics, University of Maryland, College Park, MD 20742, USA. [email protected] , § Department of Mathematics, University of Maryland, College Park, MD 20742, USA. [email protected] a r X i v : . [ m a t h . NA ] F e b R.H. Nochetto, E. Ot´arola, A.J. Salgado
The study of boundary value problems involving the fractional Laplacian is im-portant in physical applications where long range or anomalous diffusion is considered.For instance, in the flow in porous media, it is used when modeling the transport ofparticles that experience very large transitions arising from high heterogeneity andvery long spatial autocorrelation (see [10]). In the theory of stochastic processes, thefractional Laplacian is the infinitesimal generator of a stable L´evy process (see [12]).One of the main difficulties in the study of problem (1.1) is that the fractionalLaplacian is a nonlocal operator (see [46, 19, 17]). To localize it, Caffarelli andSilvestre showed in [19] that any power of the fractional Laplacian in R n can berealized as an operator that maps a Dirichlet boundary condition to a Neumann-typecondition via an extension problem on the upper half-space R n +1+ . For a boundeddomain Ω, the result by Caffarelli and Silvestre has been adapted in [20, 14, 57],thus obtaining an extension problem which is now posed on the semi-infinite cylinder C = Ω × (0 , ∞ ). This extension is the following mixed boundary value problem: div ( y α ∇ u ) = 0 , in C , u = 0 , on ∂ L C ,∂ u ∂ν α = d s f, on Ω × { } , (1.2)where ∂ L C = ∂ Ω × [0 , ∞ ) denotes the lateral boundary of C , and ∂ u ∂ν α = − lim y → + y α u y , (1.3)is the the so-called conormal exterior derivative of u with ν being the unit outernormal to C at Ω × { } . The parameter α is defined as α = 1 − s ∈ ( − , . (1.4)Finally, d s is a positive normalization constant which depends only on s ; see [19] fordetails. We will call y the extended variable and the dimension n + 1 in R n +1+ the extended dimension of problem (1.2).The limit in (1.3) must be understood in the distributional sense; see [14, 17, 19]or section 2 for more details. As noted in [19, 20, 57], the fractional Laplacian andthe Dirichlet to Neumann operator of problem (1.2) are related by d s ( − ∆) s u = ∂ u ∂ν α in Ω . Using the aforementioned ideas, we propose the following strategy to find thesolution of (1.1): given a sufficiently smooth function f we solve (1.2), thus obtaininga function u : ( x (cid:48) , y ) ∈ C (cid:55)→ u ( x (cid:48) , y ) ∈ R . Setting u : x (cid:48) ∈ Ω (cid:55)→ u ( x (cid:48) ) = u ( x (cid:48) , ∈ R ,we obtain the solution of (1.1). The purpose of this work is then to make these ideasrigorous and to analyze a discretization scheme, which consists of approximating thesolution of (1.2) via first degree tensor product finite elements. We will show sub-optimal error estimates for quasi-uniform discretizations of (1.2) in suitable weightedSobolev spaces and quasi-optimal error estimates using anisotropic elements.The main advantage of the proposed algorithm is that we solve the local problem(1.2) instead of dealing with the nonlocal operator ( − ∆) s of problem (1.1). However,this comes at the expense of incorporating one more dimension to the problem, and PDE approach to fractional diffusion ⊂ R n ( n > − ∆. However, to have a sufficiently good approximation, this re-quires the solution of a large number of eigenvalue problems which, in general, is verytime consuming. In [41, 42] the authors studied computationally problem (1.1) in theone-dimensional case and with boundary conditions of Dirichlet, Neumann and Robintype, and introduced the so-called matrix transference technique (MTT). Basically,MTT computes a spatial discretization of the fractional Laplacian by first findinga matrix approximation, A , of the Laplace operator (via finite differences or finiteelements) and then computing the s -th power of this matrix. This requires diagonal-ization of A which, again, amounts to the solution of a large number of eigenvalueproblems. For the case Ω = (0 , and s ∈ (1 / , A by writing a numerical scheme in terms of the prod-uct of a function of the matrix and a vector, f ( A ) b , where b is a suitable vector.This product is then approximated by a preconditioned Lanczos method. Under thesame setting, the work [16], makes a computational comparison of three techniquesfor the computation of f ( A ) b : the contour integral method, extended Krylov subspacemethods and the pre-assigned poles and interpolation nodes method.The outline of this paper is as follows. In § § §
3. Here we introduce a truncation of problem (1.2) and study some proper-ties of its solution. Having understood the truncation we proceed, in §
4, to study itsfinite element approximation. We prove interpolation estimates in weighted Sobolevspaces, under mild shape regularity assumptions that allow us to consider anisotropicelements in the extended variable y . Based on the regularity results of § §
5, a priori error estimates for quasi-uniform meshes which exhibit optimalregularity but suboptimal order. To restore optimal decay, we resort to the so-calledprinciple of error equidistribution and construct graded meshes in the extended vari-able y . They in turn capture the singular behavior of the solution to (1.2) and allowus to prove a quasi-optimal rate of convergence with respect to both regularity anddegrees of freedom. In §
6, to illustrate the method’s performance and theory, weprovide several numerical experiments. Finally, in §
2. Notation and preliminaries.
Throughout this work Ω is an open, boundedand connected subset of R n , n ≥
1, with Lipschitz boundary ∂ Ω, unless specifiedotherwise. We define the semi-infinite cylinder C = Ω × (0 , ∞ ) , (2.1)and its lateral boundary ∂ L C = ∂ Ω × [0 , ∞ ) . (2.2) R.H. Nochetto, E. Ot´arola, A.J. Salgado
Given Y >
0, we define the truncated cylinder C Y = Ω × (0 , Y ) . (2.3)The lateral 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, as it plays a special rˆole.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 . The upper half-space in R n +1 will be denoted by R n +1+ = { x = ( x (cid:48) , y ) : x (cid:48) ∈ R n y ∈ R , y > } . Let γ = ( γ , γ ) ∈ R and z ∈ R n +1 , the binary operation (cid:12) : R × R n +1 → R n +1 isdefined by γ (cid:12) z = ( γ z (cid:48) , γ z n +1 ) ∈ R n +1 . (2.4)The relation a (cid:46) b indicates that a ≤ Cb , with a constant C that does not dependon neither a nor b but it might depend on s and Ω. The value of C might change ateach occurrence. Given two objects X and Y in the same category, we write X (cid:44) → Y to indicate the existence of a monomorphism between them. Generally, these will beobjects in some subcategory of the topological vector spaces (metric, normed, Banach,Hilbert spaces) and, in this case, the monomorphism is nothing more than continuousembedding. If X is a vector space, we denote by X (cid:48) its dual. Let us recallsome function spaces; for details the reader is referred to [47, 49, 26, 58]. For 0 < s < | w | H s (Ω) = (cid:90) Ω (cid:90) Ω | w ( x (cid:48) ) − w ( x (cid:48) ) | | x (cid:48) − x (cid:48) | n +2 s d x (cid:48) d x (cid:48) . The Sobolev space H s (Ω) of order s is defined by H s (Ω) = (cid:8) w ∈ L (Ω) : | w | H s (Ω) < ∞ (cid:9) , (2.5)which equipped with the norm (cid:107) u (cid:107) H s (Ω) = (cid:16) (cid:107) u (cid:107) L (Ω) + | u | H s (Ω) (cid:17) , is a Hilbert space. An equivalent construction of H s (Ω) is obtained by restrictingfunctions in H s ( R n ) to Ω (cf. [58, Chapter 34]). The space H s (Ω) is defined as theclosure of C ∞ (Ω) with respect to the norm (cid:107) · (cid:107) H s (Ω) , i.e., H s (Ω) = C ∞ (Ω) H s (Ω) . (2.6)If the boundary of Ω is smooth, an equivalent approach to define fractionalSobolev spaces is given by interpolation in [47, Chapter 1]. Set H (Ω) = L (Ω), PDE approach to fractional diffusion ≤ s ≤ θ = 1 − s for the pair [ H (Ω) , L (Ω)], that is H s (Ω) = (cid:2) H (Ω) , L (Ω) (cid:3) θ . (2.7)Analogously, for s ∈ [0 , \ { } , the spaces H s (Ω) are defined as interpolation spacesof index θ = 1 − s for the pair [ H (Ω) , L (Ω)], in other words H s (Ω) = (cid:2) H (Ω) , L (Ω) (cid:3) θ , θ (cid:54) = . (2.8)The space [ H (Ω) , L (Ω)] is the so-called Lions-Magenes space, H (Ω) = (cid:2) H (Ω) , L (Ω) (cid:3) , which can be characterized as H (Ω) = (cid:26) w ∈ H (Ω) : (cid:90) Ω w ( x (cid:48) )dist( x (cid:48) , ∂ Ω) d x (cid:48) < ∞ (cid:27) , (2.9)see [47, Theorem 11.7]. Moreover, we have the strict inclusion H / (Ω) (cid:36) H / (Ω)because 1 ∈ H / (Ω) but 1 / ∈ H / (Ω). If the boundary of Ω is Lipschitz, thecharacterization (2.9) is equivalent to the definition via interpolation, and definitions(2.7) and (2.8) are also equivalent to definitions (2.5) and (2.6), respectively. To seethis, it suffices to notice that when Ω = R n these definitions yield identical spaces andequivalent norms; see [3, Chapter 7]. Consequently, using the well-known extensionresult of Stein [56] for Lipschitz domains, we obtain the asserted equivalence (see [3,Chapter 7] for details).When the boundary of Ω is Lipschitz, the space C ∞ (Ω) is dense in H s (Ω) ifand only if s ≤ , i.e., H s (Ω) = H s (Ω). If s > , we have that H s (Ω) is strictlycontained in H s (Ω); see [47, Theorem 11.1]. In particular, we have the inclusions H / (Ω) (cid:36) H / (Ω) = H / (Ω). It is important to mention thatthere is not a unique way of defining a nonlocal operator related to the fractionalLaplacian in a bounded domain. A first possibility is to suitably extend the functionsto the whole space R n and use Fourier transform F (( − ∆) s w )( ξ (cid:48) ) = | ξ (cid:48) | s F ( w )( ξ (cid:48) ) . After extension, the following point-wise formula also serves as a definition of thefractional Laplacian( − ∆) s w ( x (cid:48) ) = C n,s v.p. (cid:90) R n w ( x (cid:48) ) − w ( z (cid:48) ) | x (cid:48) − z (cid:48) | n +2 s d z (cid:48) , (2.10)where v.p. stands for the Cauchy principal value and C n,s is a positive normalizationconstant that depends only on n and s which is introduced to guarantee that thesymbol of the resulting operator is | ξ (cid:48) | s . For details we refer the reader to [17, 46, 26]and, in particular, to [46, Section 1.1] or [26, Proposition 3.3] for a proof of theequivalence of these two definitions. R.H. Nochetto, E. Ot´arola, A.J. Salgado
Even if we restrict ourselves to definitions that do not require extension, thereis more than one possibility. For instance, the so-called regional fractional Laplacian([39, 13]) is defined by restricting the Riesz integral to Ω, leading to an operatorrelated to a Neumann problem. A different operator is obtained by using the spectraldecomposition of the Dirichlet Laplace operator − ∆, see [14, 18, 20]. This approachis also different to the integral formula (2.10). Indeed, the spectral definition dependson the domain Ω considered, while the integral one at any point is independent of thedomain in which the equation is set. For more details see the discussion in [54].The definition that we shall adopt is as in [14, 18, 20] and is based on the spectraltheory of the Dirichlet Laplacian ([33, 35]) as we summarize below.We define − ∆ : L (Ω) → L (Ω) with domain Dom( − ∆) = { v ∈ H (Ω) : ∆ v ∈ L (Ω) } . This operator is unbounded, closed and, since Ω is bounded and with Lip-schitz boundary, regularity theory implies that its inverse is compact. This impliesthat the spectrum of the operator − ∆ is discrete, positive and accumulates at infinity.Moreover, there exist { λ k , ϕ k } k ∈ N ⊂ R + × H (Ω) such that { ϕ k } k ∈ N is an orthonormalbasis of L (Ω) and, for k ∈ N , (cid:40) − ∆ ϕ k = λ k ϕ k , in Ω ,ϕ k = 0 , on ∂ Ω . (2.11)Moreover, { ϕ k } k ∈ N is an orthogonal basis of H (Ω) and (cid:107)∇ x (cid:48) ϕ k (cid:107) L (Ω) = √ λ k .With this spectral decomposition at hand, fractional powers of the DirichletLaplacian ( − ∆) s can be defined for u ∈ C ∞ (Ω) by( − ∆) s u = ∞ (cid:88) k =1 u k λ sk ϕ k , (2.12)where the coefficients u k are defined by u k = (cid:82) Ω uϕ k . Therefore, if f = (cid:80) ∞ k =1 f k ϕ k ,and ( − ∆) s u = f , then u k = λ − sk f k for all k ≥ − ∆) s can be extended to the Hilbert space H s (Ω) = (cid:40) w = ∞ (cid:88) k =1 w k ϕ k ∈ L (Ω) : (cid:107) w (cid:107) H s (Ω) = ∞ (cid:88) k =1 λ sk | w k | < ∞ (cid:41) . The theory of Hilbert scales presented in [47, Chapter 1] shows that (cid:2) H (Ω) , L (Ω) (cid:3) θ = Dom( − ∆) s , where θ = 1 − s . This implies the following characterization of the space H s (Ω), H s (Ω) = H s (Ω) , s ∈ (0 , ) ,H / (Ω) , s = ,H s (Ω) , s ∈ ( , . (2.13) To exploit the Caffarelli-Silvestre extension[19], or its variants [14, 18, 20], we need to deal with a degenerate/singular ellipticequation on R n +1+ . To this end, we consider weighted Sobolev spaces (see, for instance,[34, 40, 45]), with the specific weight | y | α with α ∈ ( − , PDE approach to fractional diffusion
D ⊂ R n +1 be an open set and α ∈ ( − , L ( D , | y | α ) as the spaceof all measurable functions defined on D such that (cid:107) w (cid:107) L ( D , | y | α ) = (cid:90) D | y | α w < ∞ . Similarly we define the weighted Sobolev space H ( D , | y | α ) = (cid:8) w ∈ L ( D , | y | α ) : |∇ w | ∈ L ( D , | y | α ) (cid:9) , where ∇ w is the distributional gradient of w . We equip H ( D , | y | α ) with the norm (cid:107) w (cid:107) H ( D , | y | α ) = (cid:16) (cid:107) w (cid:107) L ( D , | y | α ) + (cid:107)∇ w (cid:107) L ( D , | y | α ) (cid:17) . (2.14)Notice that taking α = 0 in the definition above, we obtain the classical H ( D ).Properties of this weighted Sobolev space can be found in classical referenceslike [40, 45]. It is remarkable that most of the properties of classical Sobolev spaceshave a weighted counterpart and it is more so that this is not because of the specificform of the weight but rather due to the fact that the weight | y | α belongs to theso-called Muckenhoupt class A ( R n +1 ); see [34, 37, 51]. We recall the definition ofMuckenhoupt classes. Definition 2.1 (
Muckenhoupt class A p ). Let ω be a positive and measurablefunction such that ω ∈ L loc ( R N ) with N ≥ . We say ω ∈ A p ( R N ) , < p < ∞ , ifthere exists a positive constant C p,ω such that sup B (cid:18) | B | (cid:90) B ω (cid:19) (cid:18) | B | (cid:90) B ω / (1 − p ) (cid:19) p − = C p,ω < ∞ , (2.15) where the supremum is taken over all balls B in R N and | B | denotes the Lebesguemeasure of B . Since α ∈ ( − ,
1) it is immediate that | y | α ∈ A ( R n +1 ), which implies the follow-ing important result (see [37, Theorem 1]). Proposition 2.2 (
Properties of weighted Sobolev spaces ). Let
D ⊂ R n +1 bean open set and α ∈ ( − , . Then H ( D , | y | α ) , equipped with the norm (2.14) , is aHilbert space. Moreover, the set C ∞ ( D ) ∩ H ( D , | y | α ) is dense in H ( D , | y | α ) . Remark 2.3 (
Weighted L vs L ). If D is a bounded domain and α ∈ ( − , L ( D , | y | α ) ⊂ L ( D ). Indeed, since | y | − α ∈ L loc ( R n +1 ), (cid:90) D | w | = (cid:90) D | w || y | α/ | y | − α/ ≤ (cid:18)(cid:90) D | w | | y | α (cid:19) (cid:18)(cid:90) D | y | − α (cid:19) (cid:46) (cid:107) w (cid:107) L ( D , | y | α ) . The following result is given in [45, Theorem 6.3]. For completeness we presenthere a version of the proof on the truncated cylinder C Y , which will be important forthe numerical approximation of problem (1.2). Proposition 2.4 (
Embeddings in weighted Sobolev spaces ). Let Ω be a boundeddomain in R n and Y > . Then H ( C Y ) (cid:44) → H ( C Y , y α ) , for α ∈ (0 , , (2.16) and H ( C Y , y α ) (cid:44) → H ( C Y ) , for α ∈ ( − , . (2.17) R.H. Nochetto, E. Ot´arola, A.J. Salgado
Proof . Let us prove (2.16), the proof of (2.17) being similar. Since α > y α ≤ Y α , whence y α w ≤ Y α w and y α |∇ w | ≤ Y α |∇ w | a.e. on C Y for all w ∈ H ( C Y ). This implies (cid:107) w (cid:107) H ( C Y ,y α ) ≤ √ Y α/ (cid:107) w (cid:107) H ( C Y ) , which is (2.16).Define ◦ H L ( C , y α ) = (cid:8) w ∈ H ( y α ; C ) : w = 0 on ∂ L C (cid:9) . (2.18)This space can be equivalently defined as the set of measurable functions w : C → R such that w ∈ H (Ω × ( s, t )) for all 0 < s < t < ∞ , w = 0 on ∂ L C and for which thefollowing seminorm is finite (cid:107) w (cid:107) ◦ H L ( C ,y α ) = (cid:90) C y α |∇ w | ; (2.19)see [20]. As a consequence of the usual Poincar´e inequality, for any k ∈ Z and anyfunction w ∈ H (Ω × (2 k , k +1 )) with w = 0 on ∂ Ω × (2 k , k +1 ), we have (cid:90) Ω × (2 k , k +1 ) y α w ≤ C Ω (cid:90) Ω × (2 k , k +1 ) y α |∇ w | , (2.20)where C Ω denotes a positive constant that depends only on Ω. Summing up over k ∈ Z , we obtain the following weighted Poincar´e inequality : (cid:90) C y α w (cid:46) (cid:90) C y α |∇ w | . (2.21)Hence, the seminorm (2.19) is a norm on ◦ H L ( C , y α ), equivalent to (2.14).For a function w ∈ H ( C , y α ), we shall denote by tr Ω w its trace onto Ω × { } . Itis well known that tr Ω H ( C ) = H / (Ω); see [3, 58]. In the subsequent analysis weneed a characterization of the trace of functions in H ( C , y α ). For a smooth domainthis was given in [18, Proposition 1.8] for s = 1 / s ∈ (0 , \ { } . However, since the eigenvalue decomposition (2.12) of the DirichletLaplace operator holds true on a Lipschitz domain, we are able to extend this tracecharacterization to such domains. In summary, we have the following result. Proposition 2.5 (
Characterization of tr Ω ◦ H L ( C , y α ) ). Let Ω ⊂ R n be a boundedLipschitz domain. The trace operator tr Ω satisfies tr Ω ◦ H L ( C , y α ) = H s (Ω) and (cid:107) tr Ω v (cid:107) H s (Ω) (cid:46) (cid:107) v (cid:107) ◦ H L ( C ,y α ) ∀ v ∈ ◦ H L ( C , y α ) , where the space H s (Ω) is defined in (2.13) . It has been shown in [19]that any power of the fractional Laplacian in R n can be determined as an operator thatmaps a Dirichlet boundary condition to a Neumann-type condition via an extensionproblem posed on R n +1+ . For a bounded domain, an analogous result has been obtainedin [18] for s = , and in [14, 20, 57] for any s ∈ (0 , u defined on Ω. Wedefine the α -harmonic extension of u to the cylinder C , as the function u that solvesthe boundary value problem div( y α ∇ u ) = 0 , in C , u = 0 , on ∂ L C , u = u, on Ω × { } . (2.22) PDE approach to fractional diffusion u ∈ ◦ H L ( C , y α ) whenever u ∈ H s (Ω). We define the Dirichlet-to-Neumann operator Γ α, Ω : H s (Ω) → H s (Ω) (cid:48) u ∈ H s (Ω) (cid:55)−→ Γ α, Ω ( u ) = ∂ u ∂ν α ∈ H s (Ω) (cid:48) , where u solves (2.22) and ∂ u ∂ν α is given in (1.3). The space H s (Ω) (cid:48) can be character-ized as the space of distributions h = (cid:80) k h k ϕ k such that (cid:80) k | h k | λ − sk < ∞ . Thefundamental result of [19], see also [20, Lemma 2.2], is then that d s ( − ∆) s u = Γ α, Ω ( u ) , where d s is given by d s = 2 − s Γ(1 − s )Γ( s ) . (2.23)It seems remarkable that this constant does not depend on the dimension. This wasproved originally in [19] and its precise value appears in several references, for instance[14, 17].The relation between the fractional Laplacian and the extension problem is nowclear. Given f ∈ H s (Ω) (cid:48) , a function u ∈ H s (Ω) solves (1.1) if and only if its α -harmonicextension u ∈ ◦ H L ( C , y α ) solves (1.2).If u = (cid:80) k u k ϕ k , then, as shown in the proofs of [20, Proposition 2.1] and [14,Lemma 2.2], u can be expressed as u ( x ) = ∞ (cid:88) k =1 u k ϕ k ( x (cid:48) ) ψ k ( y ) , (2.24)where the functions ψ k solve ψ (cid:48)(cid:48) k + αy ψ (cid:48) k − λ k ψ k = 0 , in (0 , ∞ ) ,ψ k (0) = 1 , lim y →∞ ψ k ( y ) = 0 . (2.25)If s = , then clearly ψ k ( y ) = e −√ λ k y (see [18, Lemma 2.10]). For s ∈ (0 , \ { } instead (cf. [20, Proposition 2.1]) ψ k ( y ) = c s (cid:16)(cid:112) λ k y (cid:17) s K s ( (cid:112) λ k y ) , where K s denotes the modified Bessel function of the second kind (see [1, Chap-ter 9.6]). Using the condition ψ k (0) = 1, and formulas for small arguments of thefunction K s (see for instance § c s = 2 − s Γ( s ) . The function u ∈ ◦ H L ( C , y α ) is the unique solution of (cid:90) C y α ∇ u · ∇ φ = d s (cid:104) f, tr Ω φ (cid:105) H s (Ω) × H s (Ω) (cid:48) , ∀ φ ∈ ◦ H L ( C , y α ) , (2.26)0 R.H. Nochetto, E. Ot´arola, A.J. Salgado where (cid:104)· , ·(cid:105) H s (Ω) × H s (Ω) (cid:48) denotes the duality pairing between H s (Ω) and H s (Ω) (cid:48) which,in light of Proposition 2.5 is well defined for all f ∈ H s (Ω) (cid:48) and φ ∈ ◦ H L ( C , y α ). Thisimplies the following equalities (see [20, Proposition 2.1] for s ∈ (0 , \ { } and [18,Proposition 2.1] for s = ): (cid:107) u (cid:107) ◦ H L ( C ,y α ) = d s (cid:107) u (cid:107) H s (Ω) = d s (cid:107) f (cid:107) H s (Ω) (cid:48) . (2.27)Notice that for s = , or equivalently α = 0, problem (2.26) reduces to the weakformulation of the Laplace operator with mixed boundary conditions, which is posedon the classical Sobolev space ◦ H L ( C ). Therefore, the value s = becomes a specialcase for problem (2.26). In addition, d / = 1, and (cid:107) u (cid:107) ◦ H L ( C ) = (cid:107) u (cid:107) H / (Ω) .At this point it is important to give a precise meaning to the Dirichlet boundarycondition in (1.1). For s = , the boundary condition is interpreted in the sense ofthe Lions–Magenes space. If < s ≤
1, there is a trace operator from H s (Ω) into L ( ∂ Ω) and the boundary condition can be interpreted in this sense. For 0 < s < / f ∈ H s (Ω) (cid:48) theboundary condition does not have a clear meaning. For instance, for every s ∈ (0 , ), f = ( − ∆) s ∈ H s (Ω) (cid:48) and the solution to (1.1) for this right hand side is u = 1. If f ∈ H ζ (Ω) with ζ > − s > − s , using that ( − ∆) s is a pseudo-differential operatorof order 2 s a shift-type result is valid, i.e., u ∈ H (cid:37) (Ω) with (cid:37) = ζ + 2 s > /
2. In thiscase, the trace of u on ∂ Ω is well defined and the boundary condition is meaningful.Finally, we comment that it has been proved in [20, Lemma 2.10], that if f ∈ L ∞ (Ω)then the solution of (1.1) belongs to C , κ (Ω) with κ ∈ (0 , min { s, } ). It is important to understand the behavior of thesolution u of problem (1.2), given by (2.24). Consequently, it becomes necessary torecall some of the main properties of the modified Bessel function of the second kind K ν ( z ), ν ∈ R ; see [1, Chapter 9.6] for (i)-(iv) and [50, Theorem 5] for (v):(i) For ν > − K ν ( z ) is real and positive.(ii) For ν ∈ R , K ν ( z ) = K − ν ( z ).(iii) For ν >
0, lim z ↓ K ν ( z ) Γ( ν ) (cid:0) z (cid:1) − ν = 1 . (2.28)(iv) For k ∈ N , (cid:18) z dd z (cid:19) k ( z ν K ν ( z )) = ( − k z ν − k K ν − k ( z ) . In particular, for k = 1 and k = 2, respectively, we havedd z ( z ν K ν ( z )) = − z ν K ν − ( z ) = − z ν K − ν ( z ) , (2.29)and d d z ( z ν K ν ( z )) = z ν K − ν ( z ) − z ν − K − ν ( z ) . (2.30)(v) For z > z min { ν, / } e z K ν ( z ) is a decreasing function. PDE approach to fractional diffusion ψ k , defined in (2.25). First, for s ∈ (0 , y ↓ + y α ψ (cid:48) k ( y ) d s λ sk = − , (2.31)Property (v) provides the following asymptotic estimate for s ∈ (0 ,
1) and y ≥ | y α ψ k ( y ) ψ (cid:48) k ( y ) | ≤ C ( s ) λ sk (cid:16)(cid:112) λ k y (cid:17) (cid:12)(cid:12)(cid:12) s − (cid:12)(cid:12)(cid:12) e − √ λ k y . (2.32)Multiplying the differential equation of problem (2.25) by y α ψ k ( y ) and integrating byparts yields (cid:90) ba y α (cid:0) λ k ψ k ( y ) + ψ (cid:48) k ( y ) (cid:1) d y = y α ψ k ( y ) ψ (cid:48) k ( y ) | ba , (2.33)where a and b are real and positive constants.Let us conclude this section with some remarks on the asymptotic behavior of thefunction u that solves (2.26). Using (2.24) we obtain u ( x ) | y =0 = ∞ (cid:88) k =1 u k ϕ k ( x (cid:48) ) ψ k (0) = ∞ (cid:88) k =1 u k ϕ k ( x (cid:48) ) = u ( x (cid:48) ) . For s ∈ (0 , ∂ u ∂ν α ( x (cid:48) ,
0) = − lim y ↓ y α u y ( x (cid:48) , y ) = d s f ( x (cid:48) ) , on Ω × { } . (2.34)Notice that, if s = , then α = 0, d / = 1 and thus (2.34) reduces to ∂ u ∂ν (cid:12)(cid:12)(cid:12)(cid:12) Ω ×{ } = f ( x (cid:48) ) . For s ∈ (0 , \ { } the asymptotic behavior of the second derivative u yy as y ≈ + isa consequence of (2.30) applied to the function ψ k ( y ). For s = the behavior followsfrom ψ k ( y ) = e −√ λ k y . In conclusion, for y ≈ + , we have u yy ≈ y − α − for s ∈ (0 , \ { } , u yy ≈ s = . (2.35) Since we are interested in the approximationof the solution of problem (2.26), and this is closely related to its regularity, let usnow study the behavior of its derivatives. According to (2.34), u y ≈ y − α for y ≈ + .This clearly shows the necessity of introducing the weight, as this behavior, togetherwith the exponential decay given by (v) of § u y ∈ L ( C , y α ) \ L ( C )for s ∈ (0 , / y . Fromthe asymptotic formula (2.35), we see that, for 0 < δ (cid:28) s ∈ (0 , \ { } , (cid:90) Ω × (0 ,δ ) y α | u yy | d x (cid:48) d y ≈ (cid:90) δ y α y − − α d y = (cid:90) δ y − − α d y, (2.36)2 R.H. Nochetto, E. Ot´arola, A.J. Salgado which, since α ∈ ( − , \ { } , does not converge. However, (cid:90) Ω × (0 ,δ ) y β | u yy | d x d y ≈ (cid:90) δ y β − − α d y, converges for β > α + 1, hinting at the fact that u ∈ H ( C , y β ) \ H ( C , y α ). Thefollowing result makes these considerations rigorous. Theorem 2.6 (
Global regularity of the α -harmonic extension ). Let f ∈ H − s (Ω) ,where H − s (Ω) is defined in (2.13) for s ∈ (0 , . Let u ∈ ◦ H L ( C , y α ) solve (2.26) with f as data. Then, for s ∈ (0 , \ { } , we have (cid:107) ∆ x (cid:48) u (cid:107) L ( C ,y α ) + (cid:107) ∂ y ∇ x (cid:48) u (cid:107) L ( C ,y α ) = d s (cid:107) f (cid:107) H − s (Ω) , (2.37) and (cid:107) u yy (cid:107) L ( C ,y β ) (cid:46) (cid:107) f (cid:107) L (Ω) , with β > α + 1 . For the special case s = , we obtain (cid:107) u (cid:107) H ( C ) (cid:46) (cid:107) f (cid:107) H / (Ω) . Remark 2.7 (
Compatibility of f ). It is possible to interpret the result of The-orem 2.6 as follows. Consider s ∈ ( , α ∈ ( − , u gives us that u y ≈ − d s y − α f as y ≈ + onΩ × { } , which in turn implies that u y → y → + on Ω × { } . This is compatiblewith u = 0 on ∂ L C since this implies u y = 0 on ∂ L C . Consequently, we do not need anycompatibility condition on the data f ∈ H − s (Ω) to avoid a jump on the derivative u y . On the other hand, when α ∈ (0 , f , u y (cid:57) y → + on Ω × { } . To compensate this behavior we need the data f to vanish atthe boundary ∂ Ω at a certain rate. This condition is expressed by the requirement f ∈ H − s (Ω). Proof of Theorem 2.6.
Let us first consider s = . In this case (2.26) reduces to thePoisson problem with mixed boundary conditions. In general, the solution of a mixedboundary value problem is not smooth, even for C ∞ data. The singular behavioroccurs near the points of intersection between the Dirichlet and Neumann boundary.For instance, the solution w = √ r sin( θ/
2) of ∆ w = 0 in R , with w x = 0 for { x < , x = 0 } and w = 0 for { x ≥ , x = 0 } does not belong to H ( R ). To obtainmore regular solutions, a compatibility condition between the data, the operator andthe boundary must be imposed (see, for instance, [52]). Since in our case we havethe representation (2.24), we can explicitly compute the second derivatives and, usingthat { ϕ k } k ∈ N is an orthonormal basis of L (Ω) and { ϕ k / √ λ k } k ∈ N of H (Ω), it is notdifficult to show that f ∈ H / (Ω) implies u ∈ H ( C ), and (cid:107) u (cid:107) H ( C ) (cid:46) (cid:107) f (cid:107) H / (Ω) .In the general case s ∈ (0 , \ { } , i.e., α ∈ ( − , \ { } , using (2.33) as well asthe asymptotic properties (2.31) and (2.32), we obtain (cid:107) ∆ x (cid:48) u (cid:107) L ( C ,y α ) + (cid:107) ∂ y ∇ x (cid:48) u (cid:107) L ( C ,y α ) = ∞ (cid:88) k =1 u k λ k (cid:90) ∞ y α (cid:0) λ k ψ k ( y ) + ψ (cid:48) k ( y ) (cid:1) d y = d s ∞ (cid:88) k =1 u k λ sk = d s ∞ (cid:88) k =1 f k λ − sk = d s (cid:107) f (cid:107) H − s (Ω) , PDE approach to fractional diffusion u yy we, again, use the exact representation (2.24) and properties of Besselfunctions to conclude that any derivative with respect to the extended variable y issmooth away from the Neumann boundary Ω × { } . By virtue of (2.25) we deducethat the following partial differential equation holds in the strong sensediv( y α ∇ u ) = 0 ⇐⇒ u yy = − ∆ x (cid:48) u − αy u y . (2.38)Consider sequences { a k = 1 / √ λ k } k ≥ , { b k } k ≥ and { δ k } k ≥ with 0 < δ k ≤ a k ≤ b k .Using (2.24) we have, for k ≥ (cid:107) u yy (cid:107) L ( C ,y β ) = ∞ (cid:88) k =1 u k (cid:32) lim δ k ↓ (cid:90) a k δ k y β | ψ (cid:48)(cid:48) k ( y ) | d y + lim b k ↑∞ (cid:90) b k a k y β | ψ (cid:48)(cid:48) k ( y ) | d y (cid:33) (2.39)Let us now estimate the first integral on the right hand side of (2.39). Formulas (2.30)and (2.28) yieldlim δ k ↓ (cid:90) a k δ k y β | ψ (cid:48)(cid:48) k ( y ) | d y = c s λ − β/ − / k lim δ k ↓ (cid:90) √ λ k δ k z β (cid:12)(cid:12)(cid:12)(cid:12) d d z ( z s K s ( z )) (cid:12)(cid:12)(cid:12)(cid:12) d z (cid:46) c s λ − β/ − / k lim δ k ↓ (cid:90) √ λ k δ k z β − − α d z ≈ λ − β/ − / k (2.40)where the integral converges because β > α + 1. Let us now look at the secondintegral. Using property (v) of the modified Bessel functions, we havelim b k ↑∞ (cid:90) b k a k y β | ψ (cid:48)(cid:48) k ( y ) | d y = c s λ − β/ − / k lim b k ↑∞ (cid:90) √ λ k b k z β (cid:12)(cid:12)(cid:12)(cid:12) d d z ( z s K s ( z )) (cid:12)(cid:12)(cid:12)(cid:12) d z (cid:46) c s λ − β/ − / k . (2.41)Replacing (2.40) and (2.41) into (2.39), and using that u k = λ − sk f k , we deduce (cid:107) u yy (cid:107) L ( C ,y β ) (cid:46) ∞ (cid:88) k =1 λ − β/ − / − sk f k ≤ (cid:107) f (cid:107) L (Ω) , because 2 − s − β − = (1 + 2 α − β ) <
0. This concludes the proof.For the design of graded meshes later in § Theorem 2.8 (
Local regularity of the α -harmonic extension ). Let C ( a, b ) :=Ω × ( a, b ) for ≤ a < b ≤ . The solution u ∈ ◦ H L ( C , y α ) of (2.26) satisfies for all a, b (cid:107) ∆ x (cid:48) u (cid:107) L ( C ( a,b ) ,y α ) + (cid:107) ∂ y ∇ x (cid:48) u (cid:107) L ( C ( a,b ) ,y α ) (cid:46) ( b − a ) (cid:107) f (cid:107) H − s (Ω) , (2.42) and, with δ := β − α − > , (cid:107) u yy (cid:107) L ( C ( a,b ) ,y β ) (cid:46) (cid:0) b δ − a δ (cid:1) (cid:107) f (cid:107) L (Ω) . (2.43) Proof . To derive (2.42) we proceed as in Theorem 2.6. Since 0 ≤ a < b ≤ § | y α ψ k ( y ) ψ (cid:48) k ( y ) | (cid:46) λ sk . R.H. Nochetto, E. Ot´arola, A.J. Salgado
This, together with (2.33) and the property u k = λ − sk f k , allows us to conclude (cid:107) ∆ x (cid:48) u (cid:107) L ( C ( a,b ) ,y α ) + (cid:107) ∂ y ∇ x (cid:48) u (cid:107) L ( C ( a,b ) ,y α ) = ∞ (cid:88) k =1 u k λ k (cid:90) ba y α (cid:0) λ k ψ k ( y ) + ψ (cid:48) k ( y ) (cid:1) d y (cid:46) ( b − a ) ∞ (cid:88) k =1 u k λ sk = ( b − a ) (cid:107) f (cid:107) H − s (Ω) . To prove (2.43) we observe that the same argument used in (2.40) gives (cid:90) ba y β | ψ (cid:48)(cid:48) k ( y ) | d y (cid:46) λ − β/ − / k (cid:0) b δ − a δ (cid:1) , whence (cid:107) u yy (cid:107) L ( C ( a,b ) ,y α ) (cid:46) (cid:0) b δ − a δ (cid:1) ∞ (cid:88) k =1 f k λ − β/ − / − sk (cid:46) (cid:0) b δ − a δ (cid:1) (cid:107) f (cid:107) L (Ω) , because 2 − s − β − <
3. Truncation.
The solution u of problem (2.26) is defined on the infinite do-main C and, consequently, it cannot be directly approximated with finite element-liketechniques. In this section we will show that u decays sufficiently fast – in fact expo-nentially – in the extended direction. This suggests truncating the cylinder C to C Y ,for a suitably defined Y . The exponential decay is the content of the next result. Proposition 3.1 (
Exponential decay ). For every Y > , the solution u of (2.26) satisfies (cid:107)∇ u (cid:107) L (Ω × ( Y , ∞ ) ,y α ) (cid:46) e −√ λ Y / (cid:107) f (cid:107) H s (Ω) (cid:48) . (3.1) Proof . Recall that if u ∈ H s (Ω) has the decomposition u = (cid:80) k u k ϕ k ( x (cid:48) ), thesolution u ∈ ◦ H L ( C , y α ) to (2.26) has the representation u = (cid:80) k u k ϕ ( x (cid:48) ) ψ k ( y ), wherethe functions ψ k solve (2.25).Consider s = . In this case ψ k ( y ) = e −√ λ k y . Using the fact that { ϕ k } ∞ k =1 areeigenfunctions of Dirichlet Laplacian on Ω, orthonormal in L (Ω) and orthogonal in H (Ω) with (cid:107)∇ x (cid:48) ϕ k (cid:107) L (Ω) = √ λ k , we get (cid:90) ∞ Y (cid:90) Ω |∇ u | = (cid:90) ∞ Y (cid:90) Ω (cid:0) |∇ x (cid:48) u | + | ∂ y u | (cid:1) = ∞ (cid:88) k =1 λ k | u k | e − √ λ k Y ≤ e − √ λ Y (cid:107) u (cid:107) H / (Ω) . Since (cid:107) u (cid:107) H / (Ω) = (cid:107) f (cid:107) H / (Ω) (cid:48) , this implies (3.1).Consider now s ∈ (0 , \ { } and ψ k ( y ) = c s (cid:0) √ λ k y (cid:1) s K s ( √ λ k y ). To be able toargue as before, we need the estimates on K s and its derivative for sufficiently largearguments discussed in § (cid:90) ∞ Y (cid:90) Ω y α |∇ u | = (cid:90) ∞ Y y α (cid:90) Ω (cid:0) |∇ x (cid:48) u | + | ∂ y u | (cid:1) = ∞ (cid:88) k =1 | u k | (cid:90) ∞ Y y α (cid:0) λ k ψ k ( y ) + ψ (cid:48) k ( y ) (cid:1) d y = ∞ (cid:88) k =1 | u k | y α ψ k ( y ) ψ (cid:48) k ( y ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∞ Y (cid:46) e −√ λ Y (cid:107) u (cid:107) H s (Ω) . PDE approach to fractional diffusion (cid:107) u (cid:107) H s (Ω) = (cid:107) f (cid:107) H s (Ω) (cid:48) we get (3.1).Expression (3.1) motivates the approximation of u by a function v that solves div( y α ∇ v ) = 0 , in C Y ,v = 0 , on ∂ L C Y ∪ Ω × { Y } ,∂v∂ν α = d s f, on Ω × { } , (3.2)with Y sufficiently large. Problem (3.2) is understood in the weak sense, i.e., we definethe space ◦ H L ( C Y , y α ) = (cid:8) v ∈ H ( C , y α ) : v = 0 on ∂ L C Y ∪ Ω × { Y } (cid:9) , and seek for v ∈ ◦ H L ( C Y , y α ) such that (cid:90) C Y y α ∇ v · ∇ φ = d s (cid:104) f, tr Ω φ (cid:105) , ∀ φ ∈ ◦ H L ( C Y , y α ) . (3.3)Existence and uniqueness of v follows from the Lax-Milgram lemma. Remark 3.2 (
Zero extension ). For every Y > ◦ H L ( C Y , y α ) (cid:44) → ◦ H L ( C , y α ) . (3.4)To see this, it suffices to consider the extension by zero for y > Y .The next result shows the approximation properties of v , solution of (3.3) in C Y . Lemma 3.3 (
Exponential convergence in Y ). For any positive Y > , we have (cid:107)∇ ( u − v ) (cid:107) L ( C Y ,y α ) (cid:46) e −√ λ Y / (cid:107) f (cid:107) H s (Ω) (cid:48) . (3.5) Proof . Given φ ∈ ◦ H L ( C Y , y α ) denote by φ e its extension by zero to C . ByRemark 3.2, φ e ∈ ◦ H L ( C , y α ). Take φ e and φ as test functions in (2.26) and (3.3),respectively. Subtract the resulting expressions to obtain (cid:90) C Y y α ( ∇ u − ∇ v ) · ∇ φ = 0 ∀ φ ∈ ◦ H L ( C Y , y α ) , which implies that v is the best approximation of u in ◦ H L ( C Y , y α ), i.e., (cid:107)∇ ( u − v ) (cid:107) L ( C Y ,y α ) = inf φ ∈ ◦ H L ( C Y ,y α ) (cid:107)∇ ( u − φ ) (cid:107) L ( C Y ,y α ) . (3.6)Let us construct explicitly a function φ ∈ ◦ H L ( C Y , y α ) to use in (3.6). Define ρ ( y ) = , ≤ y ≤ Y / , Y ( Y − y ) , Y / < y < Y , , y ≥ Y . (3.7)Notice that ρ ∈ W ∞ (0 , ∞ ), | ρ ( y ) | ≤ | ρ (cid:48) ( y ) | ≤ / Y for all y >
0. Set φ ( x (cid:48) , y ) = u ( x (cid:48) , y ) ρ ( y ) for x (cid:48) ∈ Ω and y >
0. A straightforward computation shows |∇ ((1 − ρ ) u ) | ≤ (cid:0) | ρ (cid:48) | | u | + (1 − ρ ) |∇ u | (cid:1) ≤ (cid:18) Y u + |∇ u | (cid:19) , R.H. Nochetto, E. Ot´arola, A.J. Salgado so that (cid:107)∇ ( u − φ ) (cid:107) L ( C Y ,y α ) ≤ (cid:32) Y (cid:90) YY / (cid:90) Ω y α | u | + (cid:90) YY / (cid:90) Ω y α |∇ u | (cid:33) . (3.8)To estimate the first term on the right hand side of (3.8) we use the Poincar´einequality (2.20) over a dyadic partition that covers the interval [ Y / , Y ] (see thederivation of (2.21) in § (cid:90) YY / y α (cid:90) Ω | u | (cid:46) (cid:90) YY / y α (cid:90) Ω |∇ u | . (3.9)To bound the second integral in (3.8) we use (2.33) as in the proof of Proposition 3.1: (cid:90) YY / y α (cid:90) Ω |∇ u | = ∞ (cid:88) k =1 | u k | y α ψ k ( y ) ψ (cid:48) k ( y ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) YY / (cid:46) e −√ λ Y / (cid:107) f (cid:107) H s (Ω) (cid:48) . Inserting these estimates into (3.6) implies (3.5).The following result is a direct consequence of Lemma 3.3.
Remark 3.4 (
Stability ). Let Y ≥
1, then (cid:107)∇ v (cid:107) L ( C Y ,y α ) (cid:46) (cid:107) f (cid:107) H s (Ω) (cid:48) . (3.10)Indeed, by the triangle inequality (cid:107)∇ v (cid:107) L ( C Y ,y α ) ≤ (cid:107)∇ ( v − u ) (cid:107) L ( C Y ,y α ) + (cid:107)∇ u (cid:107) L ( C Y ,y α ) (cid:46) (cid:16) e −√ λ Y / + 1 (cid:17) (cid:107) f (cid:107) H s (Ω) (cid:48) . The previous two results allow us to show a full approximation estimate.
Theorem 3.5 (
Global exponential estimate ). Let Y > , then (cid:107)∇ ( u − v ) (cid:107) L ( C ,y α ) (cid:46) e −√ λ Y / (cid:107) f (cid:107) H s (Ω) (cid:48) . (3.11) In particular, for every (cid:15) > , let Y = 2 √ λ (cid:18) log C + 2 log 1 (cid:15) (cid:19) , where C depends only on s and Ω . Then, for Y ≥ max { Y , } , we have (cid:107)∇ ( u − v ) (cid:107) L ( C ,y α ) ≤ (cid:15) (cid:107) f (cid:107) H s (Ω) (cid:48) . (3.12) Proof . Extending v by zero outside of C Y we obtain (cid:107)∇ ( u − v ) (cid:107) L ( C ,y α ) = (cid:107)∇ ( u − v ) (cid:107) L ( C Y ,y α ) + (cid:107)∇ u (cid:107) L (Ω × ( Y , ∞ ) ,y α ) . Hence Lemma 3.3 and Proposition 3.1 imply (cid:107)∇ ( u − v ) (cid:107) L ( C ,y α ) ≤ Ce −√ λ Y / (cid:107) f (cid:107) H s (Ω) (cid:48) ≤ (cid:15) (cid:107) f (cid:107) H s (Ω) (cid:48) , (3.13)for all Y ≥ max { Y , } . PDE approach to fractional diffusion
4. Finite element discretization and interpolation estimates.
In this sec-tion we prove error estimates for a piecewise Q interpolation operator on anisotropicelements in the extended variable y . We consider elements of the form T = K × I ,where K ⊂ R n is an element isoparametrically equivalent to the unit cube [0 , n ,via a Q mapping and, I ⊂ R is an interval. The anisotropic character of the mesh T Y = { T } will be given by the family of intervals I .The error estimates are derived in the weighted Sobolev spaces L ( C Y , y α ) and H ( C Y , y α ), and they are valid under the condition that neighboring elements havecomparable size in the extended n + 1–dimension (see [28]). This is a mild assump-tion that includes general meshes which do not satisfy the so-called shape-regularityassumption, i.e., mesh refinements for which the quotient between outer and innerdiameter of the elements does not remain bounded (see [15, Chapter 4]).Anisotropic or narrow elements are elements with disparate sizes in each direc-tion. They arise naturally when approximating solutions of problems with a strongdirectional-dependent behavior since, using anisotropy, the local mesh size can beadapted to capture such features. Examples of this include boundary layers, shocksand edge singularities (see [28, 29]). In our problem, anisotropic elements are essen-tial in order to capture the singular/degenerate behavior of the solution u to problem(2.26) at y ≈ + given in (2.34). These elements will provide optimal error estimates,which cannot be obtained using shape-regular elements.Error estimates for weighted Sobolev spaces have been obtained in several works;see, for instance, [4, 9, 28]. The type of weight considered in [4, 9] is related to thedistance to a point or an edge, and the type of quasi-interpolators are modificationsof the well known Cl´ement [24] and Scott-Zhang [53] operators. These works aredeveloped in 3D and 2D respectively, and the analysis developed in [4] allows foranisotropy. Our approach follows the work of Dur´an and Lombardi [28], and is basedon a piecewise Q averaged interpolator on anisotropic elements. It allows us to obtainanisotropic interpolation estimates in the extended variable y and in weighted Sobolevspaces, using only that | y | α ∈ A ( R n +1 ), the Muckenhoupt class A of Definition 2.1. Let us now describe the discretizationof problem (3.2). To avoid technical difficulties we assume that the boundary of Ωis polygonal. The difficulties inherent to curved boundaries could be handled, forinstance, with the methods of [11] (see also [43, 44]). Let T Ω = { K } be a mesh of Ωmade of isoparametric quadrilaterals K in the sense of Ciarlet [22] and Ciarlet andRaviart [23]. In other words, given ˆ K = [0 , n and a family of mappings {F K ∈ Q ( ˆ K ) n } we have K = F K ( ˆ K ) (4.1)and ¯Ω = (cid:91) K ∈ T Ω K, | Ω | = (cid:88) K ∈ T Ω | K | . The collection of triangulations is denoted by T Ω .The mesh T Ω is assumed to be conforming or compatible, i.e., the intersection ofany two isoparametric elements K and K (cid:48) in T Ω is either empty or a common lowerdimensional isoparametric element.In addition, we assume that T Ω is shape regular (cf. [22, Chapter 4.3]). Thismeans that F K can be decomposed as F K = A K + B K , where A K is affine and B K R.H. Nochetto, E. Ot´arola, A.J. Salgado is a perturbation map and, if we define ˜ K = A K ( ˆ K ), h K = diam( ˜ K ), ρ K as thediameter of the largest sphere inscribed in ˜ K and the shape coefficient of K as theratio σ K = h K /ρ K , then the following two conditions are satisfied:(a) There exists a constant σ Ω > T Ω ∈ T Ω , max { σ K : K ∈ T Ω } ≤ σ Ω . (b) For all K ∈ T Ω the mapping B K is Fr´echet differentiable and (cid:107) D B K (cid:107) L ∞ ( ˆ K ) = O ( h K ) , for all K ∈ T Ω and all T Ω ∈ T Ω .As a consequence of these conditions, if h K is small enough, the mapping F K isone-to-one, its Jacobian J F K does not vanish, and J F K (cid:46) h nK , (cid:107) D F K (cid:107) L ∞ ( ˆ K ) (cid:46) h K . (4.2)The set T Ω is called quasi-uniform if for all T Ω ∈ T Ω ,max { ρ K : K ∈ T Ω } (cid:46) min { h K : K ∈ T Ω } . In this case, we define h T Ω = max K ∈ T h K .We define T Y as a triangulation of C Y into cells of the form T = K × I , where K ∈ T Ω , and I denotes an interval in the extended dimension. Notice that eachdiscretization of the truncated cylinder C Y depends on the truncation parameter Y .The set of all such triangulations is denoted by T . In order to obtain a global regu-larity assumption for T we assume the aforementioned conditions on T Ω , besides thefollowing weak regularity condition:(c) There is a constant σ such that, for all T Y ∈ T , if T = K × I , T = K × I ∈ T Y have nonempty intersection, then h I h I ≤ σ, where h I = | I | .Notice that the assumptions imposed on T are weaker than the standard shape-regularity assumptions, since they allow for anisotropy in the extended variable (cf.[28]). It is also important to notice that, given the Cartesian product structure of thecells T ∈ T Y , they are isoparametrically equivalent to ˆ T = [0 , n +1 . We will denotethe corresponding mappings by F T . Then, F T : ˆ x = (ˆ x (cid:48) , ˆ y ) ∈ ˆ T (cid:55)−→ x = ( x (cid:48) , y ) = ( F K (ˆ x (cid:48) ) , F I (ˆ y )) ∈ T = K × I, where F K is the bilinear mapping defined in (4.1) for K and, if I = ( c, d ), F I ( y ) =( y − c ) / ( d − c ). From (4.2), we immediately conclude that J F T (cid:46) h nK h I , (cid:107) D F T (cid:107) L ∞ ( ˆ T ) (cid:46) h T , (4.3)for all elements T ∈ T Y where h T = max { h K , h I } .Given T Y ∈ T , we define the finite element space V ( T Y ) by V ( T Y ) = (cid:8) W ∈ C ( C Y ) : W | T ∈ Q ( T ) ∀ T ∈ T Y , W | Γ D = 0 (cid:9) . PDE approach to fractional diffusion D = ∂ L C Y ∪ Ω × { Y } is called the Dirichlet boundary. The Galerkin approxi-mation of (3.3) is given by the unique function V T Y ∈ V ( T Y ) such that (cid:90) C Y y α ∇ V T Y · ∇ W = d s (cid:104) f, tr Ω W (cid:105) , ∀ W ∈ V ( T Y ) . (4.4)Existence and uniqueness of V T Y follows from V ( T Y ) ⊂ ◦ H L ( C Y , y α ) and the Lax-Milgram lemma.We define the space U ( T Ω ) = tr Ω V ( T Y ), which is nothing more than a Q finiteelement space over the mesh T Ω . The finite element approximation of u ∈ H s (Ω),solution of (1.1), is then given by U T Ω = tr Ω V T Y ∈ U ( T Ω ) , (4.5)and we have the following result. Theorem 4.1 (
Energy error estimate ). If V T Y ∈ V ( T Y ) solves (4.4) and U T Ω ∈ U ( T Ω ) is defined in (4.5) , then (cid:107) u − U T Ω (cid:107) H s (Ω) (cid:46) (cid:107) u − V T Y (cid:107) ◦ H L ( C ,y α ) , (4.6) and (cid:107) u − V T Y (cid:107) ◦ H L ( C ,y α ) (cid:46) (cid:15) (cid:107) f (cid:107) H s (Ω) (cid:48) + (cid:107) v − V T Y (cid:107) ◦ H L ( C Y ,y α ) . (4.7) Proof . Estimate (4.6) is just an application of the trace estimate of Proposi-tion 2.5. Inequality (4.7) is obtained by the triangle inequality and (3.12).By Galerkin orthogonality (cid:107) v − V T Y (cid:107) ◦ H L ( C Y ,y α ) = inf W ∈ V ( T Y ) (cid:107) v − W (cid:107) ◦ H L ( C Y ,y α ) . Theorem 4.1 and Galerkin orthogonality imply that the approximation estimate (4.7)depends on the regularity of u . To see this we introduce ρ ( y ) = (cid:40) , ≤ y < Y / ,p, Y / ≤ y ≤ Y , (4.8)where p is the unique cubic polynomial on [ Y / , Y ] defined by the conditions p ( Y /
2) =1, p ( Y ) = 0, p (cid:48) ( Y /
2) = 0 and p (cid:48) ( Y ) = 0. Notice that ρ ∈ W ∞ (0 , Y ), | ρ ( y ) | ≤ | ρ (cid:48) ( y ) | (cid:46) | ρ (cid:48)(cid:48) ( y ) | (cid:46)
1. Set u ( x (cid:48) , y ) = ρ ( y ) u ( x (cid:48) , y ) for x (cid:48) ∈ Ω and y ∈ [0 , Y ], andnotice that u ∈ ◦ H L ( C Y , y α ). With this construction at hand, repeating the argumentsused in the proof of Lemma 3.3, we have that (cid:107) ∆ x (cid:48) u (cid:107) L ( C Y ,y α ) (cid:46) (cid:107) ∆ x (cid:48) u (cid:107) L ( C Y ,y α ) , (cid:107) ∂ y ∇ x (cid:48) u (cid:107) L ( C Y ,y α ) (cid:46) (cid:107) ∂ y ∇ x (cid:48) u (cid:107) L ( C Y ,y α ) + (cid:107) f (cid:107) H s (Ω) (cid:48) , (cid:107) ∂ yy u (cid:107) L ( C Y ,y β ) (cid:46) (cid:107) ∂ yy u (cid:107) L ( C Y ,y β ) + (cid:107) f (cid:107) H s (Ω) (cid:48) . In addition, if we assume that there is an operatorΠ T Y : ◦ H L ( C Y , y α ) → V ( T Y ) , R.H. Nochetto, E. Ot´arola, A.J. Salgado that is stable, i.e., (cid:107) Π T Y w (cid:107) ◦ H L ( C Y ,y α ) (cid:46) (cid:107) w (cid:107) ◦ H L ( C Y ,y α ) , for all w ∈ ◦ H L ( C Y , y α ), then thefollowing estimate holds (cid:107) u − V T Y (cid:107) ◦ H L ( C ,y α ) (cid:46) (cid:15) (cid:107) f (cid:107) H s (Ω) (cid:48) + (cid:107) u − Π T Y u (cid:107) ◦ H L ( C Y ,y α ) . (4.9)To see this, we use (4.7), together with Galerkin orthogonality and the stability ofthe operator Π T Y , to obtain (cid:107) u − V T Y (cid:107) ◦ H L ( C ,y α ) (cid:46) (cid:15) (cid:107) f (cid:107) H s (Ω) (cid:48) + (cid:107) v − Π T Y v (cid:107) ◦ H L ( C Y ,y α ) (cid:46) (cid:15) (cid:107) f (cid:107) H s (Ω) (cid:48) + (cid:107) v − u (cid:107) ◦ H L ( C Y ,y α ) + (cid:107) u − Π T Y u (cid:107) ◦ H L ( C Y ,y α ) . The second term on the right hand side of the previous inequality is estimated as inLemma 3.3. We leave the details to the reader.Estimates for u − Π T Y u on weighted Sobolev spaces are derived in § u which, in light of (4.9), depends on the regularityof u . For this reason, and to lighten the notation, we shall in the sequel write u andobtain interpolation error estimates for it, even though u does not vanish at y = Y . Let us begin byintroducing some notation and terminology. Given T Y , we call N the set of its nodesand N in the set of its interior and Neumann nodes. For each vertex v ∈ N , we write v = ( v (cid:48) , v (cid:48)(cid:48) ), where v (cid:48) corresponds to a node of T Ω , and v (cid:48)(cid:48) corresponds to a node of thediscretization of the n + 1–dimension. We define h v (cid:48) = min { h K : v (cid:48) is a vertex of K } ,and h v (cid:48)(cid:48) = min { h I : v (cid:48)(cid:48) is a vertex of I } .Given v ∈ N , the star or patch around v is defined as ω v = (cid:91) v ∈ T T, and for T ∈ T Y we define its patch as ω T = (cid:91) v ∈ T ω v . Let ψ ∈ C ∞ ( R n +1 ) be such that (cid:82) ψ = 1 and D := supp ψ ⊂ B r × ( − r Y , r Y ),where B r denotes the ball in R n of radius r and centered at zero, and r ≤ /σ Ω and r Y ≤ /σ . For v ∈ N in , we rescale ψ as ψ v ( x ) = 1 h n v (cid:48) h v (cid:48)(cid:48) ψ (cid:18) v (cid:48) − x (cid:48) h v (cid:48) , v (cid:48)(cid:48) − yh v (cid:48)(cid:48) (cid:19) , and note that supp ψ v ⊂ ω v for all v ∈ N in .Given a function w ∈ L ( C Y , y α ) and a node v in N in we define, following Dur´anand Lombardi [28], the regularized Taylor polynomial of first degree of w about v as w v ( z ) = (cid:90) P ( x, z ) ψ v ( x ) d x = (cid:90) ω v P ( x, z ) ψ v ( x ) d x, (4.10)where P denotes the Taylor polynomial of degree 1 in the variable z of the function w about the point x , i.e., P ( x, z ) = w ( x ) + ∇ w ( x ) · ( z − x ) . (4.11) PDE approach to fractional diffusion L ( C Y ) (cf. [15, Proposition 4.1.12]), we conclude that P iswell defined for any function in L ( C Y , y α ).We define the average Q interpolant Π T Y w , as the unique piecewise Q functionsuch that Π T Y w ( v ) = 0 if v lies on the Dirichlet boundary Γ D and Π T Y w ( v ) = w v ( v )if v ∈ N in . If λ v denotes the Lagrange basis function associated with node v , thenΠ T Y w = (cid:88) v ∈ N in w v ( v ) λ v . There are two principal reasons to consider average interpolation. First, we areinterested in the approximation of singular functions and thus Lagrange interpolationcannot be used since point-wise values become meaningless. In fact, this motivated theintroduction of average interpolation (see [24, 53]). In addition, average interpolationhas better approximation properties when narrow elements are used (see [2]).Finally, for v ∈ N in , we define the weighted regularized average of w as Q v w = (cid:90) w ( x ) ψ v ( x ) d x = (cid:90) ω v w ( x ) ψ v ( x ) d x. (4.12) In order to obtain interpolation errorestimates in L ( C Y , y α ) and H ( C Y , y α ), it is instrumental to have a weighted Poincar´e-type inequality. Weighted Poincar´e inequalities are particularly pertinent in the studyof the nonlinear potential theory of degenerate elliptic equations, see [34, 40]. If thedomain is a ball and the weight belongs to A p , with 1 ≤ p < ∞ , this result can befound in [34, Theorem 1.3 and Theorem 1.5]. However, to the best of our knowledge,such a result is not available in the literature for more general domains. For ourspecific weight we present here a constructive proof, i.e., not based on a compactnessargument. This allows us to study the dependence of the constant on the domain. Lemma 4.2 (
Weighted Poincar´e inequality I ). Let ω ⊂ R n +1 be bounded, star-shaped with respect to a ball B , and diam ω ≈ . Let χ ∈ C (¯ ω ) with (cid:82) ω χ = 1 , and ξ α ( y ) := | a | y | + b | α for a, b ∈ R . If w ∈ H ( ω, ξ α ( y )) is such that (cid:82) ω χw = 0 , then (cid:107) w (cid:107) L ( ω,ξ α ) (cid:46) (cid:107)∇ w (cid:107) L ( ω,ξ α ) , (4.13) where the hidden constant depends only on χ , α and the radius r of B , but is inde-pendent of both a and b .Proof . The fact that α ∈ ( − ,
1) implies ξ α ∈ A ( R n +1 ) with a Muckenhouptconstant C ,ξ α in (2.15) uniform in both a and b . Define (cid:101) w = ξ α w − (cid:18)(cid:90) ω ξ α w (cid:19) χ. Clearly (cid:101) w ∈ L ( ω ) and it has vanishing mean value by construction.Since (cid:82) ω χw = 0 we obtain (cid:107) w (cid:107) L ( ω,ξ α ) = (cid:90) ω w (cid:101) w + (cid:18)(cid:90) ω ξ α w (cid:19) (cid:90) ω χw = (cid:90) ω w (cid:101) w. (4.14)Consequently, given that ω is star shaped with respect to ˆ B , and ξ α ∈ A ( R n +1 ),there exists F ∈ H ( ω, ξ α ) n +1 such that − div F = (cid:101) w , and (cid:107) F (cid:107) H ( ω,ξ − α ) n +1 (cid:46) (cid:107) (cid:101) w (cid:107) L ( ω,ξ − α ) , (4.15)2 R.H. Nochetto, E. Ot´arola, A.J. Salgado where the hidden constant in (4.15) depends on r and the constant C ,ξ α from Defi-nition 2.1 [30, Theorem 3.1].Replacing (cid:101) w by − div F in (4.14), integrating by parts and using (4.15), we get (cid:107) w (cid:107) L ( ω,ξ α ) = − (cid:90) ω w div F = (cid:90) ω ∇ w · F (cid:46) (cid:107)∇ w (cid:107) L ( ω,ξ α ) (cid:107) (cid:101) w (cid:107) L ( ω,ξ − α ) . (4.16)To estimate (cid:107) (cid:101) w (cid:107) L ( ω,ξ − α ) we use the Cauchy-Schwarz inequality and the constant C ,ξ α from Definition 2.1 as follows: (cid:107) (cid:101) w (cid:107) L ( ω,ξ − α ) ≤ (cid:18) (cid:90) ω ξ α (cid:90) ω χ ξ − α (cid:19) (cid:107) w (cid:107) L ( ω,ξ α ) (cid:46) (cid:107) w (cid:107) L ( ω,ξ α ) . Inserting the inequality above into (4.16), we obtain (4.13).We need a slightly more general form of the Poincar´e inequality for the applica-tions below. We now relax the geometric assumption on the domain ω and let thevanishing mean property hold just in a subdomain. Corollary 4.3 (
Weighted Poincar´e inequality II ). Let ω = ∪ Ni =1 ω i ⊂ R n +1 be aconnected domain and each ω i be a star-shaped domain with respect to a ball B i . Let χ i ∈ C (¯ ω i ) and ξ α be as in Lemma 4.2. If w ∈ H ( ω, ξ α ) and w i := (cid:82) ω i wχ i , then (cid:107) w − w i (cid:107) L ( ω,ξ α ) (cid:46) (cid:107)∇ w (cid:107) L ( ω,ξ α ) ∀ ≤ i ≤ N, (4.17) where the hidden constant depends on { χ i } Ni =1 , α , the radius r i of B i , and the amountof overlap between the subdomains { ω i } Ni =1 , but is independent of both a and b .Proof . This is a consequence of Lemma 4.2 and [27, Theorem 7.1]. We sketch theproof here for completeness. It suffices to deal with two subdomains, ω , ω , and theoverlapping region B = ω ∩ ω . We observe that (cid:107) w − w (cid:107) L ( ω ,ξ α ) ≤ (cid:107) w − w (cid:107) L ( ω ,ξ α ) + (cid:107) w − w (cid:107) L ( ω ,ξ α ) , together with (cid:107) w − w (cid:107) L ( ω ,ξ α ) = (cid:18) (cid:82) ω ξ α (cid:82) B ξ α (cid:19) / (cid:107) w − w (cid:107) L ( B,ξ α ) and (cid:107) w − w (cid:107) L ( B,ξ α ) (cid:46) (cid:107) w − w (cid:107) L ( ω ,ξ α ) + (cid:107) w − w (cid:107) L ( ω ,ξ α ) , imply (cid:107) w − w (cid:107) L ( ω ,ξ α ) (cid:46) (cid:107)∇ w (cid:107) L ( ω ∪ ω ,ξ α ) . This, combined with (4.13), gives (4.17)for i = 1 with a stability constant depending on the ratio (cid:82) ω ξ α (cid:82) B ξ α . L -based interpolation estimates. Owing to the weightedPoincar´e inequality of Corollary 4.3, we can adapt the proof of [28, Lemma 2.3] toobtain interpolation estimates in the weighted L -norm. These estimates allow adisparate mesh size on the extended direction, relative to the coordinate directions x i , i = 1 , . . . , n, which may in turn be graded. This is the principal difference with[28, Lemma 2.3] where, however, the domain must be a cube. Lemma 4.4 (
Weighted L -based interpolation estimates ). Let v ∈ N in . Then,for all w ∈ H ( ω v , y α ) , we have (cid:107) w − Q v w (cid:107) L ( ω v ,y α ) (cid:46) h v (cid:48) (cid:107)∇ x (cid:48) w (cid:107) L ( ω v ,y α ) + h v (cid:48)(cid:48) (cid:107) ∂ y w (cid:107) L ( ω v ,y α ) , (4.18) and, for all v ∈ H ( ω v , y α ) and j = 1 , . . . , n + 1 , we have (cid:107) ∂ x j ( w − w v ) (cid:107) L ( ω v ,y α ) (cid:46) h v (cid:48) n (cid:88) i =1 (cid:107) ∂ x j x i w (cid:107) L ( ω v ,y α ) + h v (cid:48)(cid:48) (cid:107) ∂ x j y w (cid:107) L ( ω v ,y α ) , (4.19) PDE approach to fractional diffusion where, in both inequalities, the hidden constant depends only on α , σ Ω , σ and ψ .Proof . Define by F v : ( x (cid:48) , y ) → (¯ x (cid:48) , ¯ y ) the scaling map¯ x (cid:48) = v (cid:48) − x (cid:48) h v (cid:48) , ¯ y = v (cid:48)(cid:48) − yh v (cid:48)(cid:48) , along with ω v = F v ( ω v ) and ¯ w (¯ x ) = w ( x ). Define also ¯ Q ¯ w = (cid:82) ¯ wψ, where ψ has beenintroduced in Section 4.2. Since supp ψ ⊂ ω v integration takes place only over ω v ,and (cid:82) ω v ψ = 1. Then, ¯ Q ¯ w satisfies ¯ Q ¯ w = (cid:82) ω v ¯ wψ = (cid:82) ω v wψ v = Q v w, and (cid:90) ω v ( ¯ Q ¯ w − w ) ψ d¯ x = ¯ Q ¯ w − (cid:90) ω v ¯ wψ d¯ x = 0 . (4.20)Simple scaling, using the definition of the mapping F v , yields (cid:90) ω v y α | w − Q v w | d x = h n v (cid:48) h v (cid:48)(cid:48) (cid:90) ω v ξ α | ¯ w − ¯ Q ¯ w | d¯ x, (4.21)where ξ α ( y ) := | v (cid:48)(cid:48) − ¯ yh v (cid:48)(cid:48) | α . By shape regularity, the mesh sizes h v (cid:48) , h v (cid:48)(cid:48) satisfy1 / σ ≤ h ¯ v (cid:48)(cid:48) ≤ σ and 1 / σ Ω ≤ h ¯ v (cid:48) ≤ σ Ω , respectively, and diam ω v ≈
1. In view of(4.20), we can apply Lemma 4.2 with the weight ξ α and χ = ψ , to ω = ω v to obtain (cid:107) ¯ w − ¯ Q ¯ w (cid:107) L (¯ ω v ,ξ α ) (cid:46) (cid:107) ¯ ∇ ¯ w (cid:107) L (¯ ω v ,ξ α ) , where the hidden constant depends only on α , σ Ω , σ and ψ , but not on v (cid:48)(cid:48) and h v (cid:48)(cid:48) .Applying this to (4.21), together with a change of variables with F − v , we get (4.18).The proof of (4.19) is similar. Notice that w v ( z ) = (cid:90) ω v ( w ( x ) + ∇ w ( x ) · ( z − x )) ψ v ( x ) d x = (cid:90) ω v (cid:0) ¯ w (¯ x ) + ¯ ∇ ¯ w (¯ x ) · (¯ z − ¯ x ) (cid:1) ψ (¯ x ) d¯ x =: ¯ w (¯ z ) . Since ∂ ¯ z i ¯ w (¯ z ) = (cid:82) ω v ∂ ¯ x i ¯ w (¯ x ) ψ (¯ x ) d¯ x is constant, we have the vanishing mean valueproperty (cid:90) ω v ∂ ¯ z i ( ¯ w (¯ z ) − ¯ w (¯ z )) ψ (¯ z ) d¯ z = 0 . Finally, applying Lemma 4.2 to ∂ ¯ x i ( ¯ w (¯ x ) − ¯ w (¯ x )), and scaling back via the map F v ,we obtain (4.19).By shape regularity, for all v ∈ N in and T ⊂ ω v , the quantities h v (cid:48) and h v (cid:48)(cid:48) areequivalent to h K and h I , up to a constant that depends only on σ Ω and σ , respectively.This fact leads to interpolation estimates in the weighted L -norm. Theorem 4.5 (
Stability and local interpolation estimates in the weighted L -norm ). For all T ∈ T Y and w ∈ L ( ω T , y α ) we have (cid:107) Π T Y w (cid:107) L ( T,y α ) (cid:46) (cid:107) w (cid:107) L ( ω T ,y α ) . (4.22) If, in addition, w ∈ H ( ω T , y α ) (cid:107) w − Π T Y w (cid:107) L ( T,y α ) (cid:46) h v (cid:48) (cid:107)∇ x (cid:48) w (cid:107) L ( ω T ,y α ) + h v (cid:48)(cid:48) (cid:107) ∂ y w (cid:107) L ( ω T ,y α ) . (4.23)4 R.H. Nochetto, E. Ot´arola, A.J. Salgado
The hidden constants in both inequalities depend only on σ Ω , σ , ψ and α .Proof . Let T be an element of T Y . Assume, for the moment, that Π T Y is uniformlybounded as a mapping from L ( ω T , y α ) to L ( T, y α ), i.e., (4.22).Choose an interior node v of T , i.e., a node v of T such that v ∈ N in . Since Q v w is constant, we deduce Π T Y Q v w = Q v w , whence (cid:107) w − Π T Y w (cid:107) L ( T,y α ) = (cid:107) ( I − Π T Y )( w − Q v w ) (cid:107) L ( T,y α ) (cid:46) (cid:107) w − Q v w (cid:107) L ( ω T ,y α ) , so that (4.23) follows from Corollary 4.3.It remains to show the local boundedness (4.22) of Π T Y . By definition,Π T Y w = n T (cid:88) i =1 w v i ( v i ) λ v i , where { v i } n T i =1 denotes the set of interior vertices of T . By the triangle inequality (cid:107) Π T Y w (cid:107) L ( T,y α ) ≤ n T (cid:88) i =1 (cid:107) w v i (cid:107) L ∞ ( T ) (cid:107) λ v i (cid:107) L ( T,y α ) , (4.24)so that we need to estimate (cid:107) w v i (cid:107) L ∞ ( T ) . This follows from (4.10) along with, (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) ω v i wψ v i (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:107) w (cid:107) L ( ω v i ,y α ) (cid:107) ψ v i (cid:107) L ( ω v i ,y − α ) , (4.25)and, for (cid:96) = 1 , . . . , n + 1, (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) ω v i ∂ x (cid:96) w ( x )( z (cid:96) − x (cid:96) ) ψ v i ( x ) d x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:46) (cid:107) w (cid:107) L ( ω v i ,y α ) (cid:107) ψ v i (cid:107) L ( ω v i ,y − α ) . (4.26)We get (4.26) upon integration by parts, ψ v i = 0 on ∂ω v i , and | z (cid:96) − x (cid:96) | (cid:46) h K ≈ h v (cid:48) for (cid:96) = 1 , · · · , n and | z n +1 − y | (cid:46) h I ≈ h v (cid:48)(cid:48) . Replacing (4.25) and (4.26) in (4.24), wearrive at (cid:107) Π T Y w (cid:107) L ( T,y α ) (cid:46) (cid:107) w (cid:107) L ( ω T ,y α ) n T (cid:88) i =1 (cid:107) λ v i (cid:107) L ( T,y α ) (cid:107) ψ v i (cid:107) L ( ω v i ,y − α ) (cid:46) (cid:107) w (cid:107) L ( ω T ,y α ) , where the last inequality is a consequence of λ v i and ψ being bounded in L ∞ ( ω T ), (cid:107) λ v i (cid:107) L ( T,y α ) (cid:107) ψ v i (cid:107) L ( ω v i ,y − α ) (cid:46) | ω v i | − (cid:32)(cid:90) ω v i | y | α (cid:90) ω v i | y | − α (cid:33) / , together with | y | α ∈ A ( R n +1 ); see (2.15). H -based interpolation estimates on interior elements. Here we prove interpolation estimates on the first derivatives for interior elements.The, rather technical, proof is an adaption of [28, Theorem 2.6] to our particulargeometric setting. In contrast to [28, Theorem 2.6], we do not have the symmetries of acube. However, exploiting the Cartesian product structure of the elements T = K × I , PDE approach to fractional diffusion Fig. 4.1 . A generic element T = K × I in three dimensions: a quadrilateral prism. we are capable of handling the anisotropy in the extended variable y for general shape-regular graded meshes T Y . This is the content of the following result. Theorem 4.6 (
Stability and local interpolation: interior elements ). Let T ∈ T Y be such that ∂T ∩ Γ D = ∅ . For all w ∈ H ( ω T , y α ) we have the stability bounds (cid:107)∇ x (cid:48) Π T Y w (cid:107) L ( T,y α ) (cid:46) (cid:107)∇ x (cid:48) w (cid:107) L ( ω T ,y α ) , (4.27) (cid:107) ∂ y Π T Y w (cid:107) L ( T,y α ) (cid:46) (cid:107) ∂ y w (cid:107) L ( ω T ,y α ) , (4.28) and, for all w ∈ H ( ω T , y α ) and j = 1 , . . . , n + 1 we have the error estimates (cid:107) ∂ x j ( w − Π T Y w ) (cid:107) L ( T,y α ) (cid:46) h v (cid:48) (cid:107)∇ x (cid:48) ∂ x j w (cid:107) L ( ω T ,y α ) + h v (cid:48)(cid:48) (cid:107) ∂ y ∂ x j w (cid:107) L ( ω T ,y α ) . (4.29) Proof . To exploit the particular structure of T , we label its vertices in an appro-priate way; see Figure 4.1 for the three-dimensional case. In general, if T = K × [ a, b ],we first assign a numbering { v k } k =1 ,..., n to the nodes that belong to K × { a } . If(˜ v (cid:48) , b ) is a vertex in K × { b } , then there is a v k ∈ K × { a } such that ˜ v (cid:48) = v (cid:48) k , and weset v k +2 n = ˜ v . We proceed in three steps.1 Derivative ∂ y in the extended dimension. We wish to obtain a bound for thenorm (cid:107) ∂ y ( w − Π T Y w ) (cid:107) L ( T,y α ) . Since, w − Π T Y w = ( w − w v ) + ( w v − Π T Y w ) andan estimate for the difference w − w v is given in Lemma 4.4, it suffices to consider q := w v − Π T Y w ∈ Q . Thanks to the special labeling of the nodes and the tensorproduct structure of the elements, i.e., ∂ y λ v i +2 n = − ∂ y λ v i , we get ∂ y q = n +1 (cid:88) i =1 q ( v i ) ∂ y λ v i = n (cid:88) i =1 ( q ( v i ) − q ( v i +2 n )) ∂ y λ v i , so that (cid:107) ∂ y q (cid:107) L ( T,y α ) ≤ n (cid:88) i =1 | q ( v i ) − q ( v i +2 n ) |(cid:107) ∂ y λ v i (cid:107) L ( T,y α ) . (4.30)To estimate the differences | q ( v i ) − q ( v i +2 n ) | for i = 1 , · · · , n we may, without loss ofgenerality, set i = 1. By the definitions of Π T Y and q , we have Π T Y w ( v ) = w v ( v ),whence δq ( v ) := q ( v ) − q ( v n ) = w v n ( v n ) − w v ( v n ) , R.H. Nochetto, E. Ot´arola, A.J. Salgado and by the definition (4.10) of the averaged Taylor polynomial, we have δq ( v ) = (cid:90) P ( x, v n ) ψ v n ( x ) d x − (cid:90) P ( x, v n ) ψ v ( x ) d x. (4.31)Recalling the operator (cid:12) , introduced in (2.4), we notice that, for h v = ( h v (cid:48) , h v (cid:48)(cid:48) )and z ∈ R n +1 , the vector h v (cid:12) z is uniformly equivalent to ( h K z (cid:48) , h I z (cid:48)(cid:48) ) for all T = K × I in the star ω v . Changing variables in (4.31) yields δq ( v ) = (cid:90) (cid:0) P ( v n − h v n (cid:12) z, v n ) − P ( v − h v (cid:12) z, v n ) (cid:1) ψ ( z ) d z. (4.32)To estimate this expression define θ = (0 , θ (cid:48)(cid:48) ) = (cid:16) , v (cid:48)(cid:48) n − v (cid:48)(cid:48) + ( h v (cid:48)(cid:48) − h v (cid:48)(cid:48) n ) z (cid:48)(cid:48) (cid:17) , (4.33)and F z ( t ) = P ( v − h v (cid:12) z + tθ, v n ). Using that v (cid:48) = v (cid:48) n and h v (cid:48) = h v (cid:48) n , weeasily obtain P ( v n − h v n (cid:12) z, v n ) − P ( v − h v (cid:12) z, v n ) = F z (1) − F z (0) . Consequently, δq ( v ) = (cid:90) (cid:90) F (cid:48) z ( t ) ψ ( z ) d t d z = (cid:90) (cid:90) F (cid:48) z ( t ) ψ ( z ) d z d t, (4.34)and since ψ is bounded in L ∞ and supp ψ = D ⊂ B × ( − , I ( t ) = (cid:90) D | F (cid:48) z ( t ) | d z, ≤ t ≤ . Invoking the definitions of F z and P ( x, y ), we deduce F (cid:48) z ( t ) = ∇ x P ( v − h v (cid:12) z + tθ, v n ) · θ, and ∇ x P ( x, v ) = D w ( x ) · ( v − x ) . Using these two expressions, we arrive at I ( t ) ≤ (cid:90) D (cid:0)(cid:12)(cid:12) ∂ yy w ( v − h v (cid:12) z + tθ ) (cid:12)(cid:12) (cid:12)(cid:12) v (cid:48)(cid:48) n − v (cid:48)(cid:48) + h v (cid:48)(cid:48) z (cid:48)(cid:48) − tθ (cid:48)(cid:48) (cid:12)(cid:12) + | ∂ y ∇ x (cid:48) w ( v − h v (cid:12) z + tθ ) | | v (cid:48) n − v (cid:48) + h v (cid:48) z (cid:48) | (cid:1) | θ (cid:48)(cid:48) | d z, Now, since | z (cid:48) | , | z (cid:48)(cid:48) | ≤ ≤ t ≤
1, we see that | v (cid:48) n − v (cid:48) + h v (cid:48) z (cid:48) | (cid:46) h v (cid:48) , | v (cid:48)(cid:48) n − v (cid:48)(cid:48) + h v (cid:48)(cid:48) z (cid:48)(cid:48) − tθ (cid:48)(cid:48) | (cid:46) h v (cid:48)(cid:48) . Consequently, I ( t ) (cid:46) (cid:90) D (cid:16)(cid:12)(cid:12) ∂ yy w ( v − h v (cid:12) z + tθ ) (cid:12)(cid:12) h v (cid:48)(cid:48) + | ∂ y ∇ x (cid:48) w ( v − h v (cid:12) z + tθ ) | h v (cid:48) h v (cid:48)(cid:48) (cid:1) d z. PDE approach to fractional diffusion τ = v − h v (cid:12) z + tθ , we obtain I ( t ) (cid:46) (cid:90) ω T (cid:32) h v (cid:48)(cid:48) h n v (cid:48) (cid:12)(cid:12) ∂ yy w ( τ ) (cid:12)(cid:12) + 1 h n − v (cid:48) | ∂ y ∇ x (cid:48) w ( τ ) | (cid:33) d τ, (4.35)because the support D of ψ is contained in B /σ Ω × ( − /σ Y , /σ Y ), and so is mappedinto ω v ⊂ ω T . Notice also that h v (cid:48)(cid:48) (cid:46) (1 − t ) h v (cid:48)(cid:48) + th v (cid:48)(cid:48) n . This implies I ( t ) (cid:46) (cid:32) h v (cid:48)(cid:48) h n v (cid:48) (cid:107) ∂ yy w (cid:107) L ( ω T ,y α ) + 1 h n − v (cid:48) (cid:107)∇ x (cid:48) ∂ y w (cid:107) L ( ω T ,y α ) (cid:33) (cid:107) (cid:107) L ( ω T ,y − α ) , (4.36)which, together with (4.34), yields | δq ( v ) |(cid:107) ∂ y λ v (cid:107) L ( T,y α ) (cid:46) (cid:32) h v (cid:48)(cid:48) h n v (cid:48) (cid:107) ∂ yy w (cid:107) L ( ω T ,y α ) + 1 h n − v (cid:48) (cid:107)∇ x (cid:48) ∂ y w (cid:107) L ( ω T ,y α ) (cid:33) · (cid:107) (cid:107) L ( ω T ,y − α ) (cid:107) ∂ y λ v (cid:107) L ( T,y α ) . (4.37)Since | y | α ∈ A ( R n +1 ), we have (cid:107) (cid:107) L ( ω T ,y − α ) (cid:107) ∂ y λ v (cid:107) L ( T,y α ) (cid:46) h n v (cid:48) h v (cid:48)(cid:48) (cid:18)(cid:90) I y − α (cid:19) (cid:18)(cid:90) I y α (cid:19) (cid:46) h n v (cid:48) . (4.38)Replacing this into (4.37), we obtain | δq ( v ) |(cid:107) ∂ y λ v (cid:107) L ( T,y α ) (cid:46) h (cid:48) v (cid:107)∇ x (cid:48) ∂ y w (cid:107) L ( ω T ,y α ) + h v (cid:48)(cid:48) (cid:107) ∂ yy w (cid:107) L ( ω T ,y α ) , (4.39)which, in this case, implies (4.29).2 Derivatives ∇ x (cid:48) in the domain Ω . To prove an estimate for ∇ x (cid:48) ( w − Π T Y w ) wenotice that, given a vertex v , the associated basis function λ v can be written as λ v ( x ) = Λ v (cid:48) ( x (cid:48) ) µ v (cid:48)(cid:48) ( y ), where Λ v (cid:48) is the canonical Q basis function on the variable x (cid:48) associated to the node v (cid:48) in the triangulation T Ω , and µ v (cid:48)(cid:48) corresponds to thepiecewise P basis function associated to the node v (cid:48)(cid:48) . Recall that, by construction,the basis { Λ i } n i =1 possesses the so-called partition of unity property, i.e., n (cid:88) i =1 Λ v ( x (cid:48) ) = 1 ∀ x (cid:48) ∈ K, = ⇒ n (cid:88) i =1 ∇ x (cid:48) Λ v ( x (cid:48) ) = 0 ∀ x (cid:48) ∈ K. This implies that, for every q ∈ Q ( T ), ∇ x (cid:48) q = n +1 (cid:88) i =1 q ( v i ) ∇ x (cid:48) λ v i = n (cid:88) i =1 (cid:16) q ( v i ) µ v (cid:48)(cid:48) i ( y ) + q ( v i +2 n ) µ v (cid:48)(cid:48) i +2 n ( y ) (cid:17) ∇ x (cid:48) Λ v (cid:48) i ( x (cid:48) )= n (cid:88) i =1 (cid:104) ( q ( v i ) − q ( v )) µ v (cid:48)(cid:48) i ( y ) + ( q ( v i +2 n ) − q ( v n )) µ v (cid:48)(cid:48) i +2 n ( y ) (cid:105) ∇ x (cid:48) Λ v i ( x (cid:48) ) , so that (cid:107)∇ x (cid:48) q (cid:107) L ( T,y α ) (cid:46) n (cid:88) i =1 | q ( v i ) − q ( v ) |(cid:107) µ v (cid:48)(cid:48) i ∇ x (cid:48) Λ v i (cid:107) L ( T,y α ) + n (cid:88) i =1 | q ( v n ) − q ( v i +2 n ) |(cid:107) µ v (cid:48)(cid:48) i +2 n ∇ x (cid:48) Λ v i (cid:107) L ( T,y α ) . R.H. Nochetto, E. Ot´arola, A.J. Salgado
This expression shows that the same techniques developed for the previous step allowsus to obtain (4.29).3
Stability . It remains to prove (4.27) and (4.28). By the triangle inequality, (cid:107) ∂ y Π T Y w (cid:107) L ( T,y α ) ≤ (cid:107) ∂ y ( w − Π T Y w ) (cid:107) L ( T,y α ) + (cid:107) ∂ y w (cid:107) L ( T,y α ) , so that it suffices to estimate the first term. Add and subtract w v , (cid:107) ∂ y ( w − Π T Y w ) (cid:107) L ( T,y α ) ≤ (cid:107) ∂ y ( w − w v ) (cid:107) L ( T,y α ) + (cid:107) ∂ y ( w v − Π T Y w ) (cid:107) L ( T,y α ) . (4.40)Let us estimate the first term. The definition of ψ v , together with | y | α ∈ A ( R n +1 )implies (cid:107) ψ v (cid:107) L ( ω v ,y − α ) (cid:107) (cid:107) L ( ω v ,y α ) (cid:46) , whence invoking the definition (4.10) of theregularized Taylor polynomial w v yields (cid:107) ∂ y w v (cid:107) L ( T,y α ) ≤ (cid:107) ∂ y w (cid:107) L ( ω v ,y α ) , and (cid:107) ∂ y ( w − w v ) (cid:107) L ( T,y α ) (cid:46) (cid:107) ∂ y w (cid:107) L ( T,y α ) . (4.41)To estimate the second term of the right hand side of (4.40), we repeat the stepsused to obtain (4.29), starting from (4.31). Integrating by parts and using that ψ v i = 0on ∂ω v i , we get, for (cid:96) = 1 , . . . , n + 1, (cid:90) ω v i ∂ x (cid:96) w ( x )( z (cid:96) − x (cid:96) ) ψ v i ( x ) d x = (cid:90) ω v i w ( x ) ψ v i ( x ) d x − (cid:90) ω v i w ( x )( z (cid:96) − x (cid:96) ) ∂ x (cid:96) ψ v i ( x ) d x, whence δq ( v ) = ( n + 2) (cid:18)(cid:90) w ( x ) ψ v n d x − (cid:90) w ( x ) ψ v d x (cid:19) − (cid:90) w ( x )( v n − x ) · ∇ ψ v n ( x ) d x + (cid:90) w ( x )( v − x ) · ∇ ψ v ( x ) d x = I + I . (4.42)To estimate I we consider the same change of variables used to obtain (4.32). Define G z ( t ) = ( n + 2) · w ( v − h v (cid:12) z + tθ ), with θ as in (4.33), and observe that I = (cid:90) (cid:90) G (cid:48) z ( t ) ψ ( z ) d z d t = ( n + 2) (cid:90) (cid:90) ∂ y w ( v − h v (cid:12) z + tθ ) θ (cid:48)(cid:48) ψ ( z ) d z d t. Introducing the change of variables τ = v − h v (cid:12) z + tθ , we obtain | I | (cid:46) (cid:90) ω T h n v (cid:48) | ∂ y w ( τ ) | d τ ≤ h n v (cid:48) (cid:107) ∂ y w (cid:107) L ( ω T ,y α ) (cid:107) (cid:107) L ( ω T ,y − α ) . (4.43)We now estimate I . Changing variables, I = (cid:90) (cid:0) w ( v n − h v n (cid:12) z ) − w ( v − h v (cid:12) z ) (cid:1) z (cid:48) · ∇ x (cid:48) ψ ( z ) d z + (cid:90) (cid:0) w ( v n − h v n (cid:12) z ) z (cid:48)(cid:48) − w ( v − h v (cid:12) z )( ϑ + z (cid:48)(cid:48) ) (cid:1) ∂ y ψ ( z ) d z = I , + I , , PDE approach to fractional diffusion ϑ = ( v n +11+2 n − v n +11 ) /h v (cid:48)(cid:48) . Arguing as in the derivation of (4.43) we obtain | I , | , | I , | (cid:46) (cid:90) ω T h n v (cid:48) | ∂ y w ( τ ) | d τ ≤ h n v (cid:48) (cid:107) ∂ y w (cid:107) L ( ω T ,y α ) (cid:107) (cid:107) L ( ω T ,y − α ) . (4.44)Inserting (4.43) and (4.44) in (4.42) we deduce | δq ( v ) | (cid:46) h n v (cid:48) (cid:107) ∂ y w (cid:107) L ( ω T ,y α ) (cid:107) (cid:107) L ( ω T ,y − α ) , whence | δq ( v ) |(cid:107) ∂ y λ v (cid:107) L ( T,y α ) (cid:46) (cid:107) ∂ y w (cid:107) L ( ω T ,y α ) , (4.45)because h − n v (cid:48) (cid:107) ∂ y λ v (cid:107) L ( ω T ,y α ) (cid:107) (cid:107) L ( ω T ,y − α ) ≤ C . Replacing (4.45) in (4.30), we get (cid:107) ∂ y ( w v − Π T Y w ) (cid:107) L ( T,y α ) (cid:46) (cid:107) ∂ y w (cid:107) L ( ω T ,y α ) , which, together with (4.40) and (4.41), imply the desired result (4.28). Similar argu-ments are used to prove the stability bound (4.27). H -based interpolation estimates on boundary ele-ments. Let us now extend the interpolation estimates of § T in T Y . We define the non-overlapping sets C = { T ∈ T Y : ∂T ∩ Γ D = ∅} , C = { T ∈ T Y : ∂T ∩ ∂ L C Y (cid:54) = ∅} , C = { T ∈ T Y : ∂T ∩ ( ∂ Ω × { Y } ) (cid:54) = ∅} . The elements in C are interior, so the corresponding interpolation estimate is givenin Theorem 4.6. Interpolation estimates on elements in C are a direct consequenceof [28, Theorem 3.1] and Theorem 4.7 below. This is so due to the fact that, since Y >
1, the weight y α over C is no longer singular nor degenerate. It remains only toprovide interpolation estimates for elements in C . Theorem 4.7 (
Local error interpolation estimate: elements in C ). Let T ∈ C and w ∈ H ( ω T , y α ) vanish on ∂T ∩ ∂ L C Y . Then, we have the stability bounds (cid:107)∇ x (cid:48) Π T Y w (cid:107) L ( T,y α ) (cid:46) (cid:107)∇ x (cid:48) w (cid:107) L ( ω T ,y α ) , (4.46) (cid:107) ∂ y Π T Y w (cid:107) L ( T,y α ) (cid:46) (cid:107) ∂ y w (cid:107) L ( ω T ,y α ) , (4.47) If, in addition, w ∈ H ( ω T , y α ) , then, for j = 1 , . . . , n + 1 , (cid:107) ∂ x j ( w − Π T Y w ) (cid:107) L ( T,y α ) (cid:46) h v (cid:48) (cid:107) ∂ x j ∇ x (cid:48) w (cid:107) L ( ω T ,y α ) + h v (cid:48)(cid:48) (cid:107) ∂ x j y w (cid:107) L ( ω T ,y α ) . (4.48) Proof . For simplicity we present the proof in two dimensions. Let T = (0 , a ) × (0 , b ) ∈ C . Notice that over such an element the weight becomes degenerate or singu-lar. Recall the local enumeration of vertices introduced in the proof of Theorem 4.6(see also Figure 4.1). By the definition of Π T Y we haveΠ T Y w | T = w v ( v ) λ v + w v ( v ) λ v , (4.49)0 R.H. Nochetto, E. Ot´arola, A.J. Salgado
The proofs of (4.46) and (4.47) are similar to Step 3 of Theorem 4.6. To show (4.48),we write the local difference between a function and its interpolant as ( w − Π T Y w ) | T =( w − w v ) | T + ( w v − Π T Y ) | T . Proceeding as in the proof of Lemma 4.4, we can bound ∂ x j ( w − w v ) | T for j = 1 ,
2, in the L ( T, y α )-norm, by the right hand side of (4.48)because this is independent of the trace of w . It remains then to derive a bound for( w v − Π T Y w ) | T , for which we consider two separate cases.1 Derivative in the extended direction.
We use w v ∈ Q , (4.49) and Π T Y w ( v ) =Π T Y w ( v ) = 0, to write ∂ y ( w v − Π T Y w ) | T = ( w v ( v ) − w v ( v )) ∂ y λ v + ( w v ( v ) − w v ( v )) ∂ y λ v . Since w ≡ { } × (0 , b ), then ∂ y w ≡ { } × (0 , b ). By the definition of theTaylor polynomial P , given in (4.11), and the fact that v (cid:48) = v (cid:48) , we obtain w v ( v ) − w v ( v ) = ( v (cid:48)(cid:48) − v (cid:48)(cid:48) ) (cid:90) ω T ∂ y w ( x ) ψ v ( x ) d x = ( v (cid:48)(cid:48) − v (cid:48)(cid:48) ) (cid:90) ω T (cid:90) x (cid:48) ∂ x (cid:48) y w ( σ, y ) ψ v ( x (cid:48) , y ) d σ d x (cid:48) d y. Therefore | w v ( v ) − w v ( v ) | (cid:46) h v (cid:48)(cid:48) h v (cid:48) (cid:107) ∂ x (cid:48) y w (cid:107) L ( ω T ,y α ) (cid:107) ψ v (cid:107) L ( ω T ,y − α ) (cid:46) h v (cid:48)(cid:48) h v (cid:48) h v (cid:48) h v (cid:48) h v (cid:48)(cid:48) (cid:32)(cid:90) b y − α d y (cid:33) (cid:107) ∂ x (cid:48) y w (cid:107) L ( ω T ,y α ) . Since, in view of the weak shape regularity assumption on the mesh T Y , h v (cid:48) ≈ h v (cid:48) , h v (cid:48)(cid:48) = h v (cid:48)(cid:48) , and y α ∈ A ( R n +1+ ), we conclude that | w v ( v ) − w v ( v ) |(cid:107) ∂ y λ v (cid:107) L ( T,y α ) (cid:46) h v (cid:48) h v (cid:48)(cid:48) (cid:32)(cid:90) b y − α d y (cid:90) b y α d y (cid:33) ×× (cid:107) ∂ x (cid:48) y w (cid:107) L ( ω T ,y α ) (cid:46) h v (cid:48) (cid:107) ∂ x (cid:48) y w (cid:107) L ( ω T ,y α ) . (4.50)Finally, to bound w v ( v ) − w v ( v ), we proceed as in Step 1 of the proof of Theo-rem 4.6, which is valid regardless of the trace of w , and deduce | w v ( v ) − w v ( v ) |(cid:107) ∂ y λ v (cid:107) L ( T,y α ) (cid:46) h v (cid:48) (cid:107) ∂ x (cid:48) y w (cid:107) L ( ω T ,y α ) + h v (cid:48)(cid:48) (cid:107) ∂ yy w (cid:107) L ( ω T ,y α ) . This, in conjunction with the previous estimate, yields (4.48) for the derivative in theextended direction.2
Derivative in the x (cid:48) direction. To estimate ∂ x (cid:48) ( w v − Π T Y w ) | T we proceed asin Theorem 4.6 and [28, Theorem 3.1], but we cannot exploit the symmetry of thetensor product structure now. For brevity, we shall only point out the main technicaldifferences. Using, again, that ( w v − Π T Y w ) ∈ Q , ∂ x (cid:48) ( w v − Π T Y w ) | T = w v ( v ) ∂ x (cid:48) λ v + w v ( v ) ∂ x (cid:48) λ v + ( w v ( v ) − w v ( v )) ∂ x (cid:48) λ v = w v ( v ) ∂ x (cid:48) λ v + ( w v ( v ) − w v ( v )) ∂ x (cid:48) λ v − ( w v ( v ) − w v ( v )) ∂ x (cid:48) λ v − w v ( v ) ∂ x (cid:48) λ v = J ( w v , w v ) ∂ x (cid:48) λ v + w v ( v ) ∂ x (cid:48) λ v − w v ( v ) ∂ x (cid:48) λ v , PDE approach to fractional diffusion J ( w v , w v ) = ( w v ( v ) − w v ( v )) − ( w v ( v ) − w v ( v )) . Define θ = (0 , θ (cid:48)(cid:48) ) = (0 , v − v − ( h v (cid:48)(cid:48) − h v (cid:48)(cid:48) ) z (cid:48)(cid:48) ), and rewrite J ( w v , w v ) as follows: J ( w v , w v ) = ( v (cid:48) − v (cid:48) ) (cid:90) D ( ∂ x (cid:48) w ( v − h v (cid:12) z ) − ∂ x (cid:48) w ( v − h v (cid:12) z )) ψ ( z ) d z = − ( v (cid:48) − v (cid:48) ) (cid:90) D (cid:90) ∂ x (cid:48) y w ( v − h v (cid:12) z + θt ) θ (cid:48)(cid:48) ψ ( z ) d t d z, where D = supp ψ . Denote I ( t ) = (cid:90) | ∂ x (cid:48) y w ( v − h v (cid:12) z + θt ) θ (cid:48)(cid:48) | d z. Using the change of variables z (cid:55)→ τ = v − h v (cid:12) z + θt , results in | I ( t ) | (cid:46) h v (cid:48) (cid:90) ω T | ∂ x (cid:48) y w ( τ ) | ψ ( τ ) d τ (cid:46) h v (cid:48) (cid:107) ∂ x (cid:48) y w (cid:107) L ( ω T ,y α ) (cid:107) (cid:107) L ( ω T ,y − α ) (cid:46) h − v (cid:48) (cid:107) ∂ x (cid:48) y w (cid:107) L ( ω T ,y α ) (cid:32)(cid:90) b y − α d y (cid:33) , whence | J ( w v , w v ) | (cid:46) h v (cid:48) (cid:107) ∂ x (cid:48) y w (cid:107) L ( ω T ,y α ) (cid:16)(cid:82) b y − α d y (cid:17) . This implies (cid:107) J ( w v , w v ) ∂ x (cid:48) λ v (cid:107) L ( T,y α ) (cid:46) (cid:32)(cid:90) b y − α d y (cid:33) (cid:32)(cid:90) b y α d y (cid:33) (cid:107) ∂ x (cid:48) y w (cid:107) L ( ω T ,y α ) (cid:46) h v (cid:48)(cid:48) (cid:107) ∂ x (cid:48) y w (cid:107) L ( ω T ,y α ) , which follows from the fact that y α ∈ A ( R + ), and then (4.48) holds true.The estimate of w v ( v ) ∂ x (cid:48) λ v exploits the fact that the trace of w vanishes on ∂ L C Y ; the same happens with w v ( v ) ∂ x (cid:48) λ v . In fact, we can write w v ( v ) = (cid:90) ω v (cid:90) x (cid:48) ( ∂ x (cid:48) w ( τ, y ) − ∂ x (cid:48) w ( x (cid:48) , y )) ψ v ( x (cid:48) , y ) d τ d x (cid:48) d y + (cid:90) ω v ( ∂ y w (0 , y ) − ∂ y w ( x (cid:48) , y )) yψ v ( x (cid:48) , y ) d x (cid:48) d y. To derive (4.48) we finally proceed as in the proofs of Theorem 4.6 and [28,Theorem 3.1]. We omit the details.
5. Error estimates.
The estimates of § § w ∈ H ( ω T , y α ). However, the solution u of (2.26) satisfies u yy ∈ L ( C , y β ) only when β > α + 1, according to Theorem 2.6. For this reason, inthis section we derive error estimates for both quasi-uniform and graded meshes. Theestimates of § § R.H. Nochetto, E. Ot´arola, A.J. Salgado
We start with a simple one dimensional case( n = 1) and assume that we need to approximate over the interval [0 , Y ] the function w ( y ) = y − α . Notice that w y ( y ) ≈ y − α as y ≈ + has the same behavior as thederivative in the extended direction of the α -harmonic extension u .Given M ∈ N we consider the uniform partition of the interval [0 , Y ] y k = kM Y , k = 0 , . . . , M. (5.1)and corresponding elements I k = [ y k , y k +1 ] of size h k = h = Y /M for k = 0 , . . . , M − T Y of § E k = (cid:107) ∂ y ( w − Π T Y w ) (cid:107) L ( I k ,y α ) . For k = 2 , . . . , M −
1, since y ≥ h and α < α + 1 < β , (4.29) implies E k (cid:46) h (cid:90) ω Ik y α | w yy | d y (cid:46) h α − β (cid:90) ω Ik y β | w yy | d y, (5.2)because (cid:0) yh (cid:1) α ≤ (cid:0) yh (cid:1) β . The estimate for E + E follows from from the stability ofthe operator Π T Y (4.28) and (4.47): E + E (cid:46) (cid:90) h y α | w y | (cid:46) h − α , (5.3)because w ( y ) ≈ y − α as y ≈ + . Using (5.2) and (5.3) in conjunction with 2 + α − β < − α , we obtain a global interpolation estimate (cid:107) ∂ y ( w − Π T Y w ) (cid:107) L ((0 , Y ) ,y α ) (cid:46) h (2+ α − β ) / . (5.4)These ideas can be extended to prove an error estimate for u on uniform meshes. Theorem 5.1 (
Error estimate for quasi-uniform meshes ). Let u solve (2.26) ,and V T Y be the solution of (4.4) , constructed over a quasi-uniform mesh of size h . If Y ≈ | log h | , then for all (cid:15) > (cid:107)∇ ( u − V T Y ) (cid:107) L ( C Y ,y α ) (cid:46) h s − ε (cid:107) f (cid:107) H − s (Ω) . (5.5) where the hidden constant blows up if ε tends to .Proof . Use first Theorem 3.5 and Theorem 4.1, combined with (4.9), to reduce theapproximation error to the interpolation error of u . Repeat next the steps leading to(5.2)–(5.3), but combining the interpolation estimates of Theorems 4.6 and 4.7 withthe regularity results of Theorem 2.6. Remark 5.2 (
Sharpness of (5.5) for s (cid:54) = ). According to (2.34) and (2.37), ∂ y u ≈ y − α , and this formally implies ∂ y u ∈ H s − ε ( C , y α ) for all ε > f ∈ H − s (Ω). In this sense (5.5) appears to be sharp with respect to regularity eventhough it does not exhibit the optimal rate. We verify this argument via a simplenumerical illustration for dimension n = 1. We let Ω = (0 , s = 0 .
2, right hand side f = π s sin( πx ), and note that u ( x ) = sin( πx ), and the solution u to (1.2) is u ( x, y ) = 2 − s π s Γ( s ) sin( πx ) K s ( πy ) . Figure 5.1 shows the rate of convergence for the H ( C Y , y α )-seminorm. Estimate(5.5) predicts a rate of h − . − ε . We point out that for the α -harmonic extension PDE approach to fractional diffusion −0.8 −0.6 −0.4 −0.2 Degrees of Freedom (DOFs) E rr o r DOFs − . k e k H ( y α ) Fig. 5.1 . Computational rate of convergence for quasi-uniform meshes, s = 0 . , and n = 1 . we are solving a two dimensional problem and, since the mesh T Y is quasi-uniform, T Y ≈ h − . In other words the rate of convergence, when measured in term ofdegrees of freedom, is ( T Y ) − . − ε , which is what Figure 5.1 displays. Remark 5.3 (
Case s = ). Estimate (5.5) does not hold for s = . In thiscase there is no weight and the scaling issues in (5.2) are no longer present, so that E k (cid:46) h (cid:107) v (cid:107) H ( I k ) . We thus obtain the optimal error estimate (cid:107)∇ ( u − V T Y ) (cid:107) L ( C Y ) (cid:46) h (cid:107) f (cid:107) H / (Ω) . The estimate (5.5) can be written equivalently (cid:107)∇ ( u − V T Y ) (cid:107) L ( C Y ,y α ) (cid:46) ( T Y ) − s − εn +1 (cid:107) f (cid:107) H − s (Ω) , for quasi-uniform meshes in dimension n + 1. We now show how to compensate thesingular behavior in the extended variable y by anisotropic meshes and restore theoptimal convergence rate − / ( n + 1).As in § n = 1 with the function w ( y ) = y − α over [0 , Y ]. We consider the graded partition T Y of the interval [0 , Y ] y k = (cid:18) kM (cid:19) γ Y , k = 0 , . . . , M, (5.6)where γ = γ ( α ) > / (1 − α ) >
1. If we denote by h k the length of the interval I k = [ y k , y k +1 ] = (cid:20)(cid:18) kM (cid:19) γ Y , (cid:18) k + 1 M (cid:19) γ Y (cid:21) , then h k = y k +1 − y k (cid:46) Y M γ k γ − , k = 1 , . . . , M − . We again consider the operator Π T Y of § T Y andwish to bound the local interpolation errors E k of § R.H. Nochetto, E. Ot´arola, A.J. Salgado interior elements to obtain that, for k = 2 , . . . , M − E k (cid:46) h k (cid:90) ω Ik y α | w yy | d y (cid:46) Y k γ − M γ (cid:90) ω Ik y α | w yy | d y (cid:46) Y α − β k γ − M γ (cid:18) kM (cid:19) γ ( α − β ) (cid:90) ω Ik y β | w yy | d y (cid:46) Y − α k γ (1 − α ) − M γ (1 − α ) . (5.7)because y α (cid:46) (cid:0) kM (cid:1) γ ( α − β ) Y α − β y β . Adding (5.7) over k = 2 , . . . , M −
1, and usingthat γ (1 − α ) >
3, we arrive at (cid:107) ∂ y ( w − Π T Y w ) (cid:107) L (( y , Y ) ,y α ) (cid:46) Y − α M − . (5.8)For the errors E , E we resort to the stability bounds (4.28) and (4.47) to write (cid:107) ∂ y ( w − Π T Y w ) (cid:107) L ((0 ,y ) ,y α ) (cid:46) (cid:90) ( M ) γ Y y − α d y (cid:46) Y − α M γ (1 − α ) , (5.9)where we have used (5.6). Finally, adding (5.8) and (5.9) gives (cid:107) ∂ y ( w − Π T Y w ) (cid:107) L ((0 , Y ) ,y α ) (cid:46) Y − α M − , and shows that the interpolation error exhibits optimal decay rate.We now apply this idea to the numerical solution of problem (3.3). We assume T Ω to be quasi-uniform in T Ω with T Ω ≈ M n and construct T Y ∈ T as the tensorproduct of T Ω and the partition given in (5.6), with γ > / (1 − α ). Consequently, T Y = M · T Ω ≈ M n +1 . Finally, we notice that since T Ω is shape regular andquasi-uniform, h T Ω ≈ ( T Ω ) − /n ≈ M − . Theorem 5.4 (
Error estimate for graded meshes ). Let V T ∈ V ( T Y ) solve (4.4) and U T Ω ∈ U ( T Ω ) be defined as in (4.5) . Then (cid:107) u − V T Y (cid:107) ◦ H L ( C ,y α ) (cid:46) e −√ λ Y / (cid:107) f (cid:107) H s (Ω) (cid:48) + Y (1 − α ) / ( T Y ) − / ( n +1) (cid:107) f (cid:107) H − s (Ω) , (5.10) Proof . In light of (4.9), with (cid:15) ≈ e −√ λ Y / , it suffices to bound the interpolationerror u − Π T Y u on the mesh T Y . To do so we, first of all, notice that if I and I areneighboring cells on the partition of [0 , Y ], then there is a constant σ = σ ( γ ) suchthat h I ≤ σh I , whence the weak regularity condition ( c ) holds. We can thus applythe polynomial interpolation theory of § T Y into thesets T := (cid:8) T ∈ T Y : ω T ∩ ( ¯Ω × { } ) = ∅ (cid:9) , T := (cid:8) T ∈ T Y : ω T ∩ ( ¯Ω × { } ) (cid:54) = ∅ (cid:9) . We observe that for all T = K × I k ∈ T we have k ≥ y α (cid:46) (cid:0) kM (cid:1) γ ( α − β ) Y α − β y β .Applying Theorem 4.6 and Theorem 4.7 to elements in T we obtain (cid:88) T ∈T (cid:107)∇ ( u − Π T Y u ) (cid:107) L ( T,y α ) (cid:46) (cid:88) T = K × I ∈T (cid:16) h K (cid:107)∇ x (cid:48) ∇ u (cid:107) L ( ω T ,y α ) + h I (cid:107) ∂ y ∇ x (cid:48) u (cid:107) L ( ω T ,y α ) + h I (cid:107) ∂ yy u (cid:107) L ( ω T ,y β ) (cid:17) = S + S + S . PDE approach to fractional diffusion S , which we rewrite as follows: S (cid:46) M (cid:88) k =2 Y α − β k γ − M γ (cid:18) kM (cid:19) γ ( α − β ) (cid:90) b k a k y β (cid:90) Ω | ∂ yy u | d x (cid:48) d y, with a k = (cid:0) k − M (cid:1) γ Y and b k = (cid:0) k +1 M (cid:1) γ Y . We now invoke the local estimate (2.43), aswell as the fact that b k − a k (cid:46) (cid:0) kM (cid:1) γ − Y M , to end up with S (cid:46) M (cid:88) k =2 Y − α k γ (1 − α ) − M γ (1 − α ) (cid:107) f (cid:107) L (Ω) (cid:46) Y − α M − (cid:107) f (cid:107) L (Ω) . We now handle the middle term S with the help of (2.42), which is valid for b k ≤ k ≤ k ≤ M Y − /γ , whereas for k > k we know that theestimate decays exponentially. We thus have S (cid:46) (cid:107) f (cid:107) H − s (Ω) k (cid:88) k =2 (cid:32)(cid:18) kM (cid:19) γ − Y M (cid:33) (cid:46) Y /γ M (cid:107) f (cid:107) H − s (Ω) (cid:46) Y − α M (cid:107) f (cid:107) H − s (Ω) . The first term S is easy to estimate. Since h K ≤ M − for all K ∈ T Ω , we get S (cid:46) M − (cid:107)∇ x (cid:48) ∇ v (cid:107) L ( C Y ,y α ) (cid:46) M − (cid:107) f (cid:107) H − s (Ω) (cid:46) Y − α M − (cid:107) f (cid:107) H − s (Ω) . For elements in T , we rely on the stability estimates (4.27), (4.28), (4.46) and(4.47) of Π T Y and thus repeat the arguments used to derive (5.8) and (5.9). Addingthe estimates for T and T we obtain the assertion. Remark 5.5 (
Choice of Y ). A natural choice of Y comes from equilibrating thetwo terms on the right-hand side of (5.10): (cid:15) ≈ T Y ) − n +1 ⇔ Y ≈ log( T Y )) . This implies the near-optimal estimate (cid:107) u − V T Y (cid:107) ◦ H L ( C ,y α ) (cid:46) | log( T Y ) | s · ( T Y ) − / ( n +1) (cid:107) f (cid:107) H − s (Ω) . (5.11) Remark 5.6 (
Estimate for u ). In view of (4.6), we deduce the energy estimate (cid:107) u − U T Ω (cid:107) H s (Ω) (cid:46) | log( T Y ) | s · ( T Y ) − / ( n +1) (cid:107) f (cid:107) H − s (Ω) . We can rewrite this estimate in terms of regularity u ∈ H s (Ω) and T Ω as (cid:107) u − U T Ω (cid:107) H s (Ω) (cid:46) | log( T Ω ) | s · ( T Ω ) − /n (cid:107) u (cid:107) H s (Ω) . and realize that the order is near-optimal given the regularity shift from left to right.However, our PDE approach does not allow for a larger rate ( T Ω ) (2 − s ) /n that wouldstill be compatible with piecewise bilinear polynomials but not with (5.11). Remark 5.7 (
Computational complexity ). The cost of solving the discrete prob-lem (4.4) is related to T Y , and not to T Ω , but the resulting system is sparse. Thestructure of (4.4) is so that fast multilevel solvers can be designed with complex-ity proportional to T Y . On the other hand, using an integral formulation requiressparsification of an otherwise dense matrix with associated cost ( T Ω ) . Remark 5.8 (
Fractional regularity ). The function u , solution the α -harmonicextension problem, may also have singularities in the direction of the x (cid:48) -variables and6 R.H. Nochetto, E. Ot´arola, A.J. Salgado thus exhibit fractional regularity. This depends on Ω and the right hand side f (seeRemark 2.7). The characterization of such singularities is as yet an open problemto us. The polynomial interpolation theory developed in § T Ω , which can resolve such singularities, provided wemaintain the Cartesian structure of T Y . The corresponding a posteriori error analysisis an entirely different but important direction currently under investigation.
6. Numerical experiments for the fractional Laplacian.
To illustrate theproposed techniques here we present a couple of numerical examples. The implemen-tation has been carried out with the help of the deal.II library (see [6, 7]) which,by design, is based on tensor product elements and thus is perfectly suitable for ourneeds. The main concern while developing the code was correctness and, therefore,integrals are evaluated numerically with Gaussian quadratures of sufficiently high or-der and linear systems are solved using CG with ILU preconditioner with the exitcriterion being that the (cid:96) -norm of the residual is less than 10 − . More efficienttechniques for quadrature and preconditioning are currently under investigation. Let Ω = (0 , . It is common knowledge that ϕ m,n ( x , x ) = sin( mπx ) sin( nπx ) , λ m,n = π (cid:0) m + n (cid:1) , m, n ∈ N . If f ( x , x ) = (2 π ) s sin( πx ) sin( πx ), by (2.12) we have u ( x , x ) = sin( πx ) sin( πx ) , and, by (2.24), u ( x , x , y ) = 2 − s Γ( s ) (2 π ) s/ sin( πx ) sin( πx ) y s K s ( √ πy ) . We construct a sequence of meshes { T Y k } k ≥ , where the triangulation of Ω isobtained by uniform refinement and the partition of [0 , Y k ] is as in § , Y k ]is divided with mesh points given by (5.6) with the election of the parameter γ > / (1 − α ). On the basis of Theorem 3.5, for each mesh the truncation parameter Y k is chosen so that (cid:15) ≈ ( T Y k − ) − / . This can be achieved, for instance, by setting Y k ≥ Y ,k = 2 √ λ (log C − log (cid:15) ) . With this type of meshes, (cid:107) u − U T Ω ,k (cid:107) H s (Ω) (cid:46) (cid:107) u − V T Y k (cid:107) ◦ H L ( C ,y α ) (cid:46) | log( T Y k ) | s · ( T Y k ) − / , which is near-optimal in u but suboptimal in u , since we should expect (see [15]) (cid:107) u − U T Ω ,k (cid:107) H s (Ω) (cid:46) h − s T Ω (cid:46) ( T Y k ) − (2 − s ) / . Figure 6.1 shows the rates of convergence for s = 0 . s = 0 . PDE approach to fractional diffusion −2 −1 Degrees of Freedom (DOFs) E rr o r DOFs − / k e k H ( y α ) −2 −1 Degrees of Freedom (DOFs) E rr o r DOFs − / k e k H ( y α ) Fig. 6.1 . Computational rate of convergence for the approximate solution of the fractionalLaplacian over a square with graded meshes on the extended dimension. The left panel shows therate for s = 0 . and the right one for s = 0 . . In both cases, the rate is ≈ ( T Y k ) − / in agreementwith Theorem 5.4 and Remark 5.5 Let Ω = {| x (cid:48) | ∈ R : | x (cid:48) | < } . Using polar coordi-nates it can be shown that ϕ m,n ( r, θ ) = J m ( j m,n r ) ( A m,n cos( mθ ) + B m,n sin( mθ )) , (6.1)where J m is the m -th Bessel function of the first kind; j m,n is the n -th zero of J m and A m,n , B m,n are real normalization constants that ensure (cid:107) ϕ m,n (cid:107) L (Ω) = 1 for all m, n ∈ N . It is also possible to show that λ m,n = ( j m,n ) .If f = ( λ , ) s ϕ , , then (2.12) and (2.24) show that u = ϕ , and u ( r, θ, y ) = 2 − s Γ( s ) ( λ , ) s/ ϕ , ( r, θ ) y s K s ( √ πy ) . From [1, Chapter 9], we have that j , ≈ . { T Y k } k ≥ , where the triangulation of Ω isobtained by quasi-uniform refinement and the partition of [0 , Y k ] is as in § Y k is chosen so that (cid:15) ≈ ( T Y k − ) − / . With these meshes (cid:107) u − V T Y k (cid:107) ◦ H L ( C ,y α ) (cid:46) | log( T Y k ) | s · ( T Y k ) − / , (6.2)which is near-optimal.Figure 6.2 shows the errors of (cid:107) u − V T k, Y (cid:107) H ( y α , C Y k ) for s = 0 . s = 0 .
7. Theresults, again, are in agreement with Theorem 5.4 and Remark 5.5.
7. Fractional powers of general second order elliptic operators.
Let usnow discuss how the methodology developed in previous sections extends to a gen-eral second order, symmetric and uniformly elliptic operator. This is an importantproperty of our PDE approach. Recall that, in § C . In the work of Stinga and Torrea [57],the same type of characterization has been developed for the fractional powers ofsecond order elliptic operators.Let L be a second order symmetric differential operator of the form L w = − div x (cid:48) ( A ∇ x (cid:48) w ) + cw, (7.1)8 R.H. Nochetto, E. Ot´arola, A.J. Salgado −2 −1 Degrees of Freedom (DOFs) E rr o r DOFs − / k e k H ( y α ) −1 Degrees of Freedom (DOFs) E rr o r DOFs − / k e k H ( y α ) Fig. 6.2 . Computational rate of convergence for the approximate solution of the fractionalLaplacian over a circle with graded meshes on the extended dimension. The left panel shows therate for s = 0 . and the right one for s = 0 . . In both cases, the rate is ≈ ( T Y k ) − / in agreementwith Theorem 5.4 and Remark 5.5 where c ∈ L ∞ (Ω) with c ≥ A ∈ C , (Ω , GL ( n, R )) is symmetricand positive definite, and Ω is Lipschitz. Given f ∈ L (Ω), the Lax-Milgram lemmashows that there is a unique w ∈ H (Ω) that solves L w = f in Ω , w = 0 on ∂ Ω . In addition, if Ω has a C , boundary, [38, Theorem 2.4.2.6] shows that w ∈ H (Ω) ∩ H (Ω). Since L − : L (Ω) → L (Ω) is compact and symmetric, its spectrum isdiscrete, positive and accumulates at zero. Moreover, there exists { λ k , ϕ k } k ∈ N ⊂ R + × H (Ω) such that { ϕ k } k ∈ N is an orthonormal basis of L (Ω) and for, k ∈ N , L ϕ k = λ k ϕ k in Ω , ϕ k = 0 on ∂ Ω , (7.2)and λ k → ∞ as k → ∞ . For u ∈ C ∞ (Ω) we then define the fractional powers of L as L s u = ∞ (cid:88) k =1 u k λ sk ϕ k , (7.3)where u k = (cid:82) Ω uϕ k . By density the operator L s can be extended again to H s (Ω). Thisdiscussion shows that it is legitimate to study the following problem: given s ∈ (0 , f ∈ H s (Ω) (cid:48) , find u ∈ H s (Ω) such that L s u = f in Ω . (7.4)To realize the operator L s as the Dirichlet to Neumann map of an extensionproblem we use the generalization of the result by Caffarelli and Silvestre presentedin [57]. We seek a function u : C → R that solves −L u + αy ∂ y u + ∂ yy u = 0 , in C , u = 0 , on ∂ L C ,∂ u ∂ν α = d s f, on Ω × { } , (7.5)where the constant d s is as in (2.23). In complete analogy to § d s L s u = ∂ u ∂ν α : H s (Ω) (cid:55)−→ H s (Ω) (cid:48) . PDE approach to fractional diffusion y α A ∇ u ) + y α c u , where, for all x ∈ C , A ( x ) = diag { A ( x (cid:48) ) , } ∈ GL ( n + 1 , R ).It suffices now to notice that both y α c and y α A are in A ( R n +1+ ), to concludethat, given f ∈ H s (Ω) (cid:48) , there is a unique u ∈ ◦ H L ( C , y α ) that solves (7.5), [34]. Inaddition, u = u ( · , ∈ H s (Ω) solves (7.4) and we have the stability estimate (cid:107) u (cid:107) H s (Ω) (cid:46) (cid:107)∇ u (cid:107) L ( C ,y α ) (cid:46) (cid:107) f (cid:107) H s (Ω) (cid:48) , (7.6)where the hidden constants depend on A , c , C ,y α and Ω.The representation (2.24) of u in terms of the Bessel functions is still valid. Wecan thus repeat the arguments in the proof of Theorem 3.5 to conclude that (cid:107)∇ u (cid:107) L (Ω × ( Y , ∞ ) ,y α ) (cid:46) e −√ λ Y / (cid:107) f (cid:107) H s (Ω) (cid:48) , and introduce v ∈ ◦ H L ( C Y , y α ) — solution of a truncated version of (7.5) — and showthat (cid:107)∇ ( u − v ) (cid:107) L ( C ,y α ) (cid:46) e −√ λ Y / (cid:107) f (cid:107) H s (Ω) (cid:48) . (7.7)Next, we define the finite element approximation of the solution to (7.5) as theunique function V T Y ∈ V ( T Y ) that solves (cid:90) C Y y α A ( x ) ∇ V T Y · ∇ W + y α c ( x (cid:48) ) V T Y W d x (cid:48) d y = d s (cid:104) f, tr Ω W (cid:105) , ∀ W ∈ V ( T Y ) . (7.8)We construct, as in § T Ω of Ω, which we extend to T Y ∈ T with the partition given in (5.6), with γ > / (1 − α ). Following the proof ofTheorem 5.4 we can also show the following error estimate. Theorem 7.1 (
Error estimate for general operators ). Let V T ∈ V ( T Y ) be thesolution of (7.8) and U T Ω ∈ U ( T Ω ) be defined as in (4.5) . If u , solution of (7.5) , issuch that L u , ∂ y ∇ u ∈ H ( y α , C ) , then we have (cid:107) u − U T Ω (cid:107) H s (Ω) (cid:46) (cid:107) u − V T Y (cid:107) ◦ H L ( C ,y α ) (cid:46) | log( T Y ) | s ( T Y ) − / ( n +1) (cid:107) f (cid:107) H − s (Ω) . REFERENCES[1]
M. Abramowitz and I.A. Stegun , Handbook of mathematical functions with formulas, graphs,and mathematical tables , vol. 55 of National Bureau of Standards Applied MathematicsSeries, For sale by the Superintendent of Documents, U.S. Government Printing Office,Washington, D.C., 1964.[2]
G. Acosta , Lagrange and average interpolation over 3D anisotropic elements , J. Comput.Appl. Math., 135 (2001), pp. 91–109.[3]
R.A. Adams , Sobolev spaces , Academic Press [A subsidiary of Harcourt Brace Jovanovich,Publishers], New York-London, 1975. Pure and Applied Mathematics, Vol. 65.[4]
T. Apel , Interpolation of non-smooth functions on anisotropic finite element meshes , M2ANMath. Model. Numer. Anal., 33 (1999), pp. 1149–1185.[5]
O.G. Bakunin , Turbulence and diffusion , Springer Series in Synergetics, Springer-Verlag,Berlin, 2008. Scaling versus equations.[6]
W. Bangerth, R. Hartmann, and G. Kanschat , deal.II—diferential equations analysis li-brary , Technical Reference: http//dealii.org. R.H. Nochetto, E. Ot´arola, A.J. Salgado[7] , deal.II—a general-purpose object-oriented finite element library , ACM Trans. Math.Software, 33 (2007), pp. Art. 24, 27.[8]
P.W. Bates , On some nonlocal evolution equations arising in materials science , in Nonlineardynamics and evolution equations, vol. 48 of Fields Inst. Commun., Amer. Math. Soc.,Providence, RI, 2006, pp. 13–52.[9]
Z. Belhachmi, Ch. Bernardi, and S. Deparis , Weighted Cl´ement operator and application tothe finite element discretization of the axisymmetric Stokes problem , Numer. Math., 105(2006), pp. 217–247.[10]
D.A. Benson, S.W. Wheatcraft, and M.M. Meerschaert , Application of a fractionaladvection-dispersion equation , Water Resources Res., 36 (2000), pp. 91–109.[11]
C. Bernardi , Optimal finite-element interpolation on curved domains , SIAM J. Numer. Anal.,26 (1989), pp. 1212–1240.[12]
J. Bertoin , L´evy processes , vol. 121 of Cambridge Tracts in Mathematics, Cambridge Univer-sity Press, Cambridge, 1996.[13]
K. Bogdan, K. Burdy, and Chen Z.Q. , Censored stable processes , Probab. Theory RelatedFields, 127 (2003), pp. 89–152.[14]
C. Brandle, E. Colorado, A. de Pablo, and U. Sanchez , A concave-convex elliptic probleminvolving the fractional laplacian . in press, 2012.[15]
S.C. Brenner and L.R. Scott , The mathematical theory of finite element methods , vol. 15 ofTexts in Applied Mathematics, Springer, New York, third ed., 2008.[16]
K. Burrage, K. Hale, and D Kay , An efficient implementation of an implicit fem schemefor fractional-in-space reaction-diffusion equations . OCCAM report 1358, 2012.[17]
X. Cabr´e and Y. Sire , Nonlinear equations for fractional Laplacians ii: Existence, uniquenessand qualitative properties of solutions . arXiv:1111.0796v1, 2011.[18]
X. Cabr´e and J. Tan , Positive solutions of nonlinear problems involving the square root ofthe Laplacian , Adv. Math., 224 (2010), pp. 2052–2093.[19]
L. Caffarelli and L. Silvestre , An extension problem related to the fractional Laplacian ,Comm. Partial Differential Equations, 32 (2007), pp. 1245–1260.[20]
A. Capella, J. D´avila, L. Dupaigne, and Y. Sire , Regularity of radial extremal solutionsfor some non-local semilinear equations , Comm. Partial Differential Equations, 36 (2011),pp. 1353–1384.[21]
P. Carr, H. Geman, D.B. Madan, and M. Yor , The fine structure of asset returns: Anempirical investigation , Journal of Business, 75 (2002), pp. 305–332.[22]
P.G. Ciarlet , The finite element method for elliptic problems , vol. 40 of Classics in AppliedMathematics, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA,2002. Reprint of the 1978 original [North-Holland, Amsterdam; MR0520174 (58
P.G. Ciarlet and P.-A. Raviart , Interpolation theory over curved elements, with applicationsto finite element methods , Comput. Methods Appl. Mech. Engrg., 1 (1972), pp. 217–249.[24]
P. Cl´ement , Approximation by finite element functions using local regularization , RAIROAnalyse Num´erique, 9 (1975), pp. 77–84.[25]
J. Cushman and T. Glinn , Nonlocal dispersion in media with continuously evolving scales ofheterogeneity , Trans. Porous Media, 13 (1993), pp. 123–138.[26]
E. Di Nezza, G. Palatucci, and E. Valdinoci , Hitchhiker’s guide to the fractional Sobolevspaces , Bull. Sci. Math., 136 (2012), pp. 521–573.[27]
T. Dupont and R. Scott , Polynomial approximation of functions in sobolev spaces , Math.Comp., 34 (1980), pp. 441–463.[28]
R.G. Dur´an and A.L. Lombardi , Error estimates on anisotropic Q elements for functionsin weighted Sobolev spaces , Math. Comp., 74 (2005), pp. 1679–1706 (electronic).[29] R.G. Dur´an, A.L. Lombardi, and M.I. Prieto , Superconvergence for finite element approxi-mation of a convection–diffusion equation using graded meshes , IMA Journal of NumericalAnalysis, 32 (2012), pp. 511–533.[30]
R.G. Dur´an and F. L´opez Garc´ıa , Solutions of the divergence and Korn inequalities ondomains with an external cusp , Ann. Acad. Sci. Fenn. Math., 35 (2010), pp. 421–438.[31]
G. Duvaut and J.-L. Lions , Inequalities in mechanics and physics , Springer-Verlag, Berlin,1976. Translated from the French by C. W. John, Grundlehren der Mathematischen Wis-senschaften, 219.[32]
A.C. Eringen , Nonlocal continuum field theories , Springer-Verlag, New York, 2002.[33]
L.C. Evans , Partial differential equations , vol. 19 of Graduate Studies in Mathematics, Amer-ican Mathematical Society, Providence, RI, second ed., 2010.[34]
E. B. Fabes, C.E. Kenig, and R.P. Serapioni , The local regularity of solutions of degenerateelliptic equations , Comm. Partial Differential Equations, 7 (1982), pp. 77–116.[35]
D. Gilbarg and N.S. Trudinger , Elliptic partial differential equations of second order , Clas- PDE approach to fractional diffusion sics in Mathematics, Springer-Verlag, Berlin, 2001. Reprint of the 1998 edition.[36] G. Gilboa and S. Osher , Nonlocal operators with applications to image processing , MultiscaleModel. Simul., 7 (2008), pp. 1005–1028.[37]
V. Gol (cid:48) dshtein and A. Ukhlov , Weighted Sobolev spaces and embedding theorems , Trans.Amer. Math. Soc., 361 (2009), pp. 3829–3850.[38]
P. Grisvard , Elliptic problems in nonsmooth domains , vol. 24 of Monographs and Studies inMathematics, Pitman (Advanced Publishing Program), Boston, MA, 1985.[39]
Q.Y. Guan and Z.M. Ma , Reflected symmetric α -stable processes and regional fractional lapla-cian , Probab. Theory Related Fields, 134 (2006), pp. 649–694.[40] J. Heinonen, T. Kilpel¨ainen, and O. Martio , Nonlinear potential theory of degenerate ellip-tic equations , Oxford Mathematical Monographs, The Clarendon Press Oxford UniversityPress, New York, 1993. Oxford Science Publications.[41]
M. Ilic, F. Liu, I. Turner, and V. Anh , Numerical approximation of a fractional-in-spacediffusion equation. I , Fract. Calc. Appl. Anal., 8 (2005), pp. 323–341.[42] ,
Numerical approximation of a fractional-in-space diffusion equation. II. With nonho-mogeneous boundary conditions , Fract. Calc. Appl. Anal., 9 (2006), pp. 333–349.[43]
V.G. Korneev , The construction of variational difference schemes of a high order of accuracy ,Vestnik Leningrad. Univ., 25 (1970), pp. 28–40. (In Russian).[44]
V.G. Korneev and S.E. Ponomarev , Application of curvilinear finite elements in schemesfor solution of n -order linear elliptic equations. I , ˇCisl. Metody Meh. Sploˇsn. Sredy, 5(1974), pp. 78–97. (In Russian).[45] A. Kufner , Weighted Sobolev spaces , A Wiley-Interscience Publication, John Wiley & SonsInc., New York, 1985. Translated from the Czech.[46]
N.S. Landkof , Foundations of modern potential theory , Springer-Verlag, New York, 1972.Translated from the Russian by A. P. Doohovskoy, Die Grundlehren der mathematischenWissenschaften, Band 180.[47]
J.-L. Lions and E. Magenes , Non-homogeneous boundary value problems and applications.Vol. I , Springer-Verlag, New York, 1972. Translated from the French by P. Kenneth, DieGrundlehren der mathematischen Wissenschaften, Band 181.[48]
B.M. McCay and M.N.L. Narasimhan , Theory of nonlocal electromagnetic fluids , Arch. Mech.(Arch. Mech. Stos.), 33 (1981), pp. 365–384.[49]
W. McLean , Strongly elliptic systems and boundary integral equations , Cambridge UniversityPress, Cambridge, 2000.[50]
K.S. Miller and S.G. Samko , Completely monotonic functions , Integral Transform. Spec.Funct., 12 (2001), pp. 389–402.[51]
B. Muckenhoupt , Weighted norm inequalities for the Hardy maximal function , Trans. Amer.Math. Soc., 165 (1972), pp. 207–226.[52]
G. Savar´e , Regularity and perturbation results for mixed second order elliptic problems , Comm.Partial Differential Equations, 22 (1997), pp. 869–899.[53]
L.R. Scott and S. Zhang , Finite element interpolation of nonsmooth functions satisfyingboundary conditions , Math. Comp., 54 (1990), pp. 483–493.[54]
R. Servadei and E. Valdinoci , On the spectrum of two different fractional operators . preprint,2012.[55]
S.A. Silling , Reformulation of elasticity theory for discontinuities and long-range forces , J.Mech. Phys. Solids, 48 (2000), pp. 175–209.[56]
E.M. Stein , Singular integrals and differentiability properties of functions , Princeton Mathe-matical Series, No. 30, Princeton University Press, Princeton, N.J., 1970.[57]
P.R. Stinga and J.L. Torrea , Extension problem and Harnack’s inequality for some frac-tional operators , Comm. Partial Differential Equations, 35 (2010), pp. 2092–2122.[58]
L. Tartar , An introduction to Sobolev spaces and interpolation spaces , vol. 3 of Lecture Notesof the Unione Matematica Italiana, Springer, Berlin, 2007.[59]