Numerical solution of space-fractional partial differential equations by a differential quadrature approach

TL;DR

This paper introduces a new differential quadrature method using radial basis functions to efficiently and accurately solve space-fractional PDEs, demonstrating robustness and high precision through numerical tests.

Contribution

It develops a novel RBF-based differential quadrature technique for space-fractional PDEs, enhancing accuracy and simplicity over existing methods.

Findings

Method is robust and straightforward to implement.

Achieves high accuracy with proper RBF shape parameters.

Numerical tests confirm effectiveness and validity.

Abstract

This article aims to develop a direct numerical approach to solve the space-fractional partial differential equations (PDEs) based on a new differential quadrature (DQ) technique. The fractional derivatives are approximated by the weighted linear combinations of the function values at discrete grid points on problem domain with the weights calculated via using three types of radial basis functions (RBFs) as test functions. The method in presence is robust, straight forward to apply, and highly accurate under the condition that the shape parameters of RBFs are well chosen. Numerical tests are provided to illustrate its validity and capability.

Tables

Table 1. Table 1: The numerical results when α = 1.1 𝛼 1.1 \alpha=1.1 for Example 5.1

$M$	MQ-DQM	IM-DQM	GA-DQM
5	6.4167e-02	5.7172e-02	1.9069e-01
10	6.3543e-03	8.5795e-03	2.2488e-02
15	1.3676e-03	7.2670e-04	8.5354e-04
20	4.2484e-04	4.0834e-04	1.1827e-04
25	2.2093e-04	8.2798e-05	4.6350e-06

Table 2. Table 2: The numerical results at t = 0.5 𝑡 0.5 t=0.5 when τ = 2.0 × 10 − 4 𝜏 2.0 superscript 10 4 \tau=2.0\times 10^{-4} , α = 1.8 𝛼 1.8 \alpha=1.8 for Example 5.2

$M$	MQ-DQM	IM-DQM	GA-DQM
5	2.2567e-04	2.0463e-04	2.2607e-04
10	1.5291e-05	1.2391e-05	7.3949e-06
20	4.1822e-07	1.9394e-07	1.4895e-08
25	7.6704e-08	2.9039e-08	2.0757e-09

Table 3. Table 3: A comparison of SAM and RBF-DQM at t = 1 𝑡 1 t=1 when τ = 1 / M 𝜏 1 𝑀 \tau=1/M and α = 1.5 𝛼 1.5 \alpha=1.5 .

$M$	SAM [26]	MQ-DQM	IM-DQM	GA-DQM
15	7.660e-04	1.5903e-04	1.4862e-04	1.4127e-04
20	4.493e-04	7.9355e-05	7.5052e-05	7.2001e-05
25	2.929e-04	4.6347e-05	4.5794e-05	4.5247e-05
30	2.067e-04	2.9290e-05	3.0308e-05	3.1650e-05

quadrature points		weights
-0.995332738871603	0.072181310006040	0.273056604980456	0.186173360311749
-0.958255704632838	0.262409004111262	0.270507150048226	0.165983681785832
-0.885483116532829	0.442850686515520	0.265432043865656	0.144244294127118
-0.779726461638614	0.606783228547583	0.257878671571417	0.121158215878927
-0.644926204206171	0.748098741650606	0.247917557394430	0.096941123533022
-0.486104964648995	0.861532419612587	0.235641706298296	0.071819606274043
-0.309180377381835	0.942860103381098	0.221165735833726	0.046030939495720
-0.120744600181947	0.989070420301884	0.204624806418771	0.019851626928800

Equations

$$\frac{\partial y(x,t)}{\partial t}-\kappa (x)\frac{{\partial }^{\alpha }y(x,t)}{{\partial }_{+}{x}^{\alpha }}-\upsilon (x)\frac{{\partial }^{\alpha }y(x,t)}{{\partial }_{-}{x}^{\alpha }}=f(x,t),x\in \Lambda ,0<t⩽T,$$

$$y(x,0)=\psi (x),x\in \Lambda ,$$

$$y(a,t)=g_1(t),y(b,t)=g_2(t),0<t⩽T,$$

$$\frac{{\partial }^{\alpha }y(x,t)}{{\partial }_{+}{x}^{\alpha }}$$

$$\frac{{\partial }^{\alpha }y(x,t)}{{\partial }_{-}{x}^{\alpha }}$$

$${}_a{D}_x^{\alpha }f(x)$$

$${}_x^{\ast }{D}_b^{\alpha }f(x)$$

