High-order central difference scheme for Caputo fractional derivative

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 317 (2017) 42–54 www.elsevier.com/locate/cma

High-order central difference scheme for Caputo fractional derivative Yuping Ying a,c , Yanping Lian c,1 , Shaoqiang Tang a,b,∗ , Wing Kam Liu c a HEDPS, CAPT, LTCS and College of Engineering, Peking University, Beijing 100871, China b IFSA Collaborative Innovation Center of MoE, Shanghai Jiao Tong University, Shanghai 200240, China c Department of Mechanical Engineering, Northwestern University, Evanston IL 60201, USA

Received 6 July 2016; received in revised form 1 November 2016; accepted 5 December 2016 Available online 13 December 2016

Abstract In this paper we propose a class of central difference schemes for resolving the Caputo fractional derivative. The accuracy may reach any selected integer order. More precisely, the Caputo fractional derivative operator is decomposed into symmetric and antisymmetric components. Starting from difference schemes of lower order accuracy for each component, we enhance the accuracy by a weighted average of shifted differences. The weights are calculated by matching the symbols of the scheme and the operators. We further illustrate the application of the proposed schemes to a fractional advection–diffusion equation. Together with the Crank–Nicolson algorithm, it reaches designed accuracy order, and is unconditionally stable. Numerical tests are presented to demonstrate the nice features. c 2016 Elsevier B.V. All rights reserved. ⃝

Keywords: Caputo fractional derivative; Fractional calculus; Central difference scheme; Advection–diffusion equation

1. Introduction In recent two decades, fractional partial differential equations (FPDE) have caught a lot of attentions due to their inherent ability for modelling non-local phenomena in various fields of science and engineering [1,2]. Meanwhile, the non-local feature brings about many challenges to the existing numerical methods in terms of memory storage, computational cost, accuracy, etc. In this paper, we propose a class of high order central difference schemes, and adopt them to construct unconditionally stable finite difference algorithm for solving a fractional advection–diffusion equation (FADE) as follows, which was proposed to describe anomalous diffusion [3].  α C α ∂t φ(x, t) = κ1 C 0 Dx φ(x, t) + κ2 x D L φ(x, t) − u∂x φ(x, t) + g(x, t), (1) φ(0, t) = 0; φ(L , t) = 0;  φ(x, 0) = φ0 (x), ∗ Corresponding author at: HEDPS, CAPT, LTCS and College of Engineering, Peking University, Beijing 100871, China.

E-mail address: [email protected] (S. Tang). 1 Contributing equally with Yuping Ying. http://dx.doi.org/10.1016/j.cma.2016.12.008 c 2016 Elsevier B.V. All rights reserved. 0045-7825/⃝

Y. Ying et al. / Comput. Methods Appl. Mech. Engrg. 317 (2017) 42–54

43

where κ1 ≥ 0, κ2 ≥ 0 are diffusion coefficients, and u is an advection velocity. The left-sided and right-sided Caputo α C α fractional derivative operators C 0 Dx and x D L of a function φ(x) are defined for 1 < α < 2 as follows. x

C α 0 Dx φ(x)

1 Γ (2 − α)



=

C α x D L φ(x)

1 = Γ (2 − α)



(x − η)1−α

0 L x

(η − x)1−α

d2 φ(η)dη, dη2

(2)

d2 φ(η)dη. dη2

(3)

