Second order accuracy finite difference methods for space-fractional partial differential equations

Accepted Manuscript Second order accuracy finite difference methods for space-fractional partial differential equations Yuki Takeuchi, Yoshihide Yoshimoto, Reiji Suda PII: DOI: Reference:

S0377-0427(17)30024-9 http://dx.doi.org/10.1016/j.cam.2017.01.013 CAM 10978

To appear in:

Journal of Computational and Applied Mathematics

Received date: 12 March 2016 Revised date: 27 November 2016 Please cite this article as: Y. Takeuchi, Y. Yoshimoto, R. Suda, Second order accuracy finite difference methods for space-fractional partial differential equations, Journal of Computational and Applied Mathematics (2017), http://dx.doi.org/10.1016/j.cam.2017.01.013 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Manuscript Click here to view linked References

Second order accuracy finite difference methods for space-fractional partial differential equations✩ Yuki Takeuchi, Yoshihide Yoshimoto, Reiji Suda The University of Tokyo, Graduate School of Information Science and Technology, Department of Computer Science

Abstract Space-fractional partial differential equations are used for simulations of, for example, diffusion of radioactive materials, and financial and other models, which are characterized by heavy-tailed distributions. A number of first order accuracy finite difference methods have been proposed. In the present paper, we introduce second order accuracy finite difference methods with Dirichlet boundary conditions. These methods have a parameter in these schemes, and the parameter stabilizes the schemes. This means that there exist various schemes with second order accuracy, but the stability of each scheme is different. In the present paper, we introduce the most stable scheme for any fractional calculus order by choosing the optimal parameter. In addition, we describe a phenomenon whereby the expected accuracy cannot be obtained if the analytical solution can be expanded to a series having less than second order around boundaries. This also happens in both existing methods and the proposed methods. In the present paper, we develop the stability conditions for the proposed schemes, and numerical examples of second order accuracy and accuracy decay are shown. Keywords: fractional calculus, space-fractional partial differential equations, finite difference method, difference formula, boundary problem, Dirichlet boundary condition, partial differential equations

We recieved generous support from Prof. Yuko Hatano, and this work was supported by JSPS KAKENHI Grant Number 26-6678 ✩

Preprint submitted to Computational And Applied Mathematics

November 27, 2016

1. Introduction Fractional calculus is an extended version of integer order calculus, and space-fractional differential equations have been used in many models, such as models of fractional dynamics for anomalous diffusion [Metzler and Klafter (2000, 2004)], fractional Fokker-Plank equations describing anomalous diffusion [Metzler et al. (1999); Barkai et al. (2000)], L´evy stable processes [West et al. (1997); Chechkin et al. (2006)], and porous media [de Pablo et al. (2011)]. Therefore, research on computational methods of fractional differential equations is also required. In the present paper, we treat the following one-dimensional parabolic two-sided space-fractional partial differential equation ∂u(x, t) C q = [LDxq u(x, t) + xDR u(x, t)] + f (x, t) (1) ∂t 2 for 0 < q < 2, L ≤ x ≤ R and t ≥ 0. The constant q is the fractional order. The functions f is known function, and the function u is an unknown function. The diffusion constant C is positive for 1 ≤ q < 2 and negative for 0 < q < 1. The initial condition is u(x, 0) = u0 (x) and boundary conditions are u(L, t) = a and u(R, t) = b as Dirichlet boundary conditions. In the related works, zero-Dirichlet boundary conditions are generally imposed when we consider about parabolic fractional partial differential equations. However, we treat inhomogeneous boundary conditions in addition to them. This assumption extends the range of applications and will help to consider other boundary conditions, for example, Neumann or Robin boundary conditions. This kind of two-sided space-fractional differential equations can represent L´evy stable process whose particles move with heavy-tailed distributions, and is developed by using Fourier transform[Chaves (1998); Meerschaert et al. (1999)]. The fractional derivatives in Eq. (1) are Riemann-Liouville definition, and are defined by [ ]⌈q⌉ ∫ x d 1 ϕ(ξ) q dξ, LDx ϕ(x) = dx Γ(⌈q⌉ − q) L (x − ξ)1+q−⌈q⌉ [ ]⌈q⌉ ∫ R d ϕ(ξ) 1 q dξ. xDR ϕ(x) = dx Γ(⌈q⌉ − q) x (ξ − x)1+q−⌈q⌉ The constant C is negative for 0 < q < 1 and positive for 1 < q < 2. The function u is an unknown function, and the function f is a known function. 2

For space fractional differential equations, a number of numerical solving methods have been proposed. One method is to use the Grunwald-Letnikov definition, which is given by q LDx ϕ(x)

N −1 h−q ∑ Γ(j − q) = lim ϕ(x − jh) N →∞ Γ(−q) Γ(j + 1) j=0

where h = (x − L)/N . Meerschaert and Tadjeran proposed finite difference schemes based on the Grunwald-Letnikov definition[Meerschaert and Tadjeran (2004, 2006); Meerschaert et al. (2006)]. Let hx and ht be the step size of space and time respectively, and let Ujm be an approximate solution as Ujm ≃ u(L + jhx , mht ). In addition, let fjm be fjm = f (L + jhx , mht ). Their explicit and implicit finite difference schemes are first order accuracy as O(ht ) + O(hx ). Their explicit finite difference scheme is given by

=

Ujm+1 − Ujm ht [ j+1 Γ(i − q) C ∑

hqx

i=0

Γ(−q)Γ(i + 1)

m Uj−i+1 +

Nx −j+1

∑ i=0

+fjm

Γ(i − q) m Uj+i−1 Γ(−q)Γ(i + 1)

] (2)

where Nx is defined by hx = (R − L)/Nx . Their implicit finite difference scheme is given by

=

Ujm − Ujm−1 ht [ j+1 Γ(i − q) C ∑

hqx

i=0

Γ(−q)Γ(i + 1)

m Uj−i+1 +

Nx −j+1

∑ i=0

+fjm .

Γ(i − q) Um Γ(−q)Γ(i + 1) j+i−1

] (3)

The stability conditions for the above two schemes have been analyzed using Gerschgorin’s theorem. The implicit finite difference scheme (3) is unconditionally stable, and a stability condition for the explicit finite difference scheme (2) is given by ht ≤

hqx 1 . Cq 3

(4)

This stability condition is developed by applying the matrix method which analyzes the eigenvalues of matrices in finite difference schemes. For the explicit finite difference scheme (2), Sousa obtained the following stability condition by applying von Neumann’s stability analysis[Sousa (2009)] ht ≤

hqx 1−q 2 . C

(5)