$${}_a{D}_x^{\alpha }(x-a{)}^{\beta }=\frac{\Gamma (\beta +1)}{\Gamma (\beta -\alpha +1)}(x-a{)}^{\beta },{}_x^{\ast }{D}_b^{\alpha }(b-x{)}^{\beta }=\frac{\Gamma (\beta +1)}{\Gamma (\beta -\alpha +1)}(b-x{)}^{\beta },$$

$${}_0{D}_x^{s}f(x)=\frac{{\partial }^{s}f(x)}{\partial x^s},{}_x^{\ast }{D}_b^{s}f(x)=(-1{)}^s\frac{{\partial }^{s}f(x)}{\partial x^s},$$

$$y(x,t)\cong \sum _{k=0}^M{\delta }_k(t){\varphi }_k(x),$$

$$\frac{{\partial }^{\alpha }{\varphi }_k(x_i)}{{\partial }_{+}{x}^{\alpha }}$$

$$\frac{{\partial }^{\alpha }{\varphi }_k(x_i)}{{\partial }_{-}{x}^{\alpha }}$$

$$\frac{{\partial }^{\alpha }y(x_i,t)}{{\partial }_{+}{x}^{\alpha }}\cong \sum _{k=0}^M{\delta }_k(t)\frac{{\partial }^{\alpha }{\varphi }_k(x_i)}{{\partial }_{+}{x}^{\alpha }}=\sum _{k=0}^M{\delta }_k(t)\sum _{j=0}^M{a}_{ij}^{(\alpha )}{\varphi }_k(x_j)\cong \sum _{j=0}^M{a}_{ij}^{(\alpha )}y(x_j,t),$$

$$\frac{{\partial }^{\alpha }y(x_i,t)}{{\partial }_{-}{x}^{\alpha }}\cong \sum _{k=0}^M{\delta }_k(t)\frac{{\partial }^{\alpha }{\varphi }_k(x_i)}{{\partial }_{-}{x}^{\alpha }}=\sum _{k=0}^M{\delta }_k(t)\sum _{j=0}^M{b}_{ij}^{(\alpha )}{\varphi }_k(x_j)\cong \sum _{j=0}^M{b}_{ij}^{(\alpha )}y(x_j,t),$$

$$\frac{{\partial }^{\alpha }{\phi }_k(x_i)}{{\partial }_{+}{x}^{\alpha }}$$

$$\frac{{\partial }^{\alpha }{\phi }_k(x_i)}{{\partial }_{-}{x}^{\alpha }}$$

$$\displaystyle\underbrace{\left(\begin{array}[]{cccc}\varphi_{0}(x_{0})&\varphi_{0}(x_{1})&\cdots&\varphi_{0}(x_{M})\\ \varphi_{1}(x_{0})&\varphi_{1}(x_{1})&\cdots&\varphi_{1}(x_{M})\\ \vdots&\vdots&\ddots&\vdots\\ \varphi_{M}(x_{0})&\varphi_{M}(x_{1})&\cdots&\varphi_{M}(x_{M})\end{array}\right)}_{\textbf{M}}\left(\begin{array}[]{c}\omega^{(\alpha)}_{i0}\\ \omega^{(\alpha)}_{i1}\\ \vdots\\ \omega^{(\alpha)}_{iM}\end{array}\right)=\underbrace{\left(\begin{array}[]{c}\boldsymbol{D}^{\alpha}\varphi_{0}(x_{i})\\ \boldsymbol{D}^{\alpha}\varphi_{1}(x_{i})\\ \vdots\\ \boldsymbol{D}^{\alpha}\varphi_{M}(x_{i})\end{array}\right)}_{\boldsymbol{D}^{\alpha}\boldsymbol{\varphi}(x_{i})},$$

$$\textbf{M}=\left(\begin{array}[]{cccc}\epsilon&\sqrt{(x_{1}-x_{0})^{2}+\epsilon^{2}}&\cdots&\sqrt{(x_{M}-x_{0})^{2}+\epsilon^{2}}\\ \sqrt{(x_{0}-x_{1})^{2}+\epsilon^{2}}&\epsilon&\cdots&\sqrt{(x_{M}-x_{1})^{2}+\epsilon^{2}}\\ \vdots&\vdots&\ddots&\vdots\\ \sqrt{(x_{0}-x_{M})^{2}+\epsilon^{2}}&\sqrt{(x_{1}-x_{M})^{2}+\epsilon^{2}}&\cdots&\epsilon\end{array}\right),$$