Being one of the dominant numerical tools in scientific computing and engineering applications, finite difference method (FDM) is naturally adopted to solve fractional differential equations (FDE) for its simplicity to apply. For the two most widely used fractional derivatives, the Caputo and Riemann–Liouville fractional derivatives, there exist mainly three types of FDM classified by how they are derived. The basic idea of the first type of FDM is direct numerical approximation to the expression of fractional derivatives. Among this type of methods, a certain kind of schemes called L1 formula can achieve a numerical accuracy of order O(h 2−α ) [4–7], where 0 < α < 1 is the fractional order and h stands for the grid size. This formula is established by a piecewise linear interpolation approximation to the integrand function. To improve the numerical accuracy, Gao [8] discretized the Caputo fractional derivative with a numerical accuracy of O(h 3−α ) by constructing a quadratic interpolation for the integrand function. Then Alikhanov [9] improved Gao’s scheme to prove the stability and convergence. More similar works achieving higher order up to O(h 6−α ) can be found in [10,11]. The second type of FDM is obtained by means of the Gr¨unwald–Letnikov fractional derivative, which is defined as a finite difference. Igor [12] used it as a numerical approximation for the Riemann–Liouville fractional derivative with an accuracy order of O(h). To improve the accuracy order and achieve stability for solving related differential equations, the Gr¨unwald–Letnikov fractional derivative is shifted and summed with some weights, which leads to a class of second order discretizations [13,14]. The third type of FDM is based on the fundamental work by Lubich [15], who introduced fractional linear multi-step methods to approximate fractional integrals. Later on, Zeng [16] solved a time-fractional diffusion equation with unconditional stability, where the Caputo fractional derivative was approximated by linear multi-step methods with second order accuracy. Chen [17] derived a class of fourth order approximations for the Riemann–Liouville derivative by the same methods. More finite difference schemes for fractional derivatives have been discussed in Li’s recent book [18]. In this paper, a new type of high-order central difference schemes for Caputo fractional derivative are proposed. To this end, Caputo fractional derivative is decomposed into a symmetric operator and an antisymmetric one. The symmetric part is actually related to Riesz fractional derivative up to a constant. Based on some existing low order finite difference schemes, a class of central difference schemes are designed as a weighted average of some shifted differences. The weights are obtained by matching the symbols of the schemes and the operators. We remark that these symbols are also referred to as the generating functions of finite difference schemes for Riesz fractional derivative [19–21]. Such designs can reach any desired order of accuracy, up to infinity which means exact resolution. We present the weights for finite difference schemes of up to O(h 8 ) as examples. The ‘infinite-order’ finite difference for Caputo derivative is also discussed. After this, by using the Crank–Nicolson algorithm for time integration, we present a numerical algorithm to solve (1). We also show the unconditional stability by the von Neumann analysis, and display the convergence rate of the algorithm. Algorithms of order O(h p ), where p = 2, 4, 6, 8, +∞, are shown as examples. Compared with the existing FDM in the literature, the new schemes have the following features: • The Caputo fractional derivatives of any order α > 0 can be approximated at accuracy of any desired even integer, up to infinite order. • The accuracy order of the scheme is not affected by the order of the fractional derivative, which is different from the aforementioned first type of FDM. We remark that the mechanism is similar to the Moving Least Square Reproducing Kernel interpolation [22]. • The structure of the scheme is symmetric, with all information on both sides of the reference point being used. The rest of this paper is organized as follows. In Section 2, Caputo fractional derivative is decomposed. Then the schemes for each component are obtained by matching symbols, leading to high-order schemes for Caputo fractional derivative. In Section 3, the proposed high-order schemes are used to solve (1), where unconditional stability and convergence of the algorithm are proved. In Section 4, two numerical examples are presented. Finally, the paper is concluded in Section 5.

44

Y. Ying et al. / Comput. Methods Appl. Mech. Engrg. 317 (2017) 42–54

2. High-order difference for Caputo fractional derivative In order to design a high-order finite difference scheme for the Caputo derivative, we first decompose it into a symmetric component and an antisymmetric one. Taking a weighted average of shifting second order differences, we design difference schemes of any desired order for both parts. 2.1. A decomposition for Caputo derivative Our method applies to any fractional order α > 0 without the requirement of between 1 and 2 as that in (1). Consider the left-sided Caputo derivative  x dn φ(η) 1 C α (x − η)n−α−1 dη, (4) −∞ Dx φ(x) = Γ (n − α) −∞ dηn and the right-sided one C α x D+∞ φ(x)

=

(−1)n Γ (n − α)

+∞



(η − x)n−α−1 x

dn φ(η) dη dηn

(5)

where n is a positive integer and n − 1 < α < n. Here we remark that when the integration terminal contains −∞ or +∞, the Caputo fractional derivative, Riemann–Liouville fractional derivative and Gr¨unwald–Letnikov fractional derivative share the same form [12]. We define a symmetric Caputo differential operator  1 C α α C α D+ φ(x) = (6) −∞ Dx φ(x) + x D+∞ φ(x) , 2 and an antisymmetric one  1 C α α C α D− φ(x) = D φ(x) − D φ(x) . (7) x +∞ 2 −∞ x Then the left-sided and right-sided operators are naturally evaluated from C α −∞ Dx φ(x) C α x D+∞ φ(x)

α α = D+ φ(x) + D− φ(x),

=

α D+ φ(x) −

α D− φ(x).

(8) (9)