In order to obtain high-accuracy results with finite difference methods, Meerschaert and Tadjeran applied an extrapolation method to their finite difference schemes[Tadjeran and Meerschaert (2007)]. In contrast, Chen et al. obtaind second and higher accuracy finite difference methods which is derived from Fourier transform[Chen and Deng (2014a,b); Chen et al. (2013); Tian et al. (2015)]. In addition, the authors also obtained second order accuracy finite difference method for space fractional diffusion equations[Takeuchi and Suda (2013)]. That old method is based on the second order accuracy difference formula for fractional derivatives, and this means that high-accuracy results can be obtained without any extrapolation method. However, the stability is not analyzed and that old method are not stable for some fractional order in fact. Also, that old method allows only zero Dirichlet boundary conditions. In contrast to that old method, we propose new finite difference methods using a second order accuracy scheme with stability analysis in the present paper. The proposed second order accuracy finite difference methods are stable for any fractional order and allow non-zero Dirichlet boundary conditions. In the present paper, we first introduce the proposed second order accuracy difference formula. Next, we introduce the proposed finite difference schemes and analyze their stabilities using Gerschgorin’s theorem. Finally, we conduct numerical experiments for comparison with existing methods and show that the proposed finite difference schemes are more accurate and more stable than existing first order accuracy finite difference methods. 2. Finite difference schemes 2.1. Second order accuracy difference formula The proposed second order accuracy finite difference methods are based on a second order accuracy difference formula, which shall be referred to herein as Method A. Method A, which was proposed in a previous paper 4

[Takeuchi and Suda (2012)], will be described for the convenience of the readers. Let L Dxq be the Riemann-Liouville fractional derivative operator. Let N be a grid number, and let h be the step size. Then, it holds that h = (x − L)/N , and for the fractional order q, Method A is given by the following theorem. Theorem 1. If a function f is infinitely continuously differentiable and analytic for [L, x], it holds that q LDx f (x) −1 −q N ∑

[ ] Γ(j − q) qh ′ f (x − jh) + f (x − jh) Γ(j + 1) 2 j=0 ( ) h−q 1 + q 1 −1−q + f (L)N +O Γ(−q) 2 N2

h = Γ(−q)

(6)

for N → ∞[Takeuchi and Suda (2012)]. The proof of this theorem is presented in Appendix A. Method A has second-order accuracy, and some other formulae having accuracy of the same order can be derived from Method A with appropriate approximations. For example, the following expression is a second-order accuracy formula for fractional derivative L Dxq f (x), [ ] N −1 h−q ∑ Γ(j − q) q + 2 q f (x − jh) − f (x − (j + 1)h) Γ(−q) j=0 Γ(j + 1) 2 2 +

h−q 1 + q f (L)N −1−q . Γ(−q) 2

(7)

This expression is also introduced by Tian et al.[Tian et al. (2015)] in the special case of f (L) = 0. 2.2. Second order accuracy finite difference methods for fractional differential equations In the above section, we described Method A, which is a second order accuracy difference formula. In this section, we introduce second order accuracy finite difference methods that use Method A. In order to introduce the finite difference scheme, we set the following grid, as shown in Figure 1. 5

Figure 1: Grids

Both the proposed and existing finite difference methods calculate approximate solutions at each grid point. In the figure, the red line represents the initial condition, and the green lines denote the boundary conditions. The constants Nx and Nt are the grid numbers for space and time respectively. The constant hx is the step size for space and is defined by hx = (R − L)/Nx . The constant ht is the step size for time. First, we introduce the proposed explicit and implicit schemes, followed by an analysis of these schemes. The proposed explicit finite difference scheme, i.e., Scheme A is given by

=

Ujm+1 − Ujm h [ t { j−1 C h−q ∑ Γ(i − q) 4s + q x

2

Γ(−q)

Γ(i + 1)

i=0 Nx −j−1

4

m Uj+1−i

+ (1 −

m 2s)Uj−i

4s − q m + Uj−1−i 4

}

} ∑ Γ(i − q) { 4s + q 4s − q m h−q x m m + Uj−1+i + (1 − 2s)Uj+i + Uj+1+i Γ(−q) i=0 Γ(i + 1) 4 4 +

h−q Γ(Nx − j − q) 4s + q m h−q Γ(j − q) 4s + q m x U1 + x UNx −1 Γ(−q) Γ(j + 1) 4 Γ(−q) Γ(Nx − j + 1) 4

6

} j −q Γ(j − q) Γ(j − q) 4s + q − − U0m Γ(1 − q) Γ(1 − q)Γ(j) Γ(−q)Γ(j + 1) 4 { } ] (Nx − j)−q Γ(Nx − j − q) Γ(Nx − j − q) 4s + q −q m +hx − − UNx Γ(1 − q) Γ(1 − q)Γ(Nx − j) Γ(−q)Γ(Nx − j + 1) 4 +fjm (8) +h−q x