$$\frac{{\partial }^{\alpha }{\phi }_k(x_i)}{{\partial }_{+}{x}^{\alpha }}$$

$$\frac{{\partial }^{\alpha }{\phi }_k(x_i)}{{\partial }_{-}{x}^{\alpha }}$$

$${\int }_{-1}^1(1-\zeta {)}^{\lambda }(1+\zeta {)}^{\mu }f(\zeta )d\zeta ,\lambda ,\mu >-1,$$

$$\frac{\partial y(x_i,t)}{\partial t}-\kappa (x_i)\sum _{j=0}^M{a}_{ij}^{(\alpha )}y(x_j,t)-\upsilon (x_i)\sum _{j=0}^M{b}_{ij}^{(\alpha )}y(x_j,t)=f(x_i,t),i=0,1,\cdots ,M,$$

$$\bigg{(}\textbf{I}-\tau\frac{\boldsymbol{\kappa}\textbf{A}+\boldsymbol{\upsilon}\textbf{B}}{2}\bigg{)}\textbf{Y}^{n}=\bigg{(}\textbf{I}+\tau\frac{\boldsymbol{\kappa}\textbf{A}+\boldsymbol{\upsilon}\textbf{B}}{2}\bigg{)}\textbf{Y}^{n-1}+\tau\textbf{H}^{n-1/2},$$

$$\displaystyle\textbf{A}=\left(\begin{array}[]{cccc}a^{(\alpha)}_{11}&a^{(\alpha)}_{12}&\cdots&a^{(\alpha)}_{1,M-1}\\ a^{(\alpha)}_{21}&a^{(\alpha)}_{22}&\cdots&a^{(\alpha)}_{2,M-1}\\ \vdots&\vdots&\ddots&\vdots\\ a^{(\alpha)}_{M-1,1}&a^{(\alpha)}_{M-1,2}&\cdots&a^{(\alpha)}_{M-1,M-1}\end{array}\right),\textbf{B}=\left(\begin{array}[]{cccc}b^{(\alpha)}_{11}&b^{(\alpha)}_{12}&\cdots&b^{(\alpha)}_{1,M-1}\\ b^{(\alpha)}_{21}&b^{(\alpha)}_{22}&\cdots&b^{(\alpha)}_{2,M-1}\\ \vdots&\vdots&\ddots&\vdots\\ b^{(\alpha)}_{M-1,1}&b^{(\alpha)}_{M-1,2}&\cdots&b^{(\alpha)}_{M-1,M-1}\end{array}\right),$$

$$\displaystyle\textbf{H}^{n-1/2}=\left(\begin{array}[]{c}f^{n-1/2}_{1}\\ f^{n-1/2}_{2}\\ \vdots\\ f^{n-1/2}_{M-1}\end{array}\right)+\frac{g_{0}^{n}+g_{0}^{n-1}}{2}\left(\begin{array}[]{c}\omega^{(\alpha)}_{1}\\ \omega^{(\alpha)}_{2}\\ \vdots\\ \omega^{(\alpha)}_{M-1}\end{array}\right)+\frac{g_{M}^{n}+g_{M}^{n-1}}{2}\left(\begin{array}[]{c}\tilde{\omega}^{(\alpha)}_{1}\\ \tilde{\omega}^{(\alpha)}_{2}\\ \vdots\\ \tilde{\omega}^{(\alpha)}_{M-1}\end{array}\right),$$

Full text

††journal: Elsevier

1 Introduction

In this study, we are mainly interested in an efficient method for numerically solving a class of space-fractional models in the following form

[TABLE]

subjected to the initial and boundary conditions

[TABLE]

where $1<\alpha\leqslant 2$, $\varLambda=[a,b]$, $\kappa(x)$, $\upsilon(x)$ are non-negative but do not vanish altogether. $g_{1}(t)\neq 0$ only when $\kappa(x)\equiv 0$ and $g_{2}(t)\neq 0$ only when $\upsilon(x)\equiv 0$. In Eq. (1.1), the space-fractional derivatives are defined in Caputo sense, i.e.,

[TABLE]

with the Euler’s Gamma function $\Gamma(\cdot)$.

