Numerical Methods for Fractional Diffusion
Andrea Bonito, Juan Pablo Borthagaray, Ricardo H. Nochetto, Enrique Otarola, Abner J. Salgado
NNUMERICAL METHODS FOR FRACTIONAL DIFFUSION
ANDREA BONITO, JUAN PABLO BORTHAGARAY, RICARDO H. NOCHETTO,ENRIQUE OT ´AROLA, AND ABNER J. SALGADO
Abstract.
We present three schemes for the numerical approximation of frac-tional diffusion, which build on different definitions of such a non-local process.The first method is a PDE approach that applies to the spectral definition andexploits the extension to one higher dimension. The second method is the in-tegral formulation and deals with singular non-integrable kernels. The thirdmethod is a discretization of the Dunford-Taylor formula. We discuss pros andcons of each method, error estimates, and document their performance with afew numerical experiments. Introduction
Diffusion is the tendency of a substance to evenly spread into an available space,and is one of the most common physical processes. The classical models of diffusion,which usually start from the assumption of Brownian motion [1, 57], lead to wellknown models and even better studied equations. However, in recent times, it hasbecome evident that many of the assumptions that lead to these models are notalways satisfactory or even realistic at all. For this reason, new models of diffusionhave been introduced. These, as a rule, are not based on the postulate that theunderlying stochastic process is given by Brownian motion, so that the diffusionis regarded as anomalous [87]. The evidence of anomalous diffusion processes hasbeen reported in physical and social environments, and corresponding models havebeen proposed in various areas such as electromagnetic fluids [83], ground-watersolute transport [47], biology [40], finance [41], human travel [31] and predatorsearch patterns [115].Of the many possible models of anomalous diffusion, we shall be interested herein so-called fractional diffusion, which is a nonlocal process: to evaluate fractionaldiffusion at a spatial point, information involving all spatial points is needed. Re-cently, the analysis of such operators has received a tremendous attention: fractionaldiffusion has been one of the most studied topics in the past decade [37, 38, 114].The main goal of this work is to review different techniques to approximate thesolution of problems involving fractional diffusion. To make matters precise, wewill consider the fractional powers of the Dirichlet Laplace operator ( − ∆) s , whichwe will simply call the fractional Laplacian. Given 0 < s <
1, a bounded Lipschitzdomain Ω ⊂ R d , and a function f : Ω → R , we shall be concerned with finding u : Ω → R such that(1.1) ( − ∆) s u = f in Ω , with vanishing Dirichlet boundary conditions (understood in a suitable sense). Wemust immediately comment that the efficient approximation of solutions to (1.1) Date : Draft version of July 7, 2017.AB is supported in part by NSF grant DMS-1254618.JPB has been partially supported by a CONICET doctoral fellowship.RHN has been supported in part by NSF grant DMS-1411808.EO has been supported in part by CONICYT through project FONDECYT 3160201.AJS is supported by NSF grant DMS-1418784. a r X i v : . [ m a t h . NA ] J u l A. BONITO, J.P. BORTHAGARAY, R.H. NOCHETTO, E. OT´AROLA, AND A.J. SALGADO carries two essential difficulties. The first and most important one is that ( − ∆) s is non-local . The second feature is the lack of boundary regularity, which leads toreduced rates of convergence.In what follows we review different definitions for the fractional Laplacian. Forfunctions defined over R d , there is a natural way to define the fractional Laplacianas a pseudo-differential operator with symbol | ξ | s ; namely, given a function w inthe Schwartz class S , set(1.2) ( − ∆) s w := F − (cid:0) | ξ | s F w (cid:1) , where F denotes the Fourier transform. The fractional Laplacian can be equiva-lently defined by means of the following point-wise formula (see [78, Section 1.1]and [50, Proposition 3.3])(1.3) ( − ∆) s w ( x ) = C ( d, s ) p.v. ˆ R d w ( x ) − w ( x (cid:48) ) | x − x (cid:48) | d +2 s d x (cid:48) , C ( d, s ) = 2 s s Γ( s + d ) π d/ Γ(1 − s ) , where p.v. stands for the Cauchy principal value and C ( d, s ) is a normalizationconstant chosen so that definitions (1.2) and (1.3) coincide. This clearly displaysthe non-local structure of ( − ∆) s w . We remark that, in the theory of stochasticprocesses, expression (1.3) appears as the infinitesimal generator of a 2 s -stableL´evy process [18].If Ω is a bounded domain, we consider two possible definitions of the fractionalLaplacian. For u : Ω → R , we first extend u by zero outside Ω and next usedefinition (1.3). This gives the following reinterpretation of (1.1):(1.4) ( − ∆) s ˜ u = f in Ω , ˜ u = 0 in Ω c = R d \ Ω , where the operator ( − ∆) s is understood as in (1.3) and ˜ w is the extension by zeroto R d of a function w : Ω → R in L (Ω). This definition maintains the probabilisticinterpretation of the fractional Laplacian defined over R d , that is, as the generatorof a random walk in Ω with arbitrarily long jumps, where particles are killed uponreaching Ω c ; see [32, Chapter 2]. The operator in (1.3) is well defined for smooth,compactly supported functions. Consequently, (1.3) can be extended by density to H s (Ω) which, for s ∈ [0 , / H s (Ω) := (cid:8) w | Ω : w ∈ H s ( R d ) , w | R d \ Ω = 0 (cid:9) . When ∂ Ω is Lipschitz this space is equivalent to H s (Ω) = [ L (Ω) , H (Ω)] s , the realinterpolation between L (Ω) and H (Ω) [85, 118] when s ∈ [0 ,
1] and to H s (Ω) ∩ H (Ω) when s ∈ (1 , / H − s (Ω) the dual of H s (Ω) and by (cid:104)· , ·(cid:105) the duality pairing between these two spaces.The second definition of ( − ∆) s relies on spectral theory [19]. Since − ∆ : D ( − ∆) ⊂ L (Ω) → L (Ω) is an unbounded, positive and closed operator withdense domain D ( − ∆) = H (Ω) ∩ H (Ω) and its inverse is compact, there is acountable collection of eigenpairs { λ k , ϕ k } k ∈ N ⊂ R + × H (Ω) such that { ϕ k } k ∈ N isan orthonormal basis of L (Ω) as well as an orthogonal basis of H (Ω). Fractionalpowers of the Dirichlet Laplacian can be thus defined as(1.6) ( − ∆) s w := ∞ (cid:88) k =1 λ sk w k ϕ k , w k := ˆ Ω wϕ k d x, k ∈ N , for any w ∈ C ∞ (Ω). This definition of ( − ∆) s can also be extended by densityto the space H s (Ω). We also remark that, for s ∈ [ − , H s (Ω) canequivalently be defined by H s (Ω) = (cid:40) w = ∞ (cid:88) k =1 w k ϕ k : ∞ (cid:88) k =1 λ sk w k < ∞ (cid:41) . UMERICAL METHODS FOR FRACTIONAL DIFFUSION 3
These two definitions of the fractional Laplacian, the integral one involved inproblem (1.4) and the spectral one given as in (1.6), do not coincide. In fact, asshown in [90], their difference is positive and positivity preserving, see also [35, 113].This, in particular, implies that the boundary behavior of the solutions of (1.4) and(1.6) is quite different. According to Grubb [64] the solution u of (1.4) is of theform(1.7) u ( x ) ≈ dist( x, ∂ Ω) s + v ( x ) , with v smooth; hereafter dist( x, ∂ Ω) indicates the distance from x ∈ Ω to ∂ Ω. Incontrast, Caffarelli and Stinga [35] showed that solutions of (1.6) behave like(1.8) u ( x ) ≈ dist( x, ∂ Ω) s + v ( x ) 0 < s <
12 ; u ( x ) ≈ dist( x, ∂ Ω) + v ( x ) 12 ≤ s < . This lack of boundary regularity is responsible for reduced rates of convergence.The presence of the non-integrable kernel in (1.3) is a notorious numerical diffi-culty that has hampered progress in the multidimensional case d > R d , the Caffarelli-Silvestre [37] extension converts (1.1) into the followingDirichlet-to-Neumann map formulated in the cylinder C = Ω × (0 , ∞ )(1.9) div ( y α ∇ U ) = 0 in C , ∂ ν α U = d s f on R d × { } , where y > d s := 2 − s Γ(1 − s ) / Γ( s ) is a positive nor-malization constant that depends only on s , and the parameter α is defined as α := 1 − s ∈ ( − , u = U ( · , . Cabr´e and Tan [33] and Stinga and Torrea [117] have shown that a similar extensionis valid for the spectral fractional Laplacian in a bounded domain Ω provided thata vanishing Dirichlet condition is appended on the lateral boundary ∂ L C = ∂ Ω × (0 , ∞ ); see also [39, 29]. Although (1.9) is a local problem, and thus amenable toPDE techniques, it is formulated in one higher dimension and exhibits a singularcharacter as y ↓ s ∈ (0 ,
1) and w ∈ H − s (Ω) we have(1.10) ( − ∆) − s w = 12 πi ˆ C z − s ( z + ∆) − w d z, where C is a Jordan curve oriented to have the spectrum of − ∆ to its right, and z − s is defined using the principal value of Log( z ). In addition, since the operator − ∆ is positive, one can continously deform the contour C onto the negative realaxis (around the branch cut) to obtain the so-called Balakrishnan formula(1.11) ( − ∆) − s w = sin( sπ ) π ˆ ∞ µ − s ( µ − ∆) − w d µ ;see also [122, Section IX.11] and [19, Section 10.4] for a different derivation of(1.11) using semigroup theory. This formula will be the starting point for ourDunford-Taylor approach for the spectral fractional Laplacian. These considera-tions, however, cannot be carried out for the the integral fractional Laplacian (1.3)since neither (1.10) or (1.11) are well defined quantities for this operator. There-fore we will, instead, multiply (1.4) by a test function w , integrate over R d and A. BONITO, J.P. BORTHAGARAY, R.H. NOCHETTO, E. OT´AROLA, AND A.J. SALGADO use Parseval’s equality to obtain the following weak formulation of (1.4): Given f ∈ H − s (Ω) find u ∈ H s (Ω) such that(1.12) ˆ R d | ξ | s F − (˜ u ) | ξ | s F − ( ˜ w ) d ξ = (cid:104) f, w (cid:105) , ∀ w ∈ H s (Ω) , where z denotes the complex conjugate of z ∈ C . Using again the Fourier transformand Parseval’s equality, the left hand side of the above relation can be equivalentlywritten as (see Theorem 4.5)(1.13) 2 sin( sπ ) π ˆ ∞ µ − s ˆ Ω (cid:0) ( − ∆)( I − µ ∆) − ˜ u ( x ) (cid:1) w ( x ) d x d µ. These ideas will be the starting point of the Dunford-Taylor method for the integralfractional Laplacian (1.3).The purpose of this paper is to present and briefly analyze three finite elementmethods (FEMs) to approximate (1.1). The first will consider the spectral definition(1.6), the second will deal with the integral definition (1.3), while the third approachwill be able to account for both operators by means of either (1.11) or (1.13). Wemust immediately remark that for special domain geometries, such as when d = 2and Ω is a rectangle, the use of spectral methods can be quite efficient but we donot elaborate any further as we are interested in techniques that apply to generaldomains. Our presentation is organized as follows: In Section 2, we present amethod that hinges on the extension (1.9) and uses PDE techniques. The secondmethod deals with the integral formulation and is discussed in Section 3. Finally,the third method is based on exponentially convergent quadrature approximationsof (1.11) for the spectral Laplacian and of (1.13) for the integral Laplacian. Theyare discussed in Section 4.As usual, we write a (cid:46) b to mean a ≤ Cb , with a constant C that neither de-pends on a, b or the discretization parameters and might change at each occurrence.Moreover, a ≈ b indicates a (cid:46) b and b (cid:46) a .2. The Spectral Fractional Laplacian
In this section we deal with the spectral definition of the fractional Laplacian(1.6) and, on the basis of (1.9), its discretization via PDE techniques as originallydeveloped in [93]. We must immediately remark that many of the results of thissection and section 4 extend to more general symmetric elliptic operators of theform Lw = − div( A ∇ w ) + cw , with A ∈ C , ( ¯Ω , GL ( R d )) symmetric and positivedefinite and 0 ≤ c ∈ C , ( ¯Ω , R ).Let Ω be a convex polytopal domain. Besides the semi-infinite cylinder C =Ω × (0 , ∞ ), we introduce the truncated cylinder C Y = Ω × (0 , Y ) with height Y and its lateral boundary ∂ L C Y = ∂ Ω × (0 , Y ). Since we deal with objects definedin both R d and R d +1 , it is convenient to distinguish the extended variable y . For x ∈ R d +1 , we denote x = ( x, y ) = ( x, x d +1 ) , x ∈ R d , y ∈ R + . Extension Property.
The groundbreaking extension (1.9) of Caffarelli andSilvestre [37], valid for any power s ∈ (0 , R d . Cabr´e and Tan [33]and Stinga and Torrea [117] realized that a similar extension holds for the spectralLaplacian over Ω bounded; see also [39, 29]. This problem reads(2.1) div ( y α ∇ U ) = 0 in C , U = 0 on ∂ L C , ∂ ν α U = d s f on Ω × { } , where α = 1 − s and the so-called conormal exterior derivative of U at Ω × { } is(2.2) ∂ ν α U = − lim y → + y α ∂ y U ( · , y ) . UMERICAL METHODS FOR FRACTIONAL DIFFUSION 5
The limit in (2.2) must be understood in the distributional sense [33, 37, 39, 117].With this construction at hand, the fractional Laplacian and the Dirichlet-to-Neumann operator of problem (2.1) are related by d s ( − ∆) s u = ∂ ν α U in Ω . The operator in (2.1) is in divergence form and thus amenable to variationaltechniques. However, it is nonuniformly elliptic because the weight y α either blowsup for − < α < < α < y ↓
0; the exceptional case α = 0corresponds to the regular harmonic extension for s = [33]. This entails dealingwith weighted Lebesgue and Sobolev spaces with the weight | y | α for α ∈ ( − , A ( R d +1 ), whichis the collection of weights ω so that [52, 58, 61, 89, 120] C ,ω = sup B (cid:18) B ω d x (cid:19) (cid:18) B ω − d x (cid:19) < ∞ , where the supremum is taken over all balls B in R d +1 and ffl B stands for the meanvalue over B . The Muckenhoupt characteristic C ,ω appears in all estimates involv-ing ω .If D ⊂ R d × R + , we define L ( y α , D ) as the Lebesgue space for the measure y α d x . We also define the weighted Sobolev space H ( y α , D ) = (cid:8) w ∈ L ( y α , D ) : |∇ w | ∈ L ( y α , D ) (cid:9) , where ∇ w is the distributional gradient of w . We equip H ( y α , D ) with the norm(2.3) (cid:107) w (cid:107) H ( y α ,D ) = (cid:16) (cid:107) w (cid:107) L ( y α ,D ) + (cid:107)∇ w (cid:107) L ( y α ,D ) (cid:17) . The space H ( y α , D ) is Hilbert with the norm (2.3) and C ∞ ( D ) ∩ H ( y α , D ) isdense in H ( y α , D ) because | y | α ∈ A ( R d +1 ) (cf. [61, Theorem 1], [76] and [120,Proposition 2.1.2, Corollary 2.1.6]).To analyze problem (2.1) we define the weighted Sobolev space ◦ H L ( y α , C ) = (cid:8) w ∈ H ( y α , C ) : w = 0 on ∂ L C (cid:9) . The following weighted Poincar´e inequality holds [93, inequality (2.21)] (cid:107) w (cid:107) L ( y α , C ) (cid:46) (cid:107)∇ w (cid:107) L ( y α , C ) ∀ w ∈ ◦ H L ( y α , C ) . Consequently, the seminorm on ◦ H L ( y α , C ) is equivalent to (2.3). For w ∈ H ( y α , C ),tr Ω w denotes its trace onto Ω × { } which satisfies [93, Proposition 2.5]tr Ω ◦ H L ( y α , C ) = H s (Ω) , (cid:107) tr Ω w (cid:107) H s (Ω) ≤ C tr Ω (cid:107) w (cid:107) ◦ H L ( y α , C ) . The variational formulation of (2.1) reads: find U ∈ ◦ H L ( y α , C ) such that(2.4) ˆ C y α ∇ U · ∇ w d x = d s (cid:104) f, tr Ω w (cid:105) ∀ w ∈ ◦ H L ( y α , C ) , where, we recall that, (cid:104)· , ·(cid:105) corresponds to the duality pairing between H − s (Ω) and H s (Ω). The fundamental result of Caffarelli and Silvestre [37] for Ω = R d and ofCabr´e and Tan [33, Proposition 2.2] and Stinga and Torrea [117, Theorem 1.1] forΩ bounded reads: given f ∈ H − s (Ω), if u ∈ H s (Ω) solves (1.1) and U ∈ ◦ H L ( y α , C )solves (2.4), then(2.5) u = tr Ω U , d s ( − ∆) s u = ∂ ν α U in Ω , where the first equality holds in H s (Ω), whereas the second one in H − s (Ω). A. BONITO, J.P. BORTHAGARAY, R.H. NOCHETTO, E. OT´AROLA, AND A.J. SALGADO
Regularity.
To study the finite element discretization of (2.4) we must un-derstand the regularity of U . We begin by recalling that if u = (cid:80) ∞ k =1 u k ϕ k solves(1.1), then u k = λ − sk f k , with f k being the k -th Fourier coefficient of f . The uniquesolution U of problem (1.9) thus admits the representation [93, formula (2.24)] U ( x, y ) = ∞ (cid:88) k =1 u k ϕ k ( x ) ψ k ( y ) . The functions ψ k solve the 2-point boundary value problem in R + ψ (cid:48)(cid:48) k + αy ψ (cid:48) k = λ k ψ k , in (0 , ∞ ); ψ k (0) = 1 , lim y →∞ ψ k ( y ) = 0 . Thus, if s = , we have that ψ k ( y ) = exp( −√ λ k y ) [33, Lemma 2.10]; and, if s ∈ (0 , \ { } , then [39, Proposition 2.1] ψ k ( y ) = c s ( (cid:112) λ k y ) s K s ( (cid:112) λ k y ) , where c s = 2 − s / Γ( s ) and K s denotes the modified Bessel function of the secondkind [2, Chap. 9.6]. Using asymptotic properties of K s ( ζ ) as ζ ↓ y ↓ ψ (cid:48) k ( y ) ≈ y − α , ψ (cid:48)(cid:48) k ( y ) ≈ y − α − . Exploiting these estimates, the following regularity results hold [93, Theorem 2.7].
Theorem 2.1 (global regularity of the α -harmonic extension) . Let f ∈ H − s (Ω) ,where H − s (Ω) is defined in (1.5) for s ∈ (0 , . Let U ∈ ◦ H L ( y α , C ) solve (1.9) with f as data. Then, for s ∈ (0 , \ { } , we have (2.6) (cid:107) ∆ x U (cid:107) L ( y α , C ) + (cid:107) ∂ y ∇ x U (cid:107) L ( y α , C ) = d s (cid:107) f (cid:107) H − s (Ω) , and (2.7) (cid:107) ∂ yy U (cid:107) L ( y β , C ) (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 / (Ω) . Comparing (2.6) and (2.7), we realize that the regularity of U is much worsein the extended d + 1 direction. Since Ω is convex, the following elliptic regularityestimate is valid [63]:(2.8) (cid:107) w (cid:107) H (Ω) (cid:46) (cid:107) ∆ x w (cid:107) L (Ω) ∀ w ∈ H (Ω) ∩ H (Ω) . This, combined with (2.6), yields the following estimate for the Hessian D x U inthe variable x ∈ Ω: (cid:107) D x U (cid:107) L ( y α , C ) (cid:46) (cid:107) f (cid:107) H − s (Ω) . Further regularity estimates in H¨older and Sobolev norms are derived in [35]. How-ever, we do not need them for what follows.2.3.
Truncation.
Since C is unbounded, problem (2.1) cannot be approximatedwith standard finite element techniques. However, since the solution U of problem(2.1) decays exponentially in y [93, Proposition 3.1], by truncating C to C Y :=Ω × (0 , Y ) and setting a homogeneous Dirichlet condition on y = Y , we only incurin an exponentially small error in terms of Y [93, Theorem 3.5]. If ◦ H L ( y α , C Y ) := (cid:8) w ∈ H ( y α , C Y ) : w = 0 on Γ D (cid:9) , UMERICAL METHODS FOR FRACTIONAL DIFFUSION 7 where Γ D = ∂ L C Y ∪ Ω × { Y } is the Dirichlet boundary, then the aforementionedproblem reads:(2.9) V ∈ ◦ H L ( y α , C Y ) : ˆ C Y y α ∇ V · ∇ w d x = d s (cid:104) f, tr Ω w (cid:105) ∀ w ∈ ◦ H L ( y α , C Y ) . The following exponential decay rate is proved in [93, Theorem 3.5].
Lemma 2.2 (truncation error) . If U and V denote the solutions of (2.4) and (2.9) , respectively, then (cid:107)∇ ( U − V ) (cid:107) L ( y α , C ) (cid:46) e −√ λ Y / (cid:107) f (cid:107) H − s (Ω) is valid, where λ denotes the first eigenvalue of the Dirichlet Laplace operator and Y is the truncation parameter. FEM: A Priori Error Analysis.
The first numerical work that exploits thegroundbreaking identity (2.5), designs and analyzes a FEM for (2.1) is [93]; see also[94, 101]. We briefly review now the main a priori results of [93].We introduce a conforming and shape regular mesh T = { K } of Ω, where K ⊂ R d is an element that is isoparametrically equivalent either to the unit cube[0 , d or the unit simplex in R d . Over this mesh we construct the finite elementspace(2.10) U ( T ) = (cid:8) W ∈ C ( ¯Ω) : W | K ∈ P ( K ) ∀ K ∈ T , W | ∂ Ω = 0 (cid:9) . The set P ( K ) is either the space P ( K ) of polynomials of total degree at most 1,when K is a simplex, or the space of polynomials Q ( K ) of degree not larger than1 in each variable provided K is a d -rectangle.To triangulate the truncated cylinder C Y we consider a partition I Y of the in-terval [0 , Y ] with nodes y m , m = 0 , · · · , M , and construct a mesh T Y of C Y as thetensor product of T and I Y . The set of all triangulations of C Y obtained with thisprocedure is T .For T Y ∈ T , we define the finite element space(2.11) V ( T Y ) = (cid:8) W ∈ C ( ¯ C Y ) : W | T ∈ P ( K ) ⊗ P ( I ) ∀ T ∈ T Y , W | Γ D = 0 (cid:9) , where T = K × I and observe that U ( T ) = tr Ω V ( T Y ). Note that T Y = M T ,and that T ≈ M d implies T Y ≈ M d +1 .The Galerkin approximation of (2.9) is the function V ∈ V ( T Y ) such that(2.12) ˆ C Y y α ∇ V · ∇ W d x = d s (cid:104) f, tr Ω W (cid:105) ∀ W ∈ V ( T Y ) . Existence and uniqueness of V immediately follows from V ( T Y ) ⊂ ◦ H L ( y α , C Y ) andthe Lax-Milgram Lemma. It is trivial also to obtain a best approximation result `ala C´ea, namely(2.13) (cid:107)∇ ( V − V ) (cid:107) L ( y α , C Y ) = inf W ∈ V ( T Y ) (cid:107)∇ ( V − W ) (cid:107) L ( y α , C Y ) . This reduces the numerical analysis of (2.12) to a question in approximation theory,which in turn can be answered by the study of piecewise polynomial interpolationin Muckenhoupt weighted Sobolev spaces; see [93, 96]. Exploiting the Cartesianstructure of the mesh T Y we are able to extend the anisotropic estimates of Dur´anand Lombardi [53] to our setting [93, Theorems 4.6–4.8], [96]. The following errorestimates separate in each direction. A. BONITO, J.P. BORTHAGARAY, R.H. NOCHETTO, E. OT´AROLA, AND A.J. SALGADO
Proposition 2.3 (anisotropic interpolation estimates) . There exists a quasi inter-polation operator Π T Y : L ( C Y ) → V ( T Y ) that satisfies the following anisotropicerror estimates for all j = 1 , . . . , d + 1 and all T = K × I ∈ T Y (cid:107) w − Π T Y w (cid:107) L ( y α ,T ) (cid:46) h K (cid:107)∇ x w (cid:107) L ( y α ,S T ) + h I (cid:107) ∂ y w (cid:107) L ( y α ,S T ) , (cid:107) ∂ x j ( w − Π T Y w ) (cid:107) L ( y α ,T ) (cid:46) h K (cid:107)∇ x ∂ x j w (cid:107) L ( y α ,S T ) + h I (cid:107) ∂ y ∂ x j w (cid:107) L ( y α ,S T ) , where S T stands for the patch of elements of T Y that intersect T , h K = | K | /d and h I = | I | . As a first application of Proposition 2.3 we consider a quasiuniform mesh T Y ofsize h and set w = U . We estimate U − Π T Y U for y ≥ h as follows: ˆ Y h y α (cid:107) ∂ y ( U − Π T Y U ) (cid:107) L (Ω) d y (cid:46) h ˆ Y h y α (cid:16) (cid:107) ∂ yy U (cid:107) L (Ω) + (cid:107)∇ x ∂ y U (cid:107) L (Ω) (cid:17) d y. For the first term we resort to (2.7), and recall that β > α + 1 > α , to deduce h ˆ Y h y α (cid:107) ∂ yy U (cid:107) L (Ω) d y ≤ h sup h ≤ y ≤ Y y α − β ˆ Y y β (cid:107) ∂ yy U (cid:107) L (Ω) d y ≤ h α − β (cid:107) f (cid:107) L (Ω) . For the second term we use (2.6) instead to arrive at h ˆ Y h y α (cid:107)∇ x ∂ y U (cid:107) L (Ω) d y ≤ h ˆ Y y α (cid:107)∇ x ∂ y U (cid:107) L (Ω) d y (cid:46) h (cid:107) f (cid:107) H − s (Ω) . Combining the two estimates we derive the interpolation estimate ˆ Y h y α (cid:107) ∂ y ( U − Π T Y U ) (cid:107) L (Ω) d y (cid:46) h α − β (cid:107) f (cid:107) H − s (Ω) , which is quasi-optimal in terms of regularity because 2 + α − β = 2( s − (cid:15) ) and ∂ y U ≈ y − α formally implies U ∈ H s − (cid:15) ( y α , C ) for any (cid:15) >
0; however this estimateexhibits a suboptimal rate. To restore a quasi-optimal rate, we must compensatethe behavior of ∂ yy U by a graded mesh in the extended direction, which is allowedby Proposition 2.3. Therefore, we construct a mesh I Y with nodes(2.14) y m = m γ M − γ Y , m = 0 , . . . , M, where γ > / (1 − α ) = 3 / (2 s ). Combining (2.13) with Proposition 2.3 we obtainestimates in terms of degrees of freedom [93, Theorem 5.4 and Corollary 7.11]. Theorem 2.4 (a priori error estimate) . Let T Y = T × I Y ∈ T with I Y satisfying (2.14) , and let V ( T Y ) be defined by (2.11) . If V ∈ V ( T Y ) solves (2.12) , we have (cid:107)∇ ( U − V ) (cid:107) L ( y α , C ) (cid:46) | log( T Y ) | s ( T Y ) − / ( d +1) (cid:107) f (cid:107) H − s (Ω) , where Y ≈ log( T Y ) . Alternatively, if u denotes the solution of (1.1) , then (cid:107) u − U (cid:107) H s (Ω) (cid:46) | log( T Y ) | s ( T Y ) − / ( d +1) (cid:107) f (cid:107) H − s (Ω) , where U = tr Ω V .Remark . The estimates of Theorem 2.4 hold onlyif f ∈ H − s (Ω) and the domain Ω is such that (2.8) is valid. Remark . Let V ∈ V ( T Y ) solve (2.12) over a quasi-uniform mesh T Y of C Y of size h . Combining (2.13) with the preceding discussionand accounting for the missing domain Ω × (0 , h ), which yields a better estimate[93], we obtain for Y ≈ | log h | and all ε > (cid:107)∇ ( U − V ) (cid:107) L ( y α , C Y ) (cid:46) h s − ε (cid:107) f (cid:107) H − s (Ω) , where the hidden constant blows up if ε ↓ UMERICAL METHODS FOR FRACTIONAL DIFFUSION 9
Remark . Except for a logaritmic factor, the error estimates ofTheorem 2.4 decay with a rate ( T Y ) − / ( d +1) , where d is the dimension of Ω. TheFEM (2.12) is thus sub-optimal as a method to compute in Ω. This can be improvedto an error decay ( T Y ) − /d with geometric grading; see Section 2.7. Remark s = ) . If s = , we obtain the optimal estimate (cid:107)∇ ( U − V ) (cid:107) L ( C Y ) (cid:46) h (cid:107) f (cid:107) H / (Ω) . Numerical Experiments.
We present two numerical examples for d = 2computed within the deal.II library [15, 16] using graded meshes. Integrals areevaluated with Gaussian quadratures of sufficiently high order and linear systemsare solved using CG with ILU preconditioner and the exit criterion being that the (cid:96) -norm of the residual is less than 10 − .2.5.1. Square Domain.
Let Ω = (0 , . Then ϕ 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 ), then u ( x , x ) = sin( πx ) sin( πx ) , by (1.6)and U ( x , x , y ) = 2 − s/ π s Γ( s ) − sin( πx ) sin( πx ) y s K s ( √ πy ) [93, (2.24)].We construct a sequence of meshes { T Y k } k ≥ , where T is obtained by uniformrefinement and I Y k is given by (2.14) with parameter γ > / (1 − α ). On the basisof Theorem 2.4, the truncation parameter Y k is chosen to be Y k ≥ √ λ (log C − log( T Y k − ) − / ) . With this type of meshes, (cid:107) u − tr Ω V k (cid:107) H s (Ω) (cid:46) (cid:107) U − V k (cid:107) ◦ H L ( y α , C ) (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 [30]) (cid:107) u − tr Ω V k (cid:107) H s (Ω) (cid:46) h − s (cid:46) ( T Y k ) − (2 − s ) / . Figure 2.5.1 shows the rates of convergence for s = 0 . s = 0 . −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 α ) Figure 1.
Computational rate of convergence for a square andgraded meshes. The left panel shows the rate for s = 0 . s = 0 .
8. The experimental rate of convergence is( T Y ) − / , in agreement with Theorem 2.4. −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 α ) Figure 2.
Computational rate of convergence for a circle withgraded meshes. The left panel shows the rate for s = 0 . s = 0 .
7. The experimental rate of convergence is( T Y ) − / , in agreement with Theorem 2.4.2.5.2. Circular Domain.
If Ω = { x ∈ R : | x | < } , then ϕ m,n ( r, θ ) = J m ( j m,n r ) ( A m,n cos( mθ ) + B m,n sin( mθ )) , λ m,n = j m,n , 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 normalization constants to ensure (cid:107) ϕ m,n (cid:107) L (Ω) = 1.If f = ( λ , ) s ϕ , , then (1.6) and [93, (2.24)] show that u = ϕ , and U ( r, θ, y ) = 2 − s Γ( s ) − ( λ , ) s/ ϕ , ( r, θ ) y s K s ( √ πy ) . We construct a sequence of meshes { T Y k } k ≥ as in § (cid:107) U − V k (cid:107) ◦ H L ( y α , C ) (cid:46) | log( T Y k ) | s ( T Y k ) − / , which is near-optimal. Figure 2 shows the errors of (cid:107) U − V k (cid:107) H ( y α , C Y k ) for s = 0 . s = 0 .
7. The results, again, are in agreement with Theorem 2.4.2.6.
FEM: A Posteriori Error Analysis.
A posteriori error estimation andadaptive finite element methods (AFEMs) have been the subject of intense researchsince the late 1970’s because they yield optimal performance in situations whereclassical FEM cannot. The a priori theory for (2.12) requires f ∈ H − s (Ω) and Ωconvex for (2.8) to be valid; see Remark 2.5. If either of these does not hold, then U may have singularities in the x -variables and Theorem 2.4 may not apply: aquasi-uniform refinement of Ω would not result in an efficient solution technique.An adaptive loop driven by an a posteriori error estimator is essential to recoveroptimal rates of convergence. In what follows we explore this.We first observe that we cannot rely on residual error estimators. In fact, theyhinge on the strong form of the local residual to measure the error. Let T ∈ T Y and ν be the unit outer normal to T . Integration by parts yields ˆ T y α ∇ V · ∇ W d x = ˆ ∂T y α W ∇ V · ν d σ − ˆ T div( y α ∇ V ) W d x . Since α ∈ ( − ,
1) the boundary integral is meaningless for y = 0. Even the veryfirst step in the derivation of a residual a posteriori error estimator fails! There isnothing left to do but to consider a different type of estimator.2.6.1. Local Problems over Cylindrical Stars.
Inspired by [12, 88], we deal with theanisotropy of the mesh in the extended variable y and the coefficient y α by con-sidering local problems on cylindrical stars. The solutions of these local problemsallow us to define an anisotropic a posteriori error estimator which, under certainassumptions, is equivalent to the error up to data oscillation terms. UMERICAL METHODS FOR FRACTIONAL DIFFUSION 11
Given a node v on the mesh T Y , we exploit the tensor product structure of T Y , and we write v = ( v , w ) where v and w are nodes on the meshes T and I Y respectively. For K ∈ T , we denote by N ( K ) and ◦ N ( K ) the set of nodes andinterior nodes of K , respectively. We set N ( T ) = (cid:91) K ∈ T N ( K ) , ◦ N ( T ) = (cid:91) K ∈ T ◦ N ( K ) . The star around v is S v = (cid:83) K (cid:51) v K ⊂ Ω , and the cylindrical star around v is C v := (cid:91) { T ∈ T Y : T = K × I, K (cid:51) v } = S v × (0 , Y ) ⊂ C Y . For each node v ∈ N ( T ) we define the local space W ( C v ) = (cid:8) w ∈ H ( y α , C v ) : w = 0 on ∂ C v \ Ω × { } (cid:9) , and the (ideal) estimator η v ∈ W ( C v ) to be the solution of ˆ C v y α ∇ η v · ∇ w d x = d s (cid:104) f, tr Ω w (cid:105) − ˆ C v y α ∇ V · ∇ w d x ∀ w ∈ W ( C v ) . We finally define the local indicators E v and global error estimators E T Y as follows:(2.15) E v = (cid:107)∇ η v (cid:107) L ( y α , C v ) , E T Y = (cid:88) v ∈ N ( T ) E v . We have the following key properties [42, Proposition 5.14].
Proposition 2.9 (a posteriori error estimates) . Let V ∈ ◦ H L ( y α , C Y ) and V ∈ V ( T Y ) solve (2.9) and (2.12) respectively. Then, the estimator defined in (2.15) satisfies the global bound (cid:107)∇ ( V − V ) (cid:107) L ( y α , C Y ) (cid:46) E T Y , and the local bound with constant for any v ∈ N ( T Ω ) E v ≤ (cid:107)∇ ( V − V ) (cid:107) L ( y α , C v ) . These estimates provide the best scenario for a posteriori error analysis but theyare not practical because the local space W ( C v ) is infinite dimensional. We furtherdiscretize W ( C v ) with continuous piecewise polynomials of degree > K is a quadrilateral, we use polynomials of degree ≤ K isa simplex, we employ polynomials of total degree ≤ W ( C v ) and corresponding discrete estimators instead of (2.15). Under suitableassumptions, Proposition 2.9 extends to this case [42, Section 5.4].2.6.2. Numerical Experiment.
We illustrate the performance of a practical versionof the a posteriori error estimator (2.15). We use an adaptive loop
SOLVE → ESTIMATE → MARK → REFINE with D¨orfler marking. We generate a new mesh T (cid:48) by bisecting all the elements K ∈ T contained in the marked set M based on newest-vertex bisection method;[97, 98]. We choose the truncation parameter as Y = 1 + log( T (cid:48) ) [93, Remark5.5]. We set M ≈ ( T (cid:48) ) /d and construct I (cid:48) Y by the rule (2.14). The new mesh T (cid:48) Y = REFINE ( M ) is obtained as the tensor product of T (cid:48) and I (cid:48) Y .We consider the data Ω = (0 , and f ≡
1, which is incompatible for s < / f does not have a vanishing trace whence f / ∈ H − s (Ω); see also [93, Section6.3]. Therefore, we cannot expect an optimal rate for quasi-uniform meshes T inΩ according to Theorem 2.4. Figure 3 shows that adaptive mesh refinement guidedby AFEM restores an optimal decay rate for s < / −1 Degrees of Freedom (DOFs) E rr o r DOFs − s = 0 . s = 0 . s = 0 . s = 0 . −1 Degrees of Freedom (DOFs) E s t i m a t o r DOFs − s = 0 . s = 0 . s = 0 . s = 0 . Figure 3.
Computational rate of convergence of AFEM for d = 2and s = 0 .
2, 0 .
4, 0 . s = 0 .
8. The left panel shows the errorvs. the number of degrees of freedom, the right one the total errorindicator. We recover the optimal rate T Y ) − / . For s < , theright hand side f ≡ / ∈ H − s (Ω) and a quasiuniform mesh in Ωdoes not deliver an optimal rate of convergence [93, Section 6.3].2.7. Extensions and Applications.
We conclude the discussion of this approachby mentioning several extension and applications: • Efficient solvers : Finding the solution to (2.9) entails solving a large linear systemwith a sparse matrix. In [43] the construction of efficient multilevel techniquesfor the solution of this problem was addressed. It was shown that multilevel tech-niques with line smoothers over vertical lines in the extended direction performalmost optimally; i.e., the contraction factor depends linearly on the number oflevels, and thus logarithmically on the problem size. • Time dependent problems : In [95] the time dependent problem(2.16) ∂ γt u + ( − ∆) s u = f Ω × (0 , T ) , u ( · ,
0) = u was examined. Here ∂ γt denotes the so-called Caputo derivative of order γ ∈ (0 , ∂ γt w ( x, t ) := 1Γ(1 − γ ) ˆ t t − r ) γ ∂ r w ( x, r ) d r t > , x ∈ Ω . It turns out that the solution of this problem is always singular as t ↓ u (cid:54) = 0. In fact, taking f = 0 and representing u in terms ofthe Mittag-Leffler function [62] reveals that u ( x, t ) = (cid:18) − t γ Γ(1 + γ ) ( − ∆) s + O ( t γ ) (cid:19) u ( x ) x ∈ Ω . These heuristics led to the following regularity results shown in [95](2.18) ∂ t u ∈ L log L (0 , T ; H − s (Ω)) , ∂ tt u ∈ L ( t σ ; 0 , T ; H − s (Ω)) , σ > − γ, where L log L (0 , T ) denotes the Orlicz space of functions w such that | w | log | w | ∈ L (0 , T ); see [74]. Using the extension property, problem (2.16) reduces to aquasi-stationary elliptic equation with a dynamic boundary condition. Rates ofconvergence for fully discrete schemes were derived in [95], which are consistentwith the regularity (2.18); this issue has been largely ignored in the literature.The extension of these results to a space-time fractional wave equation, i.e., γ ∈ (1 ,
2] is currently under investigation [104]. • Nonlinear problems : Elliptic and parabolic obstacle problems with the spectralfractional Laplacian were considered in [92] and [103], respectively. Rates of con-vergence for their FEM approximation were derived, and this required a carefulcombination of Sobolev regularity, as in Theorem 2.1, and H¨older regularity of
UMERICAL METHODS FOR FRACTIONAL DIFFUSION 13 the solution [36, 34]. In addition, a positivity preserving interpolant that is stablein anisotropic meshes and weighted norms had to be constructed. • PDE constrained optimization and optimal control : Optimal control problemswhere the state equation is given by either a stationary or parabolic equationwith a spectral fractional Laplacian were studied in [8] and [10], respectively.Existence and uniqueness of optimal pairs was obtained, as well as their regu-larity. In both cases fully discrete schemes were designed and their convergenceshown. The results of [8] were later improved upon considering piecewise lin-ear approximation of the optimal control variable [102] and adaptive algorithms[9]. Sparse optimal control for the spectral fractional Laplacian was studied in[105]. Finally, reference [11] presents the design and analysis of an approximationscheme for an optimal control problem where the control variable corresponds tothe order of the fractional operator [116]. • Near optimal complexity : According to Remark 2.7, the FEM on anisotropicmeshes with radical grading (2.14) is suboptimal. This deficiency can be curedby exponential grading in the extended variable y and suitable exploitation ofanalyticity properties of U ( · , y ) in y . This leads to a hybrid FEM that combinesthe hp -version in y with the h -version in Ω and exhibits an error decay ( T Y ) − /d provided f ∈ H − s (Ω) [17, 86]. Dealing with incompatible f / ∈ H − s (Ω) is openalthough this is responsible for boundary singularities governed by (1.8).3. The Integral Fractional Laplacian
Here we consider the discretization of the integral definition of the fractionalLaplacian (1.4). In view of the definition (1.5) of the space H s (Ω), and the fractionalPoincar´e inequality (cid:107) w (cid:107) L (Ω) (cid:46) | w | H s ( R d ) ∀ w ∈ H s (Ω) , we may furnish H s (Ω) with the H s ( R d )-seminorm. We also define the bilinear form (cid:74) · , · (cid:75) : H s (Ω) × H s (Ω) → R ,(3.1) (cid:74) u, w (cid:75) := C ( d, s )2 ¨ Q ( u ( x ) − u ( x (cid:48) ))( w ( x ) − w ( x (cid:48) )) | x − x (cid:48) | d +2 s d x (cid:48) d x, where Q = (Ω × R d ) ∪ ( R d × Ω) and C ( d, s ) was defined in (1.3). We denote by (cid:57) · (cid:57) the norm that (cid:74) · , · (cid:75) induces, which is just a multiple of the H s ( R d )-seminorm.The weak formulation of (1.4) is obtained upon multiplying (1.3) by a test function w ∈ H s (Ω), integrating over x ∈ Ω and exploiting symmetry to make the difference w ( x ) − w ( x (cid:48) ) appear. With the functional setting we have just described at hand,this problem is formulated as follows: find u ∈ H s (Ω) such that(3.2) (cid:74) u, w (cid:75) = (cid:104) f, w (cid:105) ∀ v ∈ H s (Ω) . Applying the Lax-Milgram lemma immediately yields well-posedness of (3.2).Although the energy norm (cid:57) · (cid:57) involves integration on Ω × R d , this norm canbe localized. In fact, due to Hardy’s inequality [44, 54], the following equivalenceshold [5, Corollary 2.6]: (cid:107) w (cid:107) H s (Ω) ≤ (cid:57) w (cid:57) (cid:46) (cid:107) w (cid:107) H s (Ω) , if s ∈ (0 , / , | w | H s (Ω) ≤ (cid:57) w (cid:57) (cid:46) | w | H s (Ω) , if s ∈ (1 / , . When s = 1 /
2, since Hardy’s inequality fails, it is not possible to bound the H ( R d )-seminorm in terms of the H (Ω)-norm for functions supported in Ω. How-ever, for the purposes we pursue in this work, it suffices to notice that the estimate (cid:57) w (cid:57) (cid:46) | w | H
12 + ε (Ω) holds for all w ∈ H + ε (Ω). From this discussion, it follows that the energy norm may be bounded in termsof fractional–order norms on Ω. Thus, in order to estimate errors in the energynorm, we may bound errors within Ω.3.1.
Regularity.
We now review some results regarding Sobolev regularity of so-lutions to problem (3.2) that are useful to deduce convergence rates of the finiteelement scheme proposed below. Regularity results for the fractional Laplacianhave been recently obtained by Grubb [64] in terms of H¨ormander µ -spaces [66].The work [5] has reinterpreted these in terms of standard Sobolev spaces. The fol-lowing result, see [27, 64], holds for domains with smooth boundaries, a conditionthat is too restrictive for a finite element analysis. Theorem 3.1 (smooth domains) . Let s ∈ (0 , , Ω be a domain with ∂ Ω ∈ C ∞ , f ∈ H r (Ω) for some r ≥ − s , u be the solution of (1.4) and α = min { s + r, / − ε } ,with ε > arbitrarily small. Then, u ∈ H s + α (Ω) and the following regularityestimate holds: (cid:107) u (cid:107) H s + α (Ω) (cid:46) (cid:107) f (cid:107) H r (Ω) , where the hidden constant depends on the domain Ω , the dimension d , s and α . As a consequence of the previous result, we see that smoothness of the righthand side f does not ensure that solutions are any smoother than H s + − ε (Ω); seealso [121]. We illustrate this phenomenon with the following example. Example 3.2 (limited regularity) . We follow [60, 106] and consider Ω = B (0 , ⊂ R d and f ≡
1. Then the solution to (3.2) is given by(3.3) u ( x ) = Γ( d )2 s Γ( d +2 s )Γ(1 + s ) (1 − | x | ) s + , where t + = max { t, } .The lack of a lifting property for the solution to (3.2) can also be explained by thefact that the eigenfunctions of this operator have reduced regularity [27, 65, 108].This is in stark contrast with the spectral fractional Laplacian (1.6), discussed inSection 2.2, whose eigenfunctions coincide with those of the Laplacian and thus aresmooth functions if the boundary of the domain is regular enough.On the other hand, H¨older regularity results for (3.2) have been obtained in[107]. They give rise to Sobolev estimates for solutions in terms of H¨older norms ofthe data, that are valid for rough domains. More precisely, we have the followingresult; see [5]. Theorem 3.3 (Lipschitz domains) . Let s ∈ (0 , and Ω be a Lipschitz domainsatisfying the exterior ball condition. If s ∈ (0 , / , let f ∈ C − s (Ω) ; if s = 1 / ,let f ∈ L ∞ (Ω) ; and if s ∈ (1 / , , let f ∈ C β (Ω) for some β > . Then, for every ε > , the solution u of (3.2) belongs to H s + − ε (Ω) , with (cid:107) u (cid:107) H s + 12 − ε (Ω) (cid:46) ε (cid:107) f (cid:107) (cid:63) , where (cid:107) ·(cid:107) (cid:63) denotes the C − s (Ω) , L ∞ (Ω) or C β (Ω) , correspondingly to whether s issmaller, equal or greater than / , and the hidden constant depends on the domain Ω , the dimension d and s . In case s > /
2, the theorem above ensures that the solution u belongs at least to H (Ω). It turns out that to prove the full H s + − ε -regularity, an intermediate stepis to ensure that the gradient of u is actually an L -function. Following [28], this UMERICAL METHODS FOR FRACTIONAL DIFFUSION 15 fact can be proved studying the behavior of the fractional seminorms | · | H − δ (Ω) ,which usually blow up as δ → δ → δ | w | H − δ (Ω) = C ( d ) | w | H (Ω) ∀ w ∈ L (Ω) . Therefore, the technique used in [5] to prove Theorem 3.3 consists of first provingthat the left-hand side of (3.4) remains bounded as δ → u of(3.2), whence u ∈ H (Ω), and next analyzing the regularity of the gradient of u .As already stated in (1.7), the solution to (3.2) behaves like dist( x, ∂ Ω) s forpoints x close to the boundary ∂ Ω. This can clearly be seen in Example 3.2 andexplains the reduced regularity obtained in Theorems 3.1 and 3.3. To capture thisbehavior, we develop estimates in fractional weighted norms, where the weight is apower of the distance to the boundary. Following [5] we introduce the notation δ ( x, x (cid:48) ) = min (cid:8) dist( x, ∂ Ω) , dist( x (cid:48) , ∂ Ω) (cid:9) , and, for (cid:96) = k + s , with k ∈ N and s ∈ (0 , κ ≥
0, we define the norm (cid:107) w (cid:107) H (cid:96)κ (Ω) := (cid:107) w (cid:107) H k (Ω) + (cid:88) | β | = k ¨ Ω × Ω | D β w ( x ) − D β w ( x (cid:48) ) | | x − x (cid:48) | d +2 s δ ( x, x (cid:48) ) κ d x (cid:48) d x. and the associated space(3.5) H (cid:96)κ (Ω) := (cid:8) w ∈ H (cid:96) (Ω) : (cid:107) w (cid:107) H (cid:96)κ (Ω) < ∞ (cid:9) . Although we are interested in the case κ ≥
0, we recall that in the definitionof weighted Sobolev spaces H kκ (Ω), with k being a nonnegative integer, arbitrarypowers of δ can be considered [75, Theorem 3.6]. On the other hand, global versions H (cid:96)κ ( R d ) are defined integrating in the space R d and taking δ as before, but somerestrictions must be taken into account to ensure their completeness. A sufficientcondition is that the weight belongs to the Muckenhoupt class A ( R d ) [73]. In thiscontext, this implies that if | κ | < / H (cid:96)κ ( R d ) are complete. Remark . In radial domains, spaces like (3.5) have been usedto characterize the mapping properties of the fractional Laplacian. In particular,when Ω is the unit ball, upon defining the weight ω ( x ) = 1 − | x | , an expliciteigendecomposition of the operator w (cid:55)→ ( − ∆) s ( ω s w ) is obtained in [6, 55]. Theeigenfunctions are products of solid harmonic polynomials and radial Jacobi poly-nomials or, in one dimension, Gegenbauer polynomials. Mapping properties of thefractional Laplacian can thus be characterized in terms of weighted Sobolev spacesdefined by means of expansions on these polynomials.The regularity in the weighted Sobolev spaces (3.5) reads as follows [5, Proposi-tion 3.12]. Theorem 3.5 (weighted Sobolev estimate) . Let Ω be a bounded, Lipschitz domainsatisfying the exterior ball condition, s ∈ (1 / , , f ∈ C − s (Ω) and u be the solutionof (3.2) . Then, for every ε > we have u ∈ H s − ε / − ε (Ω) and (cid:107) u (cid:107) H s − ε / − ε (Ω) (cid:46) ε (cid:107) f (cid:107) C − s (Ω) , where the hidden constant depends on the domain Ω , the dimension d and s . It is not straightforward to extend Theorem 3.5 to s ∈ (0 , /
2] since, in this case,we cannot invoke Theorem 3.3 to obtain that the solution u belongs to H (Ω).Circumventing this would require to introduce a weight to obtain, for some κ > u ∈ H − (cid:15)κ (Ω) for some κ >
0. However, for this is necessary to obtain aweighted version of (3.4) which, to the best of our knowledge, is not available inthe literature. In spite of this, the numerical experiments we have carried out using graded meshes and s ∈ (0 , /
2] show the same order of convergence as for s ∈ (1 / , f ∈ C β (Ω) for some β ∈ (0 , − s ), (cid:96) < min { β + 2 s, κ + s + 1 / } , then we have(3.6) | u | H (cid:96)κ (Ω) (cid:46) β + (cid:96) − s )(1 + 2( κ − s − (cid:96) )) | f | C β (Ω) , where the hidden constant depends only on Ω and the dimension d . This showsthat increasing the exponent κ of the weight allows for the differentiability order (cid:96) to increase as well. In principle, there is no restriction on κ above; however, in thenext subsection we exploit this weighted regularity by introducing approximationson a family of graded meshes. There we show that the order of convergence (withrespect to the number of degrees of freedom) is only incremented as long as κ < / FEM: A Priori Error Analysis.
The numerical approximation of the solu-tion to (3.2) presents an immediate difficulty: the kernel of the bilinear form (3.1)is singular, and consequently, special care must be taken. For this reason, mostexisting approaches restrict themselves to the one dimensional case ( d = 1). Explo-rations in this direction can be found using finite elements [49], finite differences [67]and Nystr¨om methods [6]. For several dimensions the literature is rather scarce. AMonte Carlo algorithm that avoids dealing with the singular kernel was proposedin [77].Here we present a direct finite element approximation in arbitrary dimensions.Following § T a conforming and shape regular mesh of Ω, consistingof simplicial elements K of diameter bounded by h . In order to present a unifiedapproach for the whole range s ∈ (0 ,
1) we only consider approximations of (3.2)by continuous functions. For s < / U ( T ) defined as in(2.10), the finite element approximation of (3.2) is then the unique solution to theproblem: find U ∈ U ( T ) such that(3.7) (cid:74) U, W (cid:75) = (cid:104) f, W (cid:105) ∀ W ∈ U ( T ) . From this formulation it immediately follows that U is the projection (in the energynorm) of u onto U ( T ). Consequently, we have a C´ea-like best approximation result (cid:57) u − U (cid:57) = inf W ∈ U ( T ) (cid:57) u − W (cid:57) . Thus, in order to obtain a priori rates of convergence, it just remains to bound theenergy-norm distance between the discrete spaces and the solution. One difficultaspect of dealing with fractional seminorms is that they are not additive withrespect to domain decompositions. Nevertheless, it is possible to localize thesenorms [59]: for all w ∈ H s (Ω) we have | w | H s (Ω) ≤ (cid:88) K ∈ T (cid:20) ˆ K ˆ S K | w ( x ) − w ( x (cid:48) ) | | x − x (cid:48) | d +2 s d x (cid:48) d x + C ( d, σ ) sh sK (cid:107) w (cid:107) L ( K ) (cid:21) , where S K is the patch associated with K ∈ T and σ denotes the shape-regularityparameter of the mesh T . From this inequality it follows that, to obtain a priorierror estimates, it suffices to compute interpolation errors over the set of patches { K × S K } K ∈ T . The reduced regularity of solutions implies that we need to re-sort to quasi-interpolation operators; we work with the Scott-Zhang operator Π T [112]. Local stability and approximation properties of this operator were studiedby Ciarlet Jr. in [45]. UMERICAL METHODS FOR FRACTIONAL DIFFUSION 17
Proposition 3.6 (quasi-interpolation estimate) . Let K ∈ T , max { / , s } < (cid:96) ≤ , s ∈ (0 , , and Π T be the Scott-Zhang operator. If w ∈ H (cid:96) (Ω) , then ˆ K ˆ S K | ( w − Π T w )( x ) − ( w − Π T w )( x (cid:48) ) | | x − x (cid:48) | d +2 s d x (cid:48) d x (cid:46) h (cid:96) − sK | w | H (cid:96) ( S K ) , where the hidden constant depends on d , σ , (cid:96) and blows up as s ↑ . The interpolation estimate of Proposition 3.6 shows that, if the meshsize is suf-ficiently small, we deduce an a priori error bound in the energy norm.
Theorem 3.7 (energy error estimate for quasi-uniform meshes) . Let u denote thesolution to (3.2) and denote by U ∈ U ( T ) the solution of the discrete problem (3.7) ,computed over a mesh T consisting of elements with maximum diameter h . Underthe hypotheses of Theorem 3.3 we have (cid:57) u − U (cid:57) (cid:46) h | log h |(cid:107) f (cid:107) (cid:63) , where the hidden constant depends on Ω , s and σ , and (cid:107) · (cid:107) (cid:63) denotes the C − s (Ω) , L ∞ (Ω) or C β (Ω) , correspondingly to whether s is smaller, equal or greater than / . This estimate hinges on the regularity of solutions provided by Theorem 3.3,and thus it depends on H¨older bounds for the data. We now turn our attentionto obtaining a priori error estimates in the L (Ω)-norm. Using Theorem 3.1, anAubin-Nitsche duality argument can be carried out. The proof of the followingproposition follows the steps outlined in [27, Proposition 4.3]. Proposition 3.8 ( L -error estimate) . Let u denote the solution to (3.2) and denoteby U ∈ U ( T ) the solution of the discrete problem (3.7) , computed over a mesh T consisting of elements with maximum diameter h . Under the hypotheses of Theorem3.1 we have (cid:107) u − U (cid:107) L (Ω) (cid:46) h α + β (cid:107) f (cid:107) H r (Ω) , where α = min { s + r, / − ε } , β = min { s, / − ε } , ε > may be taken arbitrarilysmall and the hidden constant depends on Ω , s , d , σ , α and blows up when ε → . Finally, for s ∈ (1 / ,
1) and d = 2, we take advantage of Theorem 3.5, fromwhich further information about the boundary behavior of solutions is available. Wepropose a standard procedure often utilized in connection with corner singularitiesor boundary layers arising in convection-dominated problems. An increased rateof convergence is achieved by resorting to a priori adapted meshes. To obtaininterpolation estimates in the weighted fractional Sobolev spaces defined by (3.5),we introduce the following Poincar´e inequality [5, Proposition 4.8]. Proposition 3.9 (weighted fractional Poincar´e inequality) . Let s ∈ (0 , , κ ∈ [0 , s ) and a domain S which is star-shaped with respect to a ball. Then, for every w ∈ H sκ ( S ) , it holds (cid:107) w − w (cid:107) L ( S ) (cid:46) d s − κS | w | H sκ ( S ) , where w = ffl S w d x , d S = diam( S ) and the hidden constant depends on the chunki-ness parameter of S and the dimension d . This inequality yields sharp quasi-interpolation estimates near the boundary ofthe domain, where the weight δ involved in (3.5) degenerates. Such bounds in turnlead to estimates for the Scott-Zhang operator in weighted fractional spaces. Weexploit them for two-dimensional problems ( d = 2) by constructing graded meshesas in [63, Section 8.4]. In addition to shape regularity, we assume that our meshes T satisfy the following property: there is a number µ ≥ h and K ∈ T , we have(3.8) h K ≤ C ( σ ) (cid:40) h µ , K ∩ ∂ Ω (cid:54) = ∅ ,h dist( K, ∂ Ω) ( µ − /µ , K ∩ ∂ Ω = ∅ . where the constant C ( σ ) depends only on the shape regularity constant σ of themesh T . The parameter µ relates the meshsize h to the number of degrees offreedom because (recall that d = 2) T ≈ (cid:40) h − , µ ≤ ,h − µ , µ > . It is now necessary to relate the parameter µ with the exponent κ of the weight δ in estimate (3.6). Increasing the parameter µ corresponds to raising κ and therebyallowing an increase of the differentiability order (cid:96) . However, if µ > µ = 2 and we have the following result. Theorem 3.10 (energy error estimates for graded meshes) . Let Ω ⊂ R and U ∈ U ( T ) be the solution to (3.7) , computed over a mesh T that satisfies (3.8) with µ = 2 . In the setting of Theorem 3.5 we have (cid:57) u − U (cid:57) (cid:46) ( T ) − | log( T ) | (cid:107) f (cid:107) C − s (Ω) , where the hidden constant depends on σ and blows up as s → / . Implementation.
Let us now discuss key details about the finite elementimplementation of (3.2) for d = 2. If { φ v } are the nodal piecewise linear basisfunctions of U ( T ), defined as in (2.10), then the entries of the stiffness matrix A = ( A vw ) are A vw = (cid:74) φ v , φ w (cid:75) = C ( d, s )2 ¨ Q ( φ v ( x ) − φ v ( x (cid:48) ))( φ w ( x ) − φ w ( x (cid:48) )) | x − x (cid:48) | s d x (cid:48) d x. Two numerical difficulties — coping with integration on unbounded domainsand handling the non-integrable singularity of the kernel — seem to discourage adirect finite element approach. However, borrowing techniques from the boundaryelement method [111], it is possible to compute accurately the entries of the matrix A . We next briefly outline the main steps of this procedure. For full details, werefer to [4] where a finite element code to solve (3.2) is documented.The integrals involved in the computation of A should be carried over R . Forthis reason it is convenient to consider a ball B containing Ω and such that thedistance from Ω to B c is an arbitrary positive number. This is needed in order toavoid difficulties caused by lack of symmetry when dealing with the integral over Ω c when Ω is not a ball. Together with B , we introduce an auxiliary triangulation T A on B \ Ω such that the complete triangulation T B over B (that is T B = T ∪ T A )remains admissible and shape-regular.We define, for 1 ≤ (cid:96), m ≤ T B and K (cid:96) , K m ∈ T B , I v , w (cid:96),m := ˆ K (cid:96) ˆ K m ( φ v ( x ) − φ v ( x (cid:48) ))( φ w ( x ) − φ w ( x (cid:48) )) | x − x (cid:48) | s d x (cid:48) d x,J v , w (cid:96) := ˆ K (cid:96) ˆ B c φ v ( x ) φ w ( x ) | x − x (cid:48) | s d x (cid:48) d x, (3.9)whence we may write A vw = C ( d, s )2 T B (cid:88) (cid:96) =1 (cid:32) T B (cid:88) m =1 I v , w (cid:96),m + 2 J v , w (cid:96) (cid:33) . UMERICAL METHODS FOR FRACTIONAL DIFFUSION 19
We reiterate that computing the integrals I v , w (cid:96),m and J v , w (cid:96) is challenging for differentreasons: the former involves a singular integrand if K (cid:96) ∩ K m (cid:54) = ∅ , while the latterneeds to be calculated in an unbounded domain.We first tackle the computation of I v , w (cid:96),m in (3.9). If K (cid:96) and K m do not touch,then the integrand is a regular function and can be integrated numerically in astandard fashion. On the other hand, if K (cid:96) ∩ K m (cid:54) = ∅ , then I v , w (cid:96),m bears some resem-blances to typical integrals appearing in the boundary element method. Indeed, thequadrature rules we employ are analogous to the ones presented in [111, Chapter5]. Basically, the scheme consists of the following steps: • Consider parametrizations χ (cid:96) : ˆ K → K (cid:96) and χ m : ˆ K → K m such that the edge orvertex shared by K (cid:96) and K m is the image of the same edge/vertex in the referenceelement ˆ K . If K (cid:96) and K m coincide, simply use the same parametrization twice. • Decompose the integration domain into certain subsimplices and then utilizeDuffy-type transformations to map these subdomains into the four-dimensionalunit hypercube. • Since the Jacobian of these Duffy transformations is regularizing, each of the inte-grals may be separated into two parts: a highly singular but explicitly integrablepart and a smooth, numerically tractable part.The second difficulty lies in the calculation of J v , w (cid:96) , namely, dealing with theunbounded domain B c . We write J v , w (cid:96) = ˆ K (cid:96) φ v ( x ) φ w ( x ) (cid:37) ( x ) d x, (cid:37) ( x ) := ˆ B c | x − x (cid:48) | s d x (cid:48) , and realize that we need to accurately compute (cid:37) ( x ) at each quadrature point x ∈ K (cid:96) ∩ ¯Ω. To do so, there are two properties we can take advantage of: theradiality of (cid:37) and the fact that (cid:37) is smooth up to the boundary of Ω because,for x ∈ ¯Ω and x (cid:48) ∈ B c | x − x (cid:48) | > dist(Ω , B c ) > . Therefore, the values of (cid:37) atquadrature nodes can be precomputed with an arbitrary degree of precision.3.4.
Numerical Experiments.
We present the outcome of two experiments posedin Ω = B (0 , ⊂ R .3.4.1. Rate of Convergence in Energy Norm.
Following Example 3.2 we have that,if f ≡
1, then the solution to (3.2) is given by (3.3). Table 1 shows computationalrates of convergence in the energy norm for several values of s , both for uniform andgraded meshes. These rates are in agreement with those predicted by Theorems3.7 and 3.10. Moreover, we observe an increased order of convergence for s ≤ / Value of s Table 1.
Example 3.4.1: Computational rates of convergence for(3.2) posed in the unit ball with right-hand side f ≡
1. Errors aremeasured in the energy norm with respect to the meshsize param-eter h . The second row corresponds to uniform meshes, while thethird to graded meshes, with µ = 2 in (3.8). Rate of Convergence in L -Norm. Remark 3.4 states that a family of explicitsolutions for (3.2) is available in the unit ball. A subclass of solutions in that familymay be expressed in terms of the Jacobi polynomials P ( α,β ) k : [ − , → R . We set f ( x ) = (cid:18) Γ(3 + s )2 − s (cid:19) P ( s, (2 | x | − , so that the solution to (3.2) is given by [55, Theorem 3] u ( x ) = (1 − | x | ) s + P ( s, (2 | x | − . We compute the orders of convergence in L (Ω) for s ∈ { . , . } ; according toProposition 3.8, it is expected to have order of convergence 0 .
75 for s = 0 .
25 and 1for s = 0 .
75 with respect to the meshsize h . The results, summarized in Table 2,agree with the predicted rates of convergence. h s = 0 . s = 0 . . . . . . . . . . . . . . . . Table 2.
Example 3.4.2: Errors in the L -norm for s = 0 .
25 and s = 0 .
75. The estimated orders of convergence with respect to themeshsize are, respectively, 0 . . FEM: A Posteriori Error Analysis.
Since solutions of (3.2) have reducedregularity, and assembling the stiffness matrix A entails a rather high computationalcost, it is of interest to devise suitable AFEMs. We now present a posteriori errorestimates of residual type and ensuing AFEM; we follow [99].We estimate the energy error (cid:57) u − U (cid:57) in terms of the residual R := f − ( − ∆) s U in H − s (Ω). To do so, we need to address two important issues: localization of thenorm in H − s (Ω) and a practical computation of R .To localize fractional norms we deviate from [59] and perform a decomposition onstars S v = supp φ v , the support of the basis functions φ v associated with node v anddiameter h v . Exploiting the partition of unity property (cid:80) v φ v = 1, and Galerkinorthogonality (cid:74) u − U, φ v (cid:75) = 0 for all v ∈ Ω, we can write, for every w ∈ H s (Ω) (cid:74) u − U, w (cid:75) = (cid:104)R , w (cid:105) = (cid:88) v (cid:104)R , wφ v (cid:105) = (cid:88) v (cid:104)R , ( w − ¯ w v ) φ v (cid:105) = (cid:88) v (cid:104) ( R− ¯ R v ) φ v , w − ¯ w v (cid:105) , where ¯ w v ∈ R are weighted mean values computed as ¯ w v = 0 provided v ∈ ∂ Ω and,otherwise, ¯ w v := (cid:104) w, φ v (cid:105)(cid:104) φ v , (cid:105) ∀ v ∈ Ω . The values of ¯ R v ∈ R are yet to be chosen. We see that R = (cid:80) v ( R − ¯ R v ) φ v whereeach term ( R − ¯ R v ) φ v ∈ H − s ( S v ) has support in S v . We have the following twoestimates for dual norms proved in [99, Lemmas 1 and 2]. Lemma 3.11 (localized upper bound of dual norm) . Let
G ∈ H − s (Ω) be decom-posed as G = (cid:80) v g v with g v ∈ H − s ( S v ) vanishing outside S v . We then have for < s < (cid:107)G(cid:107) H − s (Ω) ≤ ( d + 1) (cid:88) v (cid:107) g v (cid:107) H − s ( S v ) . UMERICAL METHODS FOR FRACTIONAL DIFFUSION 21
A key practical issue is then how to evaluate (cid:107) g v (cid:107) H − s ( S v ) for g v = ( R − ¯ R v ) φ v .For second order operators, R splits into an L -component in element interiors(provided f ∈ L (Ω)) and a singular component supported on element boundaries.In contrast, the residual R = f − ( − ∆) s U does not have a singular component andits absolutely continuous part is not always in L (Ω) for all 0 < s < f might be. This is related to singularities of ( − ∆) s φ v ( x ) as x tendsto the skeleton of S v because φ v is continuous, piecewise linear. Using that [119,Theorem XI.2.5] ( − ∆) s : (cid:102) W tp (Ω) → W t − sp (Ω) , t ∈ R , p > , is a continuous pseudo-differential operator of order 2 s , and φ v ∈ (cid:102) W p − εp (Ω) forany ε >
0, we deduce ( − ∆) s U ∈ L p (Ω) 1 p > s − . This motivates the following estimate, whose proof is given in [99, Lemma 2].
Lemma 3.12 (upper bound of local dual norm) . Let g v ∈ L p ( S v ) satisfy ´ S v g v d x =0 for each v ∈ Ω . If < s < and ≤ p < ∞ satisfies p < sd + , then (cid:107) g v (cid:107) H − s ( S v ) (cid:46) h s + d − dp v (cid:107) g v (cid:107) L p ( S v ) . However, to be able to apply Lemma 3.12 the Lebesgue exponent p must satisfy2 s − < p < sd + 12 . We note that for s < we can choose p = 2 for any dimension d . However, for ≤ s < < p <
2. For d = 1 , s <
1, but for d = 3 we have the unfortunate constraint s < .We can now choose ¯ R v as¯ R v := ´ S v R φ v d x ´ S v φ v d x , v ∈ Ωand ¯ R v = 0 otherwise, so that the local contributions satisfy ´ S v ( R − ¯ R v ) φ v d x = 0and we can then apply the bound of Lemma 3.12. The following upper a posteriorierror estimate is derived in [99, Theorem 1]. Theorem 3.13 (upper a posteriori bound) . Let f ∈ L p (Ω) and < p < ∞ ,
2. The second one is the actual computation of ( − ∆) s φ v ( x ) for d >
1, which is problematic due to its singular behavior as x tends to the skeletonon T . This is doable for d = 1 and we refer to [99, Section 8] for details andnumerical experiments. This topic is obviously open for improvement.3.6. Extensions and Applications.
We conclude the discussion by mentioningextensions and applications of this approach: • Eigenvalue problems : The eigenvalue problem for the integral fractional Lapla-cian arises, for example, in the study of fractional quantum mechanics [79]. Asalready mentioned, a major difference between the spectral fractional Laplacian and the integral one is that, for the second one, eigenfunctions have reduced reg-ularity. In [27], conforming finite element approximations were analyzed and itwas shown that the Babuˇska-Osborn theory [13] holds in this context. Regularityresults for the eigenfunctions are derived under the assumption that the domainis Lipschitz and satisfies the exterior ball condition. Numerical evidence on thenon-convex domain Ω = ( − , \ [0 , indicates that the first eigenfunction isas regular as the first one on any smooth domain. This is in contrast with theLaplacian. • Time dependent problems : In [3], problem (2.16) with γ ∈ (0 ,
2] and the integraldefinition of ( − ∆) s was considered. Regularity of solutions was studied and adiscrete scheme was proposed and analyzed. The method is based on a standardGalerkin finite element approximation in space, as described here, while in timea convolution quadrature approach was used [68, 80]. • Non-homogeneous Dirichlet conditions : An interpretation of a non-homogeneousDirichlet condition g for the integral fractional Laplacian is given by using (1.3)upon extension by g over Ω c . In [7] a mixed method for this problem is proposed;it is based on the weak enforcement of the Dirichlet condition and the incorpo-ration of a certain non-local normal derivative as a Lagrange multiplier. Thisnon-local derivative is interpreted as a non-local flux between Ω and Ω c [51]. • Non-local models for interface problems : Consider two materials with permittiv-ities/diffusivities of opposite sign, and separated by an interface with a corner.Strong singularities may appear in the classical (local) models derived from elec-tromagnetics theory. In fact, the problem under consideration is of Fredholmtype if and only if the quotient between the value of permittivities/diffusivitiestaken from both sides of the interface lies outside a so-called critical interval,which always contains the value −
1. In [26] a non-local interaction model for thematerials is proposed. Numerical evidence indicates that the non-local modelmay reduce the critical interval and that solutions are more stable than for thelocal problem.4.
Dunford-Taylor Approach for Spectral and Integral Laplacians
In this section we present an alternative approach to the ones developed in theprevious sections. It relies on the Dunford-Taylor representation (1.11) u = ( − ∆) − s f = sin( sπ ) π ˆ ∞ µ − s ( µ − ∆) − f d µ for the spectral fractional Laplacian (1.6). For the integral fractional Laplacian(1.2), instead, it hinges on the equivalent representation (1.13) of (1.12): ˆ R d | ξ | s F − (˜ u ) | ξ | s F − ( ˜ w ) d ξ = 2 sin( sπ ) π ˆ ∞ µ − s ˆ Ω (cid:0) ( − ∆)( I − µ ∆) − ˜ u ( x ) (cid:1) w ( x ) d x d µ. (4.1)In (4.1), the operators − ∆ and I − µ ∆ are defined over R d , something to be madeprecise in Theorem 4.5.In each case, the proposed method is proved to be efficient on general Lipschitzdomains Ω ⊂ R d . They rely on sinc quadratures and on finite element approxi-mations of the resulting integrands at each quadrature points. While (1.11) allowsfor a direct approximation of the solution, the approximation of (1.13) leads toa non-conforming method where the action of the stiffness matrix on a vector isapproximated.We recall that the functional spaces H s (Ω) are defined in (1.5) for s ∈ [0 , / H s (Ω) = H s (Ω) ∩ H (Ω) for s ∈ (1 , UMERICAL METHODS FOR FRACTIONAL DIFFUSION 23
Spectral Laplacian.
We follow [24, 25] and describe a method based onthe Balakrishnan representation (1.11). In order to simplify the notation, we set L := − ∆ : D ( L ) → L (Ω) and define the domain of L r , for r ∈ R , to be D ( L r ) := { v ∈ L (Ω) : L r v ∈ L (Ω) } ;this is a Banach space equipped with the norm (cid:107) v (cid:107) D ( L r ) := (cid:107) L r v (cid:107) L (Ω) . We also define the solution operator T : H − (Ω) → H (Ω) by T f := v , wherefor F ∈ H − (Ω), v ∈ H (Ω) is the unique solution of ˆ Ω ∇ v · ∇ w d x = F ( v ) , ∀ w ∈ H (Ω) . This definition directly implies that D ( L ) = range( T | L (Ω) ).4.1.1. Finite Element Discretization.
For simplicity, we assume that the domain Ωis polytopal so that it can be partitioned into a conforming subdivision T . Werecall that U ( T ) ⊂ H (Ω) stands for the subspace of globally continuous piecewiselinear polynomials with respect to T ; see Section 2.4. We denote by Π T the L (Ω)-orthogonal projection onto U ( T ) and by L T : U ( T ) → U ( T ) the finite elementapproximation of L , i.e., for V ∈ U ( T ), L T V ∈ U ( T ) solves ˆ Ω ( L T V ) W d x = ˆ Ω ∇ V · ∇ W d x, ∀ W ∈ U ( T ) . We finally denote by T T the inverse of L T , the finite element solution operator,and by h the maximum diameter of elements in T .With these notations, we are in the position to introduce the finite elementapproximation U ∈ U ( T ) of u in (1.10):(4.2) U := sin( sπ ) π ˆ ∞ µ − s ( µ + L T ) − Π T f d µ. The efficiency of the approximation of u by U depends on the efficiency of thefinite element solver ( T T Π T ) f (i.e. for the standard Laplacian), which is dictatedby the regularity of T f . This regularity aspect has been intensively discussed in theliterature [14, 20, 48, 69, 72, 91]. In this exposition, we make the following generalassumption.
Definition 4.1 (elliptic regularity) . We say that T satisfies a pick-up regularityof index 0 < α ≤ ≤ r ≤ α , the operator T is an isomorphism from H − r (Ω) to H r (Ω).Notice that α = 1 when Ω is convex, whence this definition extends (2.8) togeneral Lipschitz domains.Assuming a pick-up regularity of index α , for any r ∈ [0 , (cid:107) T w − T T Π T w (cid:107) H r (Ω) (cid:46) h α ∗ (cid:107) T w (cid:107) H α +1 (Ω) (cid:46) h α ∗ (cid:107) w (cid:107) H α − (Ω) where α ∗ := (cid:0) α + min(1 − r, α ) (cid:1) . The proof of the above estimate is classical andis based on a duality argument (Nitsche’s trick); see e.g. [25, Lemma 6.1]. Noticethat α ∗ < α when 1 − α < r , i.e, the error is measured with regularity index toolarge to take full advantage of the pick-up regularity in the duality argument.We expect that approximation (4.2) of the fractional Laplacian problem deliversthe same rate of convergence(4.3) (cid:107) u − U (cid:107) H r (Ω) (cid:46) h α ∗ (cid:107) u (cid:107) H α (Ω) (cid:46) h α ∗ (cid:107) f (cid:107) H α − s (Ω) , but the function U in (4.2) is well defined provided f ∈ L (Ω), i.e. 1 + α − s ≥ Before describing the finite element approximation result, we make the follow-ing comments. The solution u = L − s f belongs to D ( L (1+ α ) / ) provided that f ∈ D ( L (1+ α ) / − s ). Hence, estimates such as (4.3) rely on the characterizationof D ( L r/ ) for r ∈ (0 , α ]. For 0 ≤ r ≤
1, the spaces D ( L r/ ) and H r (Ω)are equivalent, as they are both scale spaces which coincide at r = 0 and r = 1.Furthermore, assuming an elliptic regularity pick-up of index 0 < α ≤
1, thischaracterization extends up to 1 + α [25, Theorem 6.4 and Remark 4.2]. Theorem 4.2 (finite element approximation) . Assume that T satisfies a pick-up regularity of index α ∈ (0 , on Ω . Given r ∈ [0 , with r ≤ s , set γ :=max { r + 2 α ∗ − s, } and α ∗ := (cid:0) α + min(1 − r, α ) (cid:1) . If f ∈ H δ (Ω) for δ ≥ γ , then (cid:107) u − U (cid:107) H r (Ω) ≤ C h h α ∗ (cid:107) f (cid:107) H δ (Ω) , where C h (cid:46) log(2 /h ) , when δ = γ and r + 2 α ∗ ≥ s, , when δ > γ and r + 2 α ∗ ≥ s, , when δ = 0 and s > r + 2 α ∗ . Sinc Quadrature.
It remains to put in place a sinc quadrature [82] to ap-proximate the integral in (4.2). We use the change of variable µ = e y so that U = sin( sπ ) π ˆ ∞−∞ e (1 − s ) y ( e y + L T ) − Π T f d y. Given k >
0, we set N + := (cid:24) π sk (cid:25) , N − := (cid:24) π − s ) k (cid:25) , and y (cid:96) := k(cid:96), and define the sinc quadrature approximation of U by(4.4) U k := sin( sπ ) π k N + (cid:88) (cid:96) = − N − e (1 − s ) y (cid:96) ( e y (cid:96) + L T ) − Π T f. The sinc quadrature consists of uniformly distributed quadrature points in the y variable, and the choice of N + and N − makes it more robust with respect to s .The decay when | z | → + ∞ and holomorphic properties of the integrand z − s ( z − L ) − in the Dunford-Taylor representation (1.10) guarantee the exponential con-vergence of the sinc quadrature [25, Theorem 7.1]. Theorem 4.3 (sinc quadrature) . For r ∈ [0 , , we have (cid:107) U − U k (cid:107) H r (Ω) (cid:46) e − π / (2 k ) (cid:107) f (cid:107) H r (Ω) . To compare with Theorem 2.4, we take r = s and assume that Ω is convex,which allows for any α in (0 , N + ≈ N − ≈ log(1 /h ) so that sinc quadrature and finite element errors are balanced.Therefore, Theorems 4.2 and 4.3 yield for the Dunford-Taylor method (cid:107) u − U k (cid:107) H s (Ω) (cid:46) ( T ) − α ∗ /d (cid:107) f (cid:107) H σ (Ω) , where σ := max(2 α ∗ − s, s ) and ( T ) − /d ≈ h for quasi-uniform subdivisions,provided we discard logarithmic terms. In contrast, the error estimate of Theorem2.4 for the extension method reads, again discarding log terms, (cid:107) u − U (cid:107) H s (Ω) (cid:46) ( T Y ) − / ( d +1) (cid:107) f (cid:107) H − s (Ω) and was derived with pick-up regularity α = 1. We first observe the presence ofthe exponent d + 1, which makes the preceding error estimate suboptimal. This UMERICAL METHODS FOR FRACTIONAL DIFFUSION 25 can be cured with geometric grading in the extended variable and hp -methodology.Section 2.7 and [17, 86] show that this approach yields the following error estimate (cid:107) u − U (cid:107) H s (Ω) (cid:46) ( T ) − /d (cid:107) f (cid:107) H − s (Ω) , with T degrees of freedom, after discarding logarithmic terms. This estimateexhibits quasi–optimal linear order for the regularity f ∈ H − s (Ω). We also see thatthe Dunford-Taylor method possesses the optimal rate of convergence 2 α ∗ = 2 − s > f ∈ H σ (Ω): σ = 2(1 − s ) when s ≤ / σ = s when s > /
3. We may also wonderwhat regularity of f would lead to the same linear order of convergence as theextension method, that requires f ∈ H − s (Ω). We argue as follows: if s ≤ ,then f ∈ H − s (Ω); otherwise, if s > , then f ∈ H s (Ω). We thus realize that theregularity of f is the same for s ≤ but it is stronger for s > .It is worth mentioning that the Dunford-Taylor algorithm seems advantageousin a multi-processor context as it appears to exhibit good strong and weak scalingproperties. The former consists of increasing the number of processors for a fixednumber of degrees of freedom, while for the latter, the number of degrees of freedomper processor is kept constant when increasing the problem size. We refer to [46]for a comparison of different methods. Remark . The method based on (4.4) requires N + + N − + 1 independent standard Laplacian finite element solves for each quadrature points y (cid:96) : V (cid:96) ∈ U ( T ) : e y (cid:96) ˆ Ω V (cid:96) W d x + ˆ Ω ∇ V (cid:96) · ∇ W d x = ˆ Ω f W d x ∀ W ∈ U ( T ) , which are then aggregated to yield U k : U k = sin( sπ ) π k N + (cid:88) (cid:96) = − N − e (1 − s ) y (cid:96) V (cid:96) . Implementation of this algorithm starting from a finite element solver for the Pois-son problem is straightforward. Numerical illustrations matching the predictedconvergence rates of Theorems 4.2 and 4.3 are provided in [24].4.2.1.
Extensions.
We now discuss several extensions. • Symmetric operators and other boundary conditions.
The operator L = − ∆ canbe replaced by any symmetric elliptic operators as long as the associated bilinearform ( u, w ) (cid:55)→ ´ Ω Lu w remains coercive and bounded in H (Ω). Different bound-ary conditions can be considered similarly as well. However, it is worth pointingout that the characterization of D ( L r ) depends on the boundary condition andmust be established.As an illustration, Figure 4 depicts the approximations using parametric sur-face finite element [56] of the solution to(4.5) ( − ∆ Γ ) s u = 1 on Γ , u = 0 on ∂ Γ , where ∆ Γ is the surface Laplacian on Γ ⊂ R , either the side boundary of acylinder or given byΓ := (cid:8) ( x + 2 sin( x ) , x + 2 cos( x ) , x ) ∈ R : ( x , x , x ) ∈ S , x ≥ (cid:9) . (4.6) • Regularly accretive operators.
The class of operators L can be extended furtherto a subclass of non-symmetric operators. They are the unbounded operatorsassociated with coercive and bounded sesquilinear forms in H (Ω) (regularlyaccretive operators [71]). In this case, fractional powers cannot be defined usinga spectral decomposition as in (1.6) but rather directly using the Dunford-Taylor Figure 4.
Dunford-Taylor Method: Numerical approximation ofthe solution to the spectral Laplacian problem (4.5) on an hyper-surface (darker = smaller values; lighter = larger values). (Left) s = 0 . s = 0 . D ( L r/ ) in termsof Sobolev regularity. It turns out that for − < r <
1, we have that D ( L r/ ) isthe same as for the symmetric operator [70] D ( L r/ ) = D (( L + L ∗ ) r/ ) = H r (Ω) , where L ∗ denotes the adjoint of L . This characterization does not generallyhold for r = 1 (Kato square root problem [71]). McKintosh [84] proved that D ( L / ) = H (Ω) = H (Ω) for sesquilinear forms of the type( ψ, ϕ ) (cid:55)→ ˆ Ω ( B ∇ ψ · ∇ ϕ + β · ∇ ψϕ + ψβ · ∇ ϕ + cψϕ ) d x, where B ∈ L ∞ (Ω , GL ( R d )), β , β ∈ L ∞ (Ω , R d ) and c ∈ L ∞ (Ω) are such thatthe form is coercive and bounded. This characterization is extended in [25,Theorem 6.4] up to r = 1 + α so that similar convergence estimates to those inTheorems 4.2 and 4.3 are established. To illustrate the method for non-symmetricoperators, we consider the following example (cid:18) − ∆ + (cid:18) (cid:19) · ∇ (cid:19) s u = 1 , in Ω , u = 0 on ∂ Ω , where Ω = ( − , \ [0 , is a L-shaped domain. Figure 5 reports the fullydiscrete approximations given by (4.4) for s = 0 . , . , . • Space and time fractional diffusion.
In [22, 21] the space-time fractional problem ∂ γt u + ( − ∆) s u = f, u (0) = u , is studied, where ∂ γt denotes the so-called Caputo derivative of order γ ∈ (0 , u ( t ) = e γ, ( − t γ L s ) u + ˆ t ξ γ − e γ,γ ( − ξ γ L s ) f ( t − ξ ) d ξ, where, again in this case, a Dunford-Taylor representation can be used to write e γ,µ ( − t γ L s ) = 12 πi ˆ C e γ,µ ( − t γ z s )( z − L ) − d z and e γ,µ defined on C is the Mittag-Leffler function. Because of the presenceof e γ,µ , the contour C cannot be deformed onto the negative real axis anymore,which prevents a representation like in (1.11). Instead, the sinc quadrature is UMERICAL METHODS FOR FRACTIONAL DIFFUSION 27
Figure 5.
Dunford-Taylor Method: Approximation of the solu-tion to fractional advection - diffusion problem on a L-shaped do-main. (Left) Solution with isolines for s = 0 .
5. (Right) Plots of u for s = 0 . , . , . s .performed directly on a hyperbolic parametrization of C in the complex plane.Nevertheless, we obtain the error estimates (cid:107) u ( t ) − U k ( t ) (cid:107) L (Ω) (cid:46) (cid:16) t − γα ∗ /s h α ∗ + t − γ e − c/k (cid:17) (cid:107) u (cid:107) L (Ω) , when f = 0 and where α ∗ := min( α, s − ) with s − denoting any number strictlysmaller than s . Here c is a constant independent of h and k . We refer the readerto [22, 21] for estimates measuring the error in higher norms or for improvedresults (in the singularity when t →
0) when u is smoother. Note that therepresentation used does not need a time-stepping method for the initial valueproblem.When u = 0 instead, a graded (a-priori known and depending on γ ) meshin time towards t = 0 is put forward. A midpoint quadrature scheme (secondorder) in time for a total of N log( N ) time steps yields the error estimates (cid:107) u ( t ) − U k, N ( t ) (cid:107) L (Ω) (cid:46) t (1 − α ∗ /s ) γ h α ∗ (cid:107) f (cid:107) L ((0 ,t ) × Ω) + max { t γ , t / γ }N − (cid:107) f (cid:107) H (0 ,t ; L (Ω)) + log( N ) e − c/k (cid:107) f (cid:107) L ∞ (0 ,t ; L (Ω)) . Notice that the method exhibits second order convergence rate (up to a loga-rithmic term) with respect to the number of time intervals. Again, we refer to[21, 22] for more details as well as additional estimates when measuring the errorin higher order norms.4.3.
Integral Laplacian.
The strategy used for the spectral Laplacian in the pre-vious section cannot be used for the integral Laplacian. In fact, formulas like (1.10)are not well defined (the integral Laplacian is not strictly positive).Instead, we rely on the following equivalent representation of the bilinear form(4.1) in the weak formulation (1.12). Recall that fractional order Sobolev spaces in R d are defined and normed by H r ( R d ) = (cid:40) w ∈ L ( R d ) : (cid:107) w (cid:107) H r ( R d ) = (cid:18) ˆ R d (1 + | ξ | ) r/ | F ( w )( ξ ) | d ξ (cid:19) / < ∞ (cid:41) for r >
0, and that the notation ˜ w stands for the zero extension of w outside Ω sothat w ∈ H r (Ω) if and only if ˜ w ∈ H r ( R d ) for r ∈ [0 , / Theorem 4.5 (equivalent representation) . Let s ∈ (0 , and ≤ r ≤ s . For η ∈ H r + s ( R d ) and θ ∈ H s − r ( R d ) , ˆ R d | ξ | r + s F ( η )( ξ ) | ξ | s − r F ( θ )( ξ ) d ξ = 2 sin( sπ ) π ˆ ∞ µ − s ˆ R d (cid:0) ( − ∆)( I − µ ∆) − η (cid:1) θ d x d µ. To prove the above theorem, it suffices to note that using Parseval’s theorem ˆ R d (cid:0) ( − ∆)( I − µ ∆) − η (cid:1) θ d x = ˆ R d | ξ | µ | ξ | F ( η )( ξ ) F ( θ )( ξ ) d ξ. and use the change of variable t = µ | ξ | together with the relation ˆ ∞ t − s t d t = π πs ) . For more details, we refer to [23, Theorem 4.1].In order to make the above representation more amenable to numerical methods,for ψ ∈ L ( R d ), we define v ( ψ, µ ) := v ( µ ) ∈ H ( R d ) to be the solution to(4.7) ˆ R d v ( µ ) φ d x + µ ˆ R d ∇ v ( µ ) · ∇ φ d x = − ˆ R d ψφ d x, ∀ φ ∈ H ( R d ) . Using this notation along with definition (1.5) of H s (Ω), we realize that for η, θ ∈ H s (Ω) with s ∈ (0 , ˆ R d | ξ | s F (˜ η )( ξ ) | ξ | s F (˜ θ )( ξ ) d ξ = 2 sin( sπ ) π ˆ ∞ µ − − s (cid:16) ˆ Ω (cid:0) η + v (˜ η, µ ) (cid:1) θ d x (cid:17) d µ ;note that v (˜ η, µ ) does not vanish outside Ω. This prompts the definition(4.8) a ( η, θ ) := 2 sin( sπ ) π ˆ ∞ µ − − s (cid:16) ˆ Ω (cid:0) η + v (˜ η, µ ) (cid:1) θ d x (cid:17) d µ, for η, θ ∈ H s (Ω). The above representation is the starting point of the proposednumerical method. The solution u ∈ H s (Ω) of the fractional Laplacian satisfies(4.9) a ( u, w ) = (cid:104) f, w (cid:105) ∀ w ∈ H s (Ω) . We discuss in Section 4.3.4 a Strang’s type argument to assess the discretizationerror from the consistency errors generated by the approximation of a ( · , · ) usingsinc quadratures, domain truncations, and finite element discretizations.4.3.1. Sinc Quadrature.
We proceed as in Section 4.2 for the spectral Laplacianand use the change of variable µ = e − y to arrive at a ( η, θ ) = sin( sπ ) π ˆ ∞−∞ e sy (cid:18) ˆ Ω (cid:0) η + v (˜ η, µ ( y )) (cid:1) θ d x (cid:19) d y. Given a quadrature spacing k > N + and N − two positive integers, the sincquadrature approximation of a ( · , · ) is given by(4.10) a k ( η, θ ) := sin( sπ ) π k N + (cid:88) (cid:96) = − N − e sy (cid:96) ˆ Ω (cid:0) η + v (˜ η, µ ( y (cid:96) )) (cid:1) θ d x. Notice that we only emphasize the dependency in k in the approximation of a ( · , · )as we will select N + and N − as a function of k .The consistency error between a k ( · , · ) and a ( · , · ) is described in the followingresult. We simply note that, as for the spectral Laplacian discussed in Section 4.1, UMERICAL METHODS FOR FRACTIONAL DIFFUSION 29 the proof of Theorem 4.6 is given in [23, Theorem 5.1 and Remark 5.1] and relieson the holomorphic property and decay as µ → ∞ of the integrand in (4.8). Theorem 4.6 (quadrature consistency) . Given θ ∈ H s (Ω) and η ∈ H δ (Ω) with s < δ ≤ min(2 − s, σ ) , where σ stands for any number strictly smaller than / .Set N + := (cid:108) π k ( δ − s ) (cid:109) and N − := (cid:108) π sk (cid:109) . Then, we have | a ( η, θ ) − a k ( η, θ ) | (cid:46) max (cid:18) δ − s , s (cid:19) e − π / (4 k ) (cid:107) η (cid:107) H δ (Ω) (cid:107) θ (cid:107) H s (Ω) . Truncated Problems.
The sinc approximation of a ( · , · ) defined by (4.10) re-quires the computation of v (˜ η, µ ( y (cid:96) )) for each quadrature point y (cid:96) (here η ∈ H β (Ω)for some s < β < / I − µ ( y (cid:96) ) ∆) − on R d . The proposed method relies on truncationsof this infinite domain problem and uses standard finite elements on the resultingbounded domains. As we shall see, the truncated domain diameter must dependon the quadrature point y (cid:96) .We let B be a convex bounded domain containing Ω and the origin of R d . With-out loss of generality, we assume that the diameter of B is 1. For a truncationparameter M , we define the dilated domains B M ( µ ) := (cid:26) { y = (1 + µ (1 + M ) x : x ∈ B } , µ ≥ , { y = (2 + M ) x : x ∈ B } , µ < , and for ψ ∈ L ( R d ), the associated functions v M ( µ ) := v M ( ψ, µ ) ∈ H ( B M ( µ ))satisfying ˆ B M ( µ ) v M ( µ ) w d x + µ ˆ B M ( µ ) ∇ v M ( µ ) · ∇ w d x = − ˆ B M ( µ ) ψw d x ∀ w ∈ H ( B M ( µ ));(4.11)compare with (4.7). The exponential decay of the function v (˜ η, µ ) yields (cid:107) v (˜ η, µ ) − v M (˜ η, µ ) (cid:107) L ( B M ( µ )) (cid:46) e − max(1 ,µ ) cM/µ (cid:107) η (cid:107) L (Ω) , where c is a constant independent of M and µ (see [23, Lemma 6.1]). As a conse-quence, the truncation consistency in using a k,M ( η, θ ) := sin( sπ ) π k N + (cid:88) (cid:96) = − N − e sy (cid:96) ˆ Ω ( η + v M (˜ η, µ ( y (cid:96) ))) θ d x instead of a k ( η, θ ) decays exponentially fast as a function of M [23, Theorem 6.2]. Theorem 4.7 (truncation consistency) . For M sufficiently large, there is positiveconstant c independent of M and k such that for all η, θ ∈ L (Ω) | a k ( η, θ ) − a k,M ( η, θ ) | (cid:46) e − cM (cid:107) η (cid:107) L (Ω) (cid:107) θ (cid:107) L (Ω) . Finite Element Discretization.
We now turn our attention to the finite ele-ment approximation of v M (˜ η, µ ) defined by (4.11). For simplicity, we assume thatthe domain Ω is polytopal so that it can be partitioned into a conforming subdivision T with elements of maximum diameter h as in Section 4.1.1. Generic constantsbelow may depend on the shape regularity and quasi-uniformity constants of T without mention of it.We need two subspaces of globally continuous piecewise linear polynomials. Thefirst one, U ( T ) ⊂ H (Ω), is defined in (2.10) relative to the partition T . Thesecond subspace, denoted U ( T M ( µ )), has a similar definition but relative to thesubdivision T M ( µ ) of B M ( µ ( y (cid:96) )). We impose that the partitions T M ( µ ) match in Ω, which implies that restrictions of functions in U ( T M ( µ )) are continuouspiecewise linears over T . We refer to [23] for details on the constructions of suchpartitions, which is the bottleneck of the proposed method.We now define for any ψ ∈ L ( R d ) the finite element approximation V M ( µ ) = V M ( ψ, µ ) ∈ U ( T M ( µ )) of the function v M ( µ ) to be ˆ B M ( µ ) V M ( µ ) W d x + µ ˆ B M ( µ ) ∇ V M ( µ ) · ∇ W d x = − ˆ B M ( µ ) ψW d x ∀ W ∈ U ( T M ( µ )) . The fully discrete approximation of the bilinear form a ( · , · ) then reads a k,M T (Ξ , Θ) := sin( sπ ) π k N + (cid:88) (cid:96) = − N − e sy (cid:96) ˆ Ω (cid:0) Ξ + V M ( (cid:101) Ξ , µ ( y (cid:96) )) (cid:1) Θ d x ∀ Ξ , Θ ∈ U ( T ) . Note that V M ( (cid:101) Ξ , µ ) is piecewise linear over T and the sum Ξ + V M ( (cid:101) Ξ , µ ) is easyto perform. The consistency error between a k,M ( · , · ) and a k,M T ( · , · ) is given next;see [23, Theorem 7.6] for a proof. Theorem 4.8 (finite element consistency) . If β ∈ ( s, / and α ∈ (0 , min( s, / ,then the following estimate is valid | a k,M (Ξ , Θ) − a k,M T (Ξ , Θ) | (cid:46) | log h | h β + α − s (cid:107) Ξ (cid:107) H β (Ω) (cid:107) Θ (cid:107) H s + α (Ω) , for all Ξ , Θ ∈ U ( T ) . Strang’s Lemma.
In addition to the three consistency estimates describedabove, Strang’s Lemma requires the U ( T )-ellipticity of the fully discrete form a k,M T ( · , · ). To show this, [23] invokes Theorem 4.6 (quadrature consistency) with δ := min { − s, β } and an inverse estimate to write for Ξ , Θ ∈ U ( T ) | a k (Ξ , Θ) − a (Ξ , Θ) | (cid:46) e − π / (4 k ) h s − δ (cid:107) Ξ (cid:107) H s (Ω) (cid:107) Θ (cid:107) H s (Ω) . This, together with the monotonicity property a k,M T (Ξ , Θ) ≥ a k (Ξ , Θ) , and the coercivity of the exact bilinear form a ( · , · ) in H s (Ω), yields the U ( T )-ellipticity of a k,M T ( · , · ) provided(4.12) e − π / (4 k ) h s − δ ≤ c for an explicit constant c .The fully discrete approximation U k,M ∈ U ( T ) of u satisfying (4.9) is given by a k,M T ( U k,M , W ) = ˆ Ω f W d x ∀ W ∈ U ( T ) . To measure the discrepancy between u and U k,M in H s (Ω), we assume the addi-tional regularity u ∈ H β (Ω) for some β ∈ ( s, / u ,solution to the integral fractional Laplacian, is discussed in Theorems 3.1 and 3.3.The theorem below, proved in [23, Theorem 7.8]), guarantees that the proposedmethod delivers an optimal rate of convergence (up to a logarithmic factor). Theorem 4.9 (error estimate) . Assume that (4.12) holds and that u ∈ H β (Ω) forsome β ∈ ( s, / . Then there is a constant c independent of h , M , and k suchthat (cid:107) u − U k,M (cid:107) H s (Ω) (cid:46) (cid:16) e − π / (4 k ) + e − cM + | log h | h β − s (cid:17) (cid:107) u (cid:107) H β (Ω) . UMERICAL METHODS FOR FRACTIONAL DIFFUSION 31
We also refer to [23, Theorem 7.8] for further discussions on mesh generation,matrix representation of the fully discrete scheme, and a preconditioned iterativemethod.Notice that the error estimate stated in Theorem 4.9 for quasi-uniform meshesis of order about h / and is similar to the one derived for the integral methodin Theorem 3.7. To see this, we choose β = s + − ε with (cid:15) > u guaranteed by Theorem 3.3, along with M ≈ log(1 /h ) and N + ≈ N − ≈ | log h | to balance the three sources of errors. Itis worth mentioning that using graded meshes for d = 2, Theorem 3.10 states anoptimal linear rate of convergence (up to a logarithmic factor) provided s ∈ (1 / , Numerical Experiment.
To illustrate the method, we depict in Figure 6 theapproximation U k,M for s = 0 . f = 1, and Ω = B (0 , R . Figure 6.
Dunford-Integral Method: Numerical approximationof the solution to the fractional integral Laplacian for s = 0 . f = 1 in the unit ball of R (darker = 0.0, whiter = 0.7). The linesrepresent the isosurfaces { u ( x ) = k/ } for k = 0 , ..., References [1] S. Abe and S. Thurner. Anomalous diffusion in view of Einstein’s 1905 theory of Brownianmotion.
Physica A: Statistical Mechanics and its Applications , 356(2–4):403 – 407, 2005.[2] M. Abramowitz and I.A. Stegun.
Handbook of mathematical functions with formulas, graphs,and mathematical tables , volume 55 of
National Bureau of Standards Applied MathematicsSeries . For sale by the Superintendent of Documents, U.S. Government Printing Office,Washington, D.C., 1964.[3] G. Acosta, F. Bersetche, and J.P. Borthagaray. Finite element approximations for fractionalevolution problems. arXiv:1705.09815 , 2017.[4] G. Acosta, F. Bersetche, and J.P. Borthagaray. A short FEM implementation for a 2d homo-geneous Dirichlet problem of a fractional Laplacian.
Comput. Math. Appl. , 2017. Acceptedfor publication.[5] G. Acosta and J.P. Borthagaray. A fractional Laplace equation: regularity of solutions andfinite element approximations.
SIAM J. Numer. Anal. , 55(2):472–495, 2017.[6] G. Acosta, J.P. Borthagaray, O. Bruno, and M. Maas. Regularity theory and high ordernumerical methods for one-dimensional fractional-Laplacian equations.
Math. Comp. , 2017.Accepted for publication. [7] G. Acosta, J.P. Borthagaray, and N. Heuer. Finite element approximations for the nonho-mogeneous fractional Dirichlet problem. In preparation.[8] H. Antil and E. Ot´arola. A FEM for an optimal control problem of fractional powers ofelliptic operators.
SIAM J. Control Optim. , 53(6):3432–3456, 2015.[9] H. Antil and E. Ot´arola. An a posteriori error analysis for an optimal control probleminvolving the fractional Laplacian.
IMA J. Numer. Anal. , 2017. (to appear).[10] H. Antil, E. Ot´arola, and A.J. Salgado. A space-time fractional optimal control problem:analysis and discretization.
SIAM J. Control Optim. , 54(3):1295–1328, 2016.[11] H. Antil, E. Ot´arola, and A.J. Salgado. Optimization with respect to order in a fractionaldiffusion model: analysis, approximation and algorithm aspects. arXiv:1612.08982, 2017.[12] I. Babuˇska and A. Miller. A feedback finite element method with a posteriori error esti-mation. I. The finite element method and some basic properties of the a posteriori errorestimator.
Comput. Methods Appl. Mech. Engrg. , 61(1):1–40, 1987.[13] I. Babuˇska and J. Osborn. Eigenvalue problems. In
Handbook of numerical analysis, Vol.II , Handb. Numer. Anal., II, pages 641–787. North-Holland, Amsterdam, 1991.[14] C. Bacuta, J.H. Bramble, and J.E. Pasciak. New interpolation results and applications tofinite element methods for elliptic boundary value problems.
East-West J. Numer. Math. ,3:179–198, 2001.[15] W. Bangerth, R. Hartmann, and G. Kanschat. deal.II—diferential equations analysis library.
Technical Reference: http//dealii.org .[16] W. Bangerth, R. Hartmann, and G. Kanschat. deal.II—a general-purpose object-orientedfinite element library.
ACM Trans. Math. Software , 33(4):Art. 24, 27, 2007.[17] L. Banjai, J.M. Melenk, R.H. Nochetto, E. Ot´arola, A.J. Salgado, and Ch. Schwab. LocalFEM for the spectral fractional Laplacian with near optimal complexity.
In preparation ,2017.[18] J. Bertoin.
L´evy processes , volume 121 of
Cambridge Tracts in Mathematics . CambridgeUniversity Press, Cambridge, 1996.[19] M.ˇS. Birman and M.Z. Solomjak.
Spektralnaya teoriya samosopryazhennykh operatorov vgilbertovom prostranstve . Leningrad. Univ., Leningrad, 1980.[20] A. Bonito, J.-L. Guermond, and F. Luddens. Regularity of the Maxwell equations in het-erogeneous media and Lipschitz domains.
J. Math. Anal. Appl. , 408(2):498–512, 2013.[21] A. Bonito, W. Lei, and J.E. Pasciak. The approximation of parabolic equations involvingfractional powers of elliptic operators.
Journal of Computational and Applied Mathematics ,315:32–48, 2017.[22] A. Bonito, W. Lei, and J.E. Pasciak. Numerical approximation of space-time fractionalparabolic equations. arXiv preprint arXiv:1704.04254 , 2017.[23] A. Bonito, W. Lei, and J.E. Pasciak. Numerical approximation of the integral fractionalLaplacian. submitted , 2017.[24] A. Bonito and J. Pasciak. Numerical approximation of fractional powers of elliptic operators.
Mathematics of Computation , 84(295):2083–2110, 2015.[25] A. Bonito and J.E. Pasciak. Numerical approximation of fractional powers of regularly ac-cretive operators.
IMA Journal of Numerical Analysis , 2016.[26] J.P. Borthagaray and P. Ciarlet, Jr. Nonlocal models for interface problems between di-electrics and metamaterials. In , 2017.[27] J.P. Borthagaray, L.M. Del Pezzo, and S. Mart´ınez. Finite element approximation for thefractional eigenvalue problem. arXiv:1603.00317 , 2016.[28] J. Bourgain, H. Brezis, and P. Mironescu. Another look at Sobolev spaces. In
OptimalControl and Partial Differential Equations , pages 439–455, 2001.[29] C. Br¨andle, E. Colorado, A. de Pablo, and U. S´anchez. A concave-convex elliptic probleminvolving the fractional Laplacian.
Proc. Roy. Soc. Edinburgh Sect. A , 143(1):39–71, 2013.[30] S.C. Brenner and L.R. Scott.
The mathematical theory of finite element methods , volume 15of
Texts in Applied Mathematics . Springer, New York, third edition, 2008.[31] D. Brockmann, L. Hufnagel, and T. Geisel. The scaling laws of human travel.
Nature ,439(7075):462–465, 2006.[32] C. Bucur and E. Valdinoci.
Nonlocal diffusion and applications , volume 20 of
Lecture Notesof the Unione Matematica Italiana . Springer; Unione Matematica Italiana, Bologna, 2016.[33] X. Cabr´e and J. Tan. Positive solutions of nonlinear problems involving the square root ofthe Laplacian.
Adv. Math. , 224(5):2052–2093, 2010.[34] L. Caffarelli and A. Figalli. Regularity of solutions to the parabolic fractional obstacle prob-lem.
J. Reine Angew. Math. , 680:191–233, 2013.[35] L. Caffarelli and Stinga P. Fractional elliptic equations, Caccioppoli estimates, and regular-ity.
Annales de l’Institut Henri Poincare (C) Non Linear Analysis , 33:767–807, 2016.
UMERICAL METHODS FOR FRACTIONAL DIFFUSION 33 [36] L. Caffarelli, S. Salsa, and L. Silvestre. Regularity estimates for the solution and the freeboundary of the obstacle problem for the fractional Laplacian.
Invent. Math. , 171(2):425–461, 2008.[37] L. Caffarelli and L. Silvestre. An extension problem related to the fractional Laplacian.
Comm. Part. Diff. Eqs. , 32(7-9):1245–1260, 2007.[38] L. Caffarelli and A. Vasseur. Drift diffusion equations with fractional diffusion and thequasi-geostrophic equation.
Ann. of Math. (2) , 171(3):1903–1930, 2010.[39] 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(8):1353–1384, 2011.[40] B. Carmichael, H. Babahosseini, S.N. Mahmoodi, and M. Agah. The fractional viscoelasticresponse of human breast tissue cells.
Physical Biology , 12(4):046001, 2015.[41] P. Carr, H. Geman, D.B. Madan, and M. Yor. The fine structure of asset returns: Anempirical investigation.
Journal of Business , 75:305–332, 2002.[42] L. Chen, R.H. Nochetto, E. Ot´arola, and A.J. Salgado. A PDE approach to fractionaldiffusion: A posteriori error analysis.
J. Comput. Phys. , 293:339–358, 2015.[43] L. Chen, R.H. Nochetto, E. Ot´arola, and A.J. Salgado. Multilevel methods for nonuniformlyelliptic operators and fractional diffusion.
Math. Comp. , 85(302):2583–2607, 2016.[44] Z.Q. Chen and R. Song. Hardy inequality for censored stable processes.
Tohoku Math. J.(2) , 55(3):439–450, 2003.[45] P. Ciarlet, Jr. Analysis of the Scott-Zhang interpolation in the fractional order Sobolevspaces.
J. Numer. Math. , 21(3):173–180, 2013.[46] R. ˇCiegis, V. Starikoviius, S. Margenov, and R. Kriauzien. Parallel solvers for fractionalpower diffusion problems.
Concurrency and Computation: Practice and Experience , pagese4216–n/a. e4216 cpe.4216.[47] J. Cushman and T. Glinn. Nonlocal dispersion in media with continuously evolving scalesof heterogeneity.
Trans. Porous Media , 13:123–138, 1993.[48] M. Dauge.
Elliptic Boundary Value Problems on Corner Domains . Lecture Notes in Math-ematics, 1341, Springer-Verlag, 1988.[49] M. D’Elia and M. Gunzburger. The fractional Laplacian operator on bounded domains asa special case of the nonlocal diffusion operator.
Comput. Math. Appl. , 66(7):1245 – 1260,2013.[50] E. Di Nezza, G. Palatucci, and E. Valdinoci. Hitchhiker’s guide to the fractional Sobolevspaces.
Bull. Sci. Math. , 136(5):521–573, 2012.[51] S. Dipierro, X. Ros-Oton, and E. Valdinoci. Nonlocal problems with Neumann boundaryconditions.
Rev. Mat. Iberoam. , 33(2):377–416, 2017.[52] J. Duoandikoetxea.
Fourier analysis , volume 29 of
Graduate Studies in Mathematics . Amer-ican Mathematical Society, Providence, RI, 2001. Translated and revised from the 1995Spanish original by David Cruz-Uribe.[53] R.G. Dur´an and A.L. Lombardi. Error estimates on anisotropic Q elements for functionsin weighted Sobolev spaces. Math. Comp. , 74(252):1679–1706 (electronic), 2005.[54] B. Dyda. A fractional order Hardy inequality.
Illinois J. Math. , 48(2):575–588, 2004.[55] B. Dyda, A. Kuznetsov, and M. Kwa´snicki. Eigenvalues of the fractional Laplace operatorin the unit ball.
J. Lond. Math. Soc. , 95(2):500–518, 2017.[56] G. Dziuk. Finite elements for the Beltrami operator on arbitrary surfaces. In
Partial differ-ential equations and calculus of variations , volume 1357 of
Lecture Notes in Math. , pages142–155. Springer, Berlin, 1988.[57] A. Einstein.
Investigations on the theory of the Brownian movement . Dover Publications,Inc., New York, 1956. Edited with notes by R. F¨urth, Translated by A. D. Cowper.[58] E.B. Fabes, C.E. Kenig, and R.P. Serapioni. The local regularity of solutions of degenerateelliptic equations.
Comm. Part. Diff. Eqs. , 7(1):77–116, 1982.[59] B. Faermann. Localization of the Aronszajn-Slobodeckij norm and application to adaptiveboundary element methods. II. The three-dimensional case.
Numer. Math. , 92(3):467–499,2002.[60] R.K. Getoor. First passage times for symmetric stable processes in space.
Trans. Amer.Math. Soc. , 101:75–90, 1961.[61] V. Gol (cid:48) dshtein and A. Ukhlov. Weighted Sobolev spaces and embedding theorems.
Trans.Amer. Math. Soc. , 361(7):3829–3850, 2009.[62] R. Gorenflo, A.A. Kilbas, F. Mainardi, and S.V. Rogosin.
Mittag-Leffler functions, relatedtopics and applications . Springer Monographs in Mathematics. Springer, Heidelberg, 2014.[63] P. Grisvard.
Elliptic problems in nonsmooth domains , volume 69 of
Classics in AppliedMathematics . Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA,2011. Reprint of the 1985 original [MR0775683], With a foreword by Susanne C. Brenner. [64] G. Grubb. Fractional Laplacians on domains, a development of H¨ormander’s theory of µ -transmission pseudodifferential operators. Adv. Math. , 268:478 – 528, 2015.[65] G. Grubb. Spectral results for mixed problems and fractional elliptic operators.
J. Math.Anal. Appl. , 421(2):1616–1634, 2015.[66] L. H¨ormander. Ch. II, Boundary problems for “classical” pseudo–differential operators.Available at , 1965.[67] Y. Huang and A.M. Oberman. Numerical methods for the fractional Laplacian: A finitedifference-quadrature approach.
SIAM J. Numer. Anal. , 52(6):3056–3084, 2014.[68] B. Jin, R. Lazarov, and Z. Zhou. Two fully discrete schemes for fractional diffusion anddiffusion-wave equations with nonsmooth data.
SIAM J. Sci. Comput. , 38(1):A146–A170,2016.[69] F. Jochmann. An H s -regularity result for the gradient of solutions to elliptic equations withmixed boundary conditions. J. Math. Anal. Appl. , 238:429–450, 1999.[70] T. Kato. Note on fractional powers of linear operators.
Proc. Japan Acad. , 36:94–96, 1960.[71] T. Kato. Fractional powers of dissipative operators.
J. Math. Soc. Japan , 13:246–274, 1961.[72] R.B. Kellogg. Interpolation between subspaces of a hilbert space. Technical report, Univ. ofMaryland,, Inst. Fluid Dynamics and App. Math., Tech. Note BN-719, 1971.[73] T. Kilpel¨ainen. Weighted Sobolev spaces and capacity.
Ann. Acad. Sci. Fenn. Ser. AI Math ,19(1):95–113, 1994.[74] M.A. Krasnosel (cid:48) ski˘ı and Ja.B. Ruticki˘ı.
Convex functions and Orlicz spaces . Translated fromthe first Russian edition by Leo F. Boron. P. Noordhoff Ltd., Groningen, 1961.[75] A. Kufner.
Weighted Sobolev spaces . A Wiley-Interscience Publication. John Wiley & SonsInc., New York, 1985. Translated from the Czech.[76] A. Kufner and B. Opic. How to define reasonably weighted Sobolev spaces.
Comment. Math.Univ. Carolin. , 25(3):537–554, 1984.[77] A. Kyprianou, A. Osojnik, and T. Shardlow. Unbiased walk-on-spheres’ Monte Carlo meth-ods for the fractional Laplacian. arXiv:1609.03127 , 2016.[78] N.S. Landkof.
Foundations of modern potential theory . Springer-Verlag, New York-Heidelberg, 1972. Translated from the Russian by A. P. Doohovskoy, Die Grundlehren dermathematischen Wissenschaften, Band 180.[79] N. Laskin. Fractional quantum mechanics and L´evy path integrals.
Physics Letters A ,268(4):298–305, 2000.[80] C. Lubich. Convolution quadrature and discretized operational calculus. I.
Numer. Math. ,52(2):129–145, 1988.[81] A. Lunardi.
Interpolation theory . Appunti. Scuola Normale Superiore di Pisa (Nuova Serie).[Lecture Notes. Scuola Normale Superiore di Pisa (New Series)]. Edizioni della Normale,Pisa, second edition, 2009.[82] J. Lund and K.L. Bowers.
Sinc methods for quadrature and differential equations . Societyfor Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1992.[83] B.M. McCay and M.N.L. Narasimhan. Theory of nonlocal electromagnetic fluids.
Arch.Mech. (Arch. Mech. Stos.) , 33(3):365–384, 1981.[84] A. McIntosh. The square root problem for elliptic operators: a survey. In
Functional-analyticmethods for partial differential equations (Tokyo, 1989) , volume 1450 of
Lecture Notes inMath. , pages 122–140. Springer, Berlin, 1990.[85] W. McLean.
Strongly elliptic systems and boundary integral equations . Cambridge Univer-sity Press, Cambridge, 2000.[86] D. Meidner, J. Pfefferer, K. Sch¨urholz, and B. Vexler. hp - finite elements for fractionaldiffusion. arXiv:1706.04066v1 , 2017.[87] R. Metzler and J. Klafter. The restaurant at the end of the random walk: recent devel-opments in the description of anomalous transport by fractional dynamics. J. Phys. A ,37(31):R161–R208, 2004.[88] P. Morin, R.H. Nochetto, and K.G. Siebert. Local problems on stars: a posteriori errorestimators, convergence, and performance.
Math. Comp. , 72(243):1067–1097 (electronic),2003.[89] B. Muckenhoupt. Weighted norm inequalities for the Hardy maximal function.
Trans. Amer.Math. Soc. , 165:207–226, 1972.[90] R. Musina and A.I. Nazarov. On fractional Laplacians.
Comm. Partial Differential Equa-tions , 39(9):1780–1790, 2014.[91] S. Nazarov and B. Plamenevsky.
Elliptic problems in domains with piecewise smooth bound-aries . De Gruyter expositions in mathematics, De Gruyter, 1994.[92] R.H. Nochetto, E. Ot´arola, and A.J. Salgado. Convergence rates for the classical, thin andfractional elliptic obstacle problems.
Philos. Trans. Roy. Soc. A , 373(2050):20140449, 14,2015.
UMERICAL METHODS FOR FRACTIONAL DIFFUSION 35 [93] R.H. Nochetto, E. Ot´arola, and A.J. Salgado. A PDE approach to fractional diffusion ingeneral domains: a priori error analysis.
Found. Comput. Math. , 15(3):733–791, 2015.[94] R.H. Nochetto, E. Ot´arola, and A.J. Salgado. A PDE approach to numerical fractionaldiffusion. In
Proceedings of the 8th International Congress on Industrial and Applied Math-ematics , pages 211–236. Higher Ed. Press, Beijing, 2015.[95] R.H. Nochetto, E. Ot´arola, and A.J. Salgado. A PDE approach to space-time fractionalparabolic problems.
SIAM J. Numer. Anal. , 54(2):848–873, 2016.[96] R.H. Nochetto, E. Ot´arola, and A.J. Salgado. Piecewise polynomial interpolation in Muck-enhoupt weighted Sobolev spaces and applications.
Numer. Math. , 132(1):85–130, 2016.[97] R.H. Nochetto, K.G. Siebert, and A. Veeser. Theory of adaptive finite element methods: anintroduction. In
Multiscale, nonlinear and adaptive approximation , pages 409–542. Springer,Berlin, 2009.[98] R.H. Nochetto and A. Veeser. Primer of adaptive finite element methods. In
Multiscale andAdaptivity: Modeling, Numerics and Applications, CIME Lectures . Springer, 2011.[99] R.H. Nochetto, T. von Petersdorff, and C.-S. Zhang. A posteriori error analysis for a classof integral equations and variational inequalities.
Numer. Math. , 116(3):519–552, 2010.[100] F.W.J. Olver.
Asymptotics and special functions . AKP Classics. A K Peters, Ltd., Wellesley,MA, 1997. Reprint of the 1974 original [Academic Press, New York; MR0435697 (55
A PDE approach to numerical fractional diffusion . ProQuest LLC, Ann Arbor,MI, 2014. Thesis (Ph.D.)–University of Maryland, College Park.[102] E. Ot´arola. A piecewise linear FEM for an optimal control problem of fractional operators:error analysis on curved domains.
ESAIM Math. Model. Numer. Anal. , 2016. (to appear).[103] E. Ot´arola and A.J. Salgado. Finite element approximation of the parabolic fractional ob-stacle problem.
SIAM J. Numer. Anal. , 54(4):2619–2639, 2016.[104] E. Ot´arola and A.J. Salgado. Regularity of solutions to space–time fractional wave equations:a PDE approach. 2017. In preparation.[105] E. Ot´arola and A.J. Salgado. Sparse optimal control for fractional diffusion.arXiv:1704.01058, 2017.[106] X. Ros-Oton. Nonlocal elliptic equations in bounded domains: a survey.
Publ. Mat. , 60(1):3–26, 2016.[107] X. Ros-Oton and J. Serra. The Dirichlet problem for the fractional Laplacian: regularity upto the boundary.
J. Math. Pures Appl. , 101(3):275 – 302, 2014.[108] X. Ros-Oton and J. Serra. Local integration by parts and Pohozaev identities for higherorder fractional Laplacians.
Discrete Contin. Dyn. Syst. , 35(5):2131–2150, 2015.[109] K. Sakamoto and M. Yamamoto. Initial value/boundary value problems for fractionaldiffusion-wave equations and applications to some inverse problems.
Journal of Mathemat-ical Analysis and Applications , 382(1):426–447, 2011.[110] S.G. Samko, A.A. Kilbas, and O.I. Marichev.
Fractional integrals and derivatives . Gordonand Breach Science Publishers, Yverdon, 1993. Theory and applications, Edited and witha foreword by S. M. Nikol (cid:48) ski˘ı, Translated from the 1987 Russian original, Revised by theauthors.[111] S.A. Sauter and C. Schwab.
Boundary element methods , volume 39 of
Springer Series inComputational Mathematics . Springer-Verlag, Berlin, 2011. Translated and expanded fromthe 2004 German original.[112] L.R. Scott and S. Zhang. Finite element interpolation of nonsmooth functions satisfyingboundary conditions.
Math. Comp. , 54(190):483–493, 1990.[113] R. Servadei and E. Valdinoci. On the spectrum of two different fractional operators.
Proc.Roy. Soc. Edinburgh Sect. A , 144(4):831–855, 2014.[114] L. Silvestre. Regularity of the obstacle problem for a fractional power of the Laplace operator.
Comm. Pure Appl. Math. , 60(1):67–112, 2007.[115] D. Sims, E. Southall, N. Humphries, G. Hays, C. Bradshaw, J. Pitchford, A. James,M. Ahmed, A. Brierley, M. Hindell, D. Morritt, M. Musyl, D. Righton, E. Shepard, V. Wear-mouth, R. Wilson, M. Witt, and J. Metcalfe. Scaling laws of marine predator search be-haviour.
Nature , 451(7182):1098–1102, 2008.[116] J. Sprekels and E. Valdinoci. A new type of identification problems: optimizing the fractionalorder in a nonlocal evolution equation.
SIAM J. Control Optim. , 55(1):70–93, 2017.[117] P.R. Stinga and J.L. Torrea. Extension problem and Harnack’s inequality for some fractionaloperators.
Comm. Partial Differential Equations , 35(11):2092–2122, 2010.[118] L. Tartar.
An introduction to Sobolev spaces and interpolation spaces , volume 3 of
LectureNotes of the Unione Matematica Italiana . Springer, Berlin, 2007.[119] M.E. Taylor. Pseudodifferential operators. Princeton Mathematical Series, vol. 34., 1981.[120] B.O. Turesson.
Nonlinear potential theory and weighted Sobolev spaces , volume 1736 of
Lecture Notes in Mathematics . Springer-Verlag, Berlin, 2000. [121] M.I. Viˇsik and G.I. `Eskin. Elliptic convolution equations in a bounded region and theirapplications.
Uspehi Mat. Nauk , 22(1 (133)):15–76, 1967.[122] K. Yosida.
Functional analysis . Second edition. Die Grundlehren der mathematischen Wis-senschaften, Band 123. Springer-Verlag New York Inc., New York, 1968.(A. Bonito)
Department of Mathematics, Texas A&M University, College Station,TX 77843, USA
E-mail address : [email protected] (J.P. Borthagaray) IMAS - CONICET and Departamento de Matem´atica, FCEyN - Uni-versidad de Buenos Aires, Buenos Aires, Argentina
E-mail address : [email protected] (R.H. Nochetto) Department of Mathematics and Institute for Physical Science andTechnology, University of Maryland, College Park, MD 20742, USA
E-mail address : [email protected] (E. Ot´arola) Departamento de Matem´atica, Universidad T´ecnica Federico Santa Mar´ıa,Valpara´ıso, Chile
E-mail address : [email protected] (A.J. Salgado) Department of Mathematics, University of Tennessee, Knoxville, TN37996, USA
E-mail address ::