May 9, 2012 - the most commonly used fractional derivatives are the GrünwaldâLetnikov derivative, the Riemannâ. Liouville derivative and the Caputo ...
May 9, 2012
10:43
WSPC/S0218-1274
1230014
International Journal of Bifurcation and Chaos, Vol. 22, No. 4 (2012) 1230014 (28 pages) c World Scientific Publishing Company DOI: 10.1142/S0218127412300145
FINITE DIFFERENCE METHODS FOR FRACTIONAL DIFFERENTIAL EQUATIONS* CHANGPIN LI† and FANHAI ZENG Department of Mathematics, Shanghai University, Shanghai 200444, P. R. China Received March 15, 2011; Revised May 6, 2011 In this review paper, the finite difference methods (FDMs) for the fractional differential equations are displayed. The considered equations mainly include the fractional kinetic equations of diffusion or dispersion with time, space and time-space derivatives. In some way, these numerical methods have similar form as the case for classical equations, some of which can be seen as the generalizations of the FDMs for the typical differential equations. And the classical tools, such as the von Neumann analysis method, the energy method and the Fourier method are extended to numerical methods for fractional differential equations accordingly. At the same time, the techniques for improving the accuracy and reducing the computation and storage are also introduced. Keywords: Fractional differential equations; finite difference method; extrapolation method; alternating direction implicit method.
1. Introduction Fractional derivatives, which are as old as the classical one, did not attract much attention until the recent decades. Nowadays, fractional differential equations (FDEs) have been found useful and applicable in science and engineering, such as physics [Barkai et al., 2000; Zaslavsky, 2002; Klages et al., 2008; Magin et al., 2008], materials [Diethelm & Freed, 1999b], control theory [Podlubny, 1999], biology [Magin, 2006] and finance [Raberto et al., 2002], and so on. In physics, fractional derivatives are adopted to model anomalous diffusion [Metzler & Klafter, 2000, 2004; Li & Zhao, 2011a], where particles spread differently than predicted by the classical Brownian motion model. Fractional kinetic equations have been proved to be particularly useful in the context of anomalous slow diffusion (subdiffusion)
[Metzler & Klafter, 2000]. Analytical methods, such as the Fourier transform method, the Laplace transform method, the Mellin transform method and the Green function method, etc. have been developed to seek the closed-form analytical solutions of the FDEs [Podlubny, 1999], while there are very few simple cases (such as some linear FDEs) in which the analytical solutions are available [Podlubny, 1999; Shen et al., 2008a]. Therefore, developing efficient and reliable numerical methods for FDEs is of great importance. Numerical methods for the classical differential equations (DEs) are abundant, while numerical methods for the FDEs are very limited. Loosely but not mathematically speaking, the fractional derivatives can be seen as a generalization of the classical derivative, so some numerical methods for the classical DEs can be extended to the FDEs in some way.
∗
The present work was partially supported by the National Natural Science Foundation of China under Grant no. 10872119 and the Key Disciplines of Shanghai Municipality under Grant no. S30104. † Author for correspondence 1230014-1
May 9, 2012
10:43
WSPC/S0218-1274
1230014
C. Li & F. Zeng
However, the differences exist in that the fractional derivatives are globally defined compared with the classical one defined in a pointwise manner, which attributes to the difficulties in the construction of numerical methods and analysis (including consistency, convergence and stability) for the FDEs, especially for the nonlinear ones. At the same time, the computation cost and storage are much expensive, since the coefficient matrices derived from the fractional systems can be full. At present, there are several numerical methods to solve FDEs, such as the finite difference method (rigorously speaking, the fractional finite difference method, here we just use it for brevity), the finite element method and the spectral method which are relatively rare [Roop, 2006, 2008; Zheng et al., 2010a, 2010b; Lin & Xu, 2007; Li & Xu, 2009] in literature, while the work on the finite difference method for FDEs are very rich and continues to be developed. This is mainly due to the following reasons. The finite difference method is a powerful tool and widely used to solve the DEs as well as the FDEs in science and engineering, which is also easy to be understood. Meanwhile, the implementation of the finite difference scheme is simple and easy to be put into practice in computer programs. In this paper, we just review almost all the existing FDMs for the FDEs of different types. If some important references have been omitted, we do apologize for those omissions. In recent publications, FDMs for the FDEs almost concentrate on the kinetic equations of diffusion (subdiffusion and superdiffusion [Li & Zhao, 2011a]), diffusion– advection, and Fokker–planck type, etc. with partial fractional derivatives that are derived asymptotically from basic random walk models, the generalized master and Langevin equations. We will outline those numerical schemes for some typical FDEs and their main results (including consistency, stability and convergence). Simultaneously, the unresolved problems are also pointed out. To some extent, some of FDMs for FDEs can be seen as the generalizations of the corresponding methods for classical equations [Yuste, 2006; Meerschaert & Tadjeran, 2006]. The tools, say von Neumann stability analysis [Yuste & Acedo, 2005; Ghazizadeh et al., 2010], Fourier analysis method [Chen et al., 2007; Sousa, 2009] and energy method [Zhuang et al., 2008; Gao & Sun, 2011], which are used for the numerical analysis of the classical PDEs, can be extended to the cases for
FPDEs. Because of the nonlocality and weak singularity of the fractional operators, the computational cost and storage of numerical methods for FDEs are more expensive compared with the corresponding methods for classical equations. In order to overcome these difficulties, some techniques were adopted to reduce the computational cost and storage of the derived methods, such as short memory principal [Diethelm & Freed, 1999a; Deng, 2007b], Richard extrapolation [Diethelm & Walz, 1997; Yuste, 2006; Tadjeran & Meerschaert, 2007] to get high order methods and ADI methods for twoand three-dimensional problems that convert highdimensional problems to separate one-dimensional ones [Meerschaert et al., 2006; Tadjeran & Meerschaert, 2007; Liu et al., 2008], and [Ford & Simpson, 2001; Wang et al., 2010]. The present paper is organized as follows. In Sec. 2, we introduce several definitions for different types of fractional derivatives and the notations used in the numerical methods for FDEs. In Sec. 3, some techniques are introduced to approximate the Riemann–Liouville and Caputo fractional derivatives. The existing FDMs for the fractional ordinary equations (FODEs) are surveyed in Sec. 4. The case with the fractional partial differential equations (FPDEs) are displayed in Sec. 5. In Sec. 6, the finite difference methods for some variable-order FPDEs are introduced. Conclusions and remarks are included in the last section.
2. Preliminaries In this section, we will introduce the definitions of fractional derivatives and notations used to describe the numerical schemes. There are several different ways to define the fractional derivatives, and the most commonly used fractional derivatives are the Gr¨ unwald–Letnikov derivative, the Riemann– Liouville derivative and the Caputo derivative [Samko et al., 1993; Podlubny, 1999; Li & Deng, 2007; Li et al., 2009; Li et al., 2011a]. We just introduce the definitions as follows. Definition 2.1. The fractional integral (or the Riemann–Liouville integral) with order α > 0 of the given function f (t) is defined as t 1 −α (t − s)α−1 f (s)ds, (1) Da,t f (t) = Γ(α) a
where Γ(·) is the Euler’s gamma function.
1230014-2
May 9, 2012
10:43
WSPC/S0218-1274
1230014
Finite Difference Methods for Fractional Differential Equations
Definition 2.2.
The left and right Gr¨ unwald– Letnikov derivatives with order α > 0 of the given function f (t) are defined as N α −α j α lim h (−1) f (t − jh), GL Da,t f (t) = h→0 j j=0 N h=t−a
(2)
Another type of the fractional derivative is the Riesz derivative [Samko et al., 1993]. The Riesz derivative has several kinds of forms [Podlubny et al., 2009], and we just introduce one case of definitions that will be used in the FDEs [Podlubny et al., 2009; Zhuang et al., 2009; Yang et al., 2009, 2010a, 2010b]. Definition 2.5. The Riesz derivative with order
and α GL Dt,b f (t)
=
N −α lim h (−1)j
h→0 N h=b−t
j=0
α > 0 of the given function f (t) is defined as
α j
(3) respectively. Definition 2.3. The left and right Riemann– Liouville derivatives with order α > 0 of the given function f (t) are defined as α RL Da,t f (t)
dn −(n−α) [a Dt f (t)] dtn dn t 1 (t − s)n−α−1 f (s)ds, = Γ(n − α) dtn a (4) =
and α RL Dt,b f (t) =
(−1)n dn Γ(n − α) dtn b × (s − t)n−α−1 f (s)ds,
−(n−α) (n) [f (t)] a Dt
1 = Γ(n − α)
t a
(t − s)n−α−1 f (n) (s)ds, (6)
and α C Dt,b f (t)
(−1)n = Γ(n − α)
t
b
Remark 2.1. Generally speaking, the above definitions of fractional derivatives are not equivalent, the differences and relations are discussed in detail in [Samko et al., 1993; Podlubny, 1999; Kilbas et al., 2006; Li & Deng, 2007; Li et al., 2009; Li et al., 2011a], and we just list one case as follows α = C Da,t f (t)
(5)
Definition 2.4. The left and right Caputo derivatives with order α > 0 of the given function f (t) are defined as
(s − t)n−α−1 f (n) (s)ds, (7)
respectively, where n is a non-negative integer and n − 1 < α < n.
(8)
In the above Definition 2.5, the parameter cα is chosen to be − 2 cos(1 πα ) (the default value in 2 this paper) in [Zhuang et al., 2009; Yang et al., 2009a, 2010a, 2010b], and 1/2 in [Podlubny et al., 2009] for the symmetric Riesz derivative, detailed information refers to [Samko et al., 1993; Podlubny et al., 2009; Zhuang et al., 2009] and the references therein.
+
respectively, where n is a non-negative integer and n − 1 ≤ α < n.
=
α α = cα (RL Da,t f (t) + RL Dt,b f (t)).
α RL Da,t f (t)
t
α C Da,t f (t)
α RZ Dt f (t)
f (t + jh),
n−1 k=0
f (k) (a)(t − a)k−α , Γ(k + 1 − α)
(9)
C n−1 [a, t]
where f ∈ and f (n) is integrable on [a, t]. α f (t) = Furthermore, if f ∈ C n [a, t], then GL Da,t α RL Da,t f (t). As to the fractional integrability and differentiability of a considered function, the reader can refer to [Li & Zhao, 2011b]. In order to better illustrate the numerical methods, some notations are introduced as follows. For single-variable function u(t) defined on [t0 , T ], define ∆t = (T − t0 )/nT to be the uniform time step, nT is a positive integer. The temporal grid points tk is defined by tk = t0 + k∆t, k = 0, . . . , nT . We denote uk to be the approximation of u(tk ). For two-variable function u(x, t) defined on [xL , xR ] ⊗ [0, T ], denote by ∆x the grid size in x-direction, ∆x = (xR − xL )/Nx and xi = xL + i∆x for i = 0, . . . , Nx , Nx is a positive integer. ∆t and tk are defined as above. Define uki as the approximation of u(xi , tk ).
1230014-3
May 9, 2012
10:43
WSPC/S0218-1274
1230014
C. Li & F. Zeng
For three-variable function u(x, y, t) defined on [xL , xR ] ⊗ [yL , yR ] ⊗ [0, T ], denote by ∆y the grid size in y-direction, ∆y = (yR − yL )/Ny and yj = yL + j∆y for j = 0, . . . , Ny , Ny is a positive integer. ∆t, ∆x, tk and xi are defined as above. Define uki,j as the approximation of u(xi , yj , tk ). For four-variable function u(x, y, z, t) defined on [xL , xR ] ⊗ [yL , yR ] ⊗ [zL , zR ] ⊗ [0, T ], denote by ∆z the grid size in z-direction, ∆z = (zR − zL )/Nz and zl = zL + l∆z for l = 0, . . . , Nz , Nz is a positive integer. ∆t, ∆x, ∆y, tk , xi and yj are defined as above. Define uki,j,l as the approximation of u(xi , yj , zl , tk ). In order to better illustrate the numerical schemes, we also introduce the following difference operators. δt uki =
uk+1 − uki i , ∆t
(10)
k+ 21
=
uk+1 + uki i , 2
(11)
δx uki =
uki+1 − uki , ∆x
(12)
ui
δxˆ uki = δx2 uki =
uki+1
uki−1
− 2∆x
,
(13)
uki+1 − 2uki + uki−1 . ∆x2
(14)
Then we define δx2 and δy2 as the central difference operators in x-direction and y-direction on uki,j respectively if it does not cause confusion for δx2 defined in (14). δx2 uki,j =
uki+1,j − 2uki,j + uki−1,j , ∆x2
(15)
δy2 uki,j =
uki,j+1 − 2uki,j + uki,j−1 . ∆y 2
(16)
The following difference operators are also introduced as follows δt uki,j,l
=
δxˆ uki,j,l = δyˆuki,j,l =
k uk+1 i,j,l − ui,j,l
∆t
,
uki+1,j,l − uki−1,j,l 2∆x uki,j+1,l − uki,j−1,l 2∆y
(17) , ,
δzˆuki,j,l δx2 uki,j,l
=
uki,j,l+1 − uki,j,l−1
=
uki+1,j,l − 2uki,j,l + uki−1,j,l
,
(21)
uki,j+1,l − 2uki,j,l + uki,j−1,l
,
(22)
uki,j,l+1 − 2uki,j,l + uki,j,l−1
.
(23)
δy2 uki,j,l = δz2 uki,j,l
=
2∆z
,
∆x2
∆y2
∆z 2
(20)
3. Approximation to the Fractional Derivatives In the present section, some numerical methods are introduced to approximate the fractional derivatives of the function u(t), t ∈ [t0 , T ]. Denote by α α RL DL,t u(tk ) and C DL,t u(tk ), the approximate soluα α u(t ) respectively, tion of RL DL,t u(tk ) and C DL,t k where ∆t, tj and nT are defined as the preceding section.
3.1. Approximation to the Riemann–Liouville derivative For a wide class of functions, both the Gr¨ unwald– Letnikov derivative and the Riemann–Liouville derivative are equivalent, especially for applications. Therefore, the Riemann–Liouville definition is suitable for the problem formulation, while the Gr¨ unwald–Letnikov definition is utilized to obtain the numerical solution [Podlubny, 1999]. In this section, some techniques and methods for approximating the Riemann–Liouville derivative are introduced and analyzed.
3.1.1. Gr¨ unwald–Letnikov approximation By the definition of the Gr¨ unwald–Letnikov derivative, it is natural to use (2) and (3) to approximate the left and right Riemann–Liouville derivatives. Denote by ωjα = (−1)j ( αj ), then the left and right Riemann–Liouville derivative operators RL Dtα0 ,t and α RL Dt,T can be approximated respectively by
(18) (19) 1230014-4
α RL Dt0 ,t u(tk )
≈ RL Dtα0 ,t u(tk ) = ∆t
−α
k j=0
(α)
ωj u(tk−j ),
(24)
May 9, 2012
10:43
WSPC/S0218-1274
1230014
Finite Difference Methods for Fractional Differential Equations (α)
w1 (z) = (1 − z)α
and α RL Dt,T u(tk )
α u(t ) ≈ RL Dt,T k −α
= ∆t
nT −k
j=0
(α) ωj u(tk+j ).
∞ = (−1)j j=0
(25) =
∞ j=0
The approximation (24) is convergent with order one for any α > 0 [Podlubny, 1999]. Both (24) and (25) are also called the standard Gr¨ unwald– Letnikov formulae, which may contribute to the unstable numerical schemes in solving FDEs [Meerschaert & Tadjeran, 2004] for 1 < α < 2, while the shifted Gr¨ unwald–Letnikov formula may be more useful for constructing the stable numerical schemes. The right and left shifted Gr¨ unwald– Letnikov formulae (one shift), which are used to approximate the left and right Riemann–Liouville derivatives, are defined respectively by α RL Dt0 ,t u(tk )
≈ RL Dtα0 ,t u(tk ) = ∆t−α
k+1 j=0
(α)
ωj u(tk−j+1 ),
(26)
−α
= ∆t
nT −k+1
j=0
(α)
ωj u(tk+j−1 ).
(α) w2 (z) (α) w3 (z)
(27)
Remark 3.1. The approximations (26) and (27) also have one-order convergence, which are useful to construct the stable numerical schemes for the FPDEs. In fact, any shift (p shifts) gives the first-order consistency, the best performance comes from minimizing |p − α| [Oldham & Spanier, 1974; Meerschaert & Tadjeran, 2004]. If 1 < α ≤ 2, the optimal choice is p = 1. The case of α = 2 reduces to the centered second difference estimator of the second derivative. In the following sections, the shifted formulae used to approximate the left and right Riemann–Liouville derivatives refer to (26) and (27) respectively.
In (24), {ωjα }kj=0 are just the first k + 1 coefficients of Taylor series of the expansion of the following function [Lubich, 1986]
zj
(α)
ωj z j .
(28)
= = =
1 3 − 2z + z 2 2 2
α ,
3 11 1 − 3z + z 2 − z 3 6 2 3
α ,
25 4 1 − 4z + 3z 2 − z 3 + z 4 12 3 4
α ,
137 10 (29) − 5z + 5z 2 − z 3 60 3 5 4 1 5 α + z − z , 4 5 15 147 20 15 (α) w6 (z) = − 6z + z 2 − z 3 + z 4 60 2 3 4 α 6 1 − z5 + z6 . 5 6
(α) w5 (z)
α u(t ) ≈ RL Dt,T k
k
Actually, higher order approximations of the form (24) exist, Lubich [1986] obtained approximations of order 2–6 in the form of (24), where the coefficients {ωjα } are just the coefficients of the Taylor series expansions of the following “generating” functions
(α) w4 (z)
and α RL Dt,T u(tk )
α
=
The coefficients {ωjα } can be easily and efficiently computed by the fast Fourier transform method [Lubich, 1986].
3.1.2. L1, L2 and L2C scheme In applications, the fractional derivative of order α(0 < α ≤ 2) is of great interest, we just list some methods to approximate the fractional derivatives in such cases. The efficient way to approximate the Riemann– Liouville derivative of order α(0 < α < 1) is the L1 scheme [Oldham & Spanier, 1974; Langlands & Henry, 2005], the derived scheme can be obtained by the following steps.
1230014-5
May 9, 2012
10:43
WSPC/S0218-1274
1230014
C. Li & F. Zeng
Using (9) leads to α RL Dt0 ,t u(tk+1 )
=
scheme (31) with t0 = 0 was derived as α 2−α α u(t u(tk+1 ) − RL D0,t , RL D0,t k+1 ) ≤ C∆t
u(t0 )(tk+1 − t0 )−α 1 + Γ(1 − α) Γ(1 − α) ×
k j=0
tj+1
tj
(32)
(tk+1 − s)−α u (s)ds. (30)
t If the integral tjj+1 (tk+1 − s)−α u (s)ds is approxit u(t )−u(t ) mated by tjj+1 j+1∆t j (tk+1 − s)−α ds, then the L1 scheme is derived
where C is dependent on α and u. The L2 scheme and its modification L2C scheme [Oldham & Spanier, 1974; Lynch et al., 2003] are suitable to discretize the Riemann– Liouville derivative with order α(1 < α < 2). Similar to the L1 scheme, one can get α RL DL,t u(tk+1 )
∆t−α u(t0 )(tk+1 − t0 )−α α + RL Dt0 ,t u(tk+1 ) = Γ(1 − α) Γ(2 − α) ×
k
=
u(t0 )(tk+1 − t0 )−α u (t0 )(tk+1 − t0 )1−α + Γ(1 − α) Γ(2 − α) k
bj,k [u(tj+1 ) − u(tj )],
1 + Γ(2 − α)
(31)
j=0
j=0
tj+1 tj
s1−α u (tk+1 − s)ds. (33)
(k −j +1)1−α −(k −j)1−α , j
= 0, . . . , k. where bj,k = In [Langlands & Henry, 2005; Sun & Wu, 2006; Lin & Xu, 2007], the error estimate of the above L1
On each subinterval [tj , tj+1 ], u (tk+1 −s) is approxu(tk−j−1 )−2u(tk−j )+u(tk+1−j ) , then the L2 imated by ∆t2 scheme is obtained as
u(t0 )(tk+1 − t0 )−α u (t0 )(tk+1 − t0 )1−α α u(t + D ) = RL t0 ,t k+1 Γ(1 − α) Γ(2 − α) k
+
∆t−α bj [u(tk−j−1 ) − 2u(tk−j ) + u(tk+1−j )] Γ(3 − α) j=0
=
u(t0 )(tk+1 − t0 )−α u (t0 )(tk+1 − t0 )1−α + Γ(1 − α) Γ(2 − α) +
k+1
Wj u(tk+1−j ),
(34)
j=−1 k where bj = (j + 1)2−α − j 2−α , and {Wj }k+1 j=−1 can be expressed by {bj }j=0 .
By letting u (tk+1 − s) ≈ scheme is obtained as α RL Dt0 ,t u(tk+1 )
=
u(tk−j−1 )−u(tk−j )+u(tk+2−j )−u(tk+1−j ) 2∆t2
on each subinterval [tj , tj+1 ], the L2C
u(t0 )(tk+1 − t0 )−α u (t0 )(tk+1 − t0 )1−α + Γ(1 − α) Γ(2 − α) k
+
∆t−α bj [u(tk−j−1 ) − u(tk−j ) + u(tk+2−j ) − u(tk+1−j )] 2Γ(3 − α) j=0
k+2 u(t0 )(tk+1 − t0 )−α u (t0 )(tk+1 − t0 )1−α ˆ j u(tk+1−j ), + + W = Γ(1 − α) Γ(2 − α) j=−1
1230014-6
(35)
May 9, 2012
10:43
WSPC/S0218-1274
1230014
Finite Difference Methods for Fractional Differential Equations
ˆ j }j=k+2 can be expressed by {bj }k . where {W j=0 j=−1 The accuracy of the L2 and L2C schemes depend on α, see [Lynch et al., 2003] for more detailed information. Remark 3.2. The above L1, L2 and L2C schemes can be extended to the approximation of the right Riemann–Liouville derivative accordingly, and we also call corresponding methods for the right Riemann–Liouville derivative L1, L2 and L2C schemes respectively.
All the above methods can be extended to approximate the Riesz fractional derivative correspondingly. For instance, the standard Riesz formula means that the left and right Riemann–Liouville derivatives contained in the Riesz derivative are approximated by the standard Gr¨ unwald–Letnikov formulae (24) and (25) respectively. The shifted Riesz formula means that the left and right Riemann–Liouville derivatives in the Riesz derivative are approximated by the right and left shifted Gr¨ unwald–Letnikov formulae (26) and (27) respectively. So do L1, L2 and L2C methods.
3.2. Approximation to the Caputo derivative Since the Riemann–Liouville derivative and the Caputo derivative have the relationship as (9), the L1, L2 and L2C method can be extended to approximate the Caputo derivative directly, we still call them the L1, L2 and L2C method if no confusion occurs. If 0 < α < 1, the Caputo derivative can be discretized by the following L1 method α C Dt0 ,t u(tk+1 )
Remark 3.3.
1 Γ(2 − α)∆t u(t
j=0
(36)
where bj = (j +1)1−α −j 1−α , and the error estimate is the same as (32). For 1 < α < 2, the L2 scheme and the L2C scheme to approximate the Caputo derivative can be written as α C Dt0 ,t u(tk+1 ) =
k+1
Wj u(tk+1−j ),
(37)
ˆ j u(tk+1−j ), W
(38)
j=−1
and α C Dt0 ,t u(tk+1 ) =
mula (24) reduces to the back forward difference method, L2 method (34) reduces to the forward difference method and L2C method (35) reduces to the centered difference method. The shifted formula (26) and L2 method (34) reduce to the centered difference method for α = 2.
=
∆t−α bj Γ(2 − α) × [u(tk+1−j ) − u(tk−j )],
Remark 3.4. If α = 1, the Gr¨ unwald–Letnikov for-
α C D0,t u(tk+1 )
k
=
k+2 j=−1
ˆ j are defined in (34) respectively, where Wj and W and (35). If 1 < α < 2, and u (0) is known, then α C D0,t u(tk+1 ) can also be approximated by the following method [Sun & Wu, 2006] k a0 δt uk+1/2 − (ak−j − ak−j−1 )δt uj−1/2 − ak u (0), (39) j=1 )−u(t )
j+1 j and where uj = u(tj ), δt uj+1/2 = ∆t 2−α ∆t 2−α 2−α −j ]. The truncation error aj = 2−α [(j + 1) α 3−α , which α u(t D0,t is |C k+1 ) − C D0,t u(tk+1 )| ≤ C∆t was presented in [Sun & Wu, 2006]. In [Odibat, 2009], a computational algorithm for approximating the Caputo derivative was also developed, and the convergence order is O(∆t2 ) for all 0 < α < 2. Another difference method of order two was derived in [Sousa, 2010] for 1 < α ≤ 2. In [Diethelm et al., 2005], numerical algorithms for solving the fractional-order integral, the Caputo derivative and the differential equations with the
Caputo derivative were provided and discussed in detail. Other methods that are used to approximate the Caputo derivative, refer to [Odibat, 2006; Murio, 2006; Schmidt & Gaul, 2006].
3.3. Matrix approach This method is based on triangular strip matrix approach to discretize the differentiation and integration operators with arbitrary order [Podlubny, 2000; Podlubny et al., 2009]. This technique can obtain all the numerical solutions at the mesh
1230014-7
May 9, 2012
10:43
WSPC/S0218-1274
1230014
C. Li & F. Zeng
grids at once, avoiding the traditional step-by-step method by moving from the previous time layer to the next one to get the numerical solutions. The matrix approach is quite simple to put into implementation. In order to develop the scheme, the matrices with the special structure such as the triangular strip type are introduced. For simplicity, we just introduce the following two kinds of triangular strip matrices: Lower triangular strip ω0 0 ω ω0 1 ω2 ω1 . . LN = . . .. .. ωN −1 . ωN −1 ωN
matrix 0 0 ω0 .. .
0 0 0 .. .
··· ··· ···
ω2 .. .
ω1
ω0
ω2
ω1
···
and upper triangular strip matrix ω0 ω1 ω2 . . . ωN −1 .. 0 ω0 ω1 . . . . 0 0 ω0 ω1 ... UN = . . 0 .. .. 0 ... . . . . . . . . . . . . ω0 0 0 0 ... 0
0 0 0
, · · · 0 ω0
ωN
(40)
α C D0,t u(t)
(j)
u(j) (0) = u0 , (41)
ω0
in (40), then approximated by (α) (α) (α) (α) T uN , uN −1 , . . . , u1 , u0 k
(α)
= BN [u(tN ), . . . , u(t1 ), u(t0 )]T , where (α) BN (α) BN
(α) BN
=
LN ∆tα ,
t ∈ [0, T ],
(43)
j = 0, 1, . . . , n − 1,
(44)
= f (t, u(t)),
with the initial condition
ωN −1 ωN −2 . ... ω1
(α) (α) (α) [uN , . . . , u1 , u0 ]T
4. Numerical Method for FODEs The fractional difference method is a useful tool to solve the FODEs. In this section, we review the difference methods to solve the following FODE
Consider the approximation of Riemann–Liouville derivative (4) (or Caputo derivative (6)). Denote (α) α u(t ) (or α by uk = RL D0,t C D0,t u(tk )), let ωk = k
(−1)k ( α )
while the fractional difference matrix is different from the Riemann–Liouville or Caputo case, refer to [Podlubny et al., 2009] for details. In [Podlubny, 2000; Podlubny et al., 2009], this approach is adopted to solve the classical differential equations such as classical diffusion equations, and the FDEs such as diffusion equations with time-fractional derivative, diffusion equations with spatial derivative in the Riesz sense, general diffusion equations with time-space-fractional derivative and fractional diffusion equations with delay. This method can also be extended to the cases of nonlinear problems, see [Podlubny, 2000; Podlubny et al., 2009] for more details.
can be
where α > 0, and n = α is the first integer not less than α. The existence and uniqueness of the initial value problem (43) can be found in [Diethelm et al., 2002] if the function f is continuous and satisfies the following Lipschitz condition |f (t, u1 ) − f (t, u2 )| ≤ L|u1 − u2 |,
where L is Lipschitz constant. The initial value problem (43) is equivalent to the following Volterra equation t n−1 tj (j) 1 u + (t − s)α−1 u(t) = j! 0 Γ(α) 0 j=0
(42)
∆t is the mesh size. We call
a kind of fractional difference matrix. If α = 1, is just the difference matrix for the classical first-order derivative corresponding to the firstorder backward difference. The right Riemann–Liouville (or Caputo) derivative can be approximated in a similar way, where the upper triangular strip matrix (41) is used accordingly. The symmetric Riesz fractional derivative can also be approximated by this approach,
(45)
× f (s, u(s))ds.
(46)
There are a number of schemes for the numerical solution of the initial value problem (43)– (44) (or (46)) in the literature, such as a class of the fractional multistep method [Lubich, 1986; Galeone & Garrappa, 2009], the fractional Adams method [Diethelm & Freed 1999a, 1999b; Diethelm et al., 2002, 2004; Li & Tao, 2009], and [Odibat & Momani, 2008b]. In [Diethelm et al., 2006], the existing methods for (43) were simply introduced, and the strengths and weaknesses of the
1230014-8
May 9, 2012
10:43
WSPC/S0218-1274
1230014
Finite Difference Methods for Fractional Differential Equations
established methods were investigated. In [Baker, 2000], some numerical methods for Volterra integral and integro-differential equations were also displayed, in which the numerical analysis of convergence and stability were discussed. The numerical solution of (46) can be obtained by the following discrete convolution quadrature u(tk ) =
n−1 tj k
j!
j=0
+ ∆tα
(j)
u0 + ∆tα s
k
ωk−j f (tj , u(tj ))
j=0
wkj f (tj , u(tj )),
(47)
j=0
where the convolution quadrature weights ωk−j are selected in order to fulfill some order conditions, and the starting weights wkj are chosen such that (47) is exact for some selected nonpolynomial terms. In [Baker, 2000; Diethelm et al., 2006], the computation of the convolution quadrature weights ωk−j , the starting weights wkj and some other related problems were discussed in detail. For more information, refer to [Baker, 2000; Diethelm et al., 2006] and the references therein. Here, we especially mention the Adams– Bashforth–Moulton method introduced in [Diethelm & Freed, 1999a, 1999b], and the derived scheme reads as
n−1 k tjk+1 (j) 1 P u + bj,k+1 f (tj , uj ), uk+1 = j! 0 Γ(α) j=0 j=0 j n−1 k t 1 k+1 (j) aj,k+1 f (tj , uj ) + ak+1,k+1 f (tk+1 , uPk+1 ), u0 + uk+1 = j! Γ(α) j=0 j=0 where aj,k+1
kα+1 − (k − α)(k + 1)α , α ∆t (k − j + 2)α+1 + (k − j)α+1 − 2(k − j + 1)α+1 , = α(α + 1) 1,
(48)
j = 0, 1 ≤ j ≤ k,
(49)
j = k + 1,
and bj,k+1
then we get the fractional Euler method
∆tα = [(k − j + 1)α − (k − j)α ], α 0 ≤ j ≤ k,
(50)
uj is the approximate solution of u(tj ). The scheme (48) can also be obtained by the following way. Considering the following integration tk+1 Ikα = (tk+1 − s)α−1 g(s)ds. (51)
uk+1 =
n−1 tjk+1 (j) u j! 0 j=0
0
≈
tk+1 0
(tk+1 − s)α−1 g˜k+1 (s)ds,
g˜k+1 (t)|[tj ,tj+1 ] =
tj+1 − t g(tj ) tj+1 − tj +
t − tj g(tj+1 ), tj+1 − tj
0 ≤ j ≤ k, (55)
(52)
then the fractional trapezoidal formula is derived n−1 k tj (j) 1 u0 + aj,k+1 f (tj , uj ) uk+1 = j! Γ(α) j=0
0 ≤ j ≤ k,
j=0
If g˜k+1 (s) is determined by
j=0
where g˜k+1 (s) is the approximation of g(s). If g˜k+1 (s) is chosen as g˜k+1 (s)|[tj ,tj+1 ) = g(tj ),
1 bj,k+1 f (tj , uj ). Γ(α) (54)
0
The above quadrature can be approximated by tk+1 (tk+1 − s)α−1 g(s)ds
k
+
+ ak+1,k+1 f (tk+1 , uk+1 ). (53) 1230014-9
(56)
May 9, 2012
10:43
WSPC/S0218-1274
1230014
C. Li & F. Zeng
The scheme (56) is implicit, which is not preferable for computing. Similar to the predictor– corrector method for the first-order ODE, we first use (54) to get uPk+1 (predictor), then we use (56) to get uk+1 (corrector) by replacing uk+1 with uPk+1 on the right-hand side of (56), which leads to (48). This method is very efficient and also suitable for longterm integration, such as the calculations of limit cycles and chaotic attractors and their synchronizations of fractional differential systems [Li & Peng, 2004; Li et al., 2006; Li & Deng, 2006; Zhou & Li, 2005]. In [Diethelm et al., 2002], the scheme (48) was investigated in a more detailed way and extended to the multiterm FODEs. The error analysis was carefully investigated in [Diethelm et al., 2004], and two conjectures (Conjectures 3.1 and 3.2 in that paper) were proposed, which respectively read as, u(T ) − uT /∆t =
k1
cj ∆t2j +
j=2
k1
dj ∆tj+α + O(hk3 ),
α > 0,
j=1
k ≥ 3,
α C D0,t u(t)
∈ C k [0, T ],
where k1 , k2 and k3 are certain constants depending only on k, k3 > max(2k1 , k2 + α); and max |u(tj ) − uj | = O(∆t1+α ),
tj ∈[,T ]
α > 1, f ∈ C 3 (G). Li and Tao [2009] further explored the error analysis of the scheme (48) and proved the following results (Theorems 3.1 and 3.4 in that article) in which they solved (negated) the above two conjectures, u(T ) − uT /∆t =
k j=1
and
∆t
j+α
2T /∆t+1
cj,i,T /∆t∆tiα ,
i=0
|u(tj ) − uj | = O(∆t1−α ), tjmax ∈[,T ]
> 0,
max |u(tj ) − uj | = O(1), tj ∈[0,T ]
under the same conditions. The detailed proofs can be referred to [Li & Tao, 2009]. Remark 4.1. In [Diethelm & Walz, 1997; Diethelm,
1997a], an algorithm based on the quadrature
α [u − formula approach was proposed for RL D0,t u0 ](t) = βu(t) + f (t), 0 < α < 1, t ∈ [0, 1], and the extrapolation technique was adopted to get high order method. Lin and Liu [2007] also discussed the consistence, convergence and stability (the derivative order α > 1) for the discrete convolution quadrature method for (43) with the fractional derivative defined in the Riemann–Liouville sense, and the numerical experiments were conducted to prove the effectiveness of their method.
Remark 4.2. The method (48) belongs to the type of the predictor–corrector method, uPk+1 is called the predictor and uk+1 the corrector. If α = 1, (48) reduces to the classical predictor–corrector method for the first-order ODEs. Choosing different predictors contributes to the different methods. In [Li et al., 2011b], gk+1 (s) is chosen as the piecewise quadratic interpolation to approximate (51), by which they constructed a higher-order predictor– corrector method to the numerical solution of (43). This idea was also adopted to solve FPDEs in [Deng, 2007a, 2007b], extended to solve the multiorder fractional differential equations describing the fractional order dynamical controlled systems in [Yang & Liu, 2006] and nonlinear dynamical systems in [Yin et al., 2007]. Remark 4.3. In [Baleanu et al., 2009], a modified Gr¨ unwald–Letnikov approach was proposed for a class of fractional optimal control problems, which was used to solve two problems of time invariant and time varying. In [Bhalekar et al., 2011], the product trapezoidal quadrature formula was extended to solve the fractional Bloch equation with delay and analyzed, which was also compared with the existing method for the same problem. Remark 4.4. Much error analysis was investigated for the scheme (48) and numerous examples were conducted to show that the Adams method was stable, while the theoretical analysis for stability has not been established for 0 < α < 1. In [Galeone & Garrappa, 2006, 2009; Garrappa, 2009], numerical methods of multistep type were proposed and stability results were derived only for the linear fracα u(t) = λu(t), 0 < α < 1, tional equation C D0,t λ ∈ C. Remark 4.5. Compared with the numerical methods for the classical differential equations, the numerical methods for the fractional differential equations
1230014-10
May 9, 2012
10:43
WSPC/S0218-1274
1230014
Finite Difference Methods for Fractional Differential Equations
need much more computational effort to get the numerical solutions because of the possible nonlocality of the fractional operators. Ford and Simpson [2001] explored the existing techniques for reducing the amount of computational effort required to solve the fractional differential equation (43) without significant loss of accuracy, detailed information refers to [Ford & Simpson, 2001]. Agrawal and Kumar [2007] compared stability, accuracy, and computational speeds of five different numerical schemes for (43) with homogeneous and nonhomogeneous initial conditions, and concluded that the choice of a numerical scheme would depend on the problem considered and the performance criteria selected. Other related works refer to [Diethelm, 1997a; Diethelm et al., 2005; Diethelm et al., 2006; Schmidt & Gaul, 2006].
5. The Difference Method for the FPDE In this section, we outline the FDMs for the FPDEs with time, space and time-space fractional derivatives. The FDMs for FPDEs in recent papers are mainly devoted to the equations of anomalous diffusion or dispersion, such as the fractional diffusion, advection–dispersion or reaction–diffusion equations, and the fractional Fokker–Planck type equations, etc. In the following, we will list the numerical schemes established by different authors for some typical equations, and find that some numerical schemes for FPDEs can be seen as the generalizations of the numerical schemes for classical PDEs. In some way, the numerical schemes for PDEs and FPDEs have the similar form, the difference depends on how the fractional derivatives of FPDEs are approximated compared with the classical one.
5.1. Finite difference methods for one-dimensional FPDEs This section concerns the existing finite difference methods for the one-dimensional time, space and time-space FPDEs.
5.1.1. Time fractional equations First, we consider the following type of onedimensional time-fractional diffusion equation (TFDE) of the form [Metzler & Klafter, 2000]
∂u ∂2u 1−γ = RL D0,t Kγ 2 + f (x, t), ∂t ∂x (x, t) ∈ (xL , xR ) ⊗ (0, T ],
(57)
with suitable initial and boundary conditions, 0 < γ ≤ 1 and Kγ > 0. Equation (57) can be rewritten into the following equivalent form [Metzler & Klafter, 2000] γ CD0,t u
= Kγ
∂2u + g(x, t), ∂x2
(x, t) ∈ (xL , xR ) ⊗ (0, T ].
(58)
There have already been several methods for the numerical solutions of (57) and (58), and it is the natural way to approximate the second spatial derivative in (57) and (58) with the second-order central difference method. In this paper, the secondorder spatial derivative is always discretized by the central difference method if we do not mention the discretization of the second-order spatial derivative in the numerical schemes. γ D0,t the approximate operator of Denote by RL γ γ γ k RL D , and RL D u is defined like RL D u(xi , t ) k 0,t γ γ k k with u(xi , tk ) replaced by ui , C D0,t and C D0,t ui are γ γ k defined similar to RL D0,t and RL D0,t ui , respectively. 0,t
0,t i
Now, we outline the numerical schemes for the numerical solution of the TFDE (57) and other anomalous equations such as the fractional Fokker– Plank equation [Zheng et al., 2010b], the fractional reaction–subdiffusion equation, and so on. Langlands and Henry [2005] proposed an implicit numerical scheme to approximate (57) with Kγ = 1 and f = 0 , which can be written as 1−γ 2 k+1 D0,t [δx ui ]. δt uki = RL
(59)
1−γ D0,t is the approximate operator defined Here RL by (31) (L1 scheme), δt and δx2 are backward Euler and central difference operators defined by (10) and (14), respectively. Langlands and Henry only numerically verified the accuracy and unconditional stability of the derived method [Langlands & Henry, 2005], while the global unconditional stability analysis for all γ(0 < γ < 1) has not been established theoretically. In fact, we can very clearly and succinctly show that the established scheme (59) is indeed unconditionally stable. We just illustrate the skeleton of
1230014-11
May 9, 2012
10:43
WSPC/S0218-1274
1230014
C. Li & F. Zeng
how to get it. Set ukj = δk eiqj∆x (i is the imaginary unit) and insert ukj into (59). By the simple computation, the following equation is derived (see Eq. (83) in [Langlands & Henry, 2005]) q∆x δk+1 = δk − 4ρ sin2 2 k γ δ0 + bk−l (δl+1 − δl ) , × (k + 1)1−γ
A class of fractional weighted average difference methods for (57) with f = 0 were established by Yuste [2006], which can be seen as the extension of the weighted average methods for ODEs. The established method is presented as follows
(60)
1−γ D0,t is the where λ is the weight factor, and RL 1−γ approximate operator of RL D0,t defined in (24) or other higher-order operators defined in the form of (24) obtained by Lubich [1986]. The accuracy of the fractional weighted average difference method (64) is given as O(∆t + ∆x2 ) except for the fractional Crank–Nicolson method (λ = 1/2), where the accuracy with respect to the timestep can be O(∆t2 ) if the secondorder approximation to the fractional derivative is used. The stability was investigated by the fractional von Neumann stability analysis which was employed in [Yuste & Acedo, 2005], and a very simple and accurate stability criterion for arbitrary γ and weighted factor λ was provided. The presented methods are unconditionally stable for 0 < λ ≤ 1/2 and conditionally stable for 1/2 < λ ≤ 1. Chen et al. [2007] proposed an implicit numerical scheme for (57). Their method can be derived D1−γ is chofrom (64) by letting λ = 0 and
1−γ δt uki = λ RL D0,t (Kγ δx2 uki ) 1−γ D0,t (Kγ δx2 uk+1 ), + (1 − λ) RL i
l=0
γ
∆t γ γ and q is a where ρ = Γ(1+γ)∆x 2 , bl = (l + 1) − l real spatial wave number. If |δk | ≤ |δ0 | for any k ≤ nT , we can assert that the scheme (59) is unconditionally stable. Denote by µ = 4ρ sin2 ( q∆x 2 ), then (60) can be written into the following equivalent form γ (1 + µ)δk+1 = δk + µ bk − δ0 (k + 1)1−γ
+µ
k (bk−l − bk−l+1 )δl .
(61)
l=1
It is easy to check that γ 1 + µ bk − 1−γ (k + 1) +µ
k
|bk−l − bk−l+1 | ≤ 1 + µ.
(62)
l=1
RL
Combining (62) and the mathematical induction method lead to |δk | ≤ |δ0 |. Therefore, the scheme (59) is unconditionally stable. In [Yuste & Acedo, 2005], the authors proposed an explicit method for (57) with f = 0 and Kγ = const. The first-order time derivative and the time fractional derivative were discretized by the forward Euler difference method and the standard Gr¨ unwald–Letnikov formula (24), respectively 1−γ 2 k D0,t (δx ui ), δt uki = Kγ RL
(63)
Yuste and Acedo investigated the stability of (63) by means of a powerful and simple new procedure close to the well-known von Neumann method for nonfractional partial differential equations, and got the truncating error as O(∆t + ∆x2 ). Meanwhile, they also checked the reliability and the stability bounds of the derived numerical algorithm.
(64)
0,t
sen to be the standard Gr¨ unwald–Letnikov derivative operator defined by (24). They successfully applied the Fourier analysis method to the convergence and stability analysis of the derived method, and obtained that the method was convergent with order O(∆t + ∆x2 ) and unconditionally stable. Cui [2009] also proposed a compact implicit difference scheme for (57) with the time fractional derivative discretized by the standard Gr¨ unwald–Letnikov formula (24), while the space derivative was approximated by a fourth-order difference method. The established implicit method was convergent with order O(∆t + ∆x4 ) and unconditionally stable. Remark 5.1. Clearly, one can recover the explicit method in [Yuste & Acedo, 2005] by letting λ = 1, and the scheme (59) is recovered by letting λ = 0 with the fractional derivative approximated by the
1230014-12
May 9, 2012
10:43
WSPC/S0218-1274
1230014
Finite Difference Methods for Fractional Differential Equations
L1 scheme. If γ = 1, the scheme (64) reduces to the classical average method. The numerical methods for Eq. (57) can be extended to other types of FPDEs. In [Chen et al., 2008b], two difference methods (explicit and implicit) were developed for the following reaction– subdiffusion equation ∂2u ∂u 1−γ Kγ 2 − Ku + f (x, t), = RL D0,t ∂t ∂x (x, t) ∈ (xL , xR ) ⊗ (0, T ],
(65)
with the initial and boundary conditions, 0 < γ ≤ 1 and Kγ , K > 0. The presented methods are similar to (64) with the fractional derivative approximated by the standard Gr¨ unwald–Letnikov formula (24), which can be written as 1−γ D0,t (Kγ δx2 uki − Kuki ) δt uki = λ RL 1−γ D0,t (Kγ δx2 uk+1 − Kuk+1 ) + (1 − λ) RL i i
+ λfik + (1 − λ)fik+1 .
(66)
If λ = 1, one gets the explicit method, and the implicit method is obtained for λ = 0. The Fourier method was taken to conduct the convergence and stability analysis of the two methods. They proved that the two methods were convergent with order O(∆t + ∆x2 ), and the implicit method was unconditionally stable, the explicit one was conditionally stable. They also used the Richardson extrapolation technique to improve the accuracy of the proposed algorithms. In [Chen et al., 2009a], the explicit and implicit numerical schemes that are similar to (66) were developed for Stokes’ first problem for a heated generalized second grade fluid with fractional derivative with a nonhomogeneous forcing term ∂2u ∂2u ∂u 1−γ = RL D0,t κ1 2 + κ2 2 + f (x, t), ∂t ∂x ∂x (x, t) ∈ (xL , xR ) ⊗ (0, T ], (67) with the initial and boundary conditions, 0 < γ ≤ 1 and κ1 , κ2 > 0. The technique for the stability and convergence analysis, the technique for improving the accuracy of the proposed methods and the results are all similar to [Chen et al., 2008b]. Chen et al. [2009b] proposed three implicit approximations for the following fractional
Fokker–Planck equation ∂u ∂ v (x) ∂2u 1−γ + Kγ 2 , = RL D0,t ∂t ∂x mηγ ∂x (x, t) ∈ (xL , xR ) ⊗ (0, T ],
(68)
with suitable initial and boundary conditions, 0 < γ ≤ 1 and Kγ > 0. They first transform (68) into an equivalent form, then three implicit methods, which were based on the Gr¨ unwald–Letnikov (time fractional derivative) and the back Euler difference method (first-order spacial derivative), L1 method and the back Euler difference method, L1 method and the central difference method, were developed for the derived equivalent equation. They proved that the three methods were unconditionally stable and convergent by the energy method, detailed information refers to [Chen et al., 2009b]. Next, we introduce a new method for the numerical solution of the FPDEs. Zhuang et al. [2008] proposed a new implicit numerical method for the following anomalous subdiffusion equation ∂2u ∂u 1−γ = RL D0,t Kγ 2 + f (x, t) , ∂t ∂x (x, t) ∈ (xL , xR ) ⊗ (0, T ],
(69)
with the suitable initial and boundary conditions, 0 < γ < 1 and Kγ > 0. Here, we just simply illustrate how the new method is obtained. First, integrating both sides of (69) on the interval [tk , tk+1 ] in the time direction, one gets u(xi , tk+1 ) = u(xi , tk ) + ×
tk+1
0
×
tk
0
1 Γ(γ)
1 Lu(xi , s) + f (xi , s) ds − 1−γ (tk+1 − s) Γ(γ)
Lu(xi , s) + f (xi , s) ds, (tk − s)1−γ
(70)
2
where Lu = Kγ ∂∂xu2 . Consider the following integration (see (52)) tk+1 (tk+1 − s)γ−1 g˜k+1 (s)ds, 0
if g˜k+1 is chosen to be as
1230014-13
g˜k+1 (s)|(tj ,tj+1 ] = g(tj+1 ),
0 ≤ j ≤ k,
(71)
May 9, 2012
10:43
WSPC/S0218-1274
1230014
C. Li & F. Zeng
then one gets
tk+1 0
Therefore,
tk+1 0
(tk+1 − s)γ−1 g˜k+1 (s)ds =
k ∆tγ bk−j g(tj+1 ), γ
bj = (j + 1)γ − j γ .
(72)
j=0
Lu(xi ,s)+f (xi ,s) (tk+1 −s)1−γ
k ∆tγ Lu(xi , s) + f (xi , s) ≈ bk−j (Lu(xi , tj+1 ) + f (xi , tj+1 )). (tk+1 − s)1−γ γ
(73)
k−1 ∆tγ Lu(xi , s) + f (xi , s) ≈ bk−j−1 (Lu(xi , tj+1 ) + f (xi , tj+1 )). (tk − s)1−γ γ
(74)
tk+1 0
Similarly, we have
j=0
tk 0
can be approximated by
j=0
Combining (70), (73) and (74), one derives the following implicit numerical scheme = uki + uk+1 i
= uki +
k
k−1
j=0
j=0
∆tγ ∆tγ j+1 bk−j (Kγ δx2 uj+1 + f ) − bk−j−1 (Kγ δx2 uj+1 + fij+1 ) i i i Γ(γ + 1) Γ(γ + 1) ∆tγ Γ(γ + 1)
b0 Kγ δx2 uk+1 + Kγ i
k−1 j=0
(bj+1 − bj )δx2 uk−j i
k−1 ∆tγ k+1 b0 fi + (bj+1 − bj )fik−j . + Γ(γ + 1)
(75)
j=0
The above scheme (75) is just the new implicit numerical method proposed by Zhuang et al. [2008] (see (2.8) in [Zhuang et al., 2008]). They proved that the derived new method was unconditionally stable and convergent with order O(∆t + ∆x2 ) by the energy method. They also improved the derived method by the extrapolation technique, at the same time, another method was proposed to get a highorder method of order O(∆t2 + ∆x2 ) by replacing g˜k+1 (s) in (71) with the piecewise interpolant (55), then repeating the procedures (72)–(75), the highorder method was obtained (see (5.2)–(5.5) in [Zhuang et al., 2008]). Wu [2009] also applied the above time discretization procedures [(70)–(75)] to solve the Stokes’ first problem (67) and got an implicit unconditionally stable method with convergence order O(∆t + ∆x2 ). Liu et al. [2009] used the same time discretization for the numerical solution of the modified anomalous subdiffusion equation of the form 2 ∂ u ∂u 1−γ1 1−γ1 = (A RL D0,t + B RL D0,t ) + f (x, t) ∂t ∂x2 + g(u, x, t),
(x, t) ∈ (xL , xR ) ⊗ (0, T ]
(76)
with suitable initial and boundary conditions, 0 < γ1 ≤ γ2 ≤ 1. In their method, the nonlinear term was treated explicitly, and they proved that the presented method was conditionally stable and convergent under suitable conditions. The time approximation of the following fractional cable equation [Liu et al., 2011] ∂2u ∂u 1−γ1 1−γ2 = RL D0,t u K 2 − µ2 RL D0,t ∂t ∂x + f (x, t),
0 < γ1 , γ2 < 1
(77)
and the fractional Fokker–Plank equation [Wu & Lu, 2010] ∂u ∂u 1−γ ∂ = RL D0,t d(x) + f (x, t) , ∂t ∂x ∂x 0 0. In his approach, the temporal discretization is similar to (48) except that the predictor is obtained by repeating procedures (52)– (54) with α = γ in (52) and tk+1 (tk+1 − s)γ−1 g(s)ds 0
≈
(82)
γ D0,t is defined by (36) (L1 scheme). Paramwhere C eter λ = 1 corresponds to the explicit method, and λ = 0 corresponds to the explicit one. Shen et al. [2006] proved that the explicit method was convergent and conditionally stable, while the implicit one was also convergent and unconditionally stable [Zhuang & Liu, 2006]. Murio [2008] also established the implicit scheme for (58) as [Zhuang & Liu, 2006], while the unconditional stability was proved by the Fourier method. In [Rihan, 2009], the similar idea was adopted to solve (58) and extended to the time fractional delay differential equations.
(83)
tk
0
(tk+1 − s)γ−1 g˜k (s)ds tk+1
+ tk
(tk+1 − s)γ−1 g(tk )ds,
(84)
where g˜k (s) is defined by (55). Both the first-order and second-order spatial derivatives in (83) were approximated by the second-order central FDMs, and he proved that the derived method was convergent with order O(∆tmax{1+2γ,2} + ∆x2 ) and conditionally stable. Scherer et al. [2008] proposed the explicit and implicit numerical methods for the fractional heat equation similar to (58) with the homogeneous initial conditions. They discretized
1230014-15
May 9, 2012
10:43
WSPC/S0218-1274
1230014
C. Li & F. Zeng
the Caputo derivative by the Gr¨ unwald–Letnikov formula, and explored the stability, convergence and error behavior of the derived schemes. Ghazizadeh et al. [2010] developed two finite difference methods, namely an explicit predictor– corrector method and totally implicit method, for fractional Cattaneo equation. They first established the explicit predictor–corrector method for the following scalar time fractional equation γ CD0,t u
+
∂f (u) = 0, ∂x
(x, t) ∈ (xL , xR ) ⊗ (0, T ], (85)
with zero initial conditions and 0 < γ < 1. The explicit predictor–corrector method is the fractional generalization of the well-known MacCormack scheme for solving hyperbolic equations. The derived predictor–corrector method is formulated as Predictor: uk+1 =− i Corrector:
k+1 j=1
ωjγ uk+1−j − i
1 = − uk+1 j 2
k+1 j=1
∆tγ k (f − fik ), ∆x i+1
(86)
ωjγ
5.1.2. Space fractional equations This section deals with the existing FDMs for onedimensional space FPDEs. Lynch et al. [2003] considered the following partial differential equations of fractional order ∂u γ = χ RL D0,x [u(x, t) − u(0, t)] + F (x, u), ∂t (x, t) ∈ (0, xR ) ⊗ (0, T ],
(90)
with suitable initial and boundary conditions, 1 ≤ γ ≤ 2. They established two schemes, namely an explicit scheme and a semi-implicit one, which can be formulated as follows Explicit scheme:
ωjγ uk+1−j j
γ D0,x [uki − uk0 ] + F (xi , uki ), δtk uki = χ RL
∆tγ k+1 , (f − f k+1 ) + uk+1 i i ∆x j+1
(87)
γ D0,x [uki − uk0 ] + F (xi , uki ) δtk uki = χ RL
− ∆t χ[w1 δtk uki−1 + w0 δtk uki + w−1 δtk uki+1 ],
(−1)j ( j ),
= and the variables with where γ bar are predicted values calculated from (86). The above predictor–corrector method (86) and (87) is called Generalized MacCormack (GMcC) scheme. The stability for GMcC scheme was analyzed by Fourier method with f (u) = cu where c = const., and the corresponding CFL number was derived. They also extended GMcC scheme to systems of time-fractional conservation equations like (85), and discussed the stability of the derived scheme. Finally, Ghazizadeh et al. [2010] established a fully implicit scheme for the following generalized Cattaneo equation ∂2u ∂u 1+γ + τ γ CD0,t u = D 2, ∂t ∂x (x, t) ∈ (xL , xR ) ⊗ (0, T ]
(88)
which can be written as 1+γ k D0,t uj = Dδx2 ukj , δ˜t ukj + τ γ C
(91)
Semi-implicit scheme:
−
3uk −4uk−1 +uk−2
1+γ k j j j where δ˜t ukj = , and C D0,t uj is 2∆x defined by (37) (L2 method). They got the experimental accuracy of the implicit scheme (89) and numerical tests suggested that the implicit method was unconditionally stable, but it has not been proven theoretically.
(89)
(92) ˆ i ) defined in (34) (or (35)), and where wi is Wi (or W γ is the L2 approximate operator defined by RL D 0,x
(34) or L2C approximate operator defined by (35) accordingly. They applied the two methods to three different types of problem to check the stability, efficiency and accuracy of the schemes, and found that the semi-implicit method was more effective than the explicit one, the L2 method was the most accurate for γ > 1.5, the L2C method was the most accurate for γ < 1.5, and around γ = 1.5 both methods had similar accuracy, which was accorded to the theoretical analysis. Liu et al. [2004] also adopted the L2 scheme to discretize the space-fractional Fokker–Plank equation. They transformed the original equation into a system of ODE which was solved by method of lines that was first introduced by Liu et al. [2004]. Shen and Liu [2005] proposed an explicit FDM based
1230014-16
May 9, 2012
10:43
WSPC/S0218-1274
1230014
Finite Difference Methods for Fractional Differential Equations
on the L2 scheme and method of lines to solve the space-fractional diffusion equation, they proved the explicit method was conditionally stable and convergent. Deng et al. [2004] discussed the numerical methods of the fractional advection–dispersion equation based on a three-point approximation for the fractional derivative and the shifted Gr¨ unwald– Letnikov formula. Meanwhile, they applied the von Neumann stability analysis to the established methods. In [Meerschaert & Tadjeran, 2004], two methods were proposed for the following space-fractional advection–dispersion equation (ADE) ∂u ∂u = −v(x) + d(x)RL DxαL ,x u(x, t) + f (x, t), ∂t ∂x (x, t) ∈ (xL , xR ) ⊗ (0, T ], (93) with suitable initial and boundary conditions, and 1 < α ≤ 2. The methods read as ! DxαL ,x uki + fik + (1 − λ) δt uki = λ −vi δx uki−1 + di RL ×
−vi δx uk+1 i−1
+
di RL DxαL ,x uk+1 i
+
fik+1
! . (94)
Meerschaert and Tadjeran proved that the explicit Euler method (λ = 1), the implicit Euler method (λ = 0) and the Crank–Nicolson method (λ = 1/2) for (93) with the spatial derivative discretized by the standard Gr¨ unwald–Letnikov formula (24) were all unconditionally unstable (in the classical ADE, the implicit Euler method and the Crank–Nicolson method are unconditionally stable), while the implicit Euler method and the Crank–Nicolson method with the spatial derivative discretized by the shifted Gr¨ unwald–Letnikov formula (26) were all unconditionally stable. Tadjeran et al. [2006] examined the Crank– Nicolson method for (93) with v(x) = 0 based on the shifted Gr¨ unwald–Letnikov formula (see (94), λ = 1/2, vi = 0). The derived Crank–Nicolson method is convergent with order O(∆t2 + ∆x) and unconditionally stable, and they obtained the asymptotic expansion of the error estimate in spatial direction which was used as the theoretical basis to construct the second-order method both in time and space. Su et al. [2010] proposed a fractional weighted average finite difference method for (93)
(v(x) and d(x) are replaced by v(x, t) and d(x, t), respectively) based on the shifted Gr¨ unwald– Letnikov formula (26) ! δt uki = λ −vik δx uki−1 + dki RL DxαL ,x uki + fik + (1 − λ) −vik+1 δx uk+1 i−1 ! k+1 α + dk+1 + fik+1 . RL DxL ,x ui i
(95)
They discussed the accuracy, consistency and stability of their method, and the results of stability analysis are similar to [Yuste, 2006]. Meerschaert and Tadjeran [2006] considered the following two-sided space-fractional partial differential equation ∂u α = c+ (x, t)RL DxαL ,x u(x, t) + c− (x, t)RL Dx,x R ∂t × u(x, t) + s(x, t),
(x, t) ∈ (xL , xR ) ⊗ (0, T ] (96)
with suitable initial and boundary conditions and 1 ≤ α ≤ 2. They established two methods, namely the explicit Euler method and the implicit Euler method, based on the shifted Gr¨ unwald–Letnikov formula (26) and (27) for (96) which can be written as DxαL ,x uki δt uki = λ + (c+ )ki RL α + (c− )ki RL Dx,x uk + ski R i
!
k+1 α + (1 − λ) (c+ )k+1 RL DxL ,x ui i
! k+1 k+1 α + (c− )k+1 D u + f . RL x,xR i i i
(97)
They proved that the explicit Euler method was conditionally stable and the implicit one was unconditionally stable. They also concluded that if the standard Gr¨ unwald–Letnikov formula was used in (97), the scheme was unstable. The stability results for the scheme (97) are a generalization and unification for the corresponding results in the classical parabolic (α = 2) and hyperbolic (α = 1) PDEs. Similar to [Meerschaert & Tadjeran, 2006], Su et al. [2009] established the fractional Crank– Nicolson scheme based on the shifted Gr¨ unwald– Letnikov formula (26) and (27) for the following
1230014-17
May 9, 2012
10:43
WSPC/S0218-1274
1230014
C. Li & F. Zeng
space-fractional advection–diffusion equation
In [Shen et al., 2008a], a discrete random walk model, the explicit and implicit FDMs were proposed for the following Riesz fractional advection– dispersion equation (RFDE)
∂u ∂u = −v(x, t) + c+ (x, t)RL DxαL ,x u(x, t) ∂t ∂x α u(x, t) + c− (x, t)RL Dx,x R
+ s(x, t),
∂u = Kα RZ Dxα u + Kβ RZ Dxβ u + f (x, t), ∂t
(x, t) ∈ (xL , xR ) ⊗ (0, T ], 1 < α ≤ 2,
(98)
with the first-order spatial derivative discretized by the backward difference method, and their method was unconditionally stable and convergent with order O(∆t2 + ∆x). Sousa [2009] also developed three difference methods (the upwind, central and Lax–Wendroff methods) based on the shifted Gr¨ unwald–Letnikov formula for (98). The derived numerical schemes can be seen as the generalizations of existing schemes for the integer-order advection diffusion equations. They proved that their methods were consistent, convergent of order one and conditionally stable by using von Neumann analysis and Fourier analysis. At the same time, the advantages and weaknesses of the three methods were analyzed. Ding et al. [2010] presented a class of weighted finite difference methods based on the shifted Gr¨ unwald–Letnikov formula for (98) with c+ (x, t) = c+ (x) and c− (x, t) = c− (x), the derived method is the same as (97) except that (c+ )ki and (c− )ki are replaced by (c+ )i and (c− )i respectively. They proved that the established method was convergent of order O(∆t + ∆x) except for the case λ = 1/2, whose accuracy in time was convergent with order O(∆t2 ). The stability results are the same as [Yuste, 2006]. Liu et al. [2007a] presented a random walk model and an explicit FDM based on the shifted Gr¨ unwald–Letnikov formula for the numerical solution of the L´evy–Feller advection–dispersion equation ∂u ∂u = aDθα u − b , ∂t ∂x
(x, t) ∈ (xL , xR ) ⊗ (0, T ], (99)
with suitable initial and boundary conditions, 1 < α ≤ 2, Dθα is the Riesz–Feller fractional derivative (in space) of order α and skewness of θ, which can be expressed by the linear combination of the left and right Riemann–Liouville differential operators. They proved that the presented method was conditionally stable and convergent of order one both in time and space.
(x, t) ∈ (xL , xR ) ⊗ (0, T ],
(100)
with suitable initial and boundary conditions, 1 < α ≤ 2, 0 < β < 1, and RZ Dxα is the Riesz space derivative operator defined by (8). The Riesz derivative operator RZ Dxα and RZ Dxβ are approximated by the shifted Gr¨ unwald–Letnikov formula and standard Gr¨ unwald–Letnikov formula respectively, the two FDMs can be written as ! Dxα uki + Kβ RZ Dxβ uki + fik + (1 − λ) δik = λ Kα RZ ! β k+1 k+1 . Dxα uk+1 + K D u + f × Kα RZ β RZ x i i i (101) They proved both methods were convergent, the explicit one (λ = 1) was conditionally stable and the implicit one (λ = 0) was unconditionally stable. In (100), if Kβ = 0, then equation reduces to the Riesz fractional diffusion equation (RFDE). Yang et al. [2010a] first established three schemes for the RFDE with the space Riesz derivative approximated by the L2 scheme, the shifted Gr¨ unwald–Letnikov formula and Matrix transformed method respectively, such that the RFDE was transformed into a system of time ODEs which were solved by using a differential/algebraic solver [Brenan et al., 1989]. The same procedures were repeated to solve the RFADE (100) with the advection term RZ Dxβ u approximated by the L1 scheme and the standard Gr¨ unwald–Letnikov formula. They conducted enough experiments to demonstrate the effectiveness and convergence of the established methods. Here, we especially mention that Wang et al. [2010] developed a fast FDM based on the shifted Gr¨ unwald–Letnikov formula for (96), which only requires storage of O(N ) and computational cost of O(N log2 N ) while retaining the accuracy and property as the regular difference method that requires storage and computational cost of O(N 2 ) and O(N 3 ) respectively.
1230014-18
May 9, 2012
10:43
WSPC/S0218-1274
1230014
Finite Difference Methods for Fractional Differential Equations
5.1.3. Time-space fractional equations
The explicit method can be derived from (105)
This section displays the difference methods for time-space fractional equations. Liu et al. [2006] presented an implicit finite difference method for the space-time fractional diffusion equation of the form
γ k+1 CD0,t ui
γ CD0,t u
α = RL D0,x u,
0 < γ ≤ 1, 1 < α ≤ 2, (x, t) ∈ (xL , xR ) ⊗ (0, T ], (102)
with the Caputo derivative discretized by (36) (L1 method), and the Riemann–Liouville derivative discretized by the shifted Gr¨ unwald–Letnikov formula γ k CD0,t ui
α uk . = RL D0,x i
(103)
They proved that the established scheme was unconditionally stable and convergent with order one. In [Yang et al., 2009b], the time and spacesymmetric fractional diffusion equation similar to (102) was considered by Yang et al. [2009a]. They proposed two methods with the Caputo time fractional derivative approximated by the FDM and the Laplace transform method, the symmetric space fractional derivative was approximated by using the matrix transform method. They conducted numerical experiments to verify their methods. For the following space-time fractional advection–dispersion equation [Liu et al., 2007b] γ CD0,t u
(x, t) ∈ (xL , xR ) ⊗ (0, T ], (104)
with suitable initial and boundary conditions, and 0 < γ ≤ 1, 0 < β ≤ 1, 1 0, κ(x, t) > 0. Liu et al. [2007b] proposed an implicit scheme and an explicit scheme for (104). The implicit scheme can be obtained with the Caputo derivative approximated by (36), the Riemann–Liouville derivatives β α RL D0,x and RL D0,x approximated by the standard Gr¨ unwald–Letnikov formula and shifted Gr¨ unwald– Letnikov formula respectively γ k+1 CD0,t ui
β = −vik+1 RL D0,x uk+1 i k+1 α + κk+1 RL D0,x ui i
+
fik+1 .
(106) They proved that the two methods were convergent with order one both in time and space, the implicit one was unconditionally stable and the explicit one was conditionally stable. A similar idea is also adopted by Zhang [2009] to establish an implicit unconditionally and convergent scheme for the space-time fractional convection–diffusion equation. Shen et al. [2011] also established the explicit and implicit difference methods similar to (105) and (106) for the space-time Riesz–Caputo fractional advection–diffusion equation. They proved that the explicit method was conditionally stable and convergent of order one, and the implicit method was unconditionally stable and convergent of order O(∆t2−γ + ∆x). At the same time, a Richardson extrapolation method was used to obtain higher-order accuracy and the short-memory principle was used to investigate the effect of the amount of computations for the established methods. For the time-space fractional Fokker–Planck equation, Yang et al. [2009a] first considered the following Fokker–Planck equation ∂ V (x) ∂u 1−γ α α = RL D0,t + Kγ RZ Dx u + f (x, t) , ∂t ∂x mηγ (x, t) ∈ (xL , xR ) ⊗ (0, T ], (107)
β α = −v(x, t)RL D0,x u + κ(x, t)RL D0,x u
+ f (x, t),
β α uk + f k+1 . = −vik RL D0,x uki + κki RL D0,x i i
with suitable initial and boundary conditions, RZ Dxα is the Riesz derivative operator defined by (8). Yang et al. [2009b] first adopted the L2 scheme, the shifted Gr¨ unwald–Letnikov formula and matrix transform method to transform (107) into a system of time FODE, then the fractional implicit trapezoidal method was used to solve the derived FODE. They checked the efficiency and accuracy of the algorithm by a numerical experiment. Yang et al. [2010a] also considered another form of the nonlinear time-space fractional Fokker– Planck equation as ∂ V (x) γ α α D u = + K D C 0,t γ RZ x u + f (u, x, t), ∂x mηγ (x, t) ∈ (xL , xR ) ⊗ (0, T ],
(105) 1230014-19
(108)
May 9, 2012
10:43
WSPC/S0218-1274
1230014
C. Li & F. Zeng
with suitable initial and boundary conditions, RZ Dxα is the Riesz derivative operator defined by (8). They proposed a numerical method with the time fractional derivative approximated by (36) (L1 scheme), the space Riesz fractional derivative approximated by the shifted Gr¨ unwald–Letnikov formula and the first-order spatial derivative approximated by the backward difference scheme if V < 0 (otherwise, the forward difference scheme can be used if V > 0) with the nonlinear term dealt with explicitly. They proved that the established method was convergent with order O(∆t + ∆x) and conditionally stable under certain conditions. In [Shen et al., 2008b], a discrete random walk model based on an explicit finite difference approximation was proposed for the following time-space fractional diffusion equation of distributed order γ CD0,t u
= Kα RZ Dxα u + Kβ RZ Dxβ u + f (x, t), (x, t) ∈ (xL , xR ) ⊗ (0, T ], (109)
with suitable initial and boundary conditions, 1 < α ≤ 2, 0 < β, γ < 1. The Riesz derivative operator was discretized as [Shen et al., 2008a], for detailed information refer to [Shen et al., 2008b].
5.2. Difference methods of high-dimensional FPDEs In this section, we review some existing FDMs for two- and three-dimensional FPDEs. Meerschaert et al. [2006] considered the following two-dimensional fractional diffusion equation ∂u = d(x, y)RL DxαL ,x u ∂t + e(x, y)RL DyβL ,y u + q(x, y, t), (x, y, t) ∈ (xL , xR ) ⊗ (yL , yR ) ⊗ (0, T ],
(110)
with the suitable initial and boundary conditions, 1 < α, β ≤ 2, d(x, y) > 0 and e(x, y) > 0. They established an implicit Euler method based on the shifted Gr¨ unwald–Letnikov formula that can be written as
k+1 = uki,j + ∆tqi,j .
(112)
Compared to (111), the derived ADI–Euler method (112) has the perturbation of order O(∆t2 ), and is unconditionally stable under certain conditions. The ADI technique can convert the two-dimensional problem into several separate onedimensional ones, which can be solved in parallel. This technique is often adopted in the numerical schemes for higher-dimensional problems. In [Tadjeran & Meerschaert, 2007], the authors presented an ADI method for (110) with the time derivative discretized by a Crank–Nicolson method (ADI–CN), and the space fractional derivative discretized by the shifted Gr¨ unwald–Letnikov formula. The derived scheme can be written as ∆t ∆t β di,j RL ei,j RL DyL ,y uk+1 DxαL ,x 1− 1− i,j 2 2 ∆t ∆t β α di,j RL ei,j RL DyL ,y = 1+ DxL ,x 1+ 2 2 ∆t k+1 k (q + qi,j ). (113) 2 i,j They proved that under certain conditions, the scheme (113) was consistent, unconditionally stable and convergent of order O(∆t2 + ∆x + ∆y). Moreover, they also proved that the truncation error had the form O(∆t2 )+ KO(∆x)+ M O(∆y)+ O(∆x2 )+ O(∆y 2 ), which is the theoretical basis for constructing the second-order method both in time and in space by the extrapolation method. Similar to (110), Chen and Liu [2008] proposed an implicit ADI–Euler method for the following two-dimensional fractional advection–dispersion equation × uki,j +
∂u = a(x)RL DxαL1 ,x u + b(y)RL DyβL1 ,y u ∂t
DxαL ,x uk+1 δt uki,j = di,j RL i,j k+1 + ei,j RL DyβL ,y uk+1 i,j + qi,j .
each time step, a very large nonsparse linear system needs to be solved, which is computationally costly. In order to efficiently solve (111) without losing accuracy, an alternating direction implicit Euler method (ADI–Euler) was developed by Meerschaert et al. [2006] which can be written as ! ! 1 − ∆tdi,j RL DxαL ,x 1 − ∆tei,j RL DyβL ,y uk+1 i,j
(111)
The above scheme (111) is unconditionally stable and convergent of order O(∆t + ∆x), while at 1230014-20
− c(x)RL DxαL2 ,x u − d(y)RL DyβL2 ,y u + q(x, y, t), (x, y, t) ∈ (xL , xR ) ⊗ (yL , yR ) ⊗ (0, T ], (114)
May 9, 2012
10:43
WSPC/S0218-1274
1230014
Finite Difference Methods for Fractional Differential Equations
with the suitable initial and boundary conditions, 1 < α1 , β1 ≤ 2, 0 < α2 , β2 ≤ 1, a(x), b(y), c(x), d(y) > 0. The derived implicit ADI–Euler method is similar to (112), while the advection terms RL Dxα2 ,x u L
and RL DyβL2 ,y u are approximated by the standard Gr¨ unwald–Letnikov formulae. They proved that their method was convergent and unconditionally stable, meanwhile, they also adopted the extrapolation technique to improve the accuracy of their method as [Tadjeran & Meerschaert, 2007]. Zhuang and Liu [2007b] constructed an implicit finite difference method for the following twodimensional time fractional diffusion equation γ CD0,t u
= a(x, y, t)
∂2u ∂2u + b(x, y, t) + f (x, y, t), ∂x2 ∂y 2
(x, y, t) ∈ (xL , xR ) ⊗ (yL , yR ) ⊗ (0, T ], (115) with the suitable initial and boundary conditions, 0 < γ < 1, a(x, y, t) > 0, b(x, y, t) > 0. The established scheme can be written as γ k+1 CD0,t ui,j
k+1 2 k+1 k+1 2 k+1 = ak+1 i,j δx ui,j + bi,j δy ui,j + fi,j ,
The time fractional derivative was approximated by the standard Gr¨ unwald–Letnikov formula, Chen et al. [2008a] established two methods for (117) which can be written as # 1−γ " 2 k+1 δt uki,j = (1 − λ) RL D0,t κ1 δx2 uk+1 i,j + κ2 δy ui,j 2 k+1 k+1 + κ3 δx2 uk+1 i,j + κ4 δy ui,j + fi,j
+λ
∂2u + κ4 2 + f (x, y, t), ∂y (x, y, t) ∈ (xL , xR ) ⊗ (yL , yR ) ⊗ (0, T ], (117) with the suitable initial and boundary conditions, and 0 < γ < 1, κ1 , κ2 , κ3 , κ4 are positive constants.
+ κ2 δy2 uki,j
k , + κ3 δx2 uki,j + κ4 δy2 uki,j + fi,j
# (118)
where λ = 1 corresponds to the explicit method and λ = 0 corresponds to the implicit one. They proved that the two methods (118) were L2 -convergent with order O(∆t + ∆x2 + ∆y 2 ), the explicit method was conditionally stable and the implicit one was unconditionally stable. Chen et al. [2010d] proposed two methods for the following two-dimensional time fractional diffusion equation which can be seen as a modification of (115) ∂2u ∂2u ∂u 1−γ = RL D0,t a(x, y, t) 2 + b(x, y, t) 2 ∂t ∂x ∂y
(116) γ D0,t is the approximate operator defined where C by (36) (L1 scheme). They proved that the derived scheme (116) was unconditionally stable and convergent of order O(∆t + ∆x2 + ∆y 2 ). Brunner et al. [2010] also considered (115) with a(x, y, t) = 1, b(x, y, t) = 1, and they proposed an algorithm that coupled an adaptive time stepping and adaptive spatial basis selection approach. Chen et al. [2008a] considered the twodimensional Rayleigh–Stokes problem for a heated generalized second grade fluid with time fractional derivative model with a nonhomogeneous term ∂2u ∂2u ∂2u ∂u 1−γ = RL D0,t κ1 2 + κ2 2 + κ3 2 ∂t ∂x ∂y ∂x
1−γ " 2 k RL D0,t κ1 δx ui,j
+ g(x, y, t), (x, y, t) ∈ (xL , xR ) ⊗ (yL , yR ) ⊗ (0, T ], (119) with the suitable initial and boundary conditions, 0 < γ < 1, a(x, y, t) > 0, b(x, y, t) > 0. The two methods, namely the implicit method and the explicit one, can be written as 1−γ " k+1 2 k+1 D0,t ai,j δx ui,j δt uki,j = (1 − λ) RL # 2 k+1 k+1 + bk+1 + gi,j i,j δy ui,j +λ
# 1−γ " k 2 k k 2 k k D δ u + b δ u a + g RL 0,t i,j x i,j i,j y i,j i,j , (120)
1−γ D0,t is defined by (31), that is, the time where RL fractional derivative is approximated by the L1 scheme. They proved that the two methods were L2 -convergent with order O(∆t + ∆x2 + ∆y 2 ), the explicit method was conditionally stable, and the implicit one was unconditionally stable. Meanwhile, the multivariate extrapolation was investigated to improve the accuracy of the derived methods.
1230014-21
May 9, 2012
10:43
WSPC/S0218-1274
1230014
C. Li & F. Zeng
For the two-dimensional time-space fractional equations, Zhuang and Liu [2007a] proposed an implicit method for the following time-space fractional diffusion equation γ CD0,t u
= a(x, y, t)RL Dxα
L
,x u
+ b(x, y, t)RL Dyβ
L
,y u
+ f (x, y, t), (x, y, t) ∈ (xL , xR ) ⊗ (yL , yR ) ⊗ (0, T ], (121) with the suitable initial and boundary conditions, 1 < α, β ≤ 2, 0 < γ ≤ 1. In their method, the time fractional derivative was approximated by (36), and the two space fractional derivatives were approximated by the shifted Gr¨ unwald–Letnikov formula. They proved that their method was unconditionally stable and convergent with order O(∆t + ∆x + ∆y). For the three-dimensional problems, Liu et al. [2008] first considered the noncontinued seepage flow in uniform media of the form ∂u " γ1 γ2 γ3 # = KxRL D0,x u + Ky RL D0,y u + Kz RL D0,z u ∂t
In real computation, another equivalent efficient factored form of (124) was presented by Liu et al. [2008, (see (2.8)–(2.10))]. They proved that the fractional ADI scheme (124) and its equivalent form were unconditionally stable, consistent and convergent of order O(∆t + ∆x + ∆y + ∆z) by the Fourier method. Then Liu et al. [2008] also established an ADI scheme based on the classical Douglas scheme for the continued seepage flow in uniform media, and the derived scheme was also unconditionally stable, consistent and convergent of order O(∆t2 + ∆x + ∆y + ∆z). They utilized the Richardson extrapolation technique to improve the space accuracy of the derived method based on the theory in [Tadjeran et al., 2006]. In [Chen & Liu, 2009], the authors established an implicit method for the following three-dimensional space Galilei invariant fractional advection–diffusion equation ∂u ∂u ∂u ∂u + κ1 + κ2 + κ3 ∂t ∂x ∂y ∂z ∂2u ∂2u ∂2u 1−γ ˜γ 2 = RL D0,t κγ 2 + κγ 2 + κ ∂x ∂y ∂z
+ f (x, y, z, t),
+ f (x, y, z, t),
(x, y, z, t) ∈ (0, xR ) ⊗ (0, yR ) ⊗ (0, zR )
(x, y, z, t) ∈ (0, xR ) ⊗ (0, yR ) ⊗ (0, zR )
⊗ (0, T ],
⊗ (0, T ], (122)
(125)
with the suitable initial and boundary conditions, 1 < γ1 , γ2 , γ3 ≤ 2, Kx , Ky , Kz > 0. The spatial fractional derivatives are approximated by the shifted Gr¨ unwald–Letnikov formula, Liu et al. [2008] proposed an implicit scheme which can be written as
with the suitable initial and boundary conditions, ˜ γ are positive 0 < γ < 1, κ1 , κ2 , κ3 , κγ , κγ , κ constants. The time fractional derivative is approximated by the standard Gr¨ unwald–Letnikov formula, Chen and Liu [2009] established the following numerical scheme
γ1 k+1 γ2 k+1 D0,x ui,j,l + Ky RL D0,y ui,j,l δt uki,j,l = Kx RL γ3 k+1 k+1 + Kz RL D0,z ui,j,l + fi,j,l .
(123)
In order to improve the efficiency of the established scheme (123), a fractional ADI scheme with the same accuracy as (123) was constructed, which can be written in the following operator form ! ! γ1 γ2 D0,x D0,x 1 − ∆tKx RL 1 − ∆tKx RL ! γ3 D0,x uk+1 × 1 − ∆tKx RL i,j,l k+1 = uki,j,l + τ fi,j,l .
(124)
k+1 k+1 δt uki,j,l + κ1 δxˆ uk+1 i,j,l + κ2 δyˆui,j,l + κ3 δzˆui,j,l
# 1−γ " 2 k+1 = RL D0,t ˜ γ δz2 uk+1 κγ δx2 uk+1 i,j,l + κγ δy ui,j,l + κ i,j,l k+1 . + fi,j,l
(126)
They proved that the above scheme (126) was solvable, L2 -convergent of order O(∆t + ∆x2 + ∆y 2 + ∆z 2 ), and unconditionally stable by using Fourier method. Meanwhile, the Richardson extrapolation technique was adopted to improve the accuracy of the derived method, which was verified by the numerical experiment.
1230014-22
May 9, 2012
10:43
WSPC/S0218-1274
1230014
Finite Difference Methods for Fractional Differential Equations
6. Finite Methods for the Variable-Order FPDEs
∂u ∂u = −κ(x, t) + ν(x, t)RZ Dxα(x,t) u(x, t) ∂t ∂x
In this section, we review some numerical methods for the variable-order FPDEs. In real world, the behavior of some diffusion processes in response to temperature changes may be better described by using variable-order equations [Lorenzo & Hartley, 2002]. Lorenzo and Hartley presented the concept of variable-order fractional integration and differentiation [Lorenzo & Hartley, 2002], for more definitions and numerical approximations to the variable-order fractional derivatives, refer to [Val´erio & Costa, 2011]. Here we just introduce one case of definitions of variable-order Gr¨ unwald–Letnikov, Riemann– Liouville and Caputo derivatives. Definitions of left variable-order Gr¨ unwald–Letnikov and Caputo fractional derivatives of f (t) with fractional order α(t) can be obtained from (2) and (6) by replacing α with α(t) in (2) and (6) respectively. The left variable-order Riemann–Liouville derivative of the given function f (t) with order α(t) is defined by 1 α(t) RL Da,t f (t) = Γ(n − α(t)) n θ d n−α(t)−1 (θ − s) f (s)ds , × dθ n a θ=t n − 1 ≤ α(t) < n. (127) The right variable-order Gr¨ unwald–Letnikov, Riemann–Liouville and Caputo derivatives can be defined accordingly, see [Zhuang et al., 2009]. The definitions for the variable-order derivatives in this paper can be seen as a direct extension of the constant-order fractional derivative, so that the relationships (such as (9)) also exist for variableorder derivatives. The left and right variable-order Riemann–Liouville derivatives can be approximated by the shifted Gr¨ unwald formulae for variable-order derivative, which are similar to (26) and (27) respectively, for detailed information refer to [Zhuang et al., 2009; Val´erio & Costa, 2011]. Now, we introduce the existing FDMs for the variable-order FPDEs, which were mainly achieved by Liu and his coauthors [Zhuang et al., 2009; Chen et al., 2010a, 2010b, 2010c]. Zhuang et al. [2009] established two methods for the following variable-order fractional advection–diffusion equation with a nonlinear source term
+ f (u, x, t),
(x, t) ∈ (xL , xR ) ⊗ (0, T ], 1 < α ≤ α(x, t) ≤ α < 2, (128)
with the suitable initial and boundary conditions, α(t) κ(x, t), ν(x, t) > 0 and RZ Dt is the variable-order Riesz derivative of order α(t) which can be derived from (8) by replacing α with α(t). The derived method was similar to (94), which can be written as ! αk δt uki = λ −νik δx uki−1 + κki RL DxLi ,x uki + (1 − λ) ! αk+1 k+1 k+1 i + κ D u × −νik+1 δx uk+1 x ,x RL L i−1 i i + fik ,
(129)
αk where RZ Dx i is the approximate operator of α(xi ,tk ) , which can be defined similarly as RZ Dx the Riesz derivative operator for constant-variable derivative based on the shifted Riesz formula. They proved that the two derived methods were both convergent of order one, the explicit method (λ = 1) was conditionally stable and the implicit one (λ = 0) was unconditionally stable. Moveover, they also presented a fractional method of lines, a matrix transfer technique, and an extrapolation method for Eq. (128). Lin et al. [2009] proposed an explicit stable finite difference method of order one for the variable-order nonlinear fractional diffusion equation of the form (128) with κ(x, t) = 0, the derived method can be seen as a special case of (129) with λ = 1, κki = 0. In 2010, Chen et al. [2010a] proposed a numerical scheme with first-order temporal accuracy and fourth-order spatial accuracy for the variable-order anomalous subdiffusion equation, which can be derived from (57) with γ replaced by γ(x, t). The derived method was unconditionally stable, which can be seen as a generalization of the numerical method proposed by Cui [2009]. Meanwhile, another improved numerical scheme with secondorder temporal accuracy and fourth-order spatial accuracy was also proposed by Chen et al. [2010a]. Chen et al. [2010c] also proposed two methods for the following two-dimensional variable-order
1230014-23
May 9, 2012
10:43
WSPC/S0218-1274
1230014
C. Li & F. Zeng
anomalous subdiffusion equation of the form ∂2u ∂2u ∂u 1−γ(x,y,t) κ1 2 + κ2 2 + f (x, y, t), = RL D0,t ∂t ∂x ∂y (x, t) ∈ (xL , xR ) ⊗ (yL , yR ) ⊗ (0, T ], (130) with the suitable initial and boundary conditions, κ1 , κ2 > 0. The variable-order fractional derivative was approximated by the Gr¨ unwald formula, the two derived numerical methods are similar to (64), which can be seen as an extension from onedimensional problem (57) to the two-dimensional case. They proved the two methods were convergent of order one, the explicit method was conditionally stable and the implicit one was unconditionally stable by Fourier method. The variable-order fractional derivative was approximated by the Gr¨ unwald formula for the variable-order derivative, the first and second spatial derivatives were approximated by the second central difference method with the nonlinear term treated explicitly, Chen et al. [2010b] proposed a numerical scheme for the following variable-order Galilei advection diffusion equation with a nonlinear source term: ∂u ∂2u ∂u 1−γ(x,t) +κ = RL D0,t κγ 2 + f (u, x, t), ∂t ∂x ∂x (x, t) ∈ (xL , xR ) ⊗ (0, T ], (131) with the suitable initial and boundary conditions, 0 < γmin ≤ γ(x, t) ≤ γmax < 1 and κ, κγ > 0. They proved the derived method was unconditionally stable, at the same time, another numerical method was developed to improve the temporal accuracy of the numerical method. Remark 6.1. Recently, the technique of the FDMs for variable-order FDEs are mainly achieved by Liu and his co-operators, and these methods are mainly based on the definition of the variable-order derivative.
7. Conclusion and Remark In this paper, we present a survey on the FDMs for the FODEs and the time, space and time-space FPDEs which characterizes anomalous diffusion or
dispersion with one, two and three dimensions. They include the fractional diffusion equations, the fractional advection–diffusion equations, the fractional advection–dispersion equations, the fractional reaction–diffusion equations and the fractional Fokker–Planck type equations, and so on. From the publications available, the FDMs and their numerical analysis (convergence and stability) for FPDEs mainly focus on the linear cases, where the classical analysis technique such as von Neumann stability analysis, Fourier method and energy method are extended to the numerical analysis of FPDEs. Some of the existing FDMs can be seen as the generalizations of the numerical methods for classical equations. At the same time, some techniques, such as Richard extrapolation method, ADI methods and the short memory principal, are adopted to try to reduce the computational cost and storage of the established schemes. However, the short memory principal is not a better choice for simulating the dynamical behaviors such as limit cycle and chaotic attractors. It should be noted that the studies for linear cases are very rich. However, very limited work has been done for nonlinear ones. To investigate the efficient numerical methods for nonlinear FDEs is a long-term duty.
References Agrawal, O. P. & Kumar, P. [2007] “Comparison of five numerical schemes for fractional differential equations,” Adv. Fract. Cal. 47, 43–60. Baker, C. T. H. [2000] “A perspective on the numerical treatment of Volterra equations,” J. Comput. Appl. Math. 125, 217–249. Baleanu, D., Defterli, O. & Agrawal, O. P. [2009] “A central difference numerical scheme for fractional optimal control problems,” J. Vib. Contr. 15, 583–597. Barkai, E., Metzler, R. & Klafter, J. [2000] “From continuous time random walks to the fractional Fokker– Planck equation,” Phys. Rev. E 61, 132–138. Bhalekar, S., Daftardar-Gejji, V., Baleanu, D. & Magin, R. [2011] “Fractional Bloch equation with delay,” Comput. Math. Appl. 61, 1355–1365. Brenan, K. E., Campbell, S. L. & Petzold, L. R. [1989] Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations (North Holland, NY). Brunner, H., Ling, L. & Yamamoto, M. [2010] “Numerical simulations of 2D fractional subdiffusion problems,” J. Comput. Phys. 229, 6613–6622. Chen, C. M., Liu, F., Turner, I. & Anh, V. [2007] “A Fourier method for the fractional diffusion
1230014-24
May 9, 2012
10:43
WSPC/S0218-1274
1230014
Finite Difference Methods for Fractional Differential Equations
equation describing sub-diffusion,” J. Comput. Phys. 227, 886–897. Chen, S. & Liu, F. [2008] “ADI–Euler and extrapolation methods for the two-dimensional fractional advectiondispersion equation,” J. Appl. Math. Comput. 26, 295–311. Chen, C. M., Liu, F. & Anh, V. [2008a] “Numerical analysis of the Rayleigh–Stokes problem for a heated generalized second grade fluid with fractional derivatives,” ANZIAM. J. 47, C168–C184. Chen, C. M., Liu, F. & Burrage, K. [2008b] “Finite difference methods and a Fourier analysis for the fractional reaction-subdiffusion equation,” Appl. Math. Comput. 198, 754–769. Chen, C. M. & Liu, F. [2009] “A numerical approximation method for solving a three-dimensional space Galilei invariant fractional advection-diffusion equation,” J. Appl. Math. Comput. 30, 219–236. Chen, C. M., Liu, F. & Anh, V. [2009a] “A Fourier method and an extrapolation technique for Stokes’ first problem for a heated generalized second grade fluid with fractional derivative,” J. Comput. Appl. Math. 223, 777–789. Chen, S., Liu, F., Zhuang, P. & Anh, V. [2009b] “Finite difference approximations for the fractional Fokker– Planck equation,” Appl. Math. Model. 33, 256–273. Chen, C. M., Liu, F., Anh, V. & Turner, I. [2010a] “Numerical schemes with high spatial accuracy for a variable-order anomalous sub-diffusion equation,” SIAM J. Sci. Comput. 32, 1740–1760. Chen, C. M., Liu, F., Anh, V. & Turner, I. [2010b] “Numerical simulation for the variable-order Galilei invariant advection diffusion equation with a nonlinear source term,” Appl. Math. Comput. 217, 5729– 5742. Chen, C. M., Liu, F., Anh, V. & Turner, I. [2010c] “Numerical methods for solving a two-dimensional variable-order anomalous subdiffusion equation,” Numer. Algor. 56, 383–403. Chen, C. M., Liu, F., Turner, I. & Anh, V. [2010d] “Numerical schemes and multivariate extrapolation of a two-dimensional anomalous sub-diffusion equation,” Numer. Algor. 54, 1–21. Cui, M. R. [2009] “Compact finite difference method for the fractional diffusion equation,” J. Comput. Phys. 228, 7792–7804. Deng, W. H. [2007a] “Numerical algorithm for the time fractional Fokker–Planck equation,” J. Comput. Phys. 227, 1510–1522. Deng, W. H. [2007b] “Short memory principle and a predictor-corrector approach for fractional differential equations,” J. Comput. Appl. Math. 206, 174–188. Deng, Z. Q., Singh, V. P. & Bengtsson, L. [2004] “Numerical solution of fractional advection-dispersion equation,” J. Hydraul. Eng. 130, 422–431.
Diethelm, K. & Walz, G. [1997] “Numerical solution of fractional order differential equations by extrapolation,” Numer. Algor. 16, 231–253. Diethelm, K. [1997] “An algorithm for the numerical solution of differential equations of fractional order,” Electron. Trans. Numer. Anal. 5, 1–6. Diethelm, K. & Freed, A. D. [1999a] “The FracPECE subroutine for the numerical solution of differential equations of fractional order,” Forschung und Wissenschaftliches Rechnen: Beitr¨ age zum Heinz-BillingPreis (Gesellschaft f¨ ur wissenschaftliche Datenverarbeitung, G¨ ottingen), pp. 57–71. Diethelm, K. & Freed, A. D. [1999b] “On the solution of nonlinear fractional differential equations used in the modeling of viscoplasticity,” Scientific Computing in Chemical Engineering II–Computational Fluid Dynamics, Reaction Engineering, and Molecular Properties (Springer, Heidelberg), pp. 217–224. Diethelm, K., Ford, N. J., & Freed, A. D. [2002] “A predictor-corrector approach for the numerical solution of fractional differential equations,” Nonlin. Dyn. 29, 3–22. Diethelm, K., Ford, N. J., & Freed, A. D. [2004] “Detailed error analysis for a fractional Adams method,” Numer. Algor. 36, 31–52. Diethelm, K., Ford, N. J., Freed, A. D. & Luchko, Y. [2005] “Algorithms for the fractional calculus: A selection of numerical methods,” Comput. Meth. Appl. Mech. Eng. 194, 43–773. Diethelm, K., Ford, N. J., Freed, A. D. & Weilbeer, M. [2006] “Pitfalls in fast numerical solvers for fractional differential equations,” J. Comput. Appl. Math. 186, 482–503. Ding, Z. Q., Xiao, A. G. & Li, M. [2010] “Weighted finite difference methods for a class of space fractional partial differential equations with variable coefficients,” J. Comput. Appl. Math. 233, 1905–1914. Du, R., Cao, W. R. & Sun, Z. Z. [2010] “A compact difference scheme for the fractional diffusion-wave equation,” J. Comput. Appl. Math. 34, 2998–3007. Ford, N. J. & Simpson, A. C. [2001] “The numerical solution of fractional differential equations: Speed versus accuracy,” Numer. Algor. 26, 333–346. Gao, G. H. & Sun, Z. Z. [2011] “A compact finite difference scheme for the fractional sub-diffusion equations,” J. Comput. Phys. 230, 586–595. Galeone, L. & Garrappa, R. [2006] “On multistep methods for differential equations of fractional order,” Mediterr. J. Math. 3, 565–580. Galeone, L. & Garrappa, R. [2009] “Explicit methods for fractional differential equations and their stability properties,” J. Comput. Appl. Math. 228, 548–560. Garrappa, R. [2009] “On some explicit Adams multistep methods for fractional differential equations,” J. Comput. Appl. Math. 229, 392–399.
1230014-25
May 9, 2012
10:43
WSPC/S0218-1274
1230014
C. Li & F. Zeng
Ghazizadeh, H. R., Maerefat, M. & Azimi, A. [2010] “Explicit and implicit finite difference schemes for fractional Cattaneo equation,” J. Comput. Phys. 229, 7042–7057. Kilbas, A. A., Srivastava, H. M. & Trujillo, J. J. [2006] Theory and Applications of Fractional Differential Equations (Elsevier, Netherlands). Klages, R., Radons, G. & Sokolov, I. M. [2008] Anomalous Transport (Wiley–VCH, Berlin). Langlands, T. A. M. & Henry, B. I. [2005] “The accuracy and stability of an implicit solution method for the fractional diffusion equation,” J. Comput. Phys. 205, 719–736. Li, C. P. & Peng, G. J. [2004] “Chaos in Chen’s system with a fractional order,” Chaos Solit. Fract. 22, 443–450. Li, C. P. & Deng, W. H. [2006] “Chaos synchronization of fractional-order differential system,” Int. J. Mod. Phys. B 20, 791–803. Li, C. P., Deng, W. H. & Xu, D. [2006] “Chaos synchronization of the Chua system with a fractional order,” Physica A 360, 171–185. Li, C. P. & Deng, W. H. [2007] “Remarks on fractional derivatives,” Appl. Math. Comput. 187, 777–784. Li, C. P. & Tao, C. X. [2009] “On the fractional Adams method,” Comput. Math. Appl. 58, 1573–1588. Li, C. P., Dao, X. H. & Guo, P. [2009] “Fractional derivatives in complex planes,” Nonlin. Anal.: TMA 71, 1857–1869. Li, X. J. & Xu, C. J. [2009] “A space-time spectral method for the time fractional diffusion equation,” SIAM J. Numer. Anal. 47, 2108–2131. Li, C. P. & Zhao, Z. G. [2011a] “Numerical approximation of nonlinear fractional differential equations with subdiffusion and superdiffusion,” Comput. Math. Appl., doi:10.1016/j.camwa.2011.02.045. Li, C. P. & Zhao, Z. G. [2011b] “ Introduction to fractional integrability and differentiability,” Eur. Phys. J. — Special Topics 193, 5–26. Li, C. P., Qian, D. L. & Chen, Y. Q. [2011a] “On Riemann–Liouville and Caputo derivatives,” Discr. Dyn. Nat. Soc. 2011, Article ID 562494. Li, C. P., Chen, A. & Ye, J. J. [2011b] “Numerical approach to fractional calculus and fractional ordinary differential equations,” J. Comput. Phys. 230, 3352–3368. Lin, R. & Liu, F. [2007] “Fractional high order methods for the nonlinear fractional ordinary differential equation,” Nonlin. Anal. 66, 856–869. Lin, Y. M. & Xu, C. J. [2007] “Finite difference/spectral approximations for the time-fractional diffusion equation,” J. Comput. Phys. 225, 1533–1552. Lin, R., Liu, F., Anh, V. & Turner, I. [2009] “Stability and convergence of a new explicit finite-difference approximation for the variable-order
nonlinear fractional diffusion equation,” Appl. Math. Comput. 212, 435–445. Liu, F., Anh, V. & Turner, I. [2004] “Numerical solution of the space fractional Fokker–Planck Equation,” J. Comput. Appl. Math. 166, 209–219. Liu, F., Zhuang, P., Anh, V. & Turner, I. [2006] “A fractional-order implicit difference approximation for the space-time fractional diffusion equation,” ANZIAM J. 47, C48–C68. Liu, Q., Liu, F., Turner, I. & Anh, V. [2007a] “Approximation of the L´evy–Feller advection-dispersion process by random walk and finite difference method,” J. Comput. Phys. 222, 57–70. Liu, F., Zhuang, P., Anh, V., Turner, I. & Burrage, K. [2007b] “Stability and convergence of the difference methods for the space-time fractional advectiondiffusion equation,” Appl. Math. Comput. 191, 12–20. Liu, Q., Liu, F., Turner, I. & Anh, V. [2008] “Numerical simulation for the 3D seepage flow with fractional derivatives in porous media,” IMA J. Appl. Math. 74, 201–229. Liu, F., Yang, C. & Burrage, K. [2009] “Numerical method and analytical technique of the modified anomalous subdiffusion equation with a nonlinear source term,” J. Comput. Appl. Math. 231, 160–176. Liu, F., Yang, Q. Q. & Turner, I. [2011] “Two new implicit numerical methods for the fractional cable equation,” J. Comput. Nonlin. Dyn. 6, 1–7. Lorenzo, C. F. & Hartley, T. T. [2002] “Variable-order and distributed order fractional operators,” Nonlin. Dyn. 29, 57–98. Lubich, C. [1986] “Discretized fractional calculus,” SIAM J. Math. Anal. 17, 704–719. Lynch, V. E., Carreras, B. A., del-Castillo-Negrete, D., Ferreira–Mejias, K. M. & Hicks, H. R. [2003] “Numerical methods for the solution of partial differential equations of fractional order,” J. Comput. Phys. 192, 406–42. Magin, R. L. [2006] Fractional Calculus in Bioengineering (Begell House Publishers). Magin, R. L., Abdullah, O., Baleanu, D. & Zhou, X. J. [2008] “Anomalous diffusion expressed through fractional order differential operators in the Bloch–Torrey equation,” J. Magn. Reson. 190, 255–270. Meerschaert, M. M. & Tadjeran, C. [2004] “Finite difference approximations for fractional advectiondispersion,” J. Comput. Appl. Math. 172, 65–77. Meerschaert, M. M., Scheffer, H. P. & Tadjeran, C. [2006] “Finite difference methods for two-dimensional fractional dispersion equation,” J. Comput. Phys. 211, 249–261. Meerschaert, M. M. & Tadjeran, C. [2006] “Finite difference approximations for two-sided space-fractional partial differential equations,” Appl. Math. Model. 56, 80–90.
1230014-26
May 9, 2012
10:43
WSPC/S0218-1274
1230014
Finite Difference Methods for Fractional Differential Equations
Metzler, R. & Klafter, J. [2000] “The random walk’s guide to anomalous diffusion: A fractional dynamics approach,” Phys. Rep. 339, 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,” J. Phys. A: Math. Gen. 37, R161–R208. Murio, D. A. [2006] “On the stable numerical evaluation of Caputo fractional derivatives,” Comput. Math. Appl. 51, 1539–1550. Murio, D. A. [2008] “Implicit finite difference approximation for time fractional diffusion equatuions,” Comput. Math. Appl. 56, 1138–1145. Odibat, Z. M. [2006] “Approximations of fractional integrals and Caputo fractional derivatives,” Appl. Math. Comput. 178, 527–533. Odibat, Z. M. [2009] “Computational algorithms for computing the fractional derivatives of functions,” Math. Comput. Simul. 79, 2013–2020. Odibat, Z. & Momani, S. [2008a] “Numerical methods for nonlinear patial equations of fractional order,” Appl. Math. Model. 32, 28–39. Odibat, Z. & Momani, S. [2008b] “An algorithm for the numerical solution of differential equations of fractional order,” JAMI 26, 15–27. Oldham, K. & Spanier, J. [1974] The Fractional Calculus (Acdemic Press, NY). Podlubny, I. [1999] Fractional Differential Euations (Acdemic Press, San Diego). Podlubny, I. [2000] “Matrix approach to discrete fractional calculus,” Fract. Cal. Appl. Anal. 4, 359–386. Podlubny, I., Chechkin, A., Skovranek, T., Chen, Y. Q. & Blas M. Vinagre Jara [2009] “Matrix approach to discrete fractional calculus II: Partial fractional differential equations,” J. Comput. Phys. 228, 3137–3153. Raberto, M., Scalas, E. & Mainardi, F. [2002] “Waitingtimes and returns in high-frequency financial data: An empirical study,” Physica A 314, 749–755. Rihan, F. A. [2009] “Computational methods for delay parabolic and time-fractional partial differential equations,” Numer. Meth. Part. D. E. 26, 1556–1571. Roop, J. P. [2006] “Computational aspects of FEM approximation of fractional advection dispersion equations on bounded domains in R2 ,” J. Comput. Appl. Math. 193, 243–268. Roop, J. P. [2008] “Numerical approximation of a one-dimensional space fractional advection-dispersion equation with boundary layer,” Comput. Math. Appl. 56, 1808–1819. Samko, S. G., Kilbas, A. A. & Marichev, O. I. [1993] Fractional Integrals and Derivatives: Theory and Applications (Gordon & Breach Science Publishers in Switzerland, Philadelphia, Pa., USA).
Scherer, R., Kalla, S. L., Boyadjiev, L. & Al–Saqabi, B. [2008] “Numerical treatment of fractional heat equations,” Appl. Numer. Math. 58, 1212–1223. Schmidt, A. & Gaul, L. [2006] “On the numerical evaluation of fractional derivatives in multi-degree-offreedom systems,” Sign. Process. 86, 2592–2601. Shen, S. & Liu, F. [2005] “Error analysis of explicit finite difference approximation for the space fractional diffusion equation with insulated ends,” ANZIAM J. 46, C871–C887. Shen, S., Liu, F., Anh, V. & Turner, I. [2006] “Detailed analysis of a conservative difference approximation for the time fractional diffusion equation,” J. Appl. Math. Comput. 3, 1–19. Shen, S., Liu, F. & Anh, V. [2008a] “Fundamental solution and discrete random walk model for a time-space fractional diffusion equation of distributed order,” J. Appl. Math. Comput. 28, 147–164. Shen, S., Liu, F., Anh, V. & Turner, I. [2008b] “The fundamental solution and numerical solution of the Riesz fractional advection-dispersion equation,” IMA J. Appl. Math. 73, 850–872. Shen, S., Liu, F. & Anh, V. [2011] “Numerical approximations and solution techniques for the space-time Riesz–Caputo fractional advection-diffusion equation,” Numer. Algor. 56, 383–403. Sousa, E. [2009] “Finite difference approximations for a fractional advection diffusion problem,” J. Comput. Phys. 228, 4038–4054. Sousa, E. [2010] “How to approximate the fractional derivative of order 1 < α ≤ 2,” The 4th IFAC Workshop Fractional Differentiation and Its Applications (Badajoz, Spain), Article no. FDA10-019, 6 pages. Su, L. J., Wang, W. Q. & Yang, Z. X. [2009] “Finite difference approximations for the fractional advectiondiffusion equation,” Phys. Lett. A 373, 4405–4408. Su, L. J., Wang, W. Q. & Xu, Q. Y. [2010] “Finite difference methods for fractional dispersion equations,” Appl. Math. Comput. 216, 3329–3334. Sun, Z. Z. & Wu, X. N. [2006] “A fully discrete difference scheme for a diffusion-wave system,” Appl. Numer. Math. 56, 193–209. Sun, H. G., Chen, W., Li, C. P. & Chen, Y. Q. [2010] “Fractional differential models for anomalous diffusion,” Physica A 389, 2719–2724. Tadjeran, C., Meerschaert, M. M. & Scheffer, H. P. [2006] “A second-order accurate numerical approximation for the fractional diffusion equation,” J. Comput. Phys. 213, 205–213. Tadjeran, C. & Meerschaert, M. M. [2007] “A secondorder accurate numerical method for the twodimensional fractional diffusion equation,” J. Comput. Phys. 220, 813–823.
1230014-27
May 9, 2012
10:43
WSPC/S0218-1274
1230014
C. Li & F. Zeng
Val´erio, S. & Costa, J. [2011] “Variable-order fractional derivatives and their numerical approximations,” Sign. Process. 91, 470–483. Wang, H., Wang, K. X. & Sircar, T. [2010] “A direct O(N log2 N ) finite difference method for fractional diffusion equations,” J. Comput. Phys. 229, 8095– 8104. Wu, C. H. [2009] “Numerical solution for Stokes’ first problem for a heated generalized second grade fluid with fractional derivative,” Appl. Numer. Math. 59, 2571–2583. Wu, C. H. & Lu, L. Z. [2010] “Implicit numerical approximation scheme for the fractional Fokker–Planck equation,” Appl. Math. Comput. 216, 1945–1955. Yang, C. & Liu, F. [2006] “A computationally effective predictor-corrector method for simulating fractional order dynamical control system,” ANZIAM J. 47, C168–C184. Yang, Q. Q., Liu, F. & Turner, I. [2009a] “Computationally efficient numerical methods for time- and space-fractional Fokker–Planck equations,” Phys. Scr. 2009, 014026. Yang, Q. Q., Turner, I. & Liu, F. [2009b] “Analytical and numerical solutions for the time and space-symmetric fractional diffusion equation,” ANZIAM J. 50, C800– C814. Yang, Q. Q., Liu, F. & Turner, I. [2010a] “Stability and convergence of an effective numerical method for the time-space fractional Fokker–Planck equation with a nonlinear source term,” Int. J. Diff. Eqs. 2010, Article ID 464321, 22 pages. Yang, Q. Q., Liu, F. & Turner, I. [2010b] “Numerical methods for fractional partial differential equations with Riesz space fractional derivatives,” Appl. Math. Model. 34, 200–218. Yin, C., Liu, F. & Anh, V. [2007] “Numerical simulation of the nonlinear fractional dynamical systems with fractional damping for the extensible and inextensible pendulum,” J. Algor. Comput. Tech. 1, 427– 447. Yuste, S. B. & Acedo, L. [2005] “An explicit finite difference method and a new von Neumann-type stability
analysis for fractional diffusion equations,” SIAM J. Numer. Anal. 42, 1862–1874. Yuste, S. B. [2006] “Weighted average finite difference methods for fractional diffusion equations,” J. Comput. Phys. 216, 264–274. Zaslavsky, G. M. [2002] “Chaos, fractional kinetics, and anomalous transport,” Phys. Rep. 371, 461–580. Zhang, Y. [2009] “A finite difference method for fractional partial differential equation,” Appl. Math. Comput. 215, 524–529. Zheng, Y. Y., Li, C. P. & Zhao, Z. G. [2010a] “A note on the finite element method for the space-fractional advection diffusion equation,” Comput. Math. Appl. 59, 1718–1726. Zheng, Y. Y., Li, C. P. & Zhao, Z. G. [2010b] “A fully discrete DG method for nonlinear fractional Fokker– Planck equation,” Math. Probl. Eng. 2010, Article ID 279038, 26 pages. Zhou, T. S. & Li, C. P. [2005] “Synchronization in fractional-order differential systems,” Physica D 212, 111–125. Zhuang, P. & Liu, F. [2006] “Implicit difference approximation for the time fractional diffusion equation,” J. Appl. Math. Comput. 33, 87–99. Zhuang, P. & Liu, F. [2007a] “Implicit difference approximation for the two-dimensional space-time fractional diffusion equation,” J. Appl. Math. Comput. 25, 269– 282. Zhuang, P. & Liu, F. [2007b] “Finite difference approximation for two-dimensional time fractional diffusion equation,” J. Algor. Comput. Tech. 1, 1–15. Zhuang, P., Liu, F., Anh, V. & Turner, I. [2008] “New solution and analytical techniques of the implicit numerical method for the anomalous subdiffusion equation,” SIAM J. Numer. Anal. 46, 1079–1095. Zhuang, P., Liu, F., Anh, V. & Turner, I. [2009] “Numerical methods for the variable-order fractional advection-diffusion equation with a nonlinear source term,” SIAM J. Numer. Anal. 47, 1760–1781.
1230014-28