The symbols  +∞of the symmetric and antisymmetric operators may be readily found by applying the Fourier transform F[φ(x)] = −∞ φ(x)e−iξ x d x to (6) and (7). They are respectively απ α |ξ | , (10) 2 απ α α F[D− ] = f − (ξ ) = i sin |ξ | sign(ξ ). (11) 2 For 1 < α < 2, the symbol for symmetric component is real and negative, hence the operator is dissipative. With a pure imaginary symbol, in contrast, the antisymmetric one is dispersive. The sum of κ1 and κ2 dictates how strong the dissipation is. The difference between them determines how dispersive the system behaves. α F[D+ ] = f + (ξ ) = cos

2.2. ‘Infinite-order’ finite difference resolution For a finite difference scheme with a uniform grid of size h, the Nyquist law for sampling dictates that the maximal resolution is to reproduce the differentiation operator exactly within the first Brillouin zone for the wave number ξ ∈ [−π/ h, π/ h], or ξ h ∈ [−π, π]. This resolution for the aforementioned two components is reached in finite differences α D+ φ(x) ∼

∞  k=−∞

(∞)

Ck

φ(x − kh),

(12)

Y. Ying et al. / Comput. Methods Appl. Mech. Engrg. 317 (2017) 42–54

45

with (∞) Ck

sin απ 2 =− π

π



z α sin(kz)dz

0

iπ α sin απ 2 = [1 F1 (α + 1; α + 2; ikπ ) − 1 F1 (α + 1; α + 2; −ikπ )] , 2(α + 1)

(13)

and α D− φ(x) ∼

∞ 

(∞)

ck

φ(x − kh),

(14)

k=−∞

with (∞) ck

cos απ 2 = π



π

z α cos(kz)dz

0

π α cos απ 2 = [1 F1 (α + 1; α + 2; ikπ ) + 1 F1 (α + 1; α + 2; −ikπ )] . 2(α + 1) ∞ 

(15)

is the hypergeometric function, where (a)k = ΓΓ(a+k) (a) = a(a + 1) · · · (a + k − 1). It is clear that the symbols for them exactly reproduce (10) and (11) in the first Brillouin zone.

Here 1 F1 (a; b; z) =

k=0

(a)k z k (b)k k!

2.3. A second order difference for the antisymmetric component It is not straightforward to obtain schemes of finite order accuracy from the infinite-order ones in the previous subsection. Instead, we take a bottom-up approach to construct it from lower order ones. For the antisymmetric operator (7), a finite difference approximation +∞ ∆α− φ(x) 1  ≡ Ck φ(x − kh) hα h α k=−∞

takes a symbol

F− (ξ h) hα

(16)

with

F[∆α− ] = F− (ξ h) =

+∞ 

Ck e−ikξ h .

(17)

k=−∞

Noticing the one-to-one correspondence between the symbol and the scheme, F− (ξ h) may be used to determine the coefficients Ck and hence referred to as the generating function of the finite difference discretization. It has a period of 2π , with the first Brillouin zone taken as [−π, π]. The design may then be performed by approximating the symbol for the Caputo derivative with that for the finite difference. In fact, we may define the accuracy for the discretization by the error order of this approximation. More precisely, the finite difference approximation (16) is of order p if   F− (ξ h) = 1 + O (ξ h) p . h α f − (ξ )

(18)

The above Eq. (18) can be comprehended as that F−h(ξα h) gives a pth-order approximation to f − (ξ ) for an arbitrary fixed wavenumber ξ . Particularly, given a function φ(x) ∈ C p+3 (R) with all its derivatives up to order p+3 belonging to L 1 , it can be proved that (18) leads to [20] ∆α− φ(x) α − D− φ(x) = O(h p ). hα

(19)

46

Y. Ying et al. / Comput. Methods Appl. Mech. Engrg. 317 (2017) 42–54

(2)

Fig. 1. Comparison of F− (z) and h α f − (ξ ) = f − (z) for α = 1.5.

It is noticed that for α = 1, (19) recovers the traditional definition ∆1− φ(x) 1 − D− φ(x) = O(h p ). h We start with a symbol for finite difference of second order accuracy. It reads (with z = ξ h)  απ  −i sin | sin(z)|α , −π ≤ z ≤ 0; (2) 2 F− (z) = i sin απ | sin(z)|α , 0 ≤ z ≤ π.  2

(20)

(21)

The second order accuracy follows from (2)

(2)