The space-fractional PDEs describe many physical phenomena such as anomalous transport, hereditary elasticity, and chaotic dynamics [14, 22, 34], while compared favorably to the integer PDEs. Since they are frequently sufficing in the absence of exact closed-form solutions, various numerical algorithms have been designed to solve them, typically including general Padé approximation [5], finite difference methods [13, 28, 29], meshless point interpolation method [12], finite element methods [7, 36], discontinuous Galerkin method [32], finite volume method [10], spline approximation method (SAM) [26], RBF Kansa method [16], polynomial and fractional spectral collocation methods [27, 35]. In [3, 6, 21, 23], a series of operational matrix methods are constructed via the approximate expansions using shifted Jacobi, Chebyshev, Legendre polynomials, and Haar wavelets functions, as elements, respectively. Some analytic techniques are referred to [15, 19, 20, 33] and references therein.

DQ method is understood as a direct numerical approach for PDEs that evaluates the derivatives via representative weighted linear combinations of function values on problem domain [1]. A group of test functions to calculate these weights can be chosen as Lagrange basis functions, RBFs, and orthogonal polynomials [2, 18, 25, 30]. DQ method is characterized by a few advantages such as high accuracy, low occupancy cost, truly mesh-free and the ease of programming.

Due to the non-locality of fractional derivatives, a great extra computational cost is usually incurred when a conventional algorithm is applied to a fractional PDE. In this work, we propose a new RBFs based DQ method (RBF-DQM) for Eqs. (1.1)-(1.3). Using three types of RBFs as test functions, the weights are successfully determined and with them, the equation under consideration degenerates to an ordinary differential system (ODS). A time-stepping RBF-DQM is derived by introducing a difference scheme in time. The presented technique inherits the features of classic DQ methods. More importantly, it is insensitive to dimensional change, so it serves as a good alternative for the high-dimensional or the other complex fractional models arising in actual applications.

The outline is as follows. In Section 2, we give a brief description of fractional derivatives. In Section 3, the weighted coefficients are calculated by commonly used RBFs, which are required to approximate the fractional derivatives. We propose a Crank-Nicolson RBF-DQM to discretize the model problem in Section 4 and test its codes on three illustrative examples in Section 5. A conclusion is drawn in the last section.

2 Fractional derivatives

At first, some basic definitions are introduced for preliminaries. Let $\alpha\in\mathbb{R}^{+}$; then the following formulas

[TABLE]

define the left and right $\alpha$-th Caputo derivatives, respectively, if $f(x)\in C^{m}(\varLambda)$, where $m=[\alpha]+1$ for $\alpha\notin\mathbb{N}$, $m=\alpha$ for $\alpha\in\mathbb{N}$, and $[\cdot]$ is the floor function.

The left and right Caputo derivatives have the properties

[TABLE]

and coincide with classic derivatives with exactness to a multiplier factor $(-1)^{s}$:

[TABLE]

where $\beta>m-1$ and $s\in\mathbb{N}$. We refer the readers to [11, 17] for more properties.

3 DQ formulations based on RBFs

In the sequel, DQ formulations for fractional derivatives based on RBFs are derived. Define a lattice on $[a,b]$ but not necessarily with equally spaced points, i.e., $a=x_{0}<x_{1}<\cdots<x_{M-1}<x_{M}=b$, $M\in\mathbb{Z}^{+}$. In general, we always approximate the exact solution of a PDE like Eqs. (1.1)-(1.3) in the form

[TABLE]

with a set of proper basis functions $\{\phi_{k}(x)\}_{k=0}^{M}$. However, if

[TABLE]

and on acting $\frac{\partial^{\alpha}}{\partial_{+}x^{\alpha}}$, $\frac{\partial^{\alpha}}{\partial_{-}x^{\alpha}}$ on both sides of Eq. (3.4), we realize that

[TABLE]

thanks to the linearity of the fractional derivatives, namely, Eqs. (3.7)-(3.8) are valid as along as Eqs. (3.5)-(3.6) are satisfied. The idea behind this is referred to as DQ [1]; $a_{ij}^{(\alpha)}$, $b_{ij}^{(\alpha)}$, $i,j=0,1,\ldots,M$, are called the weighted coefficients of fractional derivatives and will be calculated by means of typical RBFs.

3.1 Three typical RBFs

RBFs are the functions of the distance from their centers. They are popular as an effective tool to set up numerical algorithms for PDEs since the superiority of potential spectral accuracy. Here, commonly used RBFs are involved, i.e.,

Multiquadrics (MQ): $\varphi_{k}(x)=\sqrt{r_{k}^{2}+\epsilon^{2}}$

• 2.

Inverse Multiquadrics (IM): $\varphi_{k}(x)=\frac{1}{\sqrt{r_{k}^{2}+\epsilon^{2}}}$

• 3.

Gaussians (GA): $\varphi_{k}(x)=e^{-\epsilon r_{k}^{2}}$

where $r_{k}=|x-x_{k}|$, $k=0,1,\ldots,M$ and $\epsilon$ is the shape parameter. It is worthy to note that the value $\epsilon$ should be well prescribed in computation because it has a significant impact on the approximation power of a RBFs based method.

3.2 Weighted coefficients for fractional derivatives

In order to obtain the weighted coefficients of the left and right fractional derivatives, we substitute the RBFs into Eqs. (3.5)-(3.6) to get

[TABLE]

Rewriting Eqs. (3.9)-(3.10) in a matrix-vector form for each grid point $x_{i}$ yields

[TABLE]

where $\omega^{(\alpha)}_{ij}=a^{(\alpha)}_{ij}$ if $\boldsymbol{D}^{\alpha}=\frac{\partial^{\alpha}}{\partial_{+}x^{\alpha}}$ whereas $\omega_{ij}^{(\alpha)}=b^{(\alpha)}_{ij}$ if $\boldsymbol{D}^{\alpha}=\frac{\partial^{\alpha}}{\partial_{-}x^{\alpha}}$, $i,j=0,1,\ldots,M$. M is the interpolation matrix only related to the nodal distribution, being nonsingular for MQ and fully positive definite for IM, GA [4]. One has

[TABLE]

in particular, when MQ RBFs are adopted.

There are clearly no explicit expressions for $\frac{\partial^{\alpha}\varphi_{k}(x)}{\partial_{+}x^{\alpha}}$, $\frac{\partial^{\alpha}\varphi_{k}(x)}{\partial_{-}x^{\alpha}}$; fortunately, they can be approximated by numerical quadrature rules. Taking the transforms $\xi=x-\frac{(x-a)(1+\zeta)}{2}$, $\xi=x+\frac{(b-x)(1+\zeta)}{2}$ of variables, respectively, reaches to

[TABLE]

A close examination reveals that both of the two formulas are the special cases of the following weakly singular integral, i.e.,

[TABLE]

with $\lambda=0$, $\mu=1-\alpha$ that can be handled by Gauss-Jacobi quadrature rules. $a_{ij}^{(\alpha)}$, $b_{ij}^{(\alpha)}$ are then determined by solving Eqs. (3.23) for each $x_{i}$ and the fractional derivatives are removed from a fractional PDE by using Eqs. (3.7)-(3.8) as replacements, thus we obtain the required solution by solving a ODS instead.

4 A time-stepping RBF-DQM for fractional PDEs

In this section, a RBF-DQM of fully discretization is derived for the space-fractional PDEs via the above direct approximations for fractional derivatives. Define a lattice on $[0,T]$ with equally spaced points $t_{n}=n\tau$, $\tau=T/N$, $N\in\mathbb{Z}^{+}$. On substituting the weighted sums (3.7)-(3.8) into Eq. (1.1), we have

[TABLE]

actually being a first-order ODS. Also, denote $t_{n-1/2}=t_{n}-\frac{\tau}{2}$, $y^{n}_{i}=y(x_{i},t_{n})$, $f^{n-1/2}_{i}=f(x_{i},t_{n-1/2})$ for brevity. Imposing the associated constraints (1.2)-(1.3) and rewriting the ODS in matrix-vector form, a time-stepping RBF-DQM is then derived by introducing a Crank-Nicolson scheme in time, given as

[TABLE]

where I is an identity matrix, $\textbf{Y}^{n}=[y^{n}_{1},y^{n}_{2},\cdots,y^{n}_{M-1}]^{T}$, $\boldsymbol{\kappa}=\textrm{diag}(\kappa_{1},\kappa_{2},\cdots,\kappa_{M-1})$, $\boldsymbol{\upsilon}=\textrm{diag}(\upsilon_{1},\upsilon_{2},\cdots,\upsilon_{M-1})$, and A, B, $\textbf{H}^{n-1/2}$ are as follows

[TABLE]

with $\omega^{(\alpha)}_{i}=\kappa_{i}a^{(\alpha)}_{i0}+\upsilon_{i}b^{(\alpha)}_{i0}$, $\tilde{\omega}^{(\alpha)}_{i}=\kappa_{i}a^{(\alpha)}_{iM}+\upsilon_{i}b^{(\alpha)}_{iM}$. We perform the procedures on the nodal distribution $x_{i}=0.5(1-\cos\frac{i\pi}{M})\ell+a$, $\ell=b-a$, $i=0,1,\cdots,M$. An detailed implementation of RBF-DQM is summarized in the following flowchart

Input $\alpha$, $\epsilon$, $M$, $N$, and allocate $\{t_{n}\}_{n=0}^{N}$, $\{x_{i}\}_{i=0}^{M}$.

• 2.

Form M, compute $\boldsymbol{D}^{\alpha}\boldsymbol{\varphi}(x_{i})$ by Gauss-Jacobi quadrature rules, and solve Eqs. (3.23) for each $x_{i}$ so that the weighted coefficients are found.

• 3.

Do a loop from $n=1$ to $N$ to solve Eqs. (4.24) for each $t_{n}$ by forming A, B first and output the desirable approximation $\textbf{Y}^{n}$ at each time step.

5 Illustrative examples

In this part, the proposed methods, termed as MQ-DQM, IM-DQM, and GA-DQM hereinafter, are studied on three numerical examples. The shape parameter should be adjusted with the grid number $M$, so we select $\epsilon=1.25\ell/(M+1)^{0.5}$ for MQ, $\epsilon=2/(M+1)^{0.5}$ for IM, and $\epsilon=1.05(M+1)$ for GA, tentatively, by references to [8, 9, 31]. The numerical errors are all defined in $l_{\infty}$-norm and the fractional derivatives are computed with 16 quadrature points and weights, whose values corresponding to $\alpha=1.5$ are given as reference in A.

Example 5.1. Approximate ${{}_{0}}D_{x}^{1.1}\sin(x)$ by Eq. (3.7) on $[0,1]$ with the above $\epsilon$, whose explicit expression is $-x^{1.9}{{}_{1}}F_{2}(1;1.45,1.95;-0.25x^{2})/\Gamma(2.9)$, where ${{}_{1}}F_{2}$ is the hypergeometric function. The numerical results are tabulated in Table 1. As observed, the approximation improves as $M$ increases, which implies that the DQ approximations for the fractional derivatives are valid. Moreover, under the given $\epsilon$, GA-DQM seems to be more efficient than MQ-DQM and IM-DQM.

Example 5.2. Let $\kappa(x)=\upsilon(x)=1$, and $y(x,t)=t^{3}x^{2}(1-x)^{2}$; we solve Eqs. (1.1)-(1.3) on $[0,1]$ with zero initial-boundary conditions. Table 2 displays the numerical results at $t=0.5$ when $\tau=2.0\times 10^{-4}$ and $\alpha=1.8$. From the table, we find that sufficiently small errors can be achieved by MQ-DQM, IM-DQM, and GA-DQM even if a few spatial grid points are utilized and all the methods are obviously convergent by taking their own $\epsilon$, respectively.

Example 5.3. Let $\kappa(x)=\frac{x^{\alpha}\Gamma(5-\alpha)}{24}$, $\upsilon(x)=0$, and $y(x,t)=e^{-t}x^{4}$; we solve Eqs. (1.1)-(1.3) on $[0,1]$ with $\psi(x)=x^{4}$, $g_{1}(t)=0$, and $g_{2}(t)=e^{-t}$. The numerical errors of SAM [26] and our methods at $t=1$ are reported in Table 3, when $\tau=1/M$ and $\alpha=1.5$. It is observed from the table that MQ-DQM, IM-DQM, and GA-DQM outperform SAM in term of computational accuracy.

Remarks 5.1**.**

When $M$ is fixed, the value $\epsilon$ is crucial to the accuracy of a RBFs based method, so is RBF-DQM. A general trade-off principle demanding attention is that one can adjust $\epsilon$ to decrease the approximate errors, but need to pay for this by increasing the condition number of the interpolation matrix which may cause an algorithm to be instable [24], so a good $\epsilon$ that balances both the accuracy and stability is anticipated in practice. Nevertheless, how to select an optimal value is technical and is being an issue deserving to investigate.

6 Conclusion

In this research, an efficient DQ method is proposed for the space-fractional PDEs of Caputo type based on commonly used RBFs as test functions, which enjoys some properties such as high accuracy, flexibility, truly meshless, and the simplicity in implementation. Its codes are tested on three benchmark examples and the outcomes manifest that it is capable of dealing with these problems if the free parameters $\epsilon$ are well prepared. Due to its insensitivity to dimensional change, our method has potential advantages over traditional methods in finding the approximate solutions to the high-dimensional fractional equations.

Acknowledgement: This research was supported by National Natural Science Foundations of China (No.11471262 and 11501450).

Appendix A Quadrature points and weights when $\alpha=1.5$