{

for j = 1, 2, . . . , Nx −1. A significant characteristic is that Scheme A includes the parameter s. Scheme A is stable for any q by taking proper values of s. The value that the parameter s should take will be discussed in the ⃗ m+1 = next subsection. Scheme A can be represented using matrices as U m m m m m m T ⃗m m m ⃗ ⃗ ⃗ (E + A)U + f where U = (U0 , U1 , . . . , UNx ) , f = (f0 , f1 , . . . , fNmx )T and matrix E is an identity matrix. The matrix A is defined using elements ai,j as follows:   0 0 0 ... 0 0  a1,0 a1,1 a1,2 ... a1,Nx −1 a1,Nx     a2,0  a a . . . a a 2,1 2,2 2,N −1 2,N x x   A= . .. .. .. .. .. ..   . . . . . .    aNx −1,0 aNx −1,1 aNx −1,2 . . . aNx −1,Nx −1 aNx −1,Nx  0 0 0 ... 0 0 The entries ai,j for i = 1, 2, . . . , Nx − 1, j = 1, 2, . . . , Nx − 1 are defined by [ ] 4s + q Cr 2g0 (1 − 2s) + 2g1 an,n = 2 4 [ ] Cr 4s + q 4s − q 4s + q an,n+1 = an+1,n = g0 + g0 + g1 (1 − 2s) + g2 2 4 4 4 [ ] Cr 4s − q 4s + q an,n+2 = an+2,n = g1 + g2 (1 − 2s) + g3 2 4 4 [ ] Cr 4s − q 4s + q an,n+3 = an+3,n = g2 + g3 (1 − 2s) + g4 2 4 4 .. . [ ] Cr 4s + q 4s − q an,n+k = an+k,n = + gk (1 − 2s) + gk+1 gk−1 2 4 4 .. .

for k = 2, 3, . . . and n = 1, 2, . . . , Nx − 1 where r = ht /hqx and gk = Γ(k − q)/(Γ(−q)Γ(k + 1)). 7

Next, we introduce the proposed implicit scheme Scheme B, which is given by

=

Ujm − Ujm−1 h [ t { j−1 C h−q ∑ Γ(i − q) 4s + q x

2

Γ(−q)

Γ(i + 1)

i=0 Nx −j−1

4

m Uj+1−i

+ (1 −

m 2s)Uj−i

4s − q m + Uj−1−i 4

}

} ∑ Γ(i − q) { 4s + q h−q 4s − q m x m m + Uj−1+i + (1 − 2s)Uj+i + Uj+1+i Γ(−q) i=0 Γ(i + 1) 4 4

h−q Γ(j − q) 4s + q m h−q Γ(Nx − j − q) 4s + q m x U1 + x UNx −1 Γ(−q) Γ(j + 1) 4 Γ(−q) Γ(Nx − j + 1) 4 { } j −q Γ(j − q) Γ(j − q) 4s + q −q +hx − − U0m Γ(1 − q) Γ(1 − q)Γ(j) Γ(−q)Γ(j + 1) 4 { } ] (Nx − j)−q Γ(Nx − j − q) Γ(Nx − j − q) 4s + q −q m +hx − − UNx Γ(1 − q) Γ(1 − q)Γ(Nx − j) Γ(−q)Γ(Nx − j + 1) 4 +fjm . (9) +

⃗m = U ⃗ m−1 + This scheme also can be represented using matrices as (E − A)U f⃗m with the previously defined matrix A. These two proposed schemes have three features comparing to the existing schemes (2) and (3) by Meerschaert et al. One is the existence of parameter s as is mentioned above. In Schemes A and B, changing the value of s means varying the form of the Schemes or varying the weight at each grid, and the stability is also varied. In addition, the proper value of s, which makes schemes stable, must be chosen depending on the calculus order q. In the next subsection, the range that the parameter s should be included will be shown. The second feature is the symmetric Toeplitz matrix in the schemes. The two proposed Schemes A and B include the symmetric Toeplitz submatrix ai,j , 1 ≤ i, j ≤ Nx − 1. Some schemes of finite difference methods without symmetric Toeplitz submatrices were less stable than the proposed schemes for parabolic partial differential equations with constant coefficients. The schemes proposed by Meerschaert et al. also include a symmetric Toeplitz matrix. In addition, fractional differential operator is represented by Toeplitz matrices in [Tian et al. (2015); Chen and Deng (2014a)]. There are two reasons that the schemes should include the symmetric Toeplitz submatrix. One is Eq. (1) represents a diffusion process. In diffusion process without 8

any drifts, transition probabilities between two positions are the same each other. Therefore, the submatrix should be the symmetric matrix. The second reason is translational symmetry. Fractional derivatives are the extension of integer order derivatives, and the submatrix in schemes should be Toeplitz matrix for translational symmetry. The third feature is the treatment of the boundary conditions. In the abovementioned studies, the schemes proposed by Meerschaert et al. are applicable only to zero Dirichlet boundary conditions as u(L, t) = u(R, t) = 0. In contrast, the proposed schemes are applicable to both zero and non-zero Dirichlet boundary conditions. In addition, the proposed schemes exactly compute a constant function, in the sense explained below. In the present paper, this feature of the proposed schemes is referred to as error canceling. Let us consider the solution function to be a constant function, i.e., u(x) = c. Then, the error at x = jh is expressed by { } j−1 ∑ Γ(i − q) 4s + q h−q 4s − q q x c + (1 − 2s)c + c 0Dx c − Γ(−q) i=0 Γ(i + 1) 4 4 h−q Γ(j − q) 4s + q x c Γ(−q) Γ(j + 1) 4 { } j −q 4s + q Γ(j − q) Γ(j − q) −q −hx − − c Γ(1 − q) Γ(1 − q)Γ(j) Γ(−q)Γ(j + 1) 4 { } j−1 c(jh )−q −q ∑ Γ(i − q) ch−q ch Γ(j − q) x x x − − j −q − = Γ(1 − q) Γ(−q) Γ(i + 1) Γ(1 − q) Γ(j) i=0 −

= 0.

This means that the proposed schemes can be used to calculate fractional derivative of any constant function with no error. The proposed schemes are more complicated than formula (7), but provide the error canceling feature. The proposed schemes have not only better accuracy but also a wider range of application than the schemes proposed by Meerschaert et al. 2.3. Stability analysis of the proposed schemes In this subsection, we present a stability analysis of the proposed schemes and develop the stability conditions. First, we consider the proposed explicit Scheme A for 1 < q < 2. We analyze the stability by bounding eigenvalues of matrices that represents the proposed schemes. The matrix representation 9

⃗ m+1 = (E + A)U ⃗ m + f⃗m . If an arbitrary eigenvalue λ of of Scheme A is U matrix (E + A) satisfies |λ| ≤ 1, this scheme is stable. Similar to existing methods, we use Gerschgorin’s theorem to estimate eigenvalues. According to Gerschgorin’s theorem, the ∑ eigenvalues derived from the first row and the last row are 1 as |λ − 1| ≤ j̸=0 |a0,j | = 0, and the eigenvalues derived from the i-th column satisfy the following inequality for i = 1, 2, . . . , Nx − 1: ∑ |λ − (1 + ai,i )| ≤ |aj,i | . j̸=i

Note that we consider not rows, but columns, except for the first row and the last row. The above condition can be written as ∑ ∑ − |aj,i | + (1 + ai,i ) ≤ λ ≤ |aj,i | + (1 + ai,i ). j̸=i

j̸=i

Therefore, Scheme A is stable if elements ai,j satisfy the following two conditions ∑ |aj,i | + (1 + ai,i ) ≥ −1, − j̸=i

∑ j̸=i

|aj,i | + (1 + ai,i ) ≤ 1.

Here, we simplify the expressions by fixing the signs of ai,j , since we assume that the diagonal elements are negative or equal to 0 and the non-diagonal elements are positive or equal to 0. Then, we have the following two lemmas. Lemma 1. All diagonal entries are negative or equal to 0 if parameter s satisfies s≥

2−q . 4

(10)

Proof The condition whereby the diagonal entries are negative or equal to 0 is given by [ ] Cr 4s + q 2g0 (1 − 2s) + 2g1 ≤0 2 4 2−q . ⇒ s≥ 4 10

□ Lemma 2. All non-diagonal entries are positive or equal to 0 if parameter s satisfies the following three expressions: s≥

−4g1 − qg2 , 8g0 − 8g1 + 4g2

(11)

s≤

qg1 − 4g2 − qg3 , 4g1 − 8g2 + 4g3

(12)

s≥

qg2 − 4g3 − qg4 . 4g2 − 8g3 + 4g4

(13)

Proof The condition whereby the entry ai,i+1 is positive or equal to 0 is given by 2g0 s + g1 (1 − 2s) + g2 ⇒ s≥

−4g1 − qg2 . 8g0 − 8g1 + 4g2

4s + q ≥0 4

The condition whereby the entry ai,i+2 is positive or equal to 0 is given by 4s − q 4s + q + g2 (1 − 2s) + g3 ≥0 4 4 qg1 − 4g2 − qg3 . ⇒ s≤ 4g1 − 8g2 + 4g3 g1

The condition whereby the entry ai,i+k , k ≥ 3 is positive or equal to 0 is given by 4s − q 4s + q + gk (1 − 2s) + gk+1 ≥0 4 4 qgk−1 − 4gk − qgk+1 . ⇒ s≥ 4gk−1 − 8gk + 4gk+1 gk−1

Then, we have the expressions (11) and (12). Next, let fk (q) be fk (q) =

qgk−1 − 4gk − qgk+1 , 4gk−1 − 8gk + 4gk+1 11

and we prove that fk (q) > fk+1 (q) for k ≥ 3. Coefficients gk are related to the following coefficients as gk−1 = kgk /(k − 1 − q) and gk+1 = (k − q)gk /(k + 1). Then, the following holds: fk (q) =

qk g k−1−q k 4k g k−1−q k 2

− 4gk − − 8gk +

q(k−q) gk k+1 4(k−q) gk k+1

−4k + q(6 + 2q)k + (4 − q 2 )(1 + q) = . (4q + 8)(1 + q) Then, we have fk (q) − fk+1 (q) 8k + 4 − q(6 + 2q) >0 = (4q + 8)(1 + q) for 1 < q < 2 and k ≥ 3. Therefore, we have the third condition (13).

□

Under these two assumptions, the absolute value operation can be eliminated, and the two stability conditions are transformed to { ∑ x −1 − N aj,i + 2ai,i ≥ −2, ∑Nx j=1 −1 j=1 aj,i ≤ 0.

Next, let us consider a condition in which elements ai,j satisfy the above stability conditions. From lemmas 1 and 2, we have the following theorem as stability conditions of the proposed explicit Scheme A. Theorem 2. Scheme A for 1 < q < 2 is stable if the following inequalities are satisfied. −4g1 − qg2 , 8g0 − 8g1 + 4g2 qg1 − 4g2 − qg3 , s ≤ 4g1 − 8g2 + 4g3 qg2 − 4g3 − qg4 s ≥ , 4g2 − 8g3 + 4g4 4hqx ht ≤ . C((8 + 4q)s + q 2 − 4) s ≥

12

Proof First, we develop the condition for the elements ai,j to satisfy the following expression N x −1 ∑ j=1

aj,i ≤ 0.

For i = 1, it holds that N x −1 ∑ j=1

aj,1 ≤ 0

[ ] Cr 4s + q g1 + (1 − 2s)g0 ≤ 0 ⇒ 2 4 qg1 + 4g0 2−q ⇒ s≥− = 4(g1 − 2g0 ) 4

in the limit of Nx → ∞. For i = 2, Nx → ∞, it holds that N x −1 ∑ j=1

aj,k ≤ 0

[ ] Cr 4s + q Γ(3 − q) (1 − 2s) Γ(2 − q) 4s − q Γ(1 − q) ⇒ − − − ≤0 2 4q Γ(−q)Γ(3) q Γ(−q)Γ(2) 4q Γ(−q)Γ(1) qh1 − 4h2 − qh3 ⇒ s≤ 4h1 − 8h2 + 4h3

where hk is hk = Γ(k − q)/(Γ(−q)Γ(k)). For i = k, k = 3, 4, . . .,Nx → ∞, it holds that N x −1 ∑ j=1

aj,k ≤ 0

[ ] Cr 4s + q Γ(k + 1 − q) (1 − 2s) Γ(k − q) 4s − q Γ(k − 1 − q) ⇒ − − − ≤0 2 4q Γ(−q)Γ(k + 1) q Γ(−q)Γ(k) 4q Γ(−q)Γ(k − 1) qhk−1 − 4hk − qhk+1 ⇒ s≥ . 4hk−1 − 8hk + 4hk+1

In a similar mannar to the previous proof, we have

qhk−1 − 4hk − qhk+1 qhk − 4hk+1 − qhk+2 > 4hk−1 − 8hk + 4hk+1 4hk − 8hk+1 + 4hk+2 13

for k = 3, 4, . . .. Comparing the above expressions to expressions (11), (12) and (13), we have the following: −4g1 − qg2 2−q ≥ , 8g0 − 8g1 + 4g2 4 qh1 − 4h2 − qh3 qg1 − 4g2 − qg3 ≤ , s ≤ 4g1 − 8g2 + 4g3 4h1 − 8h2 + 4h3 qh2 − 4h3 − qh4 qg2 − 4g3 − qg4 ≥ s ≥ 4g2 − 8g3 + 4g4 4h2 − 8h3 + 4h4

s ≥

for 1 < q < 2. Next, we consider the following stability condition −

N x −1 ∑ j=1

aj,i + 2ai,i ≥ −2.

Since it holds that N x −1 ∑ j=1

aj,i ≤ 0,

we have the following inequality −

N x −1 ∑ j=1

aj,i + 2ai,i ≥ 2ai,i ≥ −2

[ ] Cr 4s + q ⇒ 2g0 (1 − 2s) + 2g1 ≥ −1. 2 4 Therefore, the stability condition is given by 4 (8 + 4q)s + q 2 − 4 4hqx ⇒ ht ≤ . C((8 + 4q)s + q 2 − 4) ⇒ Cr ≤

□

In a similar manner, the stability conditions for Scheme B (9) can be also developed, and we have the following theorem. 14

Theorem 3. Scheme B is stable for 1 < q < 2 if all of the following inequalities hold: −4g1 − qg2 , 8g0 − 8g1 + 4g2 qg1 − 4g2 − qg3 s ≤ , 4g1 − 8g2 + 4g3 qg2 − 4g3 − qg4 . s ≥ 4g2 − 8g3 + 4g4 s ≥

The proposed schemes have three additional stability conditions regarding parameter s. In particular, the proposed implicit Scheme B requires stability conditions for stable computation, whereas the implicit scheme (3) proposed by Meerschert et al. is unconditionally stable. The stability analysis revealed the range that parameter s should satisfy, and we have the following theorem concerning the existence of parameter s. Lemma 3. There exists a parameter s for arbitrary q with 1 < q < 2, and Scheme B can take the largest step size ht for { } −4g1 − qg2 qg2 − 4g3 − qg4 s = max , . q 8g0 − 8g1 + 4g2 4g2 − 8g3 + 4g4 In the proposed explicit method, we can take the largest step size ht by adopting the most appropriate parameter s. On the other hand, there is no constraint on ht in the proposed implicit method, and there is no most appropriate parameter s. Next, we consider the stability conditions for the proposed methods for 0 < q < 1. In a similar manner for 1 < q < 2, we have the following two theorems. Theorem 4. Scheme A is stable for 0 < q < 1 if the following are satisfied: −4g1 − qg2 , 8g0 − 8g1 + 4g2 qg1 − 4g2 − qg3 s ≥ , 4g1 − 8g2 + 4g3 hq 4 ht ≤ x . C (8 + 4q)s + q 2 − 4 s ≤

15

Theorem 5. Scheme B is stable for 0 < q < 1 if the following are satisfied: −4g1 − qg2 , 8g0 − 8g1 + 4g2 qg1 − 4g2 − qg3 s ≥ . 4g1 − 8g2 + 4g3 s ≤

Moreover, for 0 < q < 1, the proposed schemes require stability conditions for parameter s to be stable. In the same manner as for the case in which 1 < q < 2, we have the following theorem concerning the existence of parameter s. Theorem 6. There exists a parameter s for arbitrary q with 0 < q < 1, and the largest step size that the proposed explicit scheme (8) can take is ht for s=

−4g1 − qg2 . 8g0 − 8g1 + 4g2

2.4. Accuracy decay In the present paper, we state that the schemes proposed by Meerschaert et al. are first order accuracy in time and space and that the proposed Schemes A and B are second order accuracy in space. However, Method A assumes that functions are infinitely continuously differentiable and analytic in order to gain the expected accuracy. This means that the expected accuracy sometimes cannot be obtained for discontinuous or nondifferentiable functions. We herein refer to this phenomenon as accuracy decay. Accuracy decay was also reported by K. Diethelm et al.[Diethelm et al. (2004)], who revealed that the accuracy of a fractional Adams method decays for the undifferentiable function z(t) = tp , 0 < p < 1 in fractional ordinary differential equations. In the present paper, we demonstrate that accuracy decay occurs in difference formulae and finite difference schemes in fractional partial differential equations. Moreover, accuracy decay also occurs in integer order calculus. Let us consider the numerical differentiation of an undifferentiable function f (x) = x0.5 at x = 0 by forward difference with step size h. Then, we have d 1 f (x + h) − f (x) f (x) = x−0.5 ≃ = h−0.5 . dx 2 h This indicates that the numerical derivative by forward difference does not have first order accuracy. Similarly, in fractional calculus, Method A does 16

not have the expected accuracy for undifferentiable functions. However, this accuracy decay is not generally a problem in integer order calculus. If once we assume infinitely continuously differentiable functions, such as polynomials, the derivatives of these functions are also infinitely continuously differentiable, and undifferentiable functions do not occur. In contrast, fractional derivatives provide singularity to continuous and differentiable functions, as follows 0.5 0Dx x =

x0.5 . Γ(1.5)

This property allows us to treat not only fractional derivatives of continuous and differentiable functions but also those of undifferentiable functions. The degree of the accuracy decay associated with applying the difference formulae and the finite difference schemes will be demonstrated in the next section. Accuracy decay occurs in Method A and Schemes A and B. Let a function g(x) be of order O(hp ), 0 < p < 1 at x = 0. The experimental results show that the order of accuracy of Method A decays to O(h1+p ) from O(h2 ) around x = 0. In other words, this accuracy decay occurs when we numerically differentiate non-analytical functions. On the other hand, the condition under which accuracy decay occurs is different when we apply finite difference schemes for fractional differential equations. For the analytical solution function that can be expanded to a series having O(hp ) around boundaries, for example ( ) u(x, t) = exp(−t) (x − L)p + (x − L)1+p + . . . , 0 < p < 1, the order of accuracy of the schemes proposed by Meerschaert et al. decreases to O(hp ) from O(h). In addition, the order of accuracy of the proposed schemes decreases to O(hp ) from O(h2 ). Therefore, accuracy decay occurs even if the analytical solution is differentiable once or is a first order polynomial. To avoid accuracy decay, it may be effective to use adaptive finite difference methods[Oberman and Zwiers (2015); L¨otstedt et al. (2007); Jeong et al. (2014)]. Adaptive finite difference methods is the finite difference methods which utilize a graded mesh instead of a uniform mesh in space. In addition, non-uniform meshes have already applied to finite difference methods for time-fractional differential equations[Yuste and Quintana-Murillo (2012); Zhang et al. (2014)]. The author considers it is possible that the accuracy 17

will be improved by using narrower grids around boundaries. However, there is no perspective whether such a finite difference scheme is stable or not. 3. Experiments 3.1. Stability First, we conduct numerical experiments on the stability of the finite difference schemes. As mentioned above, the stability conditions for the proposed implicit scheme is given by Theorems 3 and 5. Figures 2 and 3 show the ranges in which the proposed schemes are stable for 1 < q < 2 and 0 < q < 1, respectively in x axis of q and y axis of s. In the figures,

Figure 2: Range in which parameter s should Figure 3: Range in which parameter s should satisfy 0 < q < 1 satisfy 1 < q < 2

the ’o’ symbol denotes the stability of Scheme B at that point, and the ’x’ symbol denotes instability. These results indicate that the schemes are stable if the parameter s is slightly out of the range shown in Theorems 3 and 5. The reason for this is that the stability conditions in the present paper are sufficient conditions. Since the matrix method uses Gerschgorin’s theorem to estimate eigenvalues, we obtain only sufficient conditions. Note that Scheme B is seemingly unconditionally stable in Figure 2, but there are unstable points, for example (q, s) = (1.5, 0.0), (1.8, 4.0), out of Figure 2. Next, we compare the stability condition for the proposed explicit Scheme A with that of the explicit scheme of Meerschaert et al. If the proposed scheme requires a much narrower step size ht to be stable for the same step size hx than existing schemes, the proposed scheme cannot be more accurate 18

than existing schemes. In Theorem 2, the stability condition for ht is given by ht ≤

4hqx . C((8 + 4q)s + q 2 − 4)

(14)

The stability conditions for the existing explicit scheme are obtained by two methods, the matrix method[Meerschaert and Tadjeran (2006)] and Von Neumann’s stability analysis[Sousa (2009)], and they are given by hqx 1 , Cq

(15)

hqx 1−q 2 . C

(16)

ht ≤ ht ≤

In comparing the stability conditions, the parameter s is chosen as s = −4g1 −qg2 , which is the upperbound of s used in the present paper. We can 8g0 −8g1 +4g2 choose a smaller s, but, as shown below, this s provides a looser bound for 4 ht . This means that step size ht is the narrowest and coefficient ((8+4q)s+q 2 −4) −4g1 −qg2 is the smallest for s = 8g0 −8g1 +4g2 . Figure 4 shows the stability conditions for C = 1, hx = 1/40, which indicate that we can take a larger step size ht for the same step size hx using Scheme A, as compared to the existing explicit scheme. Both accuracy and stability conditions improve by using the proposed schemes.

Figure 4: Stability conditions

19

3.2. Accuracy 3.2.1. Method A First, we conduct experiments on Method A that Schemes A and B are based on. Let us consider differentiating the following two functions. Let f1 (x) and f2 (x) be functions defined by f1 (x) = 1 − x + x2 − x3 + x4 , f2 (x) = 1 − x0.2 + x0.4 − x0.6 + x0.8 − x. We use Eq. (7) for numerical differentiation. Figures 5 and 6 show the errors

Figure 5: Errors of Method A for q = 0.3

Figure 6: Errors of Method A for q = 1.3

of Method A at x = 2 and L = 0 for q = 0.3 and q = 1.3, respectively. These results indicate that the order of accuracy decays from O(h2 ) to O(h1.2 ) for f2 , which is undifferentiable at x = 0. This means that if we differentiate a function that does not satisfy the assumption of Theorem 1 and have the order O(hp ) at x = L, the order of accuracy of Method A appears not to be O(h2 ), but rather to be O(h1+p ). 3.2.2. Space-fractional partial differential equations In this section, we conduct numerical experiments on the accuracy of the two methods, the finite difference methods of Meerschaert et al. and the proposed Schemes B. Let us assume two problems. In both problems, let constants (L, R) be (0, 1). All experiments are performed for 0 ≤ t ≤ 1 using implicit schemes. In Problem 1, let the analytical solution u(x, t) be u(x, t) = 4 exp(−t)xp (1 − x)p 20

where p is a constant to control the order of the analytical solution function around the boundaries. This means that the order of the function is O(xp ) at x = 0 and x = 1. Then, let the force term f (x, t) be f (x, t) = − exp(−t) {4xp (1 − x)p [ ] Γ(p + 1) −p, p + 1 p−q +2Cx ;x 2F1 p−q+1 Γ(p − q + 1) [ ]} Γ(p + 1) −p, p + 1 p−q +2C(1 − x) ;1 − x 2F1 p−q+1 Γ(p − q + 1) where the function 2F1 is the hypergeometric series. In Problem 1, the initial condition is given by u(x, 0) = xp (1 − x)p , and the boundary conditions are given by u(0, t) = u(1, t) = 0 as zero-Dirichlet boundary conditions. In Problem 2, let the analytical solution u(x, t) be defined by u(x, t) = exp(−t)x2 (1 − x)2 + 3(1 − x)2 − 2(1 − x)3 , and let the force term f (x, t) be f (x, t)

{ 2 (x2−q + (1 − x)2−q ) C = − exp(−t)x (1 − x) − exp(−t) 2 Γ(3 − q) } 3−q 3−q 4−q 12 (x + (1 − x) ) 24 (x + (1 − x)4−q ) − + Γ(4 − q) Γ(5 − q) { } −q 2−q x C 6 (x + (1 − x)2−q ) 12 (x3−q + (1 − x)3−q ) − − + . 2 Γ(1 − q) Γ(3 − q) Γ(4 − q) 2

2

In Problem 2, the initial condition is given by u1 (x, 0) = x2 (1 − x)2 + 3(1 − x)2 − 2(1 − x)3 , and the boundary conditions are given by u1 (0, t) = 1, u1 (1, t) = 0. Figures 7 and 8 show the absolute maximum errors defined by max |UjNt − u(jhx , Nt ht )| j

for p = 2.5, Nt = 1000 and C = 100000 in Problem 1. In the numerical experiments about accuracy, we put a large number to the constant C to emphasize the error derived from space derivatives. These results indicate that the numerical solutions of the scheme of Meerschaert et al. and 21

Figure 7: Errors of three numerical methods Figure 8: Errors of three numerical methods for q = 1.6, p = 2.5, s = 0.4 for q = 1.2, p = 2.5, s = 0.3

the proposed Scheme B are actually first order and second order accuracy respectively. Figures 9 and 10 show the errors for p = 0.5, Nt = 1000, and C = 100000 in Problem 1. Unlike the results for p = 2.5, we could not obtain the expected order of accuracy because of accuracy decay. The order of accuracy decays to O(h0.5 ) in both the existing and proposed finite difference schemes.

Figure 9: Errors of two numerical methods Figure 10: Errors of two numerical methods for q = 1.6, p = 0.5, s = 0.4 for q = 1.2, p = 0.5, s = 0.3

Next, we consider 0 < q < 1 in Problem 1. Figures 11 and 12 show the results for C = −100000, Nt = 1000, q = 0.8, 0.6, 0.4, 0.2, and s = 0.2, 0.1, 0.0, −0.1 for p = 2.5 and p = 0.5, respectively. Figure 11 indicates 22

that the accuracy of Scheme B is also higher than or equal to second order for 0 < q < 1. However, Figure 12 indicates that accuracy decay occurs for p = 0.5 in the same manner as for 1 < q < 2.

Figure 11: Errors of the proposed implicit Figure 12: Errors of the proposed implicit scheme for various q with p = 0.5 scheme for various q with p = 2.5

In Problem 2, the boundary conditions include a non-zero Dirichlet boundary condition at x = 0 as u(0, t) = 1, and the condition at x = 1 is a zero Dirichlet boundary condition because u(1, t) = 0. Figures 13 and 14 show that the numerical solutions of the existing finite difference scheme of Meerschaert et al. does not converge to the analytical solution. In contrast, we obtain almost second order accuracy results using Scheme B.

Figure 13: Errors of two schemes for Problem Figure 14: Errors of two schemes for Problem 2 with q = 1.2 2 with q = 1.6

23

4. Conclusion A number of models that use fractional calculus, which require computational methods for fractional differential equations, have been proposed previously. In the present paper, we propose new finite difference methods for space-fractional differential equations. The proposed finite difference methods provide three improvements. The first improvement is accuracy. The proposed methods are the second order accuracy due to the use of Method A, which is the second order accuracy. The second improvement is stability. The proposed finite difference schemes include parameter s, and the stability of the proposed finite difference schemes depends on parameter s. By using the appropriate parameter, the proposed explicit Scheme A is more stable than the existing explicit scheme of Meerschaert et al. As such, we can take a wider step size of time using Scheme A, as compared to the existing scheme. The third improvement is a wider scope of application. The proposed finite difference Schemes A and B can be applied for 0 < q < 2 and non-zero Dirichlet boundary conditions. The results of the experiments conducted herein reveal that the proposed schemes can actually be used for calculation in the case of 0 < q < 1 and can be used to treat non-zero Dirichlet boundary conditions. In contrast, the existing finite difference methods of Meerschaert et al. are applicable only for 1 < q < 2 and zero-Dirichlet boundary conditions. Using the proposed Schemes A and B, we can calculate space-fractional differential equations under wider conditions than the existing methods. However, there are other types of boundary conditions, for example, Neumann or Robin boundary conditions, in real applications. Therefore, in the future, we intend to develop finite difference schemes that can deal with various boundary conditions. Appendix A. Proof of Theorem 1 Theorem 1 is described by the proceeding[Takeuchi and Suda (2012)], which is not available on the Internet. Therefore, the proof is provided in the Appendix. Let 2LDxq f (x) be the second order accuracy difference formula as ] [ N −1 h−q ∑ Γ(j − q) qh ′ 2 q f (x − jh) + f (x − jh) , LDx f (x) = Γ(−q) j=0 Γ(j + 1) 2 24

and we prove the following expression q LDx f (x)

−

2 q LDx f (x)

=O

(

1 N2

)

.

Proof First, we consider the fractional integral for q < 0. The Riemann-Liouville definition for a fractional integral is given by q LDx f (x)

1 = Γ(−q)

∫

x

L

f (ξ) dξ (x − ξ)1+q

for q < 0. By applying changing variables for x − ξ = u and dividing the integral into N parts, we have 1 = Γ(−q)

∫

h

0

N −1 ∫ f (x − u) 1 ∑ (j+1)h f (x − u) du + du u1+q Γ(−q) j=1 jh u1+q

where h is h = (x − L)/N . By applying changing variables for u = jh − t and applying a Taylor expansion to the function g(y) = f (y)/(x − y)1+q , we have 1 = lim Γ(−q) ϵ→0

∫

N −1 ∫ 1 ∑ 0 g(x − jh + t)dt g(x − h + t − ϵ)dt + Γ(−q) j=1 −h −h ϵ

N −1 ∞ hn+1 1 ∑ ∑ (n) g (x − jh)(−1)n Γ(−q) j=0 n=0 (n + 1)! [ N −1 ∞ ∏ 1 ∑ h−q ∑ nk=1 (k − 1 + q) f (x − jh) q = (−1)n+1 Γ(−q) j=0 j n=1 qj n n! ] ( ) ∞ ∏n 1−q ∑ (k − 1 + q) h 1 ′ n+1 k=1 −f (x − jh) q n(−1) +O . n j n=1 qj (n + 1)! N2

=

The infinite summations in the above expressions can be transformed into finite summations using Taylor expansions, which gives q LDx f (x)

25

N −1 [ ∑

(j + 1)−q − j −q = h f (x − jh) Γ(1 − q) j=0 { }] (j + 1)−q (j + 1)1−q j 1−q ′ +hf (x − jh) − + − Γ(1 − q) Γ(2 − q) Γ(2 − q) ( ) 1 +O . N2 −q

By applying a Taylor expansion to f , we have

=

q LDx f (x) N −1 [ ∑ −q

h f (x − jh)

j=0

(j + 1)−q − j −q Γ(1 − q)

}] (j + 1)−q (j + 1)1−q j 1−q +h {f (x − jh) − f (x − (j + 1)h)} − + − Γ(1 − q) Γ(2 − q) Γ(2 − q) ( ) 1 +O N2 N −1 [ ∑ (j + 1)−q − j −q = h−q f (x − jh) Γ(1 − q) j=0 { }] (j + 1)−q − j −q (j + 1)1−q − j 1−q j 1−q − (k(j − 1))1−q −q +h f (x − jh) − + − Γ(1 − q) Γ(2 − q) Γ(2 − q) { } ( ) −q 1−q 1−q N 1 N (N − 1) −h−q f (L) − + − +O Γ(1 − q) Γ(2 − q) Γ(2 − q) N2 }] { N −1 [ ∑ (j + 1)1−q − j 1−q j 1−q − (k(j − 1))1−q −q − = h f (x − jh) Γ(2 − q) Γ(2 − q) j=0 ( ) 1 −1−q 1 −q +h f (L) N +O 2 N2 −q

where k(j − 1) =

{

{

0, j=0 j − 1, otherwise.

By comparing the coefficients of f to the following expression 2 q LDx f (x)

26

[ ] N −1 h−q ∑ Γ(j − q) qh ′ = f (x − jh) + f (x − jh) Γ(−q) j=0 Γ(j + 1) 2

h−q 1 + q f (L)N −1−q Γ(−q) 2 [ ] N −1 ∑ Γ(j − q) q(1 + q) Γ(j − 1 − q) −q = h f (x − jh) − Γ(−q)Γ(j + 1) 2 Γ(−q)Γ(j + 1) j=0 ( ) h−q 1 1 + f (L)N −1−q + O , Γ(−q) 2 N2 +

we obtain 2 q LDx f (x) N −1 ∑ −q

[{

(j + 1)1−q − j 1−q j 1−q − (k(j − 1))1−q = h f (x − jh) − Γ(2 − q) Γ(2 − q) j=0 ] Γ(j − q) q(1 + q) Γ(j − 1 − q) − + Γ(−q)Γ(j + 1) 2 Γ(−q)Γ(j + 1) ∞ ∑ (−1)n = h−q (x − L)n f (n) (x) n! n=0 )n [{ } N −1 ( ∑ j (j + 1)1−q − j 1−q j 1−q − (k(j − 1))1−q · − N Γ(2 − q) Γ(2 − q) j=0 ] Γ(j − q) q(1 + q) Γ(j − 1 − q) − + . Γ(−q)Γ(j + 1) 2 Γ(−q)Γ(j + 1)

}

By applying asymptotic expansions to the second summation, we have )n { } N −1 ( ∑ j (j + 1)1−q − j 1−q j 1−q − (k(j − 1))1−q − N Γ(2 − q) Γ(2 − q) j=0

1 1 [ n+1−q N − (N − 1)n+1−q − nN n−q − n(N − 1)n−q Γ(2 − q) N n ) ] ( n(n − 1) n−q n(n − 1) 1 n−1−q N − (N − 1) + +O n−q 2 N 2+q ] ( ) [ 1 1 (1 − q)(−q) −q q(1 − q) −1−q = N + N +O , Γ(2 − q) n−q 2 N 2+q =

27

and N −1 ( ∑ j=0

j N

)n {

Γ(j − q) q(1 + q) Γ(j − 1 − q) − Γ(−q)Γ(j + 1) 2 Γ(−q)Γ(j + 1)

}

[ ] Γ(N − q) 1 1 n(n − 1) Γ(N − q) = + Γ(−q) N n (n − q)Γ(N − n) 2 (n − 1 − q)Γ(N − n + 1) ( ) q(1 + q) 1 1 Γ(N − 1 − q) 1 − +O 2 Γ(−q) N n Γ(n − 1 − q)Γ(N − n) N 2+q [ −q ] ( ) 1 N 1 1 = − N −1−q + O . Γ(−q) n − q 2 N 2+q

Therefore, we have

q LDx f (x)

−

2 q LDx f (x)

=O

(

1 N2

)

and the theorem is proved for q < 0. Next, let us consider the case of q > 0. By taking the derivatives of both sides, we have the following expressions: [ ] d 2 q D f (x) dx L x [ ] [ ] −q N −1 ∑ h Γ(j − q) d qh ′ f (x − jh) + f (x − jh) = dx Γ(−q) j=0 Γ(j + 1) 2 [ ] −q ( ) d h 1 1+q + f (L)N −1−q + O dx Γ(−q) 2 N2 [ ] N −1 1 ∑ Γ(j − q) −q −q−1 q(1 − q) −q ′ = h f (x − jh) + h f (x − jh) Γ(−q) j=0 Γ(j + 1) N 2N [ ] N −1 1 ∑ Γ(j − q) −q ′ N − j q(N − j) −q+1 ′′ + h f (x − jh) + h f (x − jh) Γ(−q) j=0 Γ(j + 1) N 2N ( ) h−1−q −q(1 + q) 1 −2−q + f (L)N +O . Γ(−q) 2 N2 By applying the following Taylor expansions hf ′ (x − jh) = f (x − jh) − f (x − (j + 1)h) + 28

h2 ′′ f (x − jh) + O(h3 ) 2

hf ′′ (x − jh) = f ′ (x − jh) − f (x − (j + 1)h) + O(h2 ), we have

{ [ } N −1 1 ∑ Γ(j − q) −q−1 N − j − q N −j = h f (x − jh) − f (x − (j + 1)h) Γ(−q) j=0 Γ(j + 1) N N } { 1 + q −q N − j − q ′ N −j ′ + h f (x − jh) − f (x − (j + 1)h) 2 N N ] q + h−q−1 {f (x − jh) − f (x − (j + 1)h)} N ( ) h−1−q −q(1 + q) 1 −2−q + f (L)N +O Γ(−q) 2 N2 { } N −1 [ 1 ∑ −q−1 N − j − q Γ(j − q) N − j + 1 Γ(j − 1 − q) = h f (x − jh) − Γ(−q) j=0 N Γ(j + 1) N Γ(j) { } 1 + q −q ′ N − j − q Γ(j − q) N − j + 1 Γ(j − 1 − q) + h f (x − jh) − 2 N Γ(j + 1) N Γ(j) { }] ( ) Γ(j − q) Γ(j − 1 − q) 1 q + h−q−1 f (x − jh) − +O N Γ(j + 1) Γ(j) N2 h−1−q −q(1 + q) h−1−q h−1−q + f (L)N −2−q − f (L)N −2−q − qf (L)N −2−q Γ(−q) 2 Γ(−q) Γ(−q) [ ] N −1 h−1−q ∑ Γ(j − 1 − q) (1 + q)h ′ = f (x − jh) + f (x − jh) Γ(−1 − q) j=0 Γ(j + 1) 2 ( ) h−1−q 2 + q 1 −2−q + f (L)N +O Γ(−1 − q) 2 N2 ( ) 1 = 2LDx1+q f (x) + O N2

for −1 < q < 0. By repeating this operation, we have the theorem for arbitrary q. □ References Barkai, E., Metzler, R., Klafter, J., 2000. From continuous time random walks to the fractional Fokker-Planck equation. Physical Review E 61 (1), 132. 29

Chaves, A., 1998. A fractional diffusion equation to describe L´evy flights. Physics Letters A 239 (1), 13–16. Chechkin, A. V., Gonchar, V. Y., Klafter, J., Metzler, R., 2006. Fundamentals of L´evy flight processes. Advances in Chemical physics 133 (B), 439. Chen, M., Deng, W., 2014a. Fourth order difference approximations for space riemann-liouville derivatives based on weighted and shifted lubich difference operators. commun. Comput. Phys 16 (2), 516–540. Chen, M., Deng, W., 2014b. A second-order numerical method for twodimensional two-sided space fractional convection diffusion equation. Applied Mathematical Modelling 38 (13), 3244–3259. Chen, M., Deng, W., Wu, Y., 2013. Superlinearly convergent algorithms for the two-dimensional space–time caputo–riesz fractional diffusion equation. Applied Numerical Mathematics 70, 22–41. de Pablo, A., Quir´os, F., Rodriguez, A., V´azquez, J. L., 2011. A fractional porous medium equation. Advances in Mathematics 226 (2), 1378–1409. Diethelm, K., Ford, N. J., Freed, A. D., 2004. Detailed error analysis for a fractional Adams method. Numerical Algorithms 36 (1), 31–52. Jeong, D., Ha, T., Kim, M., Shin, J., Yoon, I.-H., Kim, J., 2014. An adaptive finite difference method using far-field boundary conditions for the black-scholes equation. Bulletin of the Korean Mathematical Society 51 (4), 1087–1100. L¨otstedt, P., Persson, J., von Sydow, L., Tysk, J., 2007. Space–time adaptive finite difference method for european multi-asset options. Computers & Mathematics with Applications 53 (8), 1159–1180. Meerschaert, M. M., Benson, D. A., B¨aumer, B., 1999. Multidimensional advection and fractional dispersion. Physical Review E 59 (5), 5026. Meerschaert, M. M., Scheffler, H.-P., Tadjeran, C., 2006. Finite difference methods for two-dimensional fractional dispersion equation. Journal of Computational Physics 211 (1), 249–261.

30

Meerschaert, M. M., Tadjeran, C., 2004. Finite difference approximations for fractional advection–dispersion flow equations. Journal of Computational and Applied Mathematics 172 (1), 65–77. Meerschaert, M. M., Tadjeran, C., 2006. Finite difference approximations for two-sided space-fractional partial differential equations. Applied Numerical Mathematics 56 (1), 80–90. Metzler, R., Barkai, E., Klafter, J., 1999. Anomalous diffusion and relaxation close to thermal equilibrium: a fractional Fokker-Planck equation approach. Physical Review Letters 82 (18), 3563. Metzler, R., Klafter, J., 2000. The random walk’s guide to anomalous diffusion: a fractional dynamics approach. Physics Reports 339 (1), 1–77. Metzler, R., Klafter, J., 2004. The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics. Journal of Physics A: Mathematical and General 37 (31), R161. Oberman, A. M., Zwiers, I., 2015. Adaptive finite difference methods for nonlinear elliptic and parabolic partial differential equations with free boundaries. Journal of Scientific Computing, 1–21. Sousa, E., 2009. Finite difference approximations for a fractional advection diffusion problem. Journal of Computational Physics 228 (11), 4038–4054. Tadjeran, C., Meerschaert, M. M., 2007. A second-order accurate numerical method for the two-dimensional fractional diffusion equation. Journal of Computational Physics 220 (2), 813–823. Takeuchi, Y., Suda, R., 2012. New numerical computation formula and error analysis of some existing formulae in fractional derivatives and integrals. In: The 5th IFAC Symposium on Fractional Differentiation and its Applications (FDA’12), Hohai University, Nanjing, China. Takeuchi, Y., Suda, R., 2013. Second order accuracy finite difference methods for fractional diffusion equations. In: ASME 2013 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference. American Society of Mechanical Engineers, pp. V004T08A015–V004T08A015. 31

Tian, W., Zhou, H., Deng, W., 2015. A class of second order difference approximations for solving space fractional diffusion equations. Mathematics of Computation 84 (294), 1703–1727. West, B. J., Grigolini, P., Metzler, R., Nonnenmacher, T. F., 1997. Fractional diffusion and L´evy stable processes. Physical Review E 55 (1), 99. Yuste, S. B., Quintana-Murillo, J., 2012. A finite difference method with non-uniform timesteps for fractional diffusion equations. Computer Physics Communications 183 (12), 2594–2600. Zhang, Y.-N., Sun, Z.-Z., Liao, H.-L., 2014. Finite difference methods for the time fractional diffusion equation on non-uniform meshes. Journal of Computational Physics 265, 195–210.

32