A universal solution scheme for fractional and classical PDEs
AA universal solution scheme for fractional and classical PDEs
Yixuan Wu ∗ , Yanzhi Zhang † Abstract
We propose a unified meshless method to solve classical and fractional PDE problems with( − ∆) α for α ∈ (0 , α = 2) and fractional ( α <
2) Laplacians, one local andthe other nonlocal, have distinct properties. Therefore, their numerical methods and computerimplementations are usually incompatible. We notice that for any α ≥
0, the Laplacian ( − ∆) α of generalized inverse multiquadric (GIMQ) functions can be analytically written by the Gausshypergeometric function, and thus propose a GIMQ-based method. Our method unifies thediscretization of classical and fractional Laplacians and also bypasses numerical approximationto the hypersingular integral of fractional Laplacian. These two merits distinguish our methodfrom other existing methods for the fractional Laplacian. Extensive numerical experiments arecarried out to test the performance of our method. Compared to other methods, our method canachieve high accuracy with fewer number of unknowns, which effectively reduces the storage andcomputational requirements in simulations of fractional PDEs. Moreover, the meshfree naturemakes it free of geometric constraints and enables simple implementation for any dimension d ≥
1. Additionally, two approaches of selecting shape parameters, including condition number-indicated method and random-perturbed method, are studied to avoid the ill-conditioning issueswhen large number of points.
Key words.
Fractional Laplacian, radial basis functions, generalized inverse multiquadratics,Gauss hypergeometric function, meshless method, variable shape parameters.
Over the last decade, fractional partial differential equations (PDEs) have been well recognizedfor their ability to describe anomalous diffusion phenomena in many complex systems [10, 11, 24].Mathematically, anomalous diffusion can be modeled by the fractional Laplacian ( − ∆) α (for α < ∂ xx + ∂ yy + ∂ zz for (normal) diffusion. It is well-knownthat the fractional Laplacian ( − ∆) α collapses to the classical Laplacian − ∆ as α → − , butthe properties of these two operators are essentially different. The fractional Laplacian ( − ∆) α ,representing the infinitesimal generator of a symmetric α -stable L´evy process, is a nonlocal operator,while the classical Laplacian ∆ is a local operator. Due to this distinct difference, most existingnumerical methods developed for the classical Laplacian can not be applied to solve problems withthe fractional Laplacian. So far, numerical methods for the fractional Laplacian ( − ∆) α still remainlimited. Moreover, some discretization techniques (e.g., finite element methods) for the classicaland fractional Laplacians are incompatible, and different computer codes are required for theirpractical implementation. In this work, we propose a meshless pseudospectral method based on ∗ Department of Mathematics and Statistics, Missouri University of Science and Technology, Rolla, MO 65409(Email: [email protected]) † Department of Mathematics and Statistics, Missouri University of Science and Technology, Rolla, MO 65409(Email: [email protected]; URL: http://web.mst.edu/ ∼ zhangyanz) a r X i v : . [ m a t h . NA ] J a n he generalized inverse multiquadric functions to solve normal and anomalous diffusion problems.It not only enriches the collection of numerical methods for studying the fractional Laplacian, butalso provides a unified approach to solve both classical and fractional PDEs.Let Ω ⊂ R d (for d = 1 ,
2, or 3) be an open bounded domain. Consider the following diffusionproblem with Dirichlet boundary conditions: ∂ t u ( x , t ) = − κ ( − ∆) α u + f ( x , t, u ) , for x ∈ Ω , t > , (1.1) u ( x , t ) = g ( x , t ) , for x ∈ Υ , t ≥ , (1.2)where the diffusive power α ∈ (0 , κ >
0. The notation Υ depends on α ,i.e., Υ = ∂ Ω for α = 2, while Υ = Ω c = R d \ Ω if α <
2. The Laplace operator ( − ∆) α is definedvia an α -parametric pseudo-differential form [31, 36]:( − ∆) α u ( x ) = F − (cid:2) | ξ | α F [ u ] (cid:3) , for α > , (1.3)where F represents the Fourier transform with the associated inverse transform F − . The pseudo-differential operator in (1.3) gives a unified definition to the classical and fractional Laplacians.For α = 2, it reduces to the spectral representation of the classical Laplacian − ∆. Thus, (1.1)–(1.2) collapses to a classical (normal) diffusion problem with Dirichlet boundary conditions on ∂ Ω. While α ∈ (0 , c (i.e., the complement of domain Ω).The fractional Laplacian ( − ∆) α can be also defined via a hypersingular integral [31, 36], i.e.,( − ∆) α u ( x ) = C d,α P . V . (cid:90) R d u ( x ) − u ( y ) | x − y | d + α d y , for α ∈ (0 , , (1.4)where P . V . stands for the principal value integral, and the normalization constant C d,α = 2 α − α Γ( α + d/ √ π d Γ(1 − α/ · ) being the Gamma function. This integral operator provides a pointwise definition of thefractional Laplacian. Note that the pseudo-differential operator in (1.3) unifies the definition ofclassical and fractional Laplacians so as to provide a foundation for developing compatible schemesfor these two operators, but it is challenging to incorporate non-periodic boundary conditions into(1.3). In contrast to (1.3), the pointwise definition of the fractional Laplacian in (1.4) can easilyincorporate other boundary conditions, but it is incompatible to the classical Laplacian (i.e., α (cid:54) = 2in (1.4)). As shown in [36, 30], the two definitions of the fractional Laplacian in (1.3) and (1.4) areequivalent for functions in the Schwartz space. Hence, we want to combine the advantages of bothdefinitions and develop a numerical method that is compatible for both classical and fractionalLaplacians and can easily work with Dirichlet boundary conditions as in (1.2).So far, numerical methods for the fractional Laplacian ( − ∆) α still remain scant. Most of themare based on the pointwise definition of the fractional Laplacian and thus need to approximate thehypersingular integral in (1.4); see [1, 3, 13, 2, 14, 15, 7] and references therein. On the other side,meshless RBF-based methods have been recognized for their flexiblity with complex geometries,high accuracy and efficiency, and simplicity of implementation. They have been widely applied tosolve classical PDEs [19, 18, 5]. But, the application of RBF-based methods in solving nonlocalor fractional PDEs is still very recent. To the best of our knowledge, two meshless RBF-basedmethods were recently proposed in the literature to solve problems governed by the fractional2aplacian ( − ∆) α (for α <
2) – one is the Wendland RBF method in [35], and the other is theGaussian RBF method in [8]. The Wendland RBF-based method directly approximates the Fourierintegrals in definition (1.3) and considers the extended Dirichlet boundary conditions (for α < ω ⊂ Ω c . It then solves the problem on Ω ∪ ω (instead of Ω), andconsequently requires significant amount of RBF points and demands high computational cost. Incontrast, the Gaussian RBF-based method integrates the extended boundary conditions into thescheme via the pointwise definition (1.4), and no boundary truncation is needed. Moreover, it candiscretize both classical and fractional Laplacians in a single scheme [8].The aim of this work is to develop a unified meshless pseudospectral method to solve bothclassical ( α = 2) and fractional ( α <
2) PDEs. Its unique feature – compatibility with the classicalLaplacian – makes our method distinguish from other numerical methods (e.g., in [13, 14, 3, 2, 7,1, 6, 7]) for the fractional Laplacian. Moreover, our method takes great advantage of the Laplacian( − ∆) α (for both α = 2 and α <
2) of generalized inverse multiquadric functions so as to bypassnumerical approximations to the hypersingular integral of fractional Laplacian in (1.4). Extensivenumerical experiments are carried out to test the performance of our method. Compared to othermethods, our method can achieve high accuracy with fewer number of unknowns, which effectivelyreduces the storage and computational requirements in simulations of fractional PDEs. Moreover,the meshfree nature makes it free of geometric constraints and enables simple implementation forany dimension d ≥
1. We also study two approaches, namely the condition number-indicatedapproach and random-perturbed approach, for selecting the shape parameters. Numerical studiesshow that they can effectively control the ill-conditioning issues and improve the performance ofour method.The paper is organized as follows. In Section 2, we introduce radial basis functions and theproperties of generalized (inverse) multiquadric functions. In Section 3, we propose a meshlessmethod based on generalized inverse multiquadric functions to solve the diffusion problem (1.1)–(1.2). In Section 4, we test the accuracy and compare the proposed method with the GaussianRBF-based method in [8]. Two approaches in selecting the shape parameter are also studied. InSection 5, we further test the performance of our method in solving elliptic problems and diffusionequations. Finally, discussion and conclusion are made in Section 6.
Radial basis functions (RBFs) are a family of functions that depend on the distance of point x to agiven center point y , i.e. ϕ ( x ) = ϕ ( | x − y | ) for x , y ∈ R d . RBFs have been well recognized for theirsuccess in interpolating high-dimensional scattered data. The application of RBFs in solving PDEswas first proposed by Kansa in 1990 [27, 28]. Since then, RBF-based methods have been widelyapplied to solve problems arising in various applications (see [19, 18, 5] and references therein).In the literature, RBFs are usually divided into two main categories: globally-supported func-tions (e.g., see Table 1) and compactly-supported functions (e.g., Wendland function). The Gaus-sian and multiquadric functions are two well-known globally supported RBFs, and both of themare infinitely differentiable. The multiquadric function was first proposed by Hardy to interpolatescattered data on a topographical map [25]. Later, Franke compared the multiquadric interpola-tion with other methods and concluded that the multiquadric function outperforms others in manyaspects, including numerical accuracy, computational time, and implementation simplicity [21, 22].Recently, a broad class of generalized multiquadric functions have received a lot of attention in theliterature [9, 38]. 3ame Definition ϕ ( r )Gaussian e − ( εr )2 Multiquadric (cid:112) εr ) Inverse multiquadric √ εr ) Polyharmonic spline r m ln( r ) , m ∈ N Table 1: Some examples of globally-supported RBFs with r = | x − y | . The generalized multiquadric radial basis functions take the following form [9, 38]: ϕ ( r ) = (cid:0) ε r (cid:1) β , for β ∈ R \ N , (2.1)where ε > N represents the set of non-negative integers. Thefunction in (2.1) covers a wide range of infinitely differentiable RBFs, including, e.g., • the well-known Hardy’s multiquadric function when β = [25, 26]; • the inverse multiquadric function when β = − [29, 40]; • the inverse quadratic function when β = − β > generalized multiquadric (GMQ) function, while the generalized inverse multiquadric (GIMQ) function if β <
0. It shows thatthe GIMQ functions are strictly positive definite, while the GMQ functions are strictly conditionallypositive definite of order (cid:100) β (cid:101) , where (cid:100)·(cid:101) denotes the ceiling function [9, 38]. The strictly positivedefinite is an important property to ensure the system matrix from interpolation to be invertible.Note that both GIMQ and Gaussian RBFs are strictly positive definite. The shape of Gaussianfunctions is controlled by parameter ε – the smaller the shape parameter ε , the flatter the Gaussianfunction. While the GIMQ functions are controlled by both shape parameter ε and power β ; seethe comparison of their effects in Fig. 1. It shows that the GIMQ functions monotonically decrease Figure 1: Illustration of GIMQ functions for different shape parameter ε and power β .as r → ∞ . The larger the power β (or the smaller the shape parameter ε ), the flatter the GIMQfunction. The effect of ε is more significant in determining the shape of GIMQ functions.4he discussion of GIMQ functions is recent, and thus their applications in solving PDEs remainvery limited. It shows that for β < − d/ (cid:98) ϕ ( ξ ) = 2 β ε β Γ( − β ) (cid:0) ε | ξ | (cid:1) − ( β + d ) K − ( β + d ) (cid:18) | ξ | ε (cid:19) , for ξ ∈ R d \{ } , (2.2)where K v denotes the univariate modified Bessel function of the second kind, i.e., K v ( x ) := (cid:90) ∞ e − x cosh( τ ) cosh( vτ ) dτ, for v ∈ R , x > . In this work, we will consider GIMQ functions (2.1) with a dimension-dependent power β = − ( d +1) /
2. In this case, the GIMQ function is also known as the Poisson kernel (up to a constant) [44].Substituting β = − ( d + 1) / (cid:98) ϕ ( ξ ) = 2 − d/ √ π Γ (cid:0) ( d + 1) / (cid:1) e −| ξ | /ε , ξ ∈ R d \{ } , by noticing K − ( x ) = K ( x ) = (cid:112) π/ (2 x ) e − x . Next, we will introduce the Laplacian of GIMQfunctions, which is an important building block of our method. Lemma 2.1 ( Laplacian of generalized inverse multiquadrics).
Let the GIMQ function u ( x ) = (1 + | x | ) − ( d +1) / , i.e. choosing β = − ( d + 1) / in (2.1), for dimension d ≥ . Then,the function ( − ∆) α u can be analytically given by [16]: ( − ∆) α u ( x ) = 2 − d √ π Γ (cid:0) d + α (cid:1) Γ( d/ d + 1) / F (cid:16) d + α , d + 1 + α d −| x | (cid:17) , for α ≥ , (2.3) where F denotes the Gauss hypergeometric function. Lemma 2.1 provides the analytical results of ( − ∆) α u for the GIMQ function u ( x ) = (1 + | x | ) − ( d +1) / at any point x ∈ R d . Moreover, this result holds for any power α ≥
0. In the specialcase of α = 0, the right hand side of (2.3) reduces to the GIMQ function u ( x ) = (1 + | x | ) − ( d +1) / .While α = 2, it collapses to − ∆ u ( x ) = ( d + 1) (cid:0) | x | (cid:1) − d +52 (cid:2) d − | x | (cid:3) , for x ∈ R d . Fig. 2 illustrates the results ( − ∆) α u for various α and d = 1 ,
2. Note that when d = 1 the GIMQfunction in Lemma 2.1 actually reduces to the inverse quadratic function. Fig. 2 shows that thelarger the parameter α , the faster the function ( − ∆) α u decays to zero as | x | → ∞ . For the same α , the function ( − ∆) α u decays faster in higher dimension.The uniform representation of classical and fractional Laplacian of GIMQ functions in Lemma2.1 provides a foundation for developing compatible GIMQ-based schemes for these two operators.Moreover, the analytical formulation in (2.3) allows us to avoid numerical approximation of thehypersingular integral for the fractional ( α <
2) Laplacian. The Gaussian hypergeometric function F can be calculated with the well-established methods (see [32] and references therein). If dimen-sion d is odd (i.e., d = 1, or 3), we can use the properties of the Gauss hypergeometric function F and rewrite (2.3) in terms of elementary functions. In one-dimensional ( d = 1) case, the resultin (2.3) reduces to( − ∆) α (cid:0) x (cid:1) − = Γ(1 + α ) cos (cid:2) (1 + α ) arctan | x | (cid:3) (1 + x ) − α , for x ∈ R . (2.4)5 = 0 = 0.4 = 1.5 = 2 -6 -3 0 3 6-0.50.82.13.44.76 = 0 = 0.4 = 1.5 = 2 Figure 2: Illustration of ( − ∆) α u for the GIMQ function u ( x ) = (cid:0) | x | (cid:1) − ( d +1) / .Note that cos (cid:0) arctan | x | (cid:1) = √ x , and thus (2.4) collapses to function u when α = 0. While inthree-dimensional ( d = 3) case, we have( − ∆) α (cid:0) | x | (cid:1) − = Γ(3 + α ) / , for x = , Γ(2 + α )(1 + | x | ) − α | x | sin (cid:0) (2 + α ) arctan | x | (cid:1) , otherwise , (2.5)for x ∈ R . Hence, when d = 1 or 3, one can use the results in (2.4) and (2.5) to replace (2.3) soas to avoid computing the Gauss hypergeometric function F .Additionally, the Laplace operator satisfies the following properties [8]:( − ∆) α (cid:2) u ( x − x ) (cid:3) = U ( x − x ) , for x ∈ R d , (2.6)( − ∆) α (cid:2) u ( ξ x ) (cid:3) = | ξ | α U ( ξ x ) , for ξ ∈ R , (2.7)for any α ≥
0, where we denote function U ( x ) := ( − ∆) α u ( x ). These two properties together with(2.3) play an important role in the design of our meshless methods in Section 3. In this section, we will present a new meshless method based on the GIMQ RBFs to solve thediffusion problem (1.1)–(1.2). Our method differs from other RBF-based methods in the followingaspects. First, it provides an α -parametric scheme that can solve both classical ( α = 2) andfractional ( α <
2) PDEs problems seamlessly. Second, utilizing the properties of Laplace operatorsand GIMQ functions, our method avoids numerical evaluations of fractional derivatives, which areusually approximated by quadrature rules in many other methods [33, 35]. This not only reducescomputational cost but simplifies the implementations especially in high dimensions. Last but notleast, when α < − ∆) α for α ∈ (0 , u ( x , t ) ≈ (cid:98) u ( x , t ) := ¯ N (cid:88) i =1 λ i ( t ) ϕ ε ( | x − x i | ) , for x ∈ ¯Ω , t ≥ , (3.1)6here x i are the RBF center points, and ϕ ε ( | x − x i | ) represents the GIMQ function centered at x i .In our method, we choose the GIMQ basis function as ϕ ε ( | x − x i | ) = (cid:0) ε i | x − x i | (cid:1) − d +12 , where without loss of generality, we assume that the shape parameter ε i associated with each centerpoint x i is different. Note that the power β = − ( d + 1) /
2, that is, the GIMQ basis functions changefor different dimension d ≥
1. For all α ∈ (0 , ∪ ∂ Ω. More specifically, we assume that point x i ∈ Ω if index 1 ≤ i ≤ N , while x i ∈ ∂ Ω ifindex N + 1 ≤ i ≤ ¯ N .We start with approximating the fractional ( α <
2) Laplacian subject to the Dirichlet boundaryconditions in (1.2). To this end, we consider the pointwise definition of the fractional Laplacian.Substituting the ansatz (3.1) into the fractional Laplacian (1.4) and taking the boundary conditions(1.2) into account, we obtain( − ∆) α h u ( x , t ) = C d,α (cid:18) P . V . (cid:90) Ω (cid:98) u ( x , t ) − (cid:98) u ( y , t ) | x − y | d + α d y + (cid:90) Ω c (cid:98) u ( x , t ) − g ( y , t ) | x − y | d + α d y (cid:19) = ( − ∆) α (cid:98) u ( x , t ) + C d,α (cid:90) Ω c (cid:98) u ( y , t ) − g ( y , t ) | x − y | d + α d y , for x ∈ Ω . (3.2)Note that the integral term over Ω c comes from the extended Dirichlet boundary conditions. Hence,it appears only in the fractional cases. On the other hand, we can directly apply the operator − ∆to (3.1) and obtain the approximation in the classical cases, i.e., ( − ∆) h u ( x , t ) = − ∆ (cid:98) u ( x , t ) for x ∈ Ω. Combining the classical and fractional cases, we obtain a unified approximation:( − ∆) α h u ( x , t ) = ( − ∆) α (cid:98) u ( x , t ) + ζ α C d,α (cid:90) Ω c (cid:98) u ( y , t ) − g ( y , t ) | x − y | d + α d y , for x ∈ Ω , α ∈ (0 , , (3.3)where ζ α = 1 − (cid:98) α/ (cid:99) with (cid:98)·(cid:99) being the floor function. By Lemma 2.1 and properties (2.6)–(2.7),it is easy to get( − ∆) α (cid:98) u ( x , t ) = 2 − d √ π Γ( d + α )Γ( d/ d + 1) / (cid:124) (cid:123)(cid:122) (cid:125) c d,α ¯ N (cid:88) i =1 λ i ( t ) | ε i | α F (cid:16) d + α d + 1 + α d − ε αi | x − x i | (cid:17) , (3.4)for x ∈ Ω and α ∈ (0 , c d,α for notational simplicity, which isdifferent from C d,α in the definition (1.4). It is evident that for α ∈ (0 , R d \ ¯Ω. Remark 3.1.
The scheme (3.3) – (3.4) can be straightforwardly applied to discretize the operator ( − ∆) m with m ∈ N . It can be also generalized to approximate the operator ( − ∆) α with α > but α (cid:54) = 2 m , but this generalization requires the point-wise definition of ( − ∆) α , such as the definition(1.4) for α ∈ (0 , , which is beyond the scope of this work. Choose test points x k ∈ ¯Ω for 1 ≤ k ≤ ¯ N , that is, the total number of test points is the same asthat of center points, but test points may not necessarily come from the same set of center points.More discussion on the choices of RBF center and test points can be found in [17, 38, 19]. Withoutloss of generality, we assume that test point x k ∈ Ω if index 1 ≤ k ≤ M , while x k ∈ ∂ Ω if index7 + 1 ≤ k ≤ ¯ N . For point x k ∈ Ω, the direct application of scheme (3.3)–(3.4) to diffusion problem(1.1) yields the semi-discretization form as: ¯ N (cid:88) i =1 dλ i ( t ) dt ϕ ε ( | x k − x i | ) = − κ ¯ N (cid:88) i =1 λ i ( t ) (cid:20) c d,α | ε i | α F (cid:16) d + α , d + α + 12 ; d − ε i | x k − x i | (cid:17) + ζ α C d,α (cid:90) Ω c ϕ ε ( | y − x i | ) | x k − y | d + α d y (cid:21) + ζ α C d,α (cid:90) Ω c g ( y , t ) | x k − y | d + α d y + f (cid:0) x k , t (cid:1) , t > , (3.5)for k = 1 , , . . . , M . While for test point x k ∈ ∂ Ω, we can combine (1.2) and (3.1) to obtain thediscretization of boundary conditions as: ¯ N (cid:88) i =1 λ i ( t ) ϕ ε ( | x k − x i | ) = g ( x k , t ) , t ≥ . (3.6)for k = M + 1 , M + 2 , . . . , ¯ N . It shows that boundary discretization (3.6) is only carried out forpoints on ∂ Ω (instead of on Ω c ), same for both classical and fractional cases. If α <
2, the boundaryconditions in (1.2) for x ∈ R d \ ¯Ω have been already applied to scheme (3.5) via its last integral.The initial condition at time t = 0 can be discretized as: ¯ N (cid:88) i =1 λ i (0) ϕ ε ( | x k − x i | ) = u ( x k , , for k = 1 , , . . . , ¯ N . (3.7)The scheme in (3.5)–(3.7) provides an ODE system for unknowns λ i ( t ) for 1 ≤ i ≤ ¯ N , which canbe solved by any time stepping method.In the following, we will discretize the ODE system (3.5)–(3.7) by the Crank–Nicolson method.We first denote vector λ ( t ) = (cid:0) λ ( t ) , λ ( t ) , · · · , λ ¯ N ( t ) (cid:1) T and rewrite system (3.5)–(3.6) into matrix-vector form:Φ (1: M,
1: ¯ N ) d λ ( t ) dt = − κA M × ¯ N λ ( t ) + b ( t ); Φ ( M +1: ¯ N,
1: ¯ N ) λ ( t ) = g ( t ) . (3.8)Here, Φ = (cid:8) Φ k,i (cid:9) ¯ N × ¯ N is a matrix of GIMQ basis functions with entry Φ k,i = ϕ ε (cid:0) | x k − x i | (cid:1) , andΦ ( r : r , c : c ) denotes its submatrix including row r to r and column c to c from matrix Φ. While A represents the differentiation matrix of ( − ∆) α with its entry A k,i = c d,α | ε i | α F (cid:16) d + α , d + α + 12 ; d − ε i | x k − x i | (cid:17) + ζ α C d,α (cid:90) Ω c ϕ ε ( | y − x i | ) | x k − y | d + α d y , for 1 ≤ k ≤ M and 1 ≤ i ≤ ¯ N . The vector g ( ¯ N − M ) × = (cid:0) g ( x M +1 , t ) , g ( x M +2 , t ) , . . . , g ( x ¯ N , t ) (cid:1) T ,and vector b M × has entry b k ( t ) = f ( x k , t ) + ζ α C d,α (cid:90) Ω c g ( y , t ) | x k − y | d + α d y , for k = 1 , , . . . , M. Denote time sequence t n = nτ (for n = 0 , , . . . ) with time step τ >
0. Then the Crank–Nicolsondiscretization of (3.8) yield the fully discrete scheme as (cid:32) Φ (1: M,
1: ¯ N ) + κτ A Φ ( M +1: ¯ N,
1: ¯ N ) (cid:33) λ n +1 = (cid:32) (cid:16) Φ (1: M,
1: ¯ N ) − κτ A (cid:17) λ n + τ (cid:0) b ( t n ) + b ( t n +1 ) (cid:1) g ( t n +1 ) (cid:33) , (3.9)8or n = 0 , , . . . , where λ n represents the numerical approximation of λ ( t n ). Initially, we can obtain λ by solving Φ ¯ N × ¯ N λ = U with U = (cid:0) u ( x , , u ( x , , . . . , u ( x ¯ N , (cid:1) T . Substituting λ n into(3.1) leads to the numerical solution of (1.1)–(1.2) at time t n and for any point x ∈ Ω.Since GIMQ functions are globally supported, the linear system (3.9) has a full stiffness matrixfor both classical ( α = 2) and fractional ( α <
2) problems. On the other hand, the nonlocality ofthe fractional Laplacian always leads to a linear system with full matrix even if it is discretized bylocal (e.g., finite difference/element) methods [13, 14, 1]. Hence, our method does not introduceextra computations in the fractional cases. Instead, it could save computational cost by using fewernumber of points to achieve the desired accuracy. This suggests that global methods might be morebeneficial for solving nonlocal or fractional problems.
In this section, we will test the performance of our method and compare it with the Gaussian RBF-based method recently proposed in the literature [8]. Furthermore, we will study two approachesof selecting the shape parameter ε for better performance of our method. Here, we will focus onthe one-dimensional Poisson problem:( − ∆) α u = f ( x ) , for x ∈ Ω; u ( x ) = g ( x ) , for x ∈ Υ . (4.1)Choose the domain Ω = ( − , α = 2 a two-point Dirichletboundary condition is imposed at x = ±
1, while an extended Dirichlet boundary condition isapplied on Ω c = ( −∞ , − ∪ [1 , ∞ ) if α <
2. For the purpose of easy comparison, we will choosethe RBF center points uniformly distributed on [ − ,
1] and test points from the same set of centerpoints (consequently, M = N ). Note that even though the extended Dirichlet boundary conditionsare imposed on Ω c in the fractional cases, the RBF points are only considered on ¯Ω. We remarkthat our method is flexible in choosing the center and test points, and the above choice is only anexample that we use here. Unless otherwise stated, we will use a constant shape parameter, i.e., ε i ≡ ε for 1 ≤ i ≤ ¯ N , in the following simulations.Let u j and u hj represent the exact solution and numerical approximation of u ( x ) at point x = x j ,respectively. In the following, we compute numerical errors as the root mean square (RMS) errors: (cid:107) e (cid:107) rms = (cid:18) K K (cid:88) j =1 | u j − u hj | (cid:19) / , where K (cid:29) ¯ N denotes the total number of interpolation points on domain ¯Ω. In practice, asufficiently large number K is chosen such that the numerical error (cid:107) e (cid:107) rms is insensitive to it.As we will see later, our method not only unifies the numerical discretization of the fractional(0 < α <
2) and classical ( α = 2) Laplacian in a single scheme, but also allows simple computerimplementation for any dimension d ≥ We will test the performance of our method in solving the Poisson problem (4.1) with differentboundary conditions. Most of the existing studies on fractional PDE problems focus on homoge-neous boundary conditions (e.g. g ( x ) ≡ .1.1 Homogeneous boundary conditions Here, we consider a benchmark (fractional) Poisson problem with homogeneous boundary condi-tions, i.e., g ( x ) ≡
0, and choose function f in (4.1) as f ( x ) = 2 α Γ (cid:16) α +12 (cid:17) Γ( s + 1 + α ) √ π Γ( s + 1) F (cid:16) α + 12 , − s, , x (cid:17) , for s ∈ N . For any α ∈ (0 , u ( x ) = (1 − x ) s + α + , acompact support function on [ − , s = 3.Table 2 demonstrates the RMS errors and condition numbers K of the resultant linear system,where a constant shape parameter ε is used for all RBF center points. It shows that the numericalerrors decrease with a spectral rate as the number of points ¯ N increases. Our method can achievea good accuracy even with fewer points. For instance, it has an error of O (10 − ) ∼ O (10 − )for ¯ N = 65, which is much smaller than the errors from finite difference/element methods in theliterature [13, 1]. Moreover, our method has a distinct merit – unifying the approximation scheme α = 0 . ε = 3 α = 1, ε = 3 . α = 1 . ε = 3 . α = 2, ε = 3 . N (cid:107) e (cid:107) rms K (cid:107) e (cid:107) rms K (cid:107) e (cid:107) rms K (cid:107) e (cid:107) rms K (cid:107) e (cid:107) rms and condition numbers K in solving the 1D Poisson problem (4.1),where exact solution is u ( x ) = (1 − x ) α + .for the fractional and classical Laplacians in a single scheme. Additionally, we find that as ¯ N increases, the condition number of the differentiation matrix increases quickly, which indicates thatthe system could become ill-conditioned when a large number of points is used. Different strategiescan be found in the literature to resolve the ill-conditioning issues of RBF-related methods; see[37, 38, 19] and references therein. In Section 4.3, we will study two approaches to suppress thecondition number via shape parameter ε .In Fig. 3, we further study the relation between numerical errors (cid:107) e (cid:107) rms , condition number K , total number of points ¯ N , and the constant shape parameter ε . We find that for a given ¯ N ,numerical errors generally decrease first and then increase with the shape parameter ε . There existsan optimal shape parameter ε ∗ ( ¯ N ), depending on the total number of points and solution property,at which the numerical error is minimized. The condition number could be very large if the numberof points ¯ N is large or the shape parameter ε is small. If ε is large enough, the condition numberwould monotonically decrease as ε increases. This implies that the ill-conditioning issue caused bya large number of points could be relaxed by increasing the shape parameter, but a large shapeparameter might reduce numerical accuracy (see Fig. 3 left panel).10 .5 2.5 4.5 6.5-9-7-5-3 Figure 3: Numerical errors (cid:107) e (cid:107) rms and condition numbers K versus the constant shape parameter ε in solving the 1D Poisson problem (4.1) with exact solution u ( x ) = (1 − x ) α + . In the literature, the results on fractional PDEs with nonhomogeneous Dirichlet boundary condi-tions still remain rare. In this example, we consider the Poisson problem (4.1) with f ( x ) = √ α ) √ π F (cid:16) α + 12 ; 3 + α ,
12 ; − x (cid:17) , g ( x ) = (cid:114) π sinc (cid:16) xπ (cid:17) . Its exact solution is given by u ( x ) = (cid:112) /π sinc( x/π ) for any x ∈ R and α ∈ (0 , α .We find similar observations as in Table 2 – numerical errors decrease with a spectral rate as¯ N increases. But, the condition numbers in Table 3 are generally larger because of smaller shapeparameters ε are used here. The computational cost in this example is higher due to extra efforts toevaluate the integrals of nonzero boundary conditions. Fig. 4 demonstrates the relation of numericalerrors (cid:107) e (cid:107) rms , conditional number K , total number of points ¯ N , and the shape parameter ε . Similarto observations in Fig. 3, numerical errors first decrease and then increase with the shape parameter ε , but the optimal shape parameter in this example is less sensitive to the number of points ¯ N .The minimum errors are achieved at a shape parameter ε ∗ ( ¯ N ) around 0 . ∼ . N and α (see Fig. 4 left panel). Notice that the exact solution u ∈ C , α ( ¯Ω) in Fig. 3, while u ∈ C ∞ ( ¯Ω)in Fig. 4, suggesting that solution regularity might play an important role in determining the11 = 0 . ε = 1 α = 1, ε = 1 α = 1 . ε = 1 . α = 2, ε = 1 . N (cid:107) e (cid:107) rms K (cid:107) e (cid:107) rms K (cid:107) e (cid:107) rms K (cid:107) e (cid:107) rms K (cid:107) e (cid:107) rms and condition numbers K in solving the 1D Poisson problem (4.1),where exact solution is u ( x ) = (cid:112) /π sinc( x ) for x ∈ R . Figure 4: Numerical errors (cid:107) e (cid:107) rms and condition numbers K versus the shape parameter ε in solvingthe 1D Poisson problem (4.1) with exact solution u ( x ) = (cid:112) /π sinc( x ) for x ∈ R .optimal shape parameters. The above studies also suggest that it is challenging to find a uniformoptimal shape parameter ε for different numbers of points ¯ N . In practice, it is important to selecta shape parameter that can yield good numerical accuracy and avoid ill-conditioning issues in thecomputation; see more discussion in Section 4.3.12 .2 Comparison to the Gaussian-based method Recently, a Gaussian-based method was proposed in [8] to discretize the classical and fractionalLaplacians. As noted previously, the Gaussian and GIMQ functions share many similarities, i.e.,globally supported, infinitely differentiable, and strictly positive definite. Both of them have beenwell studied in interpolation problems [21, 20]. However, very few comparison of these two functionscan be found in the literature, especially for solving PDEs. In the following, we will compare ourGIMQ-based method with the Gaussian-based method in [8] by studying their numerical errorsand condition numbers with respect to constant shape parameters.Here, we will mainly focus on the Poisson problems as studied in Section 4.1.1, and the constantshape parameters are used for all RBF center points. Our studies show that the Gaussian RBF-based method has a spectral accuracy. To compare with GIMQ results in Fig. 3, we present in Fig.5 the relation between numerical errors, condition number, number of points, and shape parametersof the Gaussian RBF-based method. It shows that numerical errors of both methods first decrease
Figure 5: Gaussian RBF-based method: Numerical errors (cid:107) e (cid:107) rms and condition numbers K versusthe shape parameter ε in solving the 1D Poisson problem with exact solution u ( x ) = (1 − x ) α + .These results are in comparison with those in Fig. 3 for the GIMQ-based method.and then increase with the constant shape parameter ε , but the errors of GIMQ method decreasewith less oscillation. There also exists an optimal constant shape parameter ε ∗ = ε ∗ ( ¯ N ) for theGaussian-based method. These two methods are different mainly in two aspects: (i) under thesame numerical setting, the optimal shape parameter ε ∗ of the Gaussian method is much larger.13lso, it is more sensitive to the number of points ¯ N . (ii) The condition numbers of both methodsmonotonically decrease after a large shape parameter (see the right panels of Figs. 3 and 5).But, the GIMQ-based method starts this monotonic decay at a smaller shape parameter, but thecondition number of the Gaussian-based method remain large for a wide range of ε . For instance,when α = 0 . N = 129, the condition number reduces to K ∼ at ε = 5 in GIMQ-basedmethod (see Fig. 3), while ε = 16 .
11 41 71 101 1310.51.62.73.84.96
11 41 71 101 13118152229
11 41 71 101 131-9-7-5-3 11 41 71 101 131-9-7-5-3
11 41 71 101 131481216
11 41 71 101 131261014
Figure 6: Comparison of the GIMQ-based method (left) and Gaussian-based method in [8] (right)for solving the 1D Poisson problem (4.1) with exact solution u ( x ) = (1 − x ) α + .14n Fig. 6, we further compare these two methods by studying their optimal shape parameter ε ∗ ( ¯ N ), and the corresponding numerical errors and condition numbers. It shows that for bothmethod, the optimal shape parameter increases almost linearly with number ¯ N , but the optimalshape parameters of Gaussian-based method are much larger and sensitive to power α . Moreprecisely, the larger the power α , the faster the parameter ε ∗ increases with respect to ¯ N . Thenumerical errors at ε ∗ are similar for both methods, but the condition numbers of Gaussian-basedmethods remain almost a constant if ε ∗ is used. These comparisons could provide some insights forpractical applications of these two methods. As we noted, the shape parameter ε plays an important role in determining the condition numberand numerical accuracy of RBF-based methods. If a small constant shape parameter (i.e., ε i ≡ ε for1 ≤ i ≤ ¯ N ) is used, the resulting linear system could easily become ill-conditioned as the numberof points ¯ N increases. While using a large shape parameter could suppress the condition number,but at the same time it might deteriorate numerical accuracy if the shape parameter is too big.It was pointed out in [34] that the optimal shape parameter depends on many factors, includingthe number of points, distribution of points, properties of functions to approximate, and evencomputer precision. More discussion can be found in [41, 38, 19] and references therein. Clearly, itis challenging to find the optimal shape parameters in practice.In the following, we study two approaches for choosing shape parameters with the goal ofbalancing numerial accuracy and condition number. In the first approach, we consider a constantshape parameter for all RBF center points, which is selected by controlling the condition number ina desired range and thus eliminates the ill-conditioning issues. In the second approach, we consider(randomly) different shape parameters for each RBF points. Similar approaches can be found in[38], but no studies have been reported in solving fractional PDEs. Condition number indicated shape parameter (denoted as ε K ). Our studies show thatwhen a constant shape parameter is used, numerical errors reduce as the number of points increases,and at the same time the condition number quickly increases (e.g. see Tables 2–3). Consequently,the system could become ill-conditioned if the number of points is too big. On the other hand, wefind when the number ¯ N is large, numerical solutions tend to be more accurate if the conditionnumber is large but before the system becomes ill-conditioned. Inspired by these observations, wewill choose a constant shape parameter such that the resulting condition number falls in a desiredrange, namely condition number indicated shape parameter . The constant shape parameter selectedfrom this approach will be denoted as ε K .In Fig. 7, we compare numerical errors from different constant shape parameters in solving thePoisson problem (4.1) with exact solution u ( x ) = (1 − x ) α + . To determine ε K , we start with asmall value of ε and then increase it gradually until the condition number of the system falls into therange of 10 ≤ K ≤ , and the corresponding shape parameter will be defined as ε K . Note thatthe shape parameter selected in this approach changes for different number ¯ N , i.e. ε K = ε K ( ¯ N ).From Fig. 7, we find that a constant shape parameter (e.g. ε = 1 or 3) eventually leads tothe growth of numerical errors when the number of points ¯ N is large enough because the systembecomes ill-conditioned. Moreover, a large shape parameter always lead to bigger numerical errorseven though the system is well-conditioned (see, e.g., ε = 5). In contrast, the K -indicated shapeparameter ε K can well balance the condition number and numerical accuracy. It can prevent thesystem from being ill-conditioned and ensure small numerical errors as ¯ N increases. This propertyis important for practical simulations. 15
50 70 90 110 130-8-6-4-2
Figure 7: Comparison of numerical errors for different shape parameters in solving the 1D Poissonproblem (4.1) with exact solution u ( x ) = (1 − x ) α + . Random-perturbed variable shape parameter (denoted as ε δ ). The condition numberindicated shape parameter ε K can effectively avoid the ill-conditioning issues in simulations withlarge ¯ N . Next, we will study another approach – random-perturbed variable shape parameter . Itperturbs a constant shape parameter with uniformly distributed random numbers at RBF centerpoints, and thus the shape parameter can be viewed as a random function of center point x i [39, 43].Let [ ε min , ε max ] denote the interval that the shape parameter is allowed. The random-perturbedshape parameter ε δ,i for center point x i is given by: ε δ,i = ε min + δ i ( ε max − ε min ) , (4.2)where δ i is a uniformly distributed random number on the interval (0,1).In Fig. 7, we test the performance of random-perturbed variable shape parameters ε δ in com-parison to constant parameters ε , where the Poisson problem (4.1) with exact solution u ( x ) =(1 − x ) α + is studied. We choose ε min = 1 and ε max = 5 in (4.2). It shows that two simula-tions with different ε δ give similar numerical errors and condition numbers. Compared to constantshape parameters, using distinct shape parameters for each RBF center point helps suppressingthe condition number and leads to good numerical errors for different number ¯ N . It can also savethe computational time in searching for a proper shape parameter (e.g., ε K ). Here, the choices of[ ε min , ε max ] is important in determining random-perturbed variable shape parameters. Our studiesin Section 4.2 suggest that this range might depend on the minimum distance between RBF centerpoints, which we will further explore in our future study. In this section, we will focus on our GIMQ-based meshless method and apply it to study bothelliptic problems and diffusion equations. In our simulations, we will use the random-perturbedshape parameters as described in (4.2), and the test points are chosen from the same set of RBFcenter points (consequently, M = N in (3.9)). 16
30 55 80 105 130-9-6-30
Figure 8: Comparison of numerical errors and condition numbers for different shape parameters insolving the 1D Poisson problem (4.1) with exact solution u ( x ) = (1 − x ) α . We consider the following elliptic problem on an L-shaped domain, i.e., Ω = { ( x, y ) ∈ ( − , \ [0 , } ,with nonhomogeneous Dirichlet boundary conditions:( − ∆) α u + 2 u = f, for x ∈ Ω; u ( x ) = e −| x | , for x ∈ Υ . (5.1)The right hand side function f in (5.1) is chosen as: f ( x ) = 2 α Γ(1 + α F (cid:16) α −| x | (cid:17) + 2 e −| x | (5.2)where F denotes the confluent hypergeometric function, and the exact solution is u = e −| x | for x ∈ R . We choose RBF center/test points as uniformly distributed grid points on domain ¯Ω, and ε min = 0 . ε max = 4 in the random-perturbed shape parameters (4.2).Table 4 shows numerical errors and condition numbers for various α . It is evident that numericalerrors decrease with a spectral rate as ¯ N increases. Compared to finite difference/element methods,our method can yield a more accurate result with fewer number of points. For example, to achievean accuracy of O (10 − ), our method only requires the number of points ¯ N = 341 (equivalently,distance between two grid points h = 0 . = 0 . α = 1 α = 1 . α = 2¯ N RMS K RMS K RMS K RMS K
21 6.238e-2 8.909e2 4.888e-2 1.147e2 3.607e-2 1.288e2 2.339e-2 2.036e265 3.795e-3 2.227e7 8.403e-3 7.030e6 4.834e-3 2.901e6 9.294e-3 1.794e5133 2.892e-4 2.420e10 3.976e-4 1.234e10 1.779e-4 9.662e9 1.702e-5 1.199e10225 3.899e-6 4.069e11 1.367e-6 7.021e11 4.178e-5 8.393e14 2.626e-6 1.236e11341 3.902e-8 1.845e15 1.336e-8 3.219e15 4.239e-8 3.178e17 2.159e-8 1.481e17Table 4: Numerical errors (cid:107) e (cid:107) rms and condition number K in solving the elliptic problem (5.1) onan L-shaped domain. We solve the diffusion problem with f ≡ κ = 0 . − , \ [0 . , . , and the initial condition u ( x ,
0) = 3 e − | x | , for x ∈ R . We then compare normal ( α = 2) and anomalous ( α <
2) diffusion by studying time evolution ofthe solution u ( x , t ). In our simulations, RBF center/test points are chosen as uniformly distributedgrid points on ¯Ω. The random-perturbed shape parameters are used with ε min = 1 and ε max = 5.Fig. 9 shows the initial solution u ( x ,
0) and dynamics of solution norm (cid:107) u ( · , t ) (cid:107) = (cid:16) (cid:90) Ω | u ( x , t ) | d x (cid:17) , t ≥ α . Here, we choose time step τ = 0 . N = 600.We have verified and ensured the convergence of solutions with smaller time step and more RBFpoints. It shows that initially the radially-symmetric solution concentrates around x = with norm Figure 9: (a) Initial condition u ( x ,
0) of the diffusion problem; (b) Dynamics of (cid:107) u ( · , t ) (cid:107) in normaland anomalous diffusion. (cid:107) u ( · , (cid:107) = 23 .
67. Then norm (cid:107) u ( · , t ) (cid:107) monotonically decreases over time owing to zero sourceterm (i.e., f = 0) and homogeneous Dirichlet boundary conditions. It shows that the smaller thepower α , the slower the decay of solution norm (cid:107) u ( · , t ) (cid:107) .18ig. 10 further compares the solution evolution at different time t . It shows that solution u ( x , t )diffuses from domain center towards boundaries over time, and it starts deforming once reachingthe inner boundaries (i.e., boundary of region [0 . , . ). It is clear that the smaller the powerFigure 10: Time evolution of solution in normal and anomalous diffusion problems. α , the slower the solution diffuses, consistent with the slower decay of norm (cid:107) u ( · , t ) (cid:107) . At time t = 1 .
5, the solution of α = 1 . x = ± y = ± α = 0 . t = 1 .
5, but is still far away fromthe outer boundaries until t = 3. In this example, we continue our study on normal and anomalous diffusion and explore the boundaryeffects in classical and fractional PDEs. Consider the problem ∂ t u ( x , t ) = − κ ( − ∆) α u + u, for x ∈ Ω , t > , (5.3)where the domain Ω = ( − , . Let the initial condition u ( x , ≡ x ∈ ¯Ω. The boundaryconditions are set as u ( x , t ) = (cid:26) sin (cid:2) π ( x − x c + ) (cid:3) sin (cid:2) π ( y + 1) (cid:3) , for x ∈ Ξ , , for x ∈ Υ \ Ξ , t ≥ , x c , x c + ] × [ − ,
1] for x c ≥
1. In other words, boundary conditionsare zero everywhere except on region Ξ, and the location of region Ξ is controlled by the value of x c . When x c = 1, region Ξ connects to domain Ω, while if x c > x c .Fig. 11 shows the dynamics of solution u ( x , t ) for x c = 1, and α = 1 . ,
2. Initially, the solution u ( x ,
0) = 0. Due to nonzero boundary conditions on Ξ = [1 , . × [ − , x = 1 quickly increases and simultaneously diffuses into domain Ω. It remains symmetric withrespect to y = 0 for any time >
0, consistent with the symmetry of boundary conditions. NoteFigure 11: Numerical solution of heat equation (5.3) at different time, where x c = 1.that when α = 2, the diffusion operator − ∆ is local, and the effective nonzero boundary conditionsoccur only at line x = 1. Fig. 12 further compares the solution evolution at y = 0 for differenttime t . It shows that the solution of classical cases diffuse fast at the initial stage, and it quickly -1 -0.5 0 0.5 100.250.50.751 -1 -0.5 0 0.5 100.250.50.751 Figure 12: Time snapshot of u ( x, , t ) for α = 1 . α = 2 (b).approaches the steady state. 20ext, we continue our study of boundary effects in Fig. 13, for x c = 1 . α = 0 . , .
4. Inthis case, the nonzero boundary conditions are imposed on region Ξ = [1 . , . × [ − ,
1] which hasno contact with the domain Ω. Hence, the solution of classical ( α = 2) problems remain u ( x , t ) ≡ t ≥
0, due to the homogeneous boundary conditions on ∂ Ω. In contrast to classical cases,Figure 13: Numerical solution of heat equation (5.3) at different time, where x c = 1 . α , the stronger the boundary effects. Com-paring the results of α = 1 . x c −
1) is smaller. These two phenomena canbe explained by the kernel function | x − y | − ( d + α ) of the fractional Laplacian in (1.4), whose valueis affected by both distance | x − y | and power α . The smaller the power α (or the shorter the dis-tance | x − y | ), the larger the kernel function, and thus the stronger the interactions from boundaryconditions. Our study might provide insights for boundary control in the study of fractional PDEs[4, 12]. We proposed a new meshless pseudospectral method based on the generalized inverse multiquadric(GIMQ) functions to solve problems with ( − ∆) α for 0 < α ≤
2. The operator ( − ∆) α representsthe local classical Laplacian − ∆ when α = 2, while the nonlocal fractional Laplacian if α < − ∆) α canbe defined via the pseudodifferential operator in (1.3) or the pointwise hypersingular integral form in211.4). We combined the advantages of these two definitions and proposed an α -parametric methodfor the Laplacian ( − ∆) α with α ∈ (0 , β = − ( d + 1) / d ≥
1. Notice that the Laplacian of GIMQ functions canbe analytically written by the hypergeometric function for any α ≥
0. This provides the keyto avoiding numerical evaluations to the hypersingular integral and unifying approximations ofclassical and fractional Laplacians. Our studies of the diffusion problems illustrated the differencesbetween the normal ( α = 2) and anomalous ( α <
2) diffusion. It showed that boundary conditionsof anomalous diffusion problems can affect the solution in domain Ω even though they may nothave direct contact, in contrast to the normal diffusion where the boundary conditions have tobe imposed on ∂ Ω. These studies could provide insights of simulations and applications with thefractional models.We compared our method with the recent Gaussian-based method in [8] and found that bothmethods have the spectral accuracy. Furthermore, we found that if uniformly distributed RBF cen-ter points are considered, the optimal shape parameter of Gaussian-based method is more sensitiveto the number of points. This suggests that different strategies might be considered in selectingthe shape parameters of Gaussian and GIMQ based methods. Two approaches in selecting shapeparameters of our GIMQ-based method were studied in detail. Numerical studies showed that thecondition-number indicated shape parameter can effectively suppress the ill-conditioning issues andthus guarantee the numerical accuracy when a larger number of points are used. This method maytake extra computational time in searching a proper shape parameter to satisfy the desired con-dition number range. While the random-perturbed shape parameters can save the time of searching.
Acknowledgements.
This work was supported by the US National Science Foundation underGrant Number DMS-1913293 and DMS-1953177.
References [1] 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.[2] G. Acosta, J. P. Borthagaray, and N. Heuer. Finite element approximations of the nonhomo-geneous fractional Dirichlet problem.
IMA J. Numer. Anal. , 39(3):1471–1501, 2019.[3] M. Ainsworth and C. Glusa.
Towards an efficient finite element method for the integral frac-tional Laplacian on polygonal domains . Contemporary computational mathematics—a cele-bration of the 80th birthday of Ian Sloan. Vol. 1, 2. Springer, Cham, 2018.[4] H. Antil, R. Khatri, and M. Warma. External optimal control of nonlocal PDEs.
InverseProblems , 35(8):084003, 35, 2019.[5] G. S. Bhatia and Arora G. Radial basis function methods for solving partial differentialequations-a review.
Indian J. Sci. Technol. , 9, 2016.[6] S. D. Bond, R. B. Lehoucq, and S. T. Rowe.
A Galerkin Radial Basis Function Method forNonlocal Diffusion . Lect. Notes Comput. Sci. Eng. Springer, 2015.227] A. Bonito, W. Lei, and J. E. Pasciak. Numerical approximation of the integral fractionalLaplacian.
Numer. Math. , 142(2):235–278, 2019.[8] J. Burkardt, Y. Wu, and Y. Zhang. A unified meshfree pseudospectral method for solvingboth classical and fractional pdes. arXiv:1201.2844v1 , 2020.[9] M. E. Chenoweth and S. Sarra. A numerical study of generalized multiquadric radial basisfunction interpolation.
SIAM Undergrad. Res. Online , 2:58–70, 01 2009.[10] A. Cusimano, N.and Bueno-Orovio, I. Turner, and K. Burrage. On the order of the fractionallaplacian in determining the spatio-temporal evolution of a space-fractional model of cardiacelectrophysiology.
PloS one , 10:1–16, 12 2015.[11] D. del Castillo-Negrete, B. A. Carreras, and V. E. Lynch. Front dynamics in reaction-diffusionsystems with l´evy flight: A fractional diffusion approach.
Phys. Rev. Lett. , 91(1):018302, 2003.[12] M. D’Elia, C. Glusa, and E Ot´arola. A priori error estimates for the optimal control of theintegral fractional laplacian.
SIAM J. Control Optim. , 54(4):2775–2798, 2019.[13] S. Duo, H. W. van Wyk, and Y. Zhang. A novel and accurate finite difference method forthe fractional Laplacian and the fractional Poisson problem.
J. Comput. Phys. , 355:233–252,2018.[14] S. Duo and Y. Zhang. Accurate numerical methods for two and three dimensional integralfractional Laplacian with applications.
Comput. Methods. Appl. Mech. Eng. , 355:639–662,2019.[15] S. Duo and Y. Zhang. Numerical approximations for the tempered fractional Laplacian: Erroranalysis and applications.
J. Sci. Comput. , 81(1):569–593, 2019.[16] B. Dyda. Fractional calculus for power functions and eigenvalues of the fractional Laplacian.
Fract. Calc. Appl. Anal. , 15(4):536–555, 2012.[17] B. Fornberg, T. A. Driscoll, G. Wright, and R. Charles. Observations on the behavior of radialbasis function approximations near boundaries.
Comput. Math. Appl. , 43(3-5):473–490, 2002.[18] B. Fornberg and N. Flyer.
A primer on radial basis functions with applications to the geo-sciences , volume 87 of
CBMS-NSF Regional Conference Series in Applied Mathematics . Societyfor Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2015.[19] B. Fornberg and N. Flyer. Solving PDEs with radial basis functions.
Acta Numer. , 24:215–258,2015.[20] B. Fornberg and C. Piret. On choosing a radial basis function and a shape parameter whensolving a convective PDE on a sphere.
J. Comput. Phys. , 227(5):2758–2780, 2008.[21] R. Franke.
A Critical Comparison of Some Methods for Interpolation of Scattered Data . Tech-nical Report NPS-53-79-03. Naval Postgraduate School, 1979.[22] R. Franke. Scattered data interpolation: Tests of some method.
Math. Comp. , 38(157):181–200, 1982.[23] K. Hamm and J. Ledford. Cardinal interpolation with general multiquadrics: convergencerates.
Adv. Comput. Math. , 44(4):1205–1233, 2018.2324] E. Hanert, E. Schumacher, and E. Deleersnijder. Front dynamics in fractional-order epidemicmodels.
J. Theor. Biol. , 279:9–16, 2011.[25] R. Hardy. Multiquadric equations of topography and other irregular surfaces.
J. Geophys.Res. , 76(3):1905–1915, 1971.[26] R.L. Hardy. Theory and applications of the multiquadric-biharmonic method.
J. Geophys.Res. , 19(8):163–208, 1990.[27] E. J. Kansa. Multiquadrics - a scattered data approximation scheme with applications to com-putational fluid-dynamics. I. Surface approximations and partial derivative estimates.
Comput.Math. Appl. , 19(8–9):127–145, 1990.[28] E. J. Kansa. Multiquadrics - a scattered data approximation scheme with applications to com-putational fluid-dynamics. II. Solutions to parabolic, hyperbolic and elliptic partial differentialequations.
Comput. Math. Appl. , 19(8–9):147–161, 1990.[29] C. Ku, C. Liu, J. Xiao, and S. Hsu. Multiquadrics without the shape parameter for solvingpartial differential equations.
Symmetry , 12(11), 2020.[30] M. Kwa´snicki. Ten equivalent definitions of the fractional Laplace operator.
Fract. Calc. Appl.Anal. , 20(1):7–51, 2017.[31] N. S. Landkof.
Foundations of Modern Potential Theory . Springer-Verlag, New York-Heidelberg, 1972.[32] J. W. Pearson, S. Olver, and M. A. Porter. Numerical methods for the computation of theconfluent and Gauss hypergeometric functions.
Numer. Algorithms , 74(3):821–866, 2017.[33] C. Piret and E. Hanert. A radial basis functions method for fractional diffusion equations.
J.Comput. Phys. , 238:71–81, 2013.[34] S. Rippa. An algorithm for selecting a good value for the parameter c in radial basis functioninterpolation. Adv. Comput. Math. , 11(2-3):193–210, 1999.[35] J. A. Rosenfeld, S. A. Rosenfeld, and W. E. Dixon. A mesh-free pseudospectral approach toestimating the fractional Laplacian via radial basis functions.
J. Comput. Phys. , 390:306–322,2019.[36] S. G. Samko, A. A. Kilbas, and O. I. Marichev.
Fractional Integrals and Derivatives . Gordonand Breach Science Publishers, Yverdon, 1993.[37] S. A. Sarra. Radial basis function approximation methods with extended precision floatingpoint arithmetic.
Eng. Anal. Bound. Elem. , 35(1):68–76, 2011.[38] S. A. Sarra and E. J. Kansa. Multiquadric radial basis function approximation methods forthe numerical solution of partial differential equations.
Advances in Computational Mechanics ,2, 01 2009.[39] S. A. Sarra and D. Sturgill. A random variable shape parameter strategy for radial basisfunction approximation methods.
Eng. Anal. Bound. Elem. , 33(11):1239–1245, 2009.2440] F. Soleymani, M. Barfeie, and F. K. Haghani. Inverse multi-quadric RBF for computing theweights of FD method: application to American options.
Commun. Nonlinear Sci. Numer.Simul. , 64:74–88, 2018.[41] C. H. Tsai, J. Kolibal, and M. Li. The golden section search algorithm for finding a goodshape parameter for meshless collocation methods.
Eng. Anal. Bound. Elem. , 34(8):738–746,2010.[42] H. Wendland.
Scattered data approximation , volume 17 of
Cambridge Monogr. Appl. Comput.
Cambridge University Press, Cambridge, 2005.[43] J. Wertz, E. J. Kansa, and L. Ling. The role of the multiquadric shape parameters in solvingelliptic partial differential equations.
Comput. Math. Appl. , 51(8):1335–1348, 2006.[44] J. Yoon. Spectral approximation orders of radial basis function interpolation on the Sobolevspace.