F− (z) F− (z) | sin(z)|α α = = = 1 − z 2 + o(z 2 ). h α f − (ξ ) f − (z) |z|α 6

(22) (2)

See Fig. 1 for an illustration with α = 1.5. By inverse Fourier transform of F− (z), the coefficients are obtained as follows.  π sin απ (2) 2 Ck = − (sin(z))α sin(kz)dz π 0 =−

kπ sin απ 2 sin 2 Γ (α + 1)  α+k   . 2α Γ 2 + 1 Γ α−k 2 +1

(23)

2.4. High order difference for the antisymmetric operator We consider a superposition of shifted second order discretization, with p an even integer. p

∆α− φ(x)

∞  2 −1  1 ( p) (2) B j Ck [φ(x − (k + j)h) + φ(x − (k − j)h)]. = 2 k=−∞ j=0

(24)

Its symbol is p

F− (z) =

2 −1 

( p)

Bj

(2)

cos( j z)F− (z).

(25)

j=0

The Taylor expansion for p   2 −1   sin z α F− (hξ ) ( p)   = B j cos( j z) ·   h α f − (ξ ) z j=0

(26)

Y. Ying et al. / Comput. Methods Appl. Mech. Engrg. 317 (2017) 42–54

47

Fig. 2. Comparison of F− (z) and h α f − (ξ ) = f − (z) for p = 2, 4, 6, 8 when α = 1.5.

Table 1 Coefficients for central difference of the antisymmetric operator. ( p)

( p)

B0

B1

p=2 p=4

1 1 + 31 α

− 13 α

p=6

9 α + 1 α2 1 + 20 12 5813 α + 4 α 2 + 5 α 3 1 + 11340 27 324

p=8

( p)

p=6 p=8

1 2 − 22 45 α − 9 α 5 2 5 3 − 2203 3780 α − 24 α − 216 α ( p)

B2

B3

7 1 2 180 α + 36 α 1 2 1 3 289 3780 α + 15 α + 108 α

71 α − 7 α 2 − 1 α 3 − 11340 1080 648

Table 2 Coefficients for central difference of symmetric operator [21]. ( p)

( p)

b0

b1

p=2 p=4

1 1 α 1 + 12

1 α − 12

p=6

17 α + 1 α 2 1 + 160 192 5297 α + 29 α 2 + 5 α 3 1 + 45360 3456 20736

p=8

( p)

p=6 p=8

41 α − 1 α 2 − 360 144 7843 α − 3 α 2 − 5 α 3 − 60480 256 13824 ( p)

b2

b3

11 1 2 1440 α + 576 α 211 7 1 2 3 15120 α + 1920 α + 6912 α

191 α − 11 α 2 − 1 α 3 − 181440 34560 41472

reads p     l  ∞ ∞    ∞ 2 −1 m j 2m z 2m m z 2m    (−1) α (−1) ( p)   1 + . Bj 1+ (2m)! l (2m + 1)! j=0 m=1 l=1 m=1

(27)

Requiring (27) to be in the form of 1 + O(z p ), we obtain 2p equations as only even powers in z appear. This ( p) determines the coefficients B j for p = 0, . . . , p/2 − 1, which are listed in Table 1 for p = 2, 4, 6, 8. As shown in Fig. 2, higher order discretizations result in better approximations to the generating function f − (z).

48

Y. Ying et al. / Comput. Methods Appl. Mech. Engrg. 317 (2017) 42–54

2.5. High order difference for the Caputo derivative α . Actually, the symmetric In the same way, high order discretization is designed for the symmetric operator D+ operator is related to the Riesz fractional derivative up to a constant

απ α ∂ α φ(x) = − sec D φ(x). α ∂|x| 2 +

(28)

The central difference schemes for Riesz fractional derivative have been constructed by [20,21]. We briefly introduce the process to obtain the schemes. We start with a second-order finite difference +∞ ∆α+ φ(x) 1  (2) ≡ c φ(x − kh), hα h α k=−∞ k

(29)

where the coefficients are (2)

ck =

(−1)k cos απ (α + 1) 2 Γ α . α Γ ( 2 − k + 1)Γ 2 + k + 1

(30)

The symbol is z α απ  2 sin  , −π ≤ z ≤ π. 2 2 Now we consider a weighted average for shift differences (2)

F+ (z) = cos

∆α+ φ(x) =

∞ 

p/2−1 

k=−∞ j=0

