Fast Fourier-like Mapped Chebyshev Spectral-Galerkin Methods for PDEs with Integral Fractional Laplacian in Unbounded Domains
Changtao Sheng, Jie Shen, Tao Tang, Li-Lian Wang, Huifang Yuan
FFAST FOURIER-LIKE MAPPED CHEBYSHEV SPECTRAL-GALERKIN METHODSFOR PDES WITH INTEGRAL FRACTIONAL LAPLACIAN IN UNBOUNDEDDOMAINS
CHANGTAO SHENG , JIE SHEN , TAO TANG , LI-LIAN WANG AND HUIFANG YUAN Abstract.
In this paper, we propose a fast spectral-Galerkin method for solving PDEs involving inte-gral fractional Laplacian in R d , which is built upon two essential components: (i) the Dunford-Taylorformulation of the fractional Laplacian; and (ii) Fourier-like bi-orthogonal mapped Chebyshev functions(MCFs) as basis functions. As a result, the fractional Laplacian can be fully diagonalised, and thecomplexity of solving an elliptic fractional PDE is quasi-optimal, i.e., O (( N log N ) d ) with N being thenumber of modes in each spatial direction. Ample numerical tests for various decaying exact solutionsshow that the convergence of the fast solver perfectly matches the order of theoretical error estimates.With a suitable time-discretisation, the fast solver can be directly applied to a large class of nonlinearfractional PDEs. As an example, we solve the fractional nonlinear Schr¨odinger equation by using thefourth-order time-splitting method together with the proposed MCF-spectral-Galerkin method. Introduction
Diffusion is the movement of a substance from an area of high concentration to an area of low concen-tration, which is a ubiquitous physical process in nature. The normal diffusion models rooted in Brownianmotion have been well-studied in years. However, numerous experimental and scientific evidences haveshown that many phenomena and complex systems involve anomalous diffusion, where the underlyingstochastic processes are non-Brownian [39, 37, 38]. Notably, the fractional models have emerged as apowerful tool in modelling anomalous diffusion in diverse fields (see, e.g., [46, 27, 36, 18, 9, 13, 48, 15]and the references therein) over the past two decades. The nonlocal operators typically involved thereininclude the Riemann-Liouville, Caputo and Riesz fractional integrals/derivatives, or the fractional Lapla-cian. They share some common and interwoven difficulties, e.g., the nonlocal and singular behaviours,so they are much more challenging and difficult to deal with than the usual local operators. The recentworks [33, 10] provide an up-to-date review in particular for numerical issues with several versions offractional Laplacian. The interested readers may also refer to [47, 21] for nonlocal modelling in manyother applications.A large volume of literature is available for numerical solutions of one-dimensional spatial and temporalfractional differential equations, which particularly include the finite difference methods/finite element
Mathematics Subject Classification.
Key words and phrases.
Integral fractional Laplacian, Dunford-Taylor formula, Mapped Chebyshev functions, bi-orthogonal, nonlocal/singular operators. Division of Mathematical Sciences, School of Physical and Mathematical Sciences, Nanyang Technological University,637371, Singapore. The research of the authors is partially supported by Singapore MOE AcRF Tier 2 Grants: MOE2018-T2-1-059 and MOE2017-T2-2-144. Emails: [email protected] (C. Sheng) and [email protected] (L. Wang). Department of Mathematics, Purdue University, West Lafayette, IN 47907-1957, USA. The work of of the author ispartially supported by NSF DMS-1620262, DMS-1720442 and AFOSR FA9550-16-1-0102. Email: [email protected] (J.Shen). Division of Science and Technology, BNU-HKBU United International College, Zhuhai, Guangdong, China, andSUSTech International Center for Mathematics, Southern University of Science and Technology, Shenzhen, China. Thework of this author is partially supported by the NSF of China (under the Grant No. 11731006) and the Science ChallengeProject (No. TZ2018001). Email: [email protected] (T. Tang). Department of Mathematics, Southern University of Science and Technology, Shenzhen, China. Email:[email protected] (H. Yuan).The last two authors would like to thank NTU and SUSTech International Center for Mathematics for hosting theirmutual visits devoted to this collaborative work. a r X i v : . [ m a t h . NA ] A ug C. SHENG, J. SHEN, T. TANG, L. WANG & H. YUAN methods (see, e.g., [19, 22, 26, 30, 31, 52] and many references therein), and spectral methods (see,e.g., [16, 28, 34]). In this work, we are mainly interested in the integral fractional Laplaican in multipledimensions, which is deemed as one of the most challenging nonlocal operators for both computation andanalysis. It is known that for s ∈ (0 , , the fractional Laplacian of u ∈ S ( R d ) (the functions of theSchwartz class) is defined by the Fourier transform:( − ∆) s u ( x ) := F − (cid:2) | ξ | s F [ u ]( ξ ) (cid:3) ( x ) , ∀ x ∈ R d . (1.1)Equivalently, it can be defined by the point-wise formula (cf. [40, Prop. 3.3]):( − ∆) s u ( x ) = C d,s p . v . (cid:90) R d u ( x ) − u ( y ) | x − y | d +2 s d y, x ∈ R d , (1.2)where “p.v.” stands for the principle value and the normalisation constant C d,s := (cid:16) (cid:90) R d − cos ξ | ξ | d +2 s d ξ (cid:17) − = 2 s s Γ( s + d/ π d/ Γ(1 − s ) . (1.3)As a result, to evaluate fractional diffusion of u at a spatial point, information involving all spatial pointsis needed. If u is defined on a bounded domain Ω, we first extend it to zero outside Ω , and then use theabove definition.As many physically motivated fractional diffusion models are naturally set in unbounded domains, thedevelopment of effective solution methods has attracted much recent attention. In general, the existingapproaches can be classified into the following two categories. • This first is to approximate the solution by the orthogonal basis functions, and fully use theanalytic properties of fractional Laplacian performing on the basis (see, e.g., [17, 35, 51, 50]).Based on some analytic fractional calculus formulas of generalised Laguerre functions, Chen et al.[17] developed an efficient spectral method for one-dimensional fractional Laplacian on the wholeline. Using the property that the Hermite functions are invariant under the Fourier transform,Mao and Shen [35] proposed the Hermite spectral-Galerkin method in the transformed domainbased on the Fourier definition (1.1). Tang et al. [51] explicitly evaluated the Hermite fractionaldifferentiation matrices and implemented the spectral-collocation methods based on some elegantanalytic tools. The idea was extended to the rational approximation in [50]. It is noteworthy thatdue to the singular and non-separable factor | ξ | s in (1.1), these methods become complicatedeven for d = 2, and computationally prohibitive for d ≥ . • The second is to use suitable equivalent formulations of the fractional Laplacian to alleviateits notorious numerical difficulties. In Caffarelli and Silvestre [14], the d -dimensional fractionalLaplacian is extended to a d + 1 dimensional elliptic operator with degenerating/singular coeffi-cients in the additional dimension. This groundbreaking extension, together with the follow-upworks for the fractional Laplacian in bounded domains, provides a viable alternative for its math-ematical and numerical treatment (see, e.g., [41, 42, 5] for finite element methods). On the otherhand, the variational form corresponding to the fractional Laplacian can be formulated as theDunford-Taylor formula (cf. [11, Thm. 4.1]): for any u, v ∈ H s ( R d ) with s ∈ (0 , (cid:0) ( − ∆) s u, ( − ∆) s v (cid:1) L ( R d ) = C s (cid:90) ∞ t − s (cid:90) R d (cid:0) ( − ∆)( I − t ∆) − u (cid:1) ( x ) v ( x ) d x d t, (1.4)where I is the identity operator and C s = 2 sin( πs ) /π. In particular, for fractional Laplacian in abounded domain Ω ⊆ R d , we have (cid:0) ( − ∆) s ˜ u, ( − ∆) s ˜ v (cid:1) L (Ω) = C s (cid:90) ∞ t − s (cid:90) Ω (cid:0) ( − ∆)( I − t ∆) − ˜ u (cid:1) ( x ) v ( x ) d x d t, (1.5)where ˜ u denotes the zero extension of u. Very recently, the finite element method with sincquadrature (in t ), was implemented and analysed in [11, 12] based on (1.5). For each quadraturenode t j , one solves the elliptic problem: − t j ∆ w j + w j = ˜ u in R d , i . e ., w j = ( I − t j ∆) − ˜ u ( x ) , (1.6) NTEGRAL FRACTIONAL LAPLACIAN 3 where the unbounded domain has to be truncated, and the side of the domain depends on t j . Infact, many sinc quadrature points should be used to resolve the singularity near t = 0 , but theproblem (1.6) becomes stiff and sharp boundary layers at ∂ Ω can occur.We also remark that direct discretization of the integral fractional Laplacian on bounded domainsbased on the definition (1.2), was discussed in some recent works (see, e.g., [29, 23] for finite differencemethods; and [2, 1, 4, 5, 20] for finite element methods).In this paper, we develop a fast spectral-Galerkin method for PDEs involving integral fractional Lapla-cian in R d . Consider, for example, the model equation:( − ∆) s u ( x ) + γu ( x ) = f ( x ) in R d ; u ( x ) = 0 as | x | → ∞ , (1.7)where s ∈ (0 ,
1) and γ > . The efficient spectral algorithm is built upon two essential components: (i) theDunford-Taylor formulation (1.4) for the fractional Laplacian; and (ii) the approximation of the solutionby the tensorial Fourier-like bi-orthogonal mapped Chebyshev functions. As a result, the complexity ofsolving (1.7) is O (( N log N ) d ) , where N is the degree of freedom along each spatial dimension. Theintegration in t (in (1.4)) can be evaluated exactly by using such a formulation and basis, so the maincomputational cost is from the MCF expansions with FFT. In fact, the framework is also applicable toHermite functions, but Hermite approximation is less compelling for at least two reasons (i) the lack ofFFT; and (ii) slow decay of the solution or the source term. As opposite to usual Laplacian, the fractionalLaplacian of a function with typical exponential or algebraic decay will decay algebraically at much slowerrate (see Propositions 4.2-4.3 of this paper). Thus, the MCF approximation is more preferable. Indeed,ample numerical results show that the fast solver for (1.7) has a convergence perfectly in agreement withthe theoretical estimate for various decaying exact solutions tested.The rest of this paper is organized as follows. In section 2, we first introduce the mapped Chebyshevfunctions and generate the Fourier-like bi-orthogonal MCFs in one dimension. In section 3, we describethe fast MCF-spectral-Galerkin method built upon the Dunford-Taylor formulation of the fractionalLaplacian. We conduct the error estimates and provide ample numerical results to show the convergenceorder of the solver is perfectly in agreement with the theoretical prediction in section 4. In the finalsection, we apply the solver to spatial discretisation of the fractional nonlinear Schr¨odinger equation, andalso conclude the paper with some final remarks.2. Fourier-like mapped Chebyshev functions
In this section, we introduce the mapped Chebyshev functions, from which we construct the Fourier-like bi-orthogonal MCFs as one of the important tools for the efficient spectral algorithms to be designedin the forthcoming section.2.1.
Mapped Chebyshev functions.
Let T n ( y ) = cos( n arccos( y )) , y ∈ Λ := ( − ,
1) be the Chebyshevpolynomial of degree n . The Chebyshev polynomials satisfy the three-term recurrence relation T n +1 ( y ) = 2 yT n ( y ) − T n − ( y ) , n ≥ , (2.1)with T ( y ) = 1 and T ( y ) = y. They form a complete orthogonal system in L ω (Λ), namely, (cid:90) Λ T n ( y ) T m ( y ) ω ( y ) d y = πc n δ nm with ω ( y ) = (1 − y ) − , (2.2)where δ nm is the Kronecker symbol, and c = 2 and c n = 1 for n ≥
1. Recall the recurrence formulas (cf.[49]): yT n ( y ) = ( T n +1 ( y ) + T n − ( y )) / , (1 − y ) T (cid:48) n ( y ) = n T n − ( y ) − T n +1 ( y )) . (2.3)We now define the mapped Chebyshev functions (MCFs) as in [25, 43, 45]. Definition 2.1.
Introduce the one-to-one algebraic mapping x = y (cid:112) − y , y = x √ x , x ∈ R , y ∈ Λ , (2.4) C. SHENG, J. SHEN, T. TANG, L. WANG & H. YUAN and define the MCFs as T n ( x ) = 1 (cid:112) c n π/ (cid:112) − y T n ( y ) = 1 (cid:112) c n π/ √ x T n (cid:16) x √ x (cid:17) , (2.5) for x ∈ R and integer n ≥ . Remark 2.1.
To enhance the resolution of MCFs, one can incorporate a scaling parameter ν > . Moreprecisely, using the mapping x = ν y (cid:112) − y , y = x √ ν + x , x ∈ R , y ∈ Λ , the scaled MFCs can be defined as T νn ( x ) := 1 (cid:112) c n π/ ν / √ ν + x T n (cid:16) x √ ν + x (cid:17) = ν − T n ( x/ν ) . In fact, it is straightforward to extend the properties and algorithms from the usual MCFs to the scaledMFCs. For clarity of presentation, we shall not carry the scaling parameter in the algorithm descriptionsand error analysis, but use it in the numerical experiments.
We have the following important properties of the MCFs, which can be shown readily by using thedefinition 2.1 and the properties of Chebyshev polynomials in (2.1)-(2.3) (also see [45, Proposition 2.4]).
Proposition 2.1.
The MCFs are orthonormal in L ( R ) , and we have S mn = S nm = (cid:90) R T (cid:48) n ( x ) T (cid:48) m ( x ) d x = c n (cid:16) (4 c n − − c n − )( n −
16 + (4 c n +1 − c n +2 )( n + 1) − c n (cid:17) , if m = n, √ c n c n +2 (cid:16) ( c n − c n +2 )( n + 1)8 − c n +1 ( n + 1) (cid:17) , if m = n + 2 , √ c n c n +4 (cid:16) c n +2 ( n + 1)( n + 3)16 (cid:17) , if m = n + 4 . (2.6)2.2. Fourier-like bi-orthogonal MCFs in one dimension.
Let P N be the set of all polynomials ofdegree at most N , and define the finite dimensional space V N := (cid:8) φ : φ ( x ) = g ( x ) ϕ ( x ) , ∀ ϕ ∈ P N (cid:9) , (2.7)where x, y are associated with the mapping (2.4) and g ( x ) := (cid:114) ω ( y ) d y d x = 1 √ x = (cid:112) − y := G ( y ) . (2.8)Note that we have V N := span (cid:8) T n ( x ) : 0 ≤ n ≤ N (cid:9) . (2.9)Following the spirit of [44], we next introduce a Fourier-like basis of V N , which is orthogonal in both L - and H -inner products. For this purpose, let S be a square matrix of order N + 1 with entries givenby (2.6), and let I the identity matrix of the same size. Note from (2.6) that S is a symmetric positivedefinite matrix with nine nonzero diagonals. Thus, all the eigenvalues are real and eigenvectors areorthonormal. To this end, let E = ( e jk ) j,k =0 , ··· ,N be the matrix formed by the orthonormal eigenvectorsof S , and Σ = diag { λ k } be the diagonal matrix of the corresponding eigenvalues. Thus, we have SE = E Σ , E t E = I . (2.10)We remark that with an even and odd separation, we can work with two symmetric positive definiteseven-diagonal sub-matrices to compute the eigenvalues and eigenvectors of S , which should be morestable for large N. NTEGRAL FRACTIONAL LAPLACIAN 5
Lemma 2.1.
Let E = ( e , e , · · · , e N ) be the matrix of the eigenvectors of S , i.e., Se p = λ p e p for ≤ p ≤ N. Define (cid:98) T p ( x ) := N (cid:88) j =0 e jp T j ( x ) , e p = ( e p , e p , · · · , e Np ) t , ≤ p ≤ N. (2.11) Then { (cid:98) T p } Np =0 form an equivalent basis of V N , and they are bi-orthogonal in the sense that ( (cid:98) T p , (cid:98) T q ) L ( R ) = δ pq , (cid:0) (cid:98) T (cid:48) p , (cid:98) T (cid:48) q (cid:1) L ( R ) = λ p δ pq , ≤ p, q ≤ N. (2.12) Proof.
In view of the definition (2.11), we infer from the orthogonality of MCFs and (2.10) that( (cid:98) T p , (cid:98) T q ) L ( R ) = N (cid:88) j =0 N (cid:88) k =0 e kp e jq ( T k , T j ) L ( R ) = N (cid:88) j =0 N (cid:88) k =0 e jq δ jk e kp = N (cid:88) k =0 e kq e kp = e tq e p = δ pq . Similarly, we can show that (cid:0)(cid:98) T (cid:48) p , (cid:98) T (cid:48) q (cid:1) L ( R ) = N (cid:88) j =0 N (cid:88) k =0 e kp e jq ( T (cid:48) k , T (cid:48) j ) L ( R ) = N (cid:88) j =0 N (cid:88) k =0 e jq S jk e kp = ( E t SE ) pq = ( Σ ) pq = λ p δ pq . This ends the proof. (cid:3) MCF-spectral-Galerkin method based on Dunford-Taylor formulation
In this section, we describe the fast MCF spectral-Galerkin algorithm for a model elliptic problemwith integral fractional Laplacian. We then apply the solver for spatial discretisation of some nonlinearfractional PDEs in the next section.3.1.
Some notation.
Denote by S ( R d ) the functions of the Schwartz class, and let S (cid:48) ( R d ) be thetopological dual of S ( R d ) . For any u ∈ S ( R d ) , its Fourier transform is given by F [ u ]( ξ ) = 1(2 π ) d (cid:90) R d u ( x ) e − i ξ · x d x, ∀ ξ ∈ R d . For real s ≥ , we define the fractional Sobolev space (cf. [40, P. 530]): H s ( R d ) = (cid:110) u ∈ L ( R d ) : (cid:107) u (cid:107) H s ( R d ) = (cid:90) R d (1 + | ξ | s ) (cid:12)(cid:12) F [ u ]( ξ ) (cid:12)(cid:12) d ξ < + ∞ (cid:111) , (3.1)and an analogous definition for the case s < H s ( R d ) = (cid:110) u ∈ S (cid:48) ( R d ) : (cid:107) u (cid:107) H s ( R d ) = (cid:90) R d (1 + | ξ | ) s (cid:12)(cid:12) F [ u ]( ξ ) (cid:12)(cid:12) d ξ < + ∞ (cid:111) , (3.2)although in this case the space H s ( R d ) is not a subset of L ( R d ) . According to [40, Prop. 3.4], we know that for s ∈ (0 , , the space H s ( R d ) can also be characterisedby the fractional Laplacian defined in (1.2), equipped with the norm (cid:107) u (cid:107) H s ( R d ) = (cid:0) (cid:107) u (cid:107) L ( R d ) + [ u ] H s ( R d ) (cid:1) , where [ u ] H s ( R d ) is so-called Gagliardo (semi)norm of u, given by[ u ] H s ( R d ) = (cid:16) (cid:90) R d (cid:90) R d | u ( x ) − u ( y ) | | x − y | d +2 s d x d y (cid:17) . (3.3)Indeed, by [40, Prop. 3.6], we have that for s ∈ (0 , , [ u ] H s ( R d ) = 2 C − d,s (cid:107) ( − ∆) s/ u (cid:107) L ( R d ) . (3.4)We have the following important space interpolation property (cf. [3, Ch. 1]), which will be used forthe error analysis later on. C. SHENG, J. SHEN, T. TANG, L. WANG & H. YUAN
Lemma 3.1.
For real r , r ≥ , let r = (1 − θ ) r + θr with θ ∈ [0 , . Then for any u ∈ H r ( R d ) ∩ H r ( R d ) , we have (cid:107) u (cid:107) H r ( R d ) ≤ (cid:107) u (cid:107) − θH r ( R d ) (cid:107) u (cid:107) θH r ( R d ) . (3.5) In particular, for s ∈ [0 , , we have (cid:107) u (cid:107) H s ( R d ) ≤ (cid:107) u (cid:107) − sL ( R d ) (cid:107) u (cid:107) sH ( R d ) . (3.6)3.2. Dunford-Taylor formulation of the fractional Laplacian.
To fix the idea, we consider( − ∆) s u ( x ) + γu ( x ) = f ( x ) in R d ; u ( x ) = 0 as | x | → ∞ , (3.7)where s ∈ (0 , , γ > f ∈ H − s ( R d ) . A weak form for (3.7) is to find u ∈ H s ( R d ) such that B ( u, v ) = (cid:0) ( − ∆) s/ u, ( − ∆) s/ v (cid:1) L ( R d ) + γ ( u, v ) L ( R d ) = [ u, v ] H s ( R d ) + γ ( u, v ) L ( R d ) = ( f, v ) L ( R d ) , ∀ v ∈ H s ( R d ) , (3.8)where [ u, v ] H s ( R d ) induces the Gagliardo (semi)norm in (3.3).By the definitions (1.1) and (3.1), we immediately obtain the continuity and coercivity of the bilinearform a ( · , · ), that is, for any u, v ∈ H s ( R d ) , |B ( u, v ) | (cid:46) (cid:107) u (cid:107) H s ( R d ) (cid:107) v (cid:107) H s ( R d ) , |B ( u, u ) | (cid:38) (cid:107) u (cid:107) H s ( R d ) . (3.9)Then, we derive from the Lax-Milgram lemma (cf. [6]) that the problem (3.8) admits a unique solutionsatisfying (cid:107) u (cid:107) H s ( R d ) (cid:46) (cid:107) f (cid:107) H − s ( R d ) . In view of the definitions in (1.1)-(1.2), we have the equivalent forms of [ u, v ] H s ( R d ) as follows[ u, v ] H s ( R d ) = (cid:90) R d (cid:90) R d ( u ( x ) − u ( y ))( v ( x ) − v ( y )) | x − y | d +2 s d x d y (3.10)= (cid:90) R d | ξ | s F [ u ]( ξ ) F [ v ]( ξ ) d ξ. (3.11)It is noteworthy that the direct implementation of a numerical scheme based on (3.10) (i.e., in physicalspace) is very difficult. Most of the existing works (see, e.g., [35, 51, 50]) are therefore mainly basedon (3.11) (i.e., in the frequency space). The Hermite function approaches can take the advantage thatthe Fourier transforms of Hermite functions are explicitly known. However, in multiple dimensions,the non-separable/singular factor | ξ | s = ( ξ + · · · + ξ d ) s makes the tensorial approach computationallyprohibitive. On the other hand, the fractional Laplacian operator may become rather complicated whena coordinate transform is applied, so the mapped Chebyshev approximation can not be applied in eitherof the above formulations.In what follows, we resort to an alternative formulation of the fractional Laplacian that can over-come these numerical difficulties. According to [11, Theorem 4.1], we have the following Dunford-Taylorformulation of the integral fractional Laplacian. Lemma 3.2.
For any u, v ∈ H s ( R d ) with s ∈ (0 , , we have (cid:0) ( − ∆) s u, ( − ∆) s v (cid:1) L ( R d ) = C s (cid:90) ∞ t − s (cid:90) R d (cid:0) ( − ∆)( I − t ∆) − u (cid:1) ( x ) v ( x ) d x d t, (3.12) where I is the identity operator and C s = 2 sin( πs ) π . (3.13)Let us denote w = w ( u, t ) := ( I − t ∆) − u ( x ) . Then there holds − t ∆ w + w = u in R d , so ( − ∆)( I − t ∆) − u = − ∆ w = t − ( u − w ) . (3.14) NTEGRAL FRACTIONAL LAPLACIAN 7
As a result, we can rewrite the weak form (3.8) as: find u ∈ H s ( R d ) such that B ( u, v ) = C s (cid:90) ∞ t − − s ( u − w, v ) L ( R d ) d t + γ ( u, v ) L ( R d ) = ( f, v ) L ( R d ) , ∀ v ∈ H s ( R d ) , (3.15)where w = w ( u, t ) solves t ( ∇ w, ∇ ψ ) L ( R d ) + ( w, ψ ) L ( R d ) = ( u, ψ ) L ( R d ) , ∀ ψ ∈ H ( R d ) . (3.16)It is evident that the wellposedness of (3.15)-(3.16) follows from its equivalence to (3.8).3.3. The MCF spectral-Galerkin scheme and its implementation.
Define V dN = V N ⊗ · · · ⊗ V N , (3.17)which is the tensor product of d copies of V N defined in (2.7). Here, V N = V N , and denote I dN : C ( R d ) → V dN the tensorial mapped Chebyshev interpolation operator. The MCF spectral-Galerkin approximationto (3.15)-(3.16) is to find u N ∈ V dN such that B N ( u N , v N ) = C s (cid:90) ∞ t − − s ( u N − w N , v N ) L ( R d ) d t + γ ( u N , v N ) L ( R d ) = ( I dN f, v N ) L ( R d ) , ∀ v N ∈ V dN , (3.18)where we find w N := w N ( u N , t ) ∈ V dN such that for any t > ,t ( ∇ w N , ∇ ψ ) L ( R d ) + ( w N , ψ ) L ( R d ) = ( u N , ψ ) L ( R d ) , ∀ ψ ∈ V dN . (3.19)Define the d -dimensional tensorial Fourier-like basis and denote the vector of the corresponding eigen-values in (2.10) by (cid:98) T n ( x ) = d (cid:89) j =1 (cid:98) T n j ( x j ) , x ∈ R d ; λ n = ( λ n , · · · , λ n d ) t . (3.20)Accordingly, we have V dN = span (cid:8)(cid:98) T n ( x ) , n ∈ Υ N (cid:9) , (3.21)where the index set Υ N := (cid:8) n = ( n , · · · , n d ) : 0 ≤ n j ≤ N, ≤ j ≤ d (cid:9) . (3.22)As an extension of (2.12), we have the following attractive property of the tensorial Fourier-like MCFs. Theorem 3.1.
For the tensorial Fourier-like MCFs, we have ( (cid:98) T p , (cid:98) T q ) L ( R d ) = δ pq ; (cid:0) ∇ (cid:98) T p , ∇ (cid:98) T q (cid:1) L ( R d ) = | λ p | δ pq , (3.23) where p, q ∈ Υ N and δ pq = d (cid:89) j =0 δ p j q j , | λ p | = λ p + · · · + λ p d . (3.24) Proof.
One verifies by using the orthogonality (2.12) and the definition (3.20) that( (cid:98) T p , (cid:98) T q ) L ( R d ) = ( (cid:98) T p , (cid:98) T q ) L ( R ) · · · ( (cid:98) T p d , (cid:98) T q d ) L ( R ) = δ p q · · · δ p d q d = δ pq , and ( ∇ (cid:98) T p , ∇ (cid:98) T q ) L ( R d ) = (cid:8) ( (cid:98) T (cid:48) p , (cid:98) T (cid:48) q ) L ( R ) ( (cid:98) T p , (cid:98) T q ) L ( R ) · · · ( (cid:98) T p d , (cid:98) T q d ) L ( R ) (cid:9) + (cid:8) ( (cid:98) T p , (cid:98) T q ) L ( R ) ( (cid:98) T (cid:48) p , (cid:98) T (cid:48) q ) L ( R ) · · · ( (cid:98) T p d , (cid:98) T q d ) L ( R ) (cid:9) + · · · + (cid:8) ( (cid:98) T p , (cid:98) T q ) L ( R ) · · · ( (cid:98) T p d − , (cid:98) T q d − ) L ( R ) ( (cid:98) T (cid:48) p d , (cid:98) T (cid:48) q d ) L ( R ) (cid:9) = λ p δ p q · · · δ p d q d + λ p δ p q · · · δ p d q d + · · · + λ p d δ p q · · · δ p d q d = ( λ p + · · · + λ p d ) δ pq = | λ p | δ pq . This ends the proof. (cid:3)
C. SHENG, J. SHEN, T. TANG, L. WANG & H. YUAN
Remarkably, the use of the Fourier-like MCF can diagonalise the integral fractional Laplacian in theDunford-Taylor formulation.
Theorem 3.2.
Using the tensorial Fourier-like MCFs as basis functions, the solution of (3.18) - (3.19) can be uniquely expressed as u N ( x ) = (cid:88) p ∈ Υ N ˆ f p γ + | λ p | s (cid:98) T p ( x ) , x ∈ R d , (3.25) where (cid:98) T p , λ p are defined in (3.20) , and ˆ f p = ( I dN f, (cid:98) T p ) L ( R d ) , p ∈ Υ N . (3.26) Proof.
Write u N = (cid:88) p ∈ Υ N ˆ u p (cid:98) T p ( x ) , w N = (cid:88) p ∈ Υ N ˆ w p (cid:98) T p ( x ) , (3.27)where w N is the unique solution of (3.19) associated with u N . For clarity, we split the proof into thefollowing steps.(i). We first show that w N can be uniquely determined by u N via w N = (cid:88) p ∈ Υ N ˆ u p t | λ p | (cid:98) T p ( x ) . (3.28)Substituting (3.27) into (3.19), and taking ψ = (cid:98) T q in (3.19), we arrive at (cid:88) p ∈ Υ N ˆ w p (cid:8) t ( ∇ (cid:98) T p , ∇ (cid:98) T q ) L ( R d ) + ( (cid:98) T p , (cid:98) T q ) L ( R d ) (cid:9) = (cid:88) p ∈ Υ N ˆ u p ( (cid:98) T p , (cid:98) T q ) L ( R d ) , ∀ q ∈ Υ N . By the orthogonality (3.23), we obtain (cid:88) p ∈ Υ N ˆ w p (cid:8) t | λ p | + 1 (cid:9) δ pq = (cid:88) p ∈ Υ N ˆ u p δ pq , ∀ q ∈ Υ N , which implies (3.28), as ˆ w p = ˆ u p t | λ p | , ∀ p ∈ Υ N . (ii). We next prove the integral identity: (cid:90) ∞ t − s | λ p | − s t | λ p | d t = π πs ) = 1 C s . (3.29)Indeed, using a change of variable y = t (cid:112) | λ p | , we find readily that (cid:90) ∞ t − s | λ p | − s t | λ p | d t = (cid:90) ∞ y − s y d y = π πs ) , where we used the known formula (3.30) below with µ = 2 − s and ν = 2. According to [24, P. 325, P.918, P.905], we have for ν ≥ µ ≥ ν (cid:54) = 0 , (cid:90) ∞ x µ − x ν d x = 1 ν B (cid:16) µν , − µν (cid:17) = 1 ν Γ (cid:16) µν (cid:17) Γ (cid:16) − µν (cid:17) = πν sin( πµ/ν ) , (3.30)where we used the properties of the Beta and Gamma functions:B( x, y ) = Γ( x )Γ( y )Γ( x + y ) , Γ(1 − x )Γ( x ) = π sin( πx ) . (iii). Finally, we can derive (3.25) with the aid of (3.28)-(3.29). It is evident that by (3.27)-(3.28), u N − w N = t (cid:88) p ∈ Υ N | λ p | t | λ p | ˆ u p (cid:98) T p ( x ) . (3.31) NTEGRAL FRACTIONAL LAPLACIAN 9
Thus, substituting (3.31) into (3.18) with v N = (cid:98) T q , we obtain from (3.23) and (3.29) that B N ( u N , (cid:98) T q ) = (cid:88) p ∈ Υ N ˆ u p (cid:26) C s (cid:90) ∞ t − s (cid:90) R d | λ p | t | λ p | (cid:98) T p ( x ) (cid:98) T q ( x )d x d t + γ ( (cid:98) T p , (cid:98) T q ) L ( R d ) (cid:27) = (cid:88) p ∈ Υ N ˆ u p (cid:26) δ pq C s (cid:90) ∞ t − s | λ p | t | λ p | d t + γδ pq (cid:27) = (cid:88) p ∈ Υ N ˆ u p (cid:0) | λ p | s + γ (cid:1) δ pq = ( I dN f, (cid:98) T q ) L ( R d ) , which implies ˆ u p = ( I dN f, (cid:98) T p ) L ( R d ) γ + | λ p | s , ∀ p ∈ Υ N . Thus, we obtain (3.25)-(3.26) immediately. (cid:3)
Remark 3.1.
It is crucial to use the Fourier-like MCFs as the basis functions for both u N and w N , sothat we can take the advantage of the bi-orthogonality and explicitly evaluate the integration in t. In otherwords, under the Fourier-like basis, the stiffness matrix of the linear system of (3.18) - (3.19) becomes adiagonal matrix of the form (cid:98) S := (cid:0) Σ ⊗ I ⊗ · · · ⊗ I + I ⊗ Σ ⊗ I ⊗ · · · ⊗ I + · · · + I ⊗ · · · ⊗ I ⊗ Σ (cid:1) s , (3.32) where Σ is defined in (2.10) and ⊗ denotes the tensor product operator as before. Remark 3.2.
The main cost of solving the system (3.3) is devoted to the evaluation of the right-handside, which can be carried out by the fast Fourier transform (FFT) related to Chebyshev polynomials. Error estimates and numerical examples
In this section, we derive some relevant MCF approximation results, which are useful for the errorestimates of the proposed MCF spectral-Galerkin scheme.4.1.
Approximation by MCFs.
We first consider d -dimensional L -orthogonal projection: π dN : L ( R d ) → V dN such that (cid:90) R d ( π dN u − u )( x ) v ( x ) d x = 0 , ∀ v ∈ V dN . (4.1)We intend to estimate the projection error in the fractional Sobolev norm, i.e., (cid:107) π dN u − u (cid:107) H s ( R d ) . Forthis purpose, we introduce some notation and spaces of functions.For notational convenience, the pairs of functions ( u, ˘ u ) and ( U, ˘ U ) associated with the mapping (2.4)have the relations u ( x ) = U ( y ( x )) , ˘ u ( x ) = u ( x ) g ( x ) = U ( y ) G ( y ) = ˘ U ( y ) , (4.2)where as in (2.8), g ( x ) = (1 + x ) − / = (cid:112) − y = G ( y ) . Define the differential operators D x j u := ∂ x j (cid:8) (1 + x j ) u (cid:9) d x j d y j = a ( x j ) ∂ x j (cid:8) (1 + x j ) u (cid:9) = ∂ y j ˘ U ,D k j x j u = a ( x j ) ∂ x j (cid:110) a ( x j ) ∂ x j (cid:110) · · · (cid:110) a ( x j ) ∂ x j (cid:110) (1 + x j ) u (cid:111)(cid:111) · · · (cid:111)(cid:111) = ∂ k j y j ˘ U , (4.3)for k j ≥ , ≤ j ≤ d, where a ( x j ) = dx j /dy j = (1 + x j ) . Correspondingly, we define the d -dimensionalSobolev space B m ( R d ) = (cid:8) u : D kx u ∈ L (cid:36) k +1 ( R d ) , ≤ | k | ≤ m (cid:9) , m = 0 , , , · · · , (4.4) where the differential operator and the weight function are D kx u = D k x · · · D k d x d u, (cid:36) k ( x ) = d (cid:89) j =1 (1 + x j ) − k j . (4.5)It is equipped with the norm and semi-norm (cid:107) u (cid:107) B m ( R d ) = (cid:16) (cid:88) ≤| k | ≤ m (cid:13)(cid:13) D kx u (cid:13)(cid:13) L (cid:36) k ( R d ) (cid:17) , | u | B m ( R d ) = (cid:16) d (cid:88) j =1 (cid:13)(cid:13) D mx j u (cid:13)(cid:13) L (cid:36) mej ( R d ) (cid:17) , (4.6)where e j = (0 , · · · , , · · · ,
0) be the j -th unit vector in R d . Theorem 4.1. If u ∈ B m ( R d ) with integer m ≥ , then we have (cid:107) π dN u − u (cid:107) H s ( R d ) ≤ cN s − m | u | B m ( R d ) , ≤ s ≤ , (4.7) where c is a positive constant independent of N and u. Proof.
In view of (3.6), it is necessary to estimate the projection errors in the L - and H -norms.According to [45, Theorems 3.1-3.2], we have (cid:107) π dN u − u (cid:107) L ( R d ) ≤ cN − m | u | B m ( R d ) , (4.8)and (cid:107)∇ ( π dN u − u ) (cid:107) L ( R d ) ≤ cN − m | u | B m ( R d ) . (4.9)Using Lemma 3.1 and (4.8)-(4.9), we arrive at (cid:107) π dN u − u (cid:107) H s ( R d ) ≤ cN s − m | u | B m ( R d ) , s ∈ [0 , . This ends the proof. (cid:3)
We now turn to the error estimate for the interpolation operator. Let { y j , ρ j } Nj =0 be the Chebyshev-Gauss quadrature nodes and weights on Λ = ( − , x j = y j (cid:113) − y j , ω j = ρ j y j , ≤ j ≤ N. (4.10)Then by the exactness of the Chebyshev-Gauss quadrature, we have (cid:90) R u ( x ) v ( x ) d x = (cid:90) Λ U ( y ) V ( y )(1 − y ) − d y = (cid:90) Λ ˘ U ( y ) ˘ V ( y )(1 − y ) − d y = N (cid:88) j =0 ˘ U ( y j ) ˘ V ( y j ) ρ j , ∀ ˘ U · ˘ V ∈ P N +1 . (4.11)which, together with (2.7), implies the exactness of quadrature (cid:90) R u ( x ) v ( x ) d x = N (cid:88) j =0 u ( x j ) v ( x j ) ω j , ∀ u · v ∈ V N +1 . (4.12)We now introduce the one-dimensional interpolation operator I N : C ( R ) → V N such that I N u ( x j ) = u ( x j ) , ≤ j ≤ N. (4.13)With a little abuse of notation, we define the d -dimensional grids by x j = ( x j , · · · , x j d ) , j ∈ Υ N , where { x j k } dk =1 are the mapped Chebyshev-Gauss nodes, and the index set Υ N is given in (3.22) as before. Wenow consider the d -dimensional MCF interpolation: C ( R d ) → V dN ,I dN u ( x j ) = I (1) N ◦ · · · ◦ I ( d ) N u ( x j ) , j ∈ Υ N , (4.14)where I ( k ) N = I N , ≤ k ≤ d, is the interpolation along x k -direction. NTEGRAL FRACTIONAL LAPLACIAN 11
In the error analysis, we also need the L -estimate of the d -dimensional MCF interpolation. For betterdescription of the error, we introduce a second semi-norm of B m ( R d ) for m ≥ u ]] B m ( R d ) := (cid:26) | u | B m ( R d ) + d (cid:88) j =1 (cid:88) k (cid:54) = j (cid:13)(cid:13) D m − x k D x j u (cid:13)(cid:13) L d(cid:36) m − ek + ej ( R d ) (cid:27) , (4.15)where the weight function (cid:36) and | u | B m ( R d ) are defined in (4.5) and (4.6) as before.We have the following L -estimates of the interpolation. Theorem 4.2.
For u ∈ B m ( R d ) with the integer m ≥ , we have (cid:107) I dN u − u (cid:107) L ( R d ) ≤ cN − m [[ u ]] B m ( R d ) , (4.16) where c is a positive constant independent of N and u. Proof.
Let I CN : C (Λ) → P N be the Chebyshev-Gauss interpolation operator. According to [43, Lemma3.6], we have that for any v ∈ L ω (Λ) and v (cid:48) ∈ L ω − (Λ) with ω ( y ) = (1 − y ) − , (cid:107) I CN v (cid:107) L ω (Λ) ≤ c (cid:0) (cid:107) v (cid:107) L ω (Λ) + N − (cid:107) (1 − y ) v (cid:48) (cid:107) L ω (Λ) (cid:1) , (4.17)where c is a positive constant independent of N and v. Moreover, by [43, Theorem 3.41], we have theone-dimensional Chebyshev-Gauss interpolation error estimates, (cid:107) ( I CN v − v ) (cid:48) (cid:107) L ω − (Λ) + N (cid:107) I CN v − v (cid:107) L ω (Λ) ≤ cN − m (cid:107) (1 − y ) m v ( m ) (cid:107) L ω (Λ) . (4.18)In view of (2.4), we have (cid:107) I dN u − u (cid:107) L ( R d ) = (cid:13)(cid:13) I C,dN ( U/G ) − ( U/G ) (cid:13)(cid:13) L ω (Λ d ) = (cid:13)(cid:13) I C,dN ˘ U − ˘ U (cid:13)(cid:13) L ω (Λ d ) , (4.19)where G ( y ) = (cid:81) dj =1 G ( y j ) and I C,dN := I C, (1) N ◦ · · · ◦ I C, ( d ) N d with I C, ( k ) N = I CN . For clarity, we only prove the results with d = 2 , as it is straightforward to extend the results to thecase with d ≥ . By virtue of the triangle inequality, (4.17) and (4.18), we obtain that for m ≥ , (cid:107) I C, N ˘ U − ˘ U (cid:107) L ω (Λ ) ≤ (cid:107) I C, (1) N ˘ U − ˘ U (cid:107) L ω (Λ ) + (cid:13)(cid:13) I C, (1) N ◦ (cid:0) I C, (2) N ˘ U − ˘ U (cid:1)(cid:13)(cid:13) L ω (Λ ) ≤ cN − m (cid:107) (1 − y ) m ∂ my ˘ U (cid:107) L ω (Λ ) + c (cid:110) (cid:107) I C, (2) N ˘ U − ˘ U (cid:107) L ω (Λ ) + N − (cid:13)(cid:13) (1 − y ) ∂ y (cid:0) I C, (2) N ˘ U − ˘ U (cid:1)(cid:13)(cid:13) L ω (Λ ) (cid:111) ≤ cN − m (cid:110)(cid:13)(cid:13) (1 − y ) m ∂ my ˘ U (cid:13)(cid:13) L ω (Λ ) + (cid:13)(cid:13) (1 − y ) m ∂ my ˘ U (cid:13)(cid:13) L ω (Λ ) + (cid:13)(cid:13) (1 − y ) (1 − y ) m − ∂ y ∂ m − y ˘ U (cid:13)(cid:13) L ω (Λ ) (cid:111) . (4.20)Note that in the above derivation, we can switch the order of I C, (1) N and I C, (2) N , so we can add the term (cid:13)(cid:13) (1 − y ) m − (1 − y ) ∂ m − y ∂ y ˘ U (cid:13)(cid:13) L ω (Λ ) in the upper bound. Then, we derive from (4.2), (4.3) and (4.20)that for m ≥ , (cid:107) I C, N ˘ U − ˘ U (cid:107) L ω (Λ ) ≤ cN − m (cid:110)(cid:13)(cid:13) (1 + x ) − m +12 (1 + x ) − D mx u (cid:13)(cid:13) L ( R ) + (cid:13)(cid:13) (1 + x ) − (1 + x ) − m +12 D mx u (cid:13)(cid:13) L ( R ) + (cid:13)(cid:13) (1 + x ) − (1 + x ) − m D x D m − x u (cid:13)(cid:13) L ( R ) + (cid:13)(cid:13) (1 + x ) − m (1 + x ) − D m − x D x u (cid:13)(cid:13) L ( R ) (cid:111) ≤ cN − m [[ u ]] B m ( R ) . It is straightforward to extend the above derivation to d ≥ . This completes the proof. (cid:3)
To conduct error analysis for the proposed MCF scheme, we assume that the error for solving theelliptic problem (3.19) is negligible (or equivalently, the quadrature errors in evaluating the fractionalLaplacian in the scheme can be ignored), so formally, we have w N = ( I − t ∆) − u N . It is noteworthy that the analysis of such an error is feasible for the finite element approximation of the fractional Laplacianin bounded domain based on the Dunford-Taylor formulation in a bounded domain, though the proofis lengthy and much involved (see [11]). However, the analysis is largely open in this situation, mainlybecause the spectrum estimate of the fractional Laplacian operator in R d appears unavailable, as oppositeto the bounded domain case (see [11]). Proposition 4.1.
Assume that the elliptic problem (3.19) in the scheme (3.18) can be solved exactly.Then we have the estimate: for u ∈ B m ( R d ) with integer m ≥ , and f ∈ B m ( R d ) with integer m ≥ , (cid:107) u − u N (cid:107) H s ( R d ) ≤ cN s − m | u | B m ( R d ) + cN − m [[ f ]] B m ( R d ) , s ∈ (0 , , (4.21) where c is a positive constant independent of u, f and N. Proof.
Under this assumption, we find from Lemma 3.2 that the scheme (3.18) can be written as: find u N ∈ V dN such that B ( u N , v N ) = ( I dN f, v N ) , ∀ v N ∈ V dN . Then by (3.8) and (3.9), we infer from a standard argument that (cid:107) u − u N (cid:107) H s ( R d ) ≤ c ( (cid:107) π dN u − u (cid:107) H s ( R d ) + (cid:107) I dN f − f (cid:107) L ( R d ) ) . Thus, the estimate (4.21) follows from Theorems 4.1 and 4.2 immediately. (cid:3)
Remark 4.1.
In what follows, we shall validate the above assumption through several numerical tests.Indeed, we shall observe that the errors of solving (3.19) are insignificant, and the order of the numericalerrors perfectly agrees with the estimated order.
Useful analytic formulas.
We first derive analytical formulas for fractional Laplacian of somefunctions with typical exponential or algebraic decay, upon which we construct the exact solutions to testthe accuracy of the proposed method, and to validate the assumption in Proposition 4.1. Moreover, wereveal that the fractional Laplacian has a very different property from the usual Laplacian. For example,the image of an exponential function decays algebraically (see Proposition 4.2 below), as opposite to theusual one.We have the following exact formulas for the Gaussian function and rational function, whose derivationsare sketched in Appendix A and Appendix B, respectively.
Proposition 4.2.
For real s > and integer d ≥ , we have that ( − ∆) s (cid:8) e −| x | (cid:9) = 2 s Γ( s + d/ d/ F (cid:16) s + d d −| x | (cid:17) . (4.22) Moreover for non-integer s > and | x | → ∞ , we have the asymptotic behaviour ( − ∆) s (cid:8) e −| x | (cid:9) = − s sin( πs ) π Γ( s + d/ s ) | x | d +2 s (cid:8) O ( | x | − ) (cid:9) . (4.23) Proposition 4.3.
For real s, r > , and integer d ≥ , we have ( − ∆) s (cid:110) | x | ) r (cid:111) = 2 s Γ( s + γ )Γ( s + d/ γ )Γ( d/ F (cid:16) s + r, s + d d −| x | (cid:17) . (4.24) Moreover for non-integer s > and | x | → ∞ , we have the asymptotic properties: for r (cid:54) = d/ , ( − ∆) s (cid:110) | x | ) r (cid:111) ∼ | x | ) s + µ , (4.25) where µ = min { r, d/ } . if r = d/ , then ( − ∆) s (cid:110) | x | ) r (cid:111) ∼ ln(1 + | x | )(1 + | x | ) s + d/ ; (4.26) NTEGRAL FRACTIONAL LAPLACIAN 13
Numerical results.
We apply the MCF spectral-Galerkin method to solve the model problem (3.7)in various situations.
Example 4.1. (Accuracy test).
We first consider (3.7) with the following exact solutions: u e ( x ) = e −| x | , u a ( x ) = (1 + | x | ) − r , r > , x ∈ R d . (4.27)In view of (4.22) and (4.24), the source terms f e ( x ) and f a ( x ) are respectively given by f e ( x ) = γe −| x | + 2 s Γ( s + d/ d/ F (cid:16) s + d d −| x | (cid:17) ,f a ( x ) = γ (1 + | x | ) − r + 2 s Γ( s + r )Γ( s + d/ r )Γ( d/ F (cid:16) s + r, s + d d −| x | (cid:17) . (4.28)Now, we intend to use the error estimates in Propositions 4.2 and 4.3 to analytically calculate theexpected order of convergence by the MCF scheme, and then verify the convergence order numerically.For this purpose, we consider a generic function of algebraic decay as follows w ( x ) = 1(1 + | x | ) µ , x ∈ R d , µ > , (4.29)Using the definitions (4.3) and (4.15), we obtain from direct calculation that D x j w ( x ) = (1 + x j ) ∂ x j (cid:8) (1 + x j ) w ( x ) (cid:9) = (1 + x j )(1 + | x | ) − µ − ( − µx j (1 + x j )+ x j (1 + | x | )) ∼ | x j | − µ +3 (cid:89) l (cid:54) = j | x l | − µ , (4.30)and similarly, D x k D x j w ( x ) ∼ | x k | − µ +3 | x j | − µ +3 (cid:89) l (cid:54) = j | x l | − µ , D x k D x j w ( x ) ∼ | x k | − µ +5 | x j | − µ +3 (cid:89) l (cid:54) = j,k | x l | − µ . By an induction argument, we can show D m − x k D x j w ( x ) ∼ | x k | − µ +2 m − | x j | − µ +3 (cid:89) l (cid:54) = j,k | x l | − µ , j (cid:54) = k,D mx k w ( x ) ∼ | x k | − µ +2 m +1 (cid:89) l (cid:54) = k | x l | − µ . (4.31)Thus, for j (cid:54) = k , we have I kj ( x ) := | D m − x k D x j w ( x ) | (cid:36) m − e k + e j ( x ) ∼ | x k | − µ +2 m − | x j | − µ +4 (cid:89) l (cid:54) = j,k | x l | − µ − , (4.32)and for 1 ≤ k ≤ d, I kk ( x ) := | D mx k w ( x ) | (cid:36) me k ( x ) ∼ | x k | − µ +2 m (cid:89) l (cid:54) = k | x l | − µ − . (4.33)Then by (4.6), (4.15) and (4.33), we find that if m < µ − , then | w | B m ( R d ) = d (cid:88) k =1 (cid:90) R d I kk ( x )d x < ∞ , [[ w ]] B m ( R d ) = | w | B m ( R d ) + d (cid:88) j =1 (cid:88) k (cid:54) = j (cid:90) R d I kj ( x )d x < ∞ . (4.34)For the exact solution u e ( x ) = e −| x | , we have from (4.23) and (4.28) that f e ( x ) ∼ (1 + | x | ) s + d/ . Therefore, in this case, the error is dominated by the MCF interpolation approximation of f e ( x ) . There-fore, using (4.34) with µ = s + d/ , we conclude from Proposition 4.2 that the expected convergence is O ( N − (2 s + d )+1 / ε ) for small ε > . Remarkably, the numerical results in Figure 4.1 (a), (c), (e) perfectlyagree with the theoretical prediction (see the dashed reference lines). Indeed, it is very different from theusual Laplacian (see Proposition 4.2), we do not expect the exponential convergence, but algebraic decayof the errors. log (N) l og ( E rr o r) -6.5-6-5.5-5-4.5-4-3.5-3 s=0.3s=0.7 O(N -1.1 )O(N -1.9 ) (a) d = 1 and u e ( x ) = e − x log (N) l og ( E rr o r) -6.5-6-5.5-5-4.5-4-3.5-3 s=0.3s=0.7 O(N -1.1 )O(N -1.9 ) (b) d = 1 and u a ( x ) = (1 + | x | ) − . log (N) l og ( E rr o r) -9-8.5-8-7.5-7-6.5-6-5.5 s=0.3s=0.7 O(N -2.1 )O(N -2.9 ) (c) d = 2 and u e ( x ) = e −| x | log (N) l og ( E rr o r) -9-8.5-8-7.5-7-6.5-6-5.5-5-4.5-4 s=0.3s=0.7 O(N -2.9 )O(N -2.1 ) (d) d = 2 and u a ( x ) = (1 + | x | ) − . log (N) l og ( E rr o r) -11-10.5-10-9.5-9-8.5-8-7.5-7-6.5 s=0.3s=0.7 O(N -3.1 )O(N -3.9 ) (e) d = 3 and u e ( x ) = e −| x | log (N) l og ( E rr o r) -11-10-9-8-7-6-5-4-3-2 s=0.3s=0.7 O(N -3.1 )O(N -3.4 ) (f) d = 3 and u a ( x ) = (1 + | x | ) − . Figure 4.1:
Decay of H s -errors of the MCF scheme with γ = 1 and the scaling factor ν = 2 . s = 0 . , . r = 2 .
3. The dashed reference lines are expected orderspredicted by Proposition 4.1.
We now turn to the second case with the exact solution u a ( x ) = (1 + | x | ) − r , where we take r = 2 . r > d/ , we derive from Proposition 4.3 that f a ( x ) ∼ (1 + | x | ) s + d/ . Then by4.34 and Proposition 4.2, we have the convergence behaviour (cid:107) u a − u N (cid:107) H s ( R d ) = O ( N s − m ) + O ( N − m ) , m < r − , m < s + d − . This implies the convergence order: O ( N − min { r − s, s + d } + + ε ) . Indeed, we observe from Figure 4.1 (b),(d), (f) a perfect agreement again. For example, for d = 3 and r = 2 . , we have the rate O ( N − . ε ) for NTEGRAL FRACTIONAL LAPLACIAN 15 s = 0 . , while the rate O ( N − . ε ) for s = 0 . . Interestingly, there is a pre-asymptotic range where oneobserves a sub-geometric convergence from Figure 4.1 (d), (f) for d = 2 ,
3, but after the pre-asymptoticrange, the convergence rates become algebraic as predicted. Such a phenomenon has been also observedfor the Laguerre function approximation (cf. [43, P. 278]).We highlight that the above numerical tests validate log (N) l og ( E rr o r) -9-8-7-6-5-4-3 s=0.3, ν =1s=0.7, ν =1s=0.3, ν =4s=0.7, ν =4 (a) Different scaling factors log (N) l og ( E rr o r) -7-6.5-6-5.5-5-4.5-4-3.5-3-2.5 Hermite s=0.3Hermite s=0.7MCFs s=0.3MCFs s=0.7 (b) Comparison with Hermite approach Figure 4.2: (a). Maximum error for the exact solution u ( x ) = (1 + x ) − . with different scaling factor ν , and s = 0 . , .
7; (b). A comparison of maximum error between our method and Hermite-Galerkin method [35] forexact solution u ( x ) = (1 + x ) − . , with s = 0 . , . ν = 1. Example 4.2. (Effect of the scaling factor).
In this example, we first show the influence of thescaling factor ν to the accuracy. It is known that, with a proper choice of the scaling parameter, theaccuracy of spectral method on the unbounded domain can all be improved. Here, we plot in Figure4.2 (a) the maximum error in log-log scale of our MCF algorithm with different scaling parameter ν .We observe that for any fixed s , the two error curves are nearly parallel, which implies a proper scalingwill improve the accuracy, but will not change the convergence rate. In Figure 4.2 (b), we compare themaximum errors of our algorithm using MCFs as basis functions with the Hermite spectral method in[35], for which we take r = 2 .
3. As we can see from Figure 4.2 (b) that the convergence rates of ourapproach are faster than that of the Hermite spectral method in [35].
Example 4.3. (Accuracy for given source term f ( x ) ). Here, we further compare our MCF methodwith the Hermite function approximation in [35], where the tests were provided for given source termswith unknown solutions. We therefore compute the reference “exact” solutions with large N = 600. InFigure 4.3 (a)-(c), we compare the L -errors of our algorithm with the Hermite spectral method in [35]in one and two dimensions. It is noteworthy that the algorithm in [35] is computationally prohibitive for d = 3 . In all cases, our approach outperforms the Hermite method in both accuracy and efficiency. Wereport in Figure 4.3 (d)-(f) the maximum point-wise errors against various N with d = 2 , . The MCFmethod performs consistently well.We also tabulate in Table 4.1 the L -errors and the convergence orders of two methods (see Table 2and 3 of [35] for the data of the Hermite method). Here, f ( x ) = (1 + x ) e − x , s = 0 . , . ν = 2 . Example 4.4. (Multi-term fractional equations).
Consider the three-dimensional multi-term frac-tional Laplacian equation J (cid:88) j =1 ρ j ( − ∆) s j u ( x ) = f ( x ) , in R ; u ( x ) = 0 as | x | → ∞ . (4.35) log (N) l og ( E rr o r) -7-6.5-6-5.5-5-4.5-4-3.5-3 Hermite s=0.6Hermite s=0.9MCFs s=0.6MCFs s=0.9 (a) d = 1 and f ( x ) = (1 + x ) e − x log (N) l og ( E rr o r) -6.5-6-5.5-5-4.5-4-3.5-3 Hermite s=0.6Hermite s=0.9MCFs s=0.6MCFs s=0.9 (b) d = 1 and f ( x ) = (1 + x ) − log (N) l og ( E rr o r) -10.5-10-9.5-9-8.5-8-7.5-7-6.5-6 Hermite s=0.6Hermite s=0.9MCFs s=0.6MCFs s=0.9 (c) f ( x , x ) = (1 + x )(1 + x ) e − | x | log (N) l og ( E rr o r) -10-9.5-9-8.5-8-7.5-7-6.5-6-5.5-5 s=0.3s=0.7 (d) f ( x , x ) = (1 + 2 x + 5 x ) − log (N) l og ( E rr o r) -11.5-11-10.5-10-9.5-9-8.5-8-7.5-7-6.5 s=0.3s=0.7 (e) f ( x , x , x ) = (1 + x + 2 x + 3 x ) e − | x | log (N) l og ( E rr o r) -14-12-10-8-6-4-2 s=0.3s=0.7 (f) f ( x , x , x ) = (1 + x + 2 x + 3 x ) − . Figure 4.3: (a)-(c): A comparison of L -errors between our method and Hermite-Galerkin method in [35] fordifferent source function f ( x ); (d)-(f): The maximum errors for different source function f ( x ) with d = 2 ,
3. Inthe tests, we take and γ = 1 , ν = 2 . In Figure 4.4 (a), we plot in log-log scale the maximum errors of (4.35) against various N , where we take u ( x ) = (1 + | x | ) − π , J = 4 and s = 0 . , s = 0 . , s = 0 . , s = 0 , ρ = 1 , ρ = 2 , ρ = √ , ρ = 1 . (4.36) NTEGRAL FRACTIONAL LAPLACIAN 17
Table 4.1:
A comparison of L -error for f ( x ) = (1 + x ) e − x . s = 0 . s = 0 . N Hermite [35] Order MCF Order Hermite [35] Order MCF Order80 2.77e-04 6.29e-05 2.06e-05 2.21e-06100 2.21e-04 1.01 4.38e-05 1.61 1.53e-05 1.32 1.35e-06 2.20120 1.84e-04 1.02 3.26e-05 1.61 1.20e-05 1.33 9.10e-07 2.19140 1.57e-04 1.03 2.56e-05 1.57 9.80e-06 1.33 6.51e-07 2.18160 1.36e-04 1.03 2.08e-05 1.52 8.20e-06 1.34 4.86e-07 2.18180 1.21e-04 1.04 1.75e-05 1.49 7.00e-06 1.34 3.75e-07 2.20200 1.08e-04 1.04 1.49e-05 1.49 6.08e-06 1.34 2.96e-07 2.24220 9.79e-05 1.05 1.29e-05 1.52 5.35e-06 1.34 2.38e-07 2.29240 8.93e-05 1.06 1.12e-05 1.58 4.76e-06 1.35 1.94e-07 2.35
In Figure 4.4 (b), we plot in log-log scale the maximum errors of (4.35) against various N , where we take f ( x ) = (1 + x + 2 x + 3 x ) e − | x | , J = 4 and s = 0 . , s = 0 . , s = 0 . , s = 0 , ρ = 2 , ρ = 1 , ρ = 0 . , ρ = 1 . (4.37)We observe the algebraic decay of the errors, and the method is as accurate and efficient as the previouscases. log (N) l og ( E rr o r) -10-9-8-7-6-5-4-3 s =0.77, s =0.33, s =0.21, s =0 (a) d = 3 with given exact solution log (N) l og ( E rr o r) -10-9.5-9-8.5-8-7.5-7-6.5-6 s =0.76, s =0.41, s =0.23, s =0 (b) d = 3 with given source term Figure 4.4: (a). The maximum error for (4.35) with u ( x ) = (1 + | x | ) − π and ν = 2 .
5; (b). The maximum errorfor (4.35) with f ( x ) = (1 + x + 2 x + 3 x ) e − | x | and ν = 2 . MCF approximation of nonlinear fractional Schr¨odinger equations
In this section, we apply the fast algorithm to some nonlinear PDEs involving fractional Laplacian.As an example, we consider nonlinear fractional Schr¨odinger equation (fNLS) (cf. [32]):i ψ t = 12 ( − ∆) s ψ + γ | ψ | p ψ, x ∈ R d , t ∈ (0 , T ] ,ψ ( x,
0) = ψ ( x ) , x ∈ R d , | ψ | → , | x | → ∞ , (5.1)where i = − ψ ( x, t ) is a complex-valued wave function, the parameters γ and p are real constants, and ψ is given. It is noteworthy that the mass is conserved (cf. [7, 32]): M ( t ) = (cid:90) R d | ψ ( x, t ) | d x = M (0) , t > . (5.2) The scheme.
We adopt the time-splitting technique, and start with rewriting the fNLS (5.1) asfollows i ψ t = Aψ + Bψ, (5.3)where Aψ = γ | ψ ( x, t ) | ψ ( x, t ) , Bψ = 12 ( − ∆) s ψ ( x, t ) . The notion of time-splitting is to solve the following two subproblems:i ∂ψ ( x, t ) ∂t = Aψ ( x, t ) = γ | ψ ( x, t ) | ψ ( x, t ) , x ∈ R d , (5.4)and i ∂ψ ( x, t ) ∂t = Bψ ( x, t ) = 12 ( − ∆) s ψ ( x, t ) , x ∈ R d . (5.5)The essence of the splitting method is to solve the two sub-problems iteratively at each time step.(i). We first consider the subproblem (5.4). Multiplying (5.4) by ¯ ψ ( x, t ), we find from the resultedequation that | ψ ( x, t ) | invariant in t (see e.g., [7]). More precisely, for t ≥ t s ( t s is any given time), (5.4)becomes i ∂ψ ( x, t ) ∂t = γ | ψ ( x, t s ) | ψ ( x, t ) , t ≥ t s , x ∈ R d , (5.6)which can be integrated exactly, i.e., ψ ( x, t ) = e − i γ | ψ ( x,t s ) | ( t − t s ) ψ ( x, t s ) , t ≥ t s , x ∈ R d . (5.7)(ii). We now turn to the subproblem (5.5). Remarkably, the Fourier-like basis can diagonalise theoperator B so that e − i B ∆ t ψ can be efficiently evaluated (which is crucial for the final scheme to be timereversible and time transverse invariant). More precisely, we seek ψ N ( x, t ) ∈ V dN as an approximatesolution to (5.5), such thati (cid:0) ∂ t ψ N , v (cid:1) L ( R d ) = (cid:0) Bψ N , v (cid:1) L ( R d ) = 12 (cid:0) ( − ∆) s ψ N , v (cid:1) L ( R d ) , ∀ v ∈ V dN . (5.8)Using the Fourier-like MCF basis, we write ψ N ( x, t ) = (cid:88) k ∈ Υ N ˆ ψ k ( t ) (cid:98) T k ( x ) , x ∈ R d . (5.9)Substituting it into (5.8), and taking the inner product with (cid:98) T m ( x ), we deduce from (3.2) thati ∂ ˆ ψ m ( t ) ∂t = 12 | λ m | s ˆ ψ m ( t ) , m ∈ Υ N . (5.10)Then, we derive from (5.10) that the solution for (5.8), i.e., the numerical solution of (5.5), is given by ψ N ( x, t ) = e − i B ( t − t s ) ψ N ( x, t s ) = (cid:88) k ∈ Υ N e − i2 | λ k | s ( t − t s ) ˆ ψ k ( t s ) (cid:98) T k ( x ) , t ≥ t s . (5.11)With the exact solution (5.7) and the approximate solution (5.11) for two subproblems (5.4) and (5.5),respectively, we now describe the implementation of the fourth-order time splitting method (TS4) forsolving (5.1). Let { x p } p ∈ Υ N be tensorial grids as in (4.14), and t n = n ∆ t be the time-stepping grids. Let ψ np be the approximation of ψ ( x p , t n ) , and denote by ψ n the solution vector with components { ψ np } p ∈ Υ N .For notational convenience, we define the solution map related to (5.11): T N [ ω ; Ψ p ]( x ) = (cid:88) k ∈ Υ N e − i ω | λ k | s ∆ t ˆ Ψ k (cid:98) T k ( x ) , (5.12)where { ˆ Ψ k } are the MCF expansion coefficients computed from the sampling of Ψ ∈ V dN on the grids { x p } , and ω > NTEGRAL FRACTIONAL LAPLACIAN 19
Following [7], we carry out the fourth-order time-splitting method for the fNLS (5.1), from time t = t n to t = t n +1 , as follows ψ (1) p = e − ω γ ∆ t | ψ np | ψ np , ψ (2) p = T N [ ω ; ψ (1) p ]( x p ) ,ψ (3) p = e − ω γ ∆ t | ψ (2) p | ψ (2) p , ψ (4) p = T N [ ω ; ψ (3) p ]( x p ) ,ψ (5) p = e − ω γ ∆ t | ψ (4) p | ψ (4) p , ψ (6) p = T N [ ω , ψ (5) p ]( x p ) ,ψ n +1 p = e − ω γ ∆ t | φ (6) p | ψ (6) p , ∀ p ∈ Υ N , (5.13)where the weights are given by ω = 0 . , ω = 0 . ,ω = − . , ω = − . . (5.14)To show the stability of fourth-order splitting method, we further define (cid:107) ψ n (cid:107) N = (cid:88) j ∈ Υ N | ψ nj | ω j := N (cid:88) j =0 · · · N d (cid:88) j d =0 ψ ( x j , · · · , x j d ) ω j · · · ω j d , (5.15)where ψ nj = ψ n ( x j ), and { x j , ω j } j ∈ Υ N are the corresponding tensorial nodes and weights as in (4.10).Following [7, Lemma 3.1], we can show the property stated below. Theorem 5.1.
The
T S has the normalisation conservation, i.e., (cid:107) ψ n (cid:107) N = (cid:88) j ∈ Υ N | ψ nj | ω j = (cid:88) j ∈ Υ N | ψ ( x j ) | ω j = (cid:107) ψ (cid:107) N , n ≥ . (5.16)5.2. Numerical results.
In the computation, we take d = 2, and the initial condition to be ψ ( x , x ) = sech( x )sech( x ) exp(i( x + x )) , ( x , x ) ∈ R . (5.17)In order to test the fourth-order accuracy in time of the TS4 method, we compute a numerical solutionwith focusing case γ = − s = 0 .
7, a very fine mesh, e.g., N = 300, and a very small time step∆ t = 0 . ψ . Let ψ ∆ t be the numerical solution with N = 300 and time stepside ∆ t . Table 5.1 lists the maximum error and L -error at T = 2 for different time step size ∆ t . Theresults in Table 5.1 demonstrate the fourth-order accuracy in time of the TS4 method (5.13).In Figure 5.1, we plot the maximum errors and L -error versus space discretization N and timediscretization ∆ t . They indicate that the numerical errors decay algebraically as N increases/or ∆ t decreases. Table 5.1: Time discretization errors for the TS4 method (5.13) at T = 2 with N = 300.∆ t − L -error 2.557e-03 2.235e-04 1.616e-05 1.084e-06 6.553e-08 4.435e-09order − In Figure 5.2 (a)-(d), we depict the modulus squared of the numerical solution with defocusing case( γ = 1) obtained by TS4. Here, we take N = 200, T = 1 ,
2, and different values of fractional order s = 0 . , .
7. We observe that the solution diffused as expected. On the other hand, the blow-up of thesolution might happen for focusing case γ = − T = 1 with N = 200 and s = 0 . , . . We can observethe expected blow-up phenomenon. log (N) l og ( E rr o r) -5.5-5-4.5-4-3.5-3-2.5-2-1.5 T=2, γ =-1, s=0.7, L -ErrorT=2, γ =-1, s=0.7, L ∞ -Error (a) Errors vs. N Δ t10 -3 -2 -1 E rr o r -10 -8 -6 -4 -2 T=2, γ =-1, s=0.7, L ∞ -ErrorT=2, γ =-1, s=0.7, L -Error (b) Errors vs. ∆ t Figure 5.1: (a). The numerical error of (5.17) with s = 0 . , γ = − , T = 2; (b). The numerical error of (5.17)with s = 0 . , γ = − , T = 2. (a) T = 1, s = 0 . γ = 1 (b) T = 2, s = 0 . γ = 1(c) T = 1, s = 0 . γ = − T = 1, s = 0 . γ = − Figure 5.2:
Profiles of the modulus square of the numerical solutions at different time and with different fractionalorders.
Concluding Remarks.
We developed a fast MCF-spectral-Galerkin method for PDEs involvingintegral fractional Laplacian in R d . The fast solver is integrated with two critical components: (i) theDunford-Taylor formulation for the fractional Laplacian; and (ii) Fourier-like bi-orthogonal MCFs as NTEGRAL FRACTIONAL LAPLACIAN 21 basis functions. The fast spectral algorithm could achieve a quasi-optimal computational cost. Differentfrom the existing works on bounded domains (cf. [11, 12]), the integration in t is evaluated explicitly,and the fractional Laplacian can be fully diagonalised under (i) and (ii). Indeed, the existing approachesfor fractional Laplacian in unbounded domains are either too complicated or computational prohibitiveeven for d = 2 . However, the fast solver works for any dimension, and can be easily incorporated withe.g., the hyperbolic cross and sparse grids (cf. [45]) when the dimension is high.The proposed method can be extended to invert the operator D s := ( − ∆ + γ I ) s with s ∈ (0 ,
1) and γ > . In fact, one can verify readily that the Dunford-Taylor formulation in Lemma 3.2 takes the form (cid:0) D s u, D s v (cid:1) L ( R d ) = C s (cid:90) ∞ t − s (cid:90) R d (cid:0) ( − ∆ + γ I ) (cid:0) I + t ( − ∆ + γ I ) (cid:1) − u (cid:1) ( x ) v ( x ) d x d t. (5.18)Then the fast algorithm in Theorem 3.2 is extendable to this case straightforwardly. Appendix A. Proof of Proposition 4.2
The results with d = 1 were derived in [51], so it suffices to prove them for integer d ≥
2. Note that F (cid:8) e −| x | (cid:9) ( ξ ) = 1(2 π ) d/ (cid:90) R d e −| x | e − i x · ξ d x = 1(2 π ) d/ (cid:90) R e − x e − i x ξ d x · · · (cid:90) R e − x d e − i x d ξ d d x d = 12 d/ e − | ξ | , where we used the identity (cf. [24, P. 339]): (cid:90) R e − x e − i xξ d x = √ πe − ξ . Thus from the definition (1.1), we obtain( − ∆) s (cid:8) e −| x | (cid:9) ( x ) = F − (cid:110) | ξ | s F (cid:8) e −| x | (cid:9) ( ξ ) (cid:111) = 12 d/ (2 π ) d/ (cid:90) R d | ξ | s e − | ξ | e i x · ξ d ξ = 2 d d/ (2 π ) d/ (cid:90) R d + | ξ | s e − | ξ | cos( x ξ ) cos( x ξ ) · · · cos( x d ξ d )d ξ. (A.1)We proceed with the calculation by using the d -dimensional spherical coordinates: ξ = r cos θ ; ξ = r sin θ cos θ ; · · · · · · ; ξ d − = r sin θ · · · sin θ d − cos θ d − ; ξ d = r sin θ · · · sin θ d − sin θ d − , r = | ξ | , (A.2)so we can write ( − ∆) s (cid:8) e −| x | (cid:9) ( x ) = 1 π d/ (cid:90) ∞ r s + d − e − r I ( r ; x )d r, (A.3)where I ( r ; x ) = (cid:90) [0 , π ] d − cos (cid:0) rx cos θ (cid:1) cos (cid:0) rx sin θ cos θ (cid:1) · · · cos (cid:0) rx d − sin θ · · · sin θ d − cos θ d − (cid:1) cos (cid:0) rx d sin θ · · · sin θ d − sin θ d − (cid:1) (sin θ ) d − (sin θ ) d − · · · (sin θ d − ) d θ d θ · · · d θ d − . We first integrate I ( r ; x ) with respect to θ d − . To do this, we recall the integral formula involving theBessel functions (cf. [24, P. 732]): for real µ, ν > − a, b > (cid:90) π J ν ( a sin θ ) J µ ( b cos θ ) sin ν +1 θ cos µ +1 θ d θ = a ν b µ J ν + µ +1 (cid:0) √ a + b (cid:1) ( a + b ) ( ν + µ +1) / , (A.4) Then using the identity cos z = (cid:112) πz/ J − / ( z ) and (A.4) (with a = rx d − sin θ · · · sin θ d − , b = rx d sin θ · · · sin θ d − and µ = ν = − / (cid:90) π cos (cid:0) rx d − sin θ · · · sin θ d − cos θ d − (cid:1) cos (cid:0) rx d sin θ · · · sin θ d − sin θ d − (cid:1) d θ d − = π J (cid:0) r sin θ · · · sin θ d − (cid:113) x d − + x d (cid:1) . Substituting the above into I ( r, x ), and applying the same argument to θ d − , θ d − , · · · , θ iteratively d − I ( r ; x ) = (cid:16) π (cid:17) d ( r | x | ) − d J d − ( r | x | ) . (A.5)We proceed with the integral identity (cf. [24, P. 713]): for real µ + ν > − p > (cid:90) R + J µ ( bt ) e − p t t ν − d t = b µ Γ(( µ + ν ) / µ +1 p ν + µ Γ( µ + 1) F (cid:16) µ + ν µ + 1; − b p (cid:17) . (A.6)Then, substituting (A.5) into (A.3) and using (A.6) (with µ = d/ − ν = 2 s + d/ − ∆) s (cid:8) e −| x | (cid:9) = | x | − d d/ (cid:90) ∞ r s + d e − r J d − ( r | x | )d r = 2 s Γ( s + d/ d/ F (cid:16) s + d d −| x | (cid:17) . This yields (4.22). The asymptotic behaviour (4.23) follows from the property (cf. [8, P. 278]): F ( a ; b ; z ) = Γ( b )Γ( b − a ) ( − z ) − a (cid:8) O ( | z | − ) (cid:9) . (A.7)Then (4.23) follows. This completes the proof. Appendix B. Proof of Proposition 4.3
The identity with d = 1 can be found in [50], so we assume that d ≥
2. Using the d -spherical coordinatesystem in (A.2), we obtain from (A.5) that F (cid:110) | x | ) γ (cid:111) ( ξ ) = 1(2 π ) d/ (cid:90) R d e − i x · ξ (1 + | x | ) γ d x = 2 d (2 π ) d/ (cid:90) R d + cos( x ξ ) cos( x ξ ) · · · cos( x d ξ d )(1 + | x | ) γ d x = (cid:16) π (cid:17) d (cid:90) ∞ r d − (1 + r ) γ I ( r ; ξ ) d r = | ξ | − d (cid:90) ∞ r d (1 + r ) γ J d − ( r | ξ | )d r. Recall the integral formula (cf. [24, P. 686]): for − < ν < µ + and a, b > (cid:90) ∞ x ν +1 ( x + a ) µ +1 J ν ( bx )d x = a ν − µ b µ µ Γ( µ + 1) K ν − µ ( ab ) , (B.1)where K ν ( x ) is the modified Bessel functions of the second kind. Note that K − ν ( x ) = K ν ( x ). Thenletting µ = γ − ν = d/ − F (cid:110) | x | ) γ (cid:111) ( ξ ) = | ξ | γ − d γ − Γ( γ ) K γ − d ( | ξ | ) . We also use the integral formula (cf. [24, P. 692]): for real a >
0, real b , and ν − λ + 1 > | µ | , (cid:90) ∞ x − λ K µ ( ax ) J ν ( bx )d x = b ν Γ (cid:0) ( ν − λ + µ + 1) / (cid:1) Γ (cid:0) ( ν − λ − µ + 1) / λ +1 a ν − λ +1 Γ( ν + 1) × F (cid:16) ν − λ + µ + 12 , ν − λ − µ + 12 ; ν + 1; − b a (cid:17) . (B.2)Once again, using the d -spherical coordinate system (A.2), (A.5) and (B.2) (with λ = − s − γ , µ = γ − d/ ν = d/ − − ∆) s (cid:110) | x | ) γ (cid:111) = 1(2 π ) d γ − Γ( γ ) (cid:90) R d e i x · ξ | ξ | s + γ − d K γ − d ( | ξ | ) d ξ NTEGRAL FRACTIONAL LAPLACIAN 23 = 2 d (2 π ) d γ − Γ( γ ) (cid:90) R d + cos( x ξ ) cos( x ξ ) · · · cos( x d ξ d ) | ξ | s + γ − d K γ − d ( | ξ | ) d ξ = 2 d − γ +1 π d Γ( γ ) (cid:90) ∞ r s + γ + d − K γ − d ( r ) I ( r, x ) d r = 2 − γ +1 Γ( γ ) | x | − d (cid:90) ∞ r s + γ K γ − d ( r ) J d − ( r | x | )d r = 2 s Γ( s + γ )Γ( s + d )Γ( γ )Γ( d ) F (cid:16) s + γ, s + d d −| x | (cid:17) . This completes the derivation of (4.24).According to [8, P. 76], the asymptotic behaviour of the hypergeometric function for large | x | (unless a − b is an integer) is F ( a, b ; c ; x ) = λ | x | − a + λ | x | − b + O ( | x | − a − ) + O ( | x | − b − ) . (B.3)where λ and λ are constants; if a − b is an integer, z − a or z − b has to be multiplied by a factor ln( x ).Then we have the asymptotic behaviour of ( − ∆) s (cid:8) | x | ) γ (cid:9) as | x | → ∞ in (4.25)-(4.26). References [1]
G. Acosta, F. M. Bersetche, and J. P. Borthagaray , A short FE implementation for a 2d homogeneous Dirichletproblem of a fractional Laplacian , Comput. Math. Appl., 74 (2017), pp. 784–816.[2]
G. Acosta and J. P. Borthagaray , A fractional Laplace equation: regularity of solutions and finite element approx-imations , SIAM J. Numer. Anal., 55 (2017), pp. 472–495.[3]
M. Agranovich , Sobolev spaces, their generalizations and elliptic problems in smooth and Lipschitz domains , Springer,2015.[4]
M. Ainsworth and C. Glusa , Aspects of an adaptive finite element method for the fractional laplacian: a priori anda posteriori error estimates, efficient implementation and multigrid solver , Comput. Methods Appl. Mech. Engrg., 327(2017), pp. 4–35.[5] ,
Hybrid finite element–spectral method for the fractional Laplacian: approximation theory and efficient solver ,SIAM J. Sci. Comput., 40 (2018), pp. A2383–A2405.[6]
I. Babuska , Survey lectures on the mathematical foundations of the finite element method , The Mathematical Foun-dations of the Finite Element Method with Applicaions to Partial Differential Equations, (1972), pp. 3–359.[7]
W. Bao and J. Shen , A fourth-order time-splitting Laguerre–Hermite pseudospectral method for Bose–Einstein con-densates , SIAM J. Sci. Comput., 26 (2005), pp. 2010–2028.[8]
H. Bateman , Higher transcendental functions [volumes i-iii] , 1953.[9]
D. A. Benson, S. W. Wheatcraft, and M. M. Meerschaert , Application of a fractional advection-dispersionequation , Water Resour. Res., 36 (2000), pp. 1403–1412.[10]
A. Bonito, J. P. Borthagaray, R. H. Nochetto, E. Ot´arola, and A. J. Salgado , Numerical methods for fractionaldiffusion , Comput. Vis. Sci., 19 (2018), pp. 19–46.[11]
A. Bonito, W. Lei, and J. E. Pasciak , Numerical approximation of the integral fractional Laplacian , Numer. Math.,142 (2019), pp. 235–278.[12] ,
On sinc quadrature approximations of fractional powers of regularly accretive operators , J. Numer. Math., 27(2019), pp. 57–68.[13]
D. Brockmann, L. Hufnagel, and T. Geisel , The scaling laws of human travel , Nature, 439 (2006), p. 462.[14]
L. Caffarelli and L. Silvestre , An extension problem related to the fractional Laplacian , Comm. Partial DifferentialEquations, 32 (2007), pp. 1245–1260.[15]
B. Carmichael, H. Babahosseini, S. Mahmoodi, and M. Agah , The fractional viscoelastic response of human breasttissue cells , Phys. Biol., 12 (2015), p. 046001.[16]
L. Chen, Z. Mao, and H. Li , Jacobi-Galerkin spectral method for eigenvalue problems of Riesz fractional differentialequations , arXiv preprint arXiv:1803.03556, (2018).[17]
S. Chen, J. Shen, and L.-L. Wang , Laguerre functions and their applications to tempered fractional differentialequations on infinite intervals , J. Sci. Comput., 74 (2018), pp. 1286–1313.[18]
J. H. Cushman and T. Ginn , Nonlocal dispersion in media with continuously evolving scales of heterogeneity , Transp.Porous Media, 13 (1993), pp. 123–138.[19]
W. Deng , Finite element method for the space and time fractional Fokker–Planck equation , SIAM J. Numer. Anal.,47 (2008), pp. 204–226.[20]
W. Deng, B. Li, Z. Qian, and H. Wang , Time discretization of a tempered fractional Feynman–Kac equation withmeasure data , SIAM J. Numer. Anal., 56 (2018), pp. 3249–3275.[21]
Q. Du , Nonlocal modeling, analysis, and computation , vol. 94, CBMS-NSF Regional Conference Series in AppliedMathematics, SIAM, 2019. [22]
S. Duo and Y. Zhang , Computing the ground and first excited states of the fractional Schr¨odinger equation in aninfinite potential well , Commun. Comput. Phys., 18 (2015), pp. 321–350.[23] ,
Finite difference methods for two and three dimensional fractional Laplacian with applications to solve thefractional reaction-diffusion equations , arXiv preprint arXiv:1804.02718, (2018).[24]
I. S. Gradshteyn and I. M. Ryzhik , Table of Integrals, Series, and Products , Elsevier/Academic Press, Amsterdam,eighth ed., 2015. Translated from the Russian, Translation edited and with a preface by Daniel Zwillinger and VictorMoll, Revised from the seventh edition [MR2360010].[25]
B. Guo and Z. Wang , Modified Chebyshev rational spectral method for the whole line , in Proceedings of the fourthinternational conference on dynamical systems and differential equations, 2002, pp. 365–374.[26]
X. Guo, Y. Li, and H. Wang , A high order finite difference method for tempered fractional diffusion equations withapplications to the CGMY model , SIAM J. Sci. Comput., 40 (2018), pp. A3322–A3343.[27]
Y. Hatano and N. Hatano , Dispersive transport of ions in column experiments: An explanation of long-tailed profiles ,Water Resour. Res., 34 (1998), pp. 1027–1033.[28]
D. Hou and C. Xu , A fractional spectral method with applications to some singular problems , Adv. Comput. Math.,43 (2017), pp. 911–944.[29]
Y. Huang and A. Oberman , Numerical methods for the fractional Laplacian: a finite difference-quadrature approach ,SIAM J. Numer. Anal., 52 (2014), pp. 3056–3084.[30]
B. Jin, R. Lazarov, and Z. Zhou , Error estimates for a semidiscrete finite element method for fractional orderparabolic equations , SIAM J. Numer. Anal., 51 (2013), pp. 445–466.[31]
B. Jin, B. Li, and Z. Zhou , Numerical analysis of nonlinear subdiffusion equations , SIAM J. Numer. Anal., 56 (2018),pp. 1–23.[32]
C. Klein, C. Sparber, and P. Markowich , Numerical study of fractional nonlinear Schr¨odinger equations , Proc.Ser. A Math. Phys. Eng. Sci., 470 (2014), p. 20140364.[33]
A. Lischke, G. Pang, M. Gulian, F. Song, C. Glusa, X. Zheng, Z. Mao, W. Cai, M. M. Meerschaert,M. Ainsworth, et al. , What is the fractional Laplacian? , arXiv preprint arXiv:1801.09767, (2018).[34]
Z. Mao, S. Chen, and J. Shen , Efficient and accurate spectral method using generalized Jacobi functions for solvingRiesz fractional differential equations , Appl. Numer. Math., 106 (2016), pp. 165–181.[35]
Z. Mao and J. Shen , Hermite spectral methods for fractional PDEs in unbounded domains , SIAM J. Sci. Comput.,39 (2017), pp. A1928–A1950.[36]
B. McCay and M. N. L. Narasimhan , Theory of nonlocal electromagnetic fluids , Arch. Mech., 33 (1981), pp. 365–384.[37]
R. Metzler and J. Klafter , The random walk’s guide to anomalous diffusion: a fractional dynamics approach ,Phys. Rep., 339 (2000), pp. 1–77.[38] ,
The restaurant at the end of the random walk: recent developments in the description of anomalous transportby fractional dynamics , J. Phys. A, 37 (2004), p. R161.[39]
E. W. Montroll and G. H. Weiss , Random walks on lattices. II , J. Math. Phys., 6 (1965), pp. 167–181.[40]
E. D. Nezza, G. Palatucci, and E. Valdinoci , Hitchhiker’s guide to the fractional Sobolev spaces , Bull. Sci. Math.,136 (2012), pp. 521–573.[41]
R. H. Nochetto, E. Ot´arola, and A. J. Salgado , A PDE approach to fractional diffusion in general domains: apriori error analysis , Found. Comput. Math., 15 (2014), pp. 733–791.[42]
R. H. Nochetto, E. Otarola, and A. J. Salgado , A PDE approach to space-time fractional parabolic problems ,SIAM J. Numer. Anal., 54 (2016), pp. 848–873.[43]
J. Shen, T. Tang, and L.-L. Wang , Spectral methods: algorithms, analysis and applications , vol. 41, Springer Science& Business Media, 2011.[44]
J. Shen and L. Wang , Some recent advances on spectral methods for unbounded domains , Commun. Comput. Phys.,5 (2009), pp. 195–241.[45]
J. Shen, L.-L. Wang, and H. Yu , Approximations by orthonormal mapped Chebyshev functions for higher-dimensionalproblems in unbounded domains , J. Comput. Appl. Math., 265 (2014), pp. 264–275.[46]
M. Shlesinger, B. West, and J. Klafter , L´evy dynamics of enhanced diffusion: Application to turbulence , Phys.Rev. Lett., 58 (1987), p. 1100.[47]
S. A. Silling , Reformulation of elasticity theory for discontinuities and long-range forces , J. Mech. Phys. Solids, 48(2000), pp. 175–209.[48]
D. W. Sims, E. J. Southall, N. E. Humphries, G. C. Hays, C. J. Bradshaw, J. W. Pitchford, A. James, M. Z.Ahmed, A. S. Brierley, M. A. Hindell, et al. , Scaling laws of marine predator search behaviour , Nature, 451(2008), p. 1098.[49]
G. Szeg¨o , Orthogonal polynomials , vol. 23, American Mathematical Soc., 1939.[50]
T. Tang, L.-L. Wang, H. Yuan, and T. Zhou , Rational spectral methods for PDEs involving fractional Laplacian inunbounded domains , arXiv preprint arXiv:1905.02476, (2019).[51]
T. Tang, H. Yuan, and T. Zhou , Hermite spectral collocation methods for fractional PDEs in unbounded domains ,Commun. Comput. Phys., 24 (2018), pp. 1143–1168.[52]