1 ( p) (2) b c [φ(x − (k + j)h) + φ(x − (k − j)h)] . 2 j k

(31)

(32)

Its symbol reads p

F+ (z) =

2 −1 

( p)

(2)

b j cos( j z)F+ (z).

(33)

j=0

The Taylor expansion for p   2 −1   2 sin 2z α F+ (hξ ) ( p)   = b cos( j z) · j  z  h α f + (ξ ) j=0

(34)

p       l  ∞    ∞ ∞ m z 2m 2 −1 m j 2m z 2m    (−1) α (−1) ( p) 2   1 + . bj 1+ (2m)! l (2m + 1)! j=0 m=1 l=1 m=1

(35)

is

Let (35) to be of the form 1 + O(z p ) and we calculate the weights b j . In Table 2, the coefficients for p = 2, 4, 6, 8 are listed as examples. Fig. 3 shows the approximation of F+ (z) to f + (z) for α = 1.5. By all the analyses above, the high order difference schemes for both the left-sided and right-sided Caputo derivatives are linear combinations of the discretized symmetric and antisymmetric operators, as suggested by (8)–(9). 3. Application to FADE To solve (1) in [0, L] with the proposed schemes constructed for (4) and (5), we extend the function φ(x, t) and its initial data φ0 (x) to be 0 outside of [0, L]. Then the fractional derivatives in (1) may be put as the left-sided and right-sided ones, and discretized in the above-mentioned manner.

Y. Ying et al. / Comput. Methods Appl. Mech. Engrg. 317 (2017) 42–54

49

Fig. 3. Comparison of F+ (z) and h α f + (ξ ) = f + (z) for p = 2, 4, 6, 8 when α = 1.5.

Consider an equidistance space grid with h = L/N , and denote xn = nh, φn (t) = φ(xn , t). φ˙ n (t) =

1 2h α

p/2−1  

n−1 

 ( p) ( p) (2) (2) (2) (2) (κ1 + κ2 )b j (cm− j + cm+ j ) + (κ1 − κ2 )B j (Cm− j + Cm+ j ) φn−m

m=n−N +1 j=0 ( p)

− u D1 φn + g(xn , t). ( p)

Here for an even integer p, D1 φn =

(36) 1 h

p/2 

ak (φn+k − φn−k ) is the standard pth order central difference for the first

k=1

order derivative. The coefficients ak are solved from 2

p/2 

p/2 

ak k = 1,

( j = 1, . . . , p/2 − 1).

ak k 2 j+1 = 0,

(37)

k=1

k=1

Its symbol is

F1 (ξ h) h

=

1 h

p/2 

2iak sin(khξ ).

k=1

The linear system (36) may be studied by Fourier analysis. Noticing the homogeneous Dirichlet boundary condition, we take a sine-transform in space φn (t) =

N −1 

φ˜ m (t) sin(ξm n),

m=1

ξm =

mπ . N

The amplitude φ˜ m (t) then solves   (κ1 + κ2 )F+ (ξm h) (κ1 − κ2 )F− (ξm h) F1 (ξm h) ˙ ˜ φ m (t) = φ˜ m (t) + g˜ m . + −u hα hα h

(38)

(39)

Here it can be readily verified that F+ (z) is real and non-positive for p = 2, 4, 6, 8, +∞ and α between 1 and 2 (see Appendix). Besides, both F− (z) and F1 (z) are pure imaginary. Hence we know the semi-discrete dynamical system (36) is stable. For larger order p, Fig. 3 shows the trend of F+ (z) to approach f + (z) with p increasing. Then the stability of (36) is also expected. Further, (36) may be recast to a matrix form ˙ Φ(t) = AΦ(t) + G(t),

(40)

where Φ = (φ1 , . . . , φ N −1 )T , G = (g1 , . . . , g N −1 )T , and A represents the constant coefficient matrix for the derivatives. For this stable system, we adopt a Crank–Nicolson algorithm for time integration. 1 s+1 1 (Φ − Φ s ) = (A(Φ s+1 + Φ s ) + G s+1 + G s ). τ 2

(41)

50

Y. Ying et al. / Comput. Methods Appl. Mech. Engrg. 317 (2017) 42–54

By von Neumann analysis, the amplification matrix is −1    1 1 I − τA I + τA . 2 2

(42)

Its eigenvectors are the Fourier modes while the eigenvalues are 2 + τ λm , 2 − τ λm

m = 1, . . . , N − 1,

(43)

with F1 (ξm h) (κ1 + κ2 )F+ (ξm h) (κ1 − κ2 )F− (ξm h) + −u . (44) α α h h h They are readily found to be with modulus no greater than 1. Therefore the Crank–Nicolson discretization with the designed high-order difference in space provides an unconditionally stable scheme with an accuracy of pth order in space and second order in time. This unconditional stability can be demonstrated by confirming F+ (z) is non-positive for a required p. We also remark that this stability analysis can be applied to other types of FDE discretized by the proposed central differences. λm =

4. Numerical examples Two numerical examples are presented to illustrate the proposed high order central difference schemes. Example 1.  C 1.5 C 1.5  ∂t φ(x, t) = 0.8 0 Dx φ(x, t) + 0.2 x D L φ(x, t) − ∂x φ(x, t) + f (x, t), φ(0, t) = 0; φ(2, t) = 0, 0 ≤ t ≤ 0.5,   φ(x, 0) = x 8 (2 − x)8 , 0 ≤ x ≤ 2,

(45)

where f (x, t) = −x 8 (2 − x)8 sin t −

8  k=0

cos t

  8 (−1)k 28−k k

  Γ (9 + k) Γ (9 + k) 6.5+k 6.5+k 7+k x + 0.2 (2 − x) − (8 + k)x . 0.8 Γ (7.5 + k) Γ (7.5 + k)

(46)

The solution is φ(x, t) = x 8 (2 − x)8 cos t.

(47)

For this example, the scheme (41) is applied with various temporal and spatial mesh sizes. The maximal absolute error at t = 0.5 is used to calculate the numerical convergence rate. As listed in Table 3, the numerical results demonstrate that the convergence rate are consistent with the designed accuracy order in Section 2. For the ‘infinite-order’ finite difference scheme, the maximal absolute error at t = 0.5 is plotted in terms of N = h2 which is the number of grid points. In Fig. 4, we observe an exponential convergence, which is so for a spectral method. For another example, we consider long-time computation for a pure fractional diffusion equation. Example 2.  C 1.5 C 1.5  ∂t φ(x, t) = κ1 0 Dx φ(x, t) + κ2 x D L φ(x, t), φ(0, t) = 0; φ(20, t) = 0, 0 ≤ t ≤ 10,  2  φ(x, 0) = e−(x−30) /4 , 0 < x < 60.

(48)

Y. Ying et al. / Comput. Methods Appl. Mech. Engrg. 317 (2017) 42–54

51

Fig. 4. For the ‘infinite-order’ finite difference scheme, the maximal absolute error decreases with N = 2/ h exponentially. Table 3 The maximal absolute error and convergence rate of scheme (36) for solving (45). O(τ 2 + h 6 )

τ

Error at t = 0.5

Rate in time

h = 1/48

1/20 1/40 1/80 1/160

1.96e−4 4.89e−5 1.22e−5 3.06e−6

2.00 2.00 1.99

O(τ 2 + h 2 )

h

Error at t = 0.5

Rate in space

τ = 1/10000

1/8 1/16 1/32 1/64

1.70e−2 4.30e−3 1.10e−3 2.75e−4

1.98 1.97 2.00

O(τ 2 + h 4 )

h

Error at t = 0.5

Rate in space

τ = 1/10000

1/8 1/16 1/32 1/64

2.37e−3 1.70e−4 1.10e−5 6.97e−7

3.80 3.95 3.98

O(τ 2 + h 6 )

h

Error at t = 0.5

Rate in space

τ = 1/10000

1/6 1/12 1/24 1/48

2.20e−3 5.34e−5 9.51e−7 1.77e−8

5.36 5.81 5.75

O(τ 2 + h 8 )

h

Error at t = 0.5

Rate in space

τ = 1/10000

1/8 1/10 1/12 1/14

1.29e−4 2.64e−5 6.92e−6 2.18e−6

7.11 7.34 7.49

This equation is computed by (41) of the order O(τ 2 + h 8 ) with τ = 0.1 and h = 0.5. We take two different sets of diffusion coefficients to study the effects for the fractional derivative terms. When (κ1 , κ2 ) = (0.5, 0.5), (48) is 1 ∂ 1.5 φ(x) ∂t φ(x, t) = √ . 2 ∂|x|1.5 The numerical results at t = 0, 10 are displayed in the left subplot of Fig. 5, showing a typical diffusion process.

(49)

52

Y. Ying et al. / Comput. Methods Appl. Mech. Engrg. 317 (2017) 42–54

Fig. 5. Numerical results for (48) at t = 0, 10. Left: (κ1 , κ2 ) = (0.5, 0.5); right: (κ1 , κ2 ) = (0.3, −0.3).

Fig. 6. Numerical results of (48) at t = 0, 10 for (κ1 , κ2 ) = (0.8, 0.2).

On the other hand, when (κ1 , κ2 ) = (0.3, −0.3), (despite lack of physical meaning for κ2 < 0,) the dispersion during wave propagation is clearly observed in the right subplot of Fig. 5. We further consider a generic case with (κ1 , κ2 ) = (0.8, 0.2) = (0.5 + 0.3, 0.5 − 0.3). Fig. 6 presents the profile of φ(x, t) at t = 0, 10. The anomalous diffusion can be understood as a combination of a symmetric diffusion and a dispersive convection. 5. Conclusion High order central difference schemes for Caputo fractional derivative are proposed. After a decomposition to Caputo fractional derivative, schemes of any desired accuracy order for the resulting components are obtained by matching the symbols of the operators and an overall scheme with weighted average of shifted differences. We use the schemes and the Crank–Nicolson algorithm to solve FADE. Unconditional stability of the scheme at order O(τ 2 +h p ) are shown. We prove the cases of p = 2, 4, 6, 8, +∞ as examples. Following the same way, cases of higher p may be confirmed. It is commonly conceived that anomalous diffusion contains components of diffusion and inherent convection. The decomposition of Caputo fractional derivative offers an alternative interpretation. The diffusion is symmetric and governed by the symmetric component (equivalent to Riesz fractional derivative), whereas the convection is dispersive which deforms the solution in an antisymmetrical manner. Though the accuracy order of the schemes can be arbitrarily high, how to choose an optimal one for a practical problem is still an open question. Besides, it should be noted that φ(x, t) in (1) is extended to be 0 outside of [0, L] to exploit the proposed schemes. Sufficient smoothness of φ(x, t) at x = 0, L is required to give a good approximation. However, this condition may not be satisfied for the cases with inhomogeneous boundary conditions. Then, there are much further work to be done for applications to these boundary conditions, such as Dirichlet boundary condition, Neumann boundary condition, and Robin boundary condition. As is well-known in the nonlinear parabolic equation community, singularities, blow-ups and boundary singular effects can combine to give a complex scenario. Further explorations are challenging and desirable.

Y. Ying et al. / Comput. Methods Appl. Mech. Engrg. 317 (2017) 42–54

53

Acknowledgements This research is partially supported by NSFC under Grant Nos. 11272009 and 11521202. Yuping Ying is grateful to the support from Chinese Scholarship Council. Yanping Lian and Wing Kam Liu are grateful to the support from ARO under Grant No. W911NF-15-1-0569. Appendix In this appendix, we prove F+ (z) is real and non-positive for p = 2, 4, 6, 8, ∞ when 1 < α < 2. Proof. When p = +∞, it is obvious that απ α |z| ≤ 0 − π ≤ z ≤ π. 2 When p = 2, it can be readily checked that z α απ  (2) F+ (z) = F+ (z) = cos 2 sin  ≤ 0 − π ≤ z ≤ π. 2 2 For p > 2, we have p  2 −1  ( p) (2) F+ (z) =  b j cos( j z) F+ (z). (∞)

F+ (z) = F+ (z) = cos

(50)

(51)

(52)

j=0

For p = 4, it holds that p  2 −1  1 1 ( p)  b j cos( j z) = 1 + α − α cos z ≥ 1. 12 12 j=0

(53)

For p = 6, it holds that  p     2 −1  17 1 2 41 1 2 11 1 2 ( p)  b j cos( j z) = 1 + α+ α + − α− α cos z + α+ α cos(2z) 160 192 360 144 1440 576 j=0     17 41 11 1 1 1 ≥ 1+ − − − − α+ α2 160 360 1440 192 144 576 > 0. For p = 8, it holds that  p   2 −1  5297 29 2 5 7843 3 2 5 ( p) 3 3   b j cos( j z) = 1 + α+ α + α + − α− α − α cos z 45360 3456 20736 60480 256 13824 j=0   211 7 2 1 3 + α+ α + α cos(2z) 15120 1920 6912   191 11 2 1 + − α− α − α 3 cos(3z) 181440 34560 41472   5297 7843 211 191 ≥ 1+ − − − α 45360 60480 15120 181440   29 3 7 11 + − − − α2 3456 256 1920 34560   5 5 1 1 + − − − α3 20736 13824 6912 41472

(54)

54

Y. Ying et al. / Comput. Methods Appl. Mech. Engrg. 317 (2017) 42–54

> 0.

(55)

We conclude that F+ (z) ≤ 0 for p = 2, 4, 6, 8. References [1] Y.A. Rossikhin, M.V. Shitikova, Application of fractional calculus for dynamic problems of solid mechanics: novel trends and recent results, Appl. Mech. Rev. 63 (1) (2010) 010801. [2] B.J. West, Colloquium: Fractional calculus view of complexity: A tutorial, Rev. Modern Phys. 86 (4) (2014) 1169. [3] R. Metzler, J. Klafter, The random walk’s guide to anomalous diffusion: a fractional dynamics approach, Phys. Rep. 339 (1) (2000) 1–77. [4] K. Oldham, J. Spanier, The Fractional Calculus: Theory and Applications of Differentiation and Integration to Arbitrary Order, Vol. 111, Academic Press, New York, London, 1974. [5] T.A.M. Langlands, B.I. Henry, The accuracy and stability of an implicit solution method for the fractional diffusion equation, J. Comput. Phys. 205 (2) (2005) 719–736. [6] Y. Lin, C. Xu, Finite difference/spectral approximations for the time-fractional diffusion equation, J. Comput. Phys. 225 (2) (2007) 1533–1552. [7] G.-h. Gao, Z.-z. Sun, A compact finite difference scheme for the fractional sub-diffusion equations, J. Comput. Phys. 230 (3) (2011) 586–595. [8] G.-H. Gao, Z.-Z. Sun, Hong-Wei Zhang, A new fractional numerical differentiation formula to approximate the Caputo fractional derivative and its applications, J. Comput. Phys. 259 (2014) 33–50. [9] A.A. Alikhanov, A new difference scheme for the time fractional diffusion equation, J. Comput. Phys. 280 (2015) 424–438. [10] J. Cao, C. Li, Y.Q. Chen, High-order approximation to Caputo derivatives and Caputo-type advection-diffusion equations (ii), Fract. Calc. Appl. Anal. 18 (3) (2015) 735–761. [11] H. Li, J. Cao, C. Li, High-order approximation to Caputo derivatives and Caputo-type advection–diffusion equations (III), J. Comput. Appl. Math. 299 (2016) 159–175. [12] I. Podlubny, Fractional Differential Equations, Academic press, New York, 1999. [13] C. Tadjeran, M.M. Meerschaert, H.-P. Scheffler, A second-order accurate numerical approximation for the fractional diffusion equation, J. Comput. Phys. 213 (1) (2006) 205–213. [14] W. Tian, H. Zhou, W. Deng, A class of second order difference approximations for solving space fractional diffusion equations, Math. Comp. 84 (294) (2015) 1703–1727. [15] Ch. Lubich, Discretized fractional calculus, SIAM J. Math. Anal. 17 (3) (1986) 704–719. [16] F. Zeng, C. Li, F. Liu, I. Turner, Numerical algorithms for time-fractional subdiffusion equation with second-order accuracy, SIAM J. Sci. Comput. 37 (1) (2015) A55–A78. [17] M. Chen, W. Deng, Fourth order accurate scheme for the space fractional diffusion equations, SIAM J. Numer. Anal. 52 (3) (2014) 1418–1438. [18] C. Li, F. Zeng, Numerical Methods for Fractional Calculus, Vol. 24, Chapman and Hall/CRC, Boca Raton, USA, 2015. [19] M.D. Ortigueira, Riesz potential operators and inverses via fractional centred derivatives, Int. J. Math. Math. Sci. 2006 (2006) 1–12. [20] C. C ¸ elik, M. Duman, Crank–Nicolson method for the fractional diffusion equation with the Riesz fractional derivative, J. Comput. Phys. 231 (4) (2012) 1743–1750. [21] H. Ding, C. Li, Y. Chen, High-order algorithms for Riesz derivative and their applications (II), J. Comput. Phys. 293 (2015) 218–237. [22] S. Li, W.K. Liu, Moving least-square reproducing kernel method Part II: Fourier analysis, Comput. Methods Appl. Mech. Engrg. 139 (1) (1996) 159–193.