- Email: [email protected]
Implicit numerical approximation scheme for the fractional Fokker–Planck equation
Applied Mathematics and Computation 216 (2010) 1945–1955
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Implicit numerical approximation scheme for the fractional Fokker–Planck equation q Chunhong Wu a,*, Linzhang Lu b a b
Department of Mathematics and Physics, Xiamen University of Technology, Xiamen 361024, PR China School of Mathematical Sciences, Xiamen University, Xiamen 361005, PR China
a r t i c l e
i n f o
a b s t r a c t In this paper, we consider an anomalous subdiffusion process, governed by fractional Fokker–Planck equation. An effective numerical method for approximating Fokker–Planck equation in a bounded domain is presented. The stability and convergence of the numerical method are analyzed. Some numerical examples are presented to show the application of the present technique. The numerical results exhibit the good performance of our theoretical analysis. Ó 2010 Elsevier Inc. All rights reserved.
Keywords: Fractional Fokker–Planck equation Subdiffusion Numerical method Stability Convergence
1. Introduction Fractional kinetic equations of the diffusion, diffusion–advection, and Fokker–Planck type have received much attention because of their practical importance [1,2]. Since Metzler and Klafter [2] firstly reviewed the fractional Fokker–Planck equation, a large number of research papers developing fractional dynamics further, or applying it to various systems have appeared, such as, Metzler and Klafter [3] discussed the theoretical modelling of sub- and superdiffusive processes and their applications, and introduced extensively the occurrence of anomalous dynamics in various fields ranging from nanoscale over biological to geophysical and environmental systems. Sokolov and Klafter [4] presented a continuous-time random walk model for approximating the subdiffusive processes, and derived two different but equivalent forms of kinetic equations, which reduce to known fractional diffusion or Fokker–Planck equations for waiting-time distributions following a power law. Yuste and Acedo [5] obtained an explicit finite difference method and a new von Neumann-type stability analysis for the fractional subdiffusion equation. Langlands and Henry [6] investigated the accuracy and stability of an implicit numerical scheme for solving the fractional diffusion equation. Zahran et al. [7] derived the counterpart of fractional Fokker–Planck equation using comb-like model as fractal background medium. The resulting equation which is the fractional Fokker–Planck equation essentially differs from the usual Fokker–Planck equation (i.e., instead of first-order derivative of time evolution operator in Fokker–Planck equation, the derivatives of fractional order b arises in fractional Fokker–Planck equation)
uðx; tÞ ¼ Db t LFP uðx; tÞ;
0 < b 6 1;
ð1Þ
where
LFP ¼
@ @ dðxÞ v ðxÞ ; @x @x
q Research supported by National Natural Science Foundation of PR China (10531080) and the Scientific Research Foundation of Xiamen University of Technology (YKJ09016R). * Corresponding author. E-mail addresses: [email protected] (C. Wu), [email protected] (L. Lu).
0096-3003/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2010.03.024
1946
C. Wu, L. Lu / Applied Mathematics and Computation 216 (2010) 1945–1955
and
Db t uðx; tÞ ¼
1 CðbÞ
Z
t
0
uðx; gÞ ðt gÞ1b
dg:
The above expression is defined as the Riemann–Liouville fractional integral operator [8], CðÞ is the gamma function, dðxÞ > 0 and v ðxÞ are the diffusion and drift coefficients related to the fractal system under consideration, respectively. In order to construct a suitable numerical method, we differentiate both sides of (1) and let v ðxÞ ¼ 0, then the following equation is obtained
@uðx; tÞ @ @ ¼ D1b dðxÞ uðx; tÞ ; t @t @x @x where the symbol D1b is the Riemann–Liouville fractional derivative operator, defined by t
Dt1b uðx; tÞ ¼
1 @ CðbÞ @t
Z 0
t
uðx; gÞ ðt gÞ1b
dg:
Without loss of generality, a source term is added, then we consider the fractional Fokker–Planck equation (FFPE) without the drift term
@uðx; tÞ @ @ ¼ D1b dðxÞ uðx; tÞ þ f ðx; tÞ ; t @t @x @x
0 6 x 6 a; 0 < t 6 T
ð2Þ
with boundary conditions
uð0; tÞ ¼ u1 ðtÞ;
0 6 t 6 T;
uða; tÞ ¼ u2 ðtÞ;
0 6 t 6 T;
ð3Þ
and initial condition
uðx; 0Þ ¼ xðxÞ;
0 6 x 6 a;
ð4Þ
where 0 < b 6 1. When the fractional differential equation is derived from various fields, it is urgent to research that how to solve it. Some researchers are interested in the analytical solutions, for example, Mainardi et al. [9] discussed the Cauchy problem for the space–time fractional diffusion-wave equation, investigated the fundamental solution by Fourier–Laplace representation and provided a general representation of the solution in terms of Mellin-Barnes integrals in the complex plane by using the Mellin transform. Qi and Xu [10] studied the unsteady flow of viscoelastic fluid with the fractional derivative Maxwell model (FDMM) in a channel, the exact solutions were obtained for an arbitrary pressure gradient by means of the finite Fourier cosine transform and the Laplace transform. Wang and Xu [11] considered the unsteady axial Couette flow of fractional second grade fluid (FSGF) and fractional Maxwell fluid (FMF) between two infinitely long concentric circular cylinders, and obtained the analytical solutions by using the theory of Laplace transform and Weber transform. As analytic solutions of most fractional differential equations are not usually expressed explicitly, many authors resort to numerical solution strategies based on convergence and stability analysis. Liu et al. [12,13] considered the space fractional partial differential equation and transformed the equation into a system of ordinary differential equations that was then solved using backward differentiation formulas. Meerschaert and Tadjeran [14] presented practical numerical methods to solve the one-dimensional fractional advection–dispersion equations with variable coefficients on a finite domain. Shen and Liu [15] proposed an explicit finite difference approximation for the space fractional diffusion equation and gave an error analysis. Meerschaert et al. [16] developed practical numerical methods to solve two-dimensional fractional dispersion equations with variable coefficients on a finite domain and obtained the first-order accuracy in space and time. And they proposed the stability analysis by using matrix analysis method. Roop [17] investigated the computational aspects for the Galerkin approximation using continuous piecewise polynomial basis functions on a regular triangulation of a bounded domain in R2 . Liu et al. [18] discussed an approximation of the Lévy–Feller advection–dispersion process by a random walk and finite difference method. Chen et al. [19,20] proposed a Fourier method for the fractional diffusion equation describing subdiffusion and the Galilei invariant fractional advection diffusion equation. Zhuang et al. [21] proposed a new implicit numerical method and two solution techniques for improving its convergence order for solving an anomalous subdiffusion equation. The rest of this paper is organized as follows. In Section 2, we construct an effective implicit numerical approximation scheme for the FFPE (2) with boundary and initial conditions (3) and (4). In Sections 3 and 4, the stability and convergence analyses are discussed. Finally, in Section 5, some numerical results are given to testify our method and theoretical analyses.
2. An implicit numerical approximation scheme for the FFPE In this section, an implicit numerical approximation scheme for the FFPE (2) with boundary and initial conditions (3), (4) is proposed.
C. Wu, L. Lu / Applied Mathematics and Computation 216 (2010) 1945–1955
1947
Let K ¼ ½0; a ½0; T, then define the function space associated with a special function d(x)
(
@
@
@4
@
)
dðxÞ uðx; tÞ 2 C 2 ðKÞ; 3 UðKÞ ¼ uðx; tÞ dðxÞ uðx; tÞ 2 CðKÞ : @x @x @x @x @t And we suppose the FFPE (2)–(4) has a solution uðx; tÞ 2 Ud ðKÞ for the special function d(x) in (2). We now discretize space and time by grid points and time instants as follows:
a ; m T k ¼ 0; 1; 2; . . . ; K; s ¼ ; K
xi ¼ ih;
i ¼ 0; 1; 2; . . . ; m; h ¼
t k ¼ k s; where h and
xiþ1 2
s are the space and time steps, respectively. And we introduce the dual subdivision points
1 h; ¼ iþ 2
i ¼ 0; 1; 2; . . . ; m 1:
Then, we can discretize the variable uðxi ; tk Þ. For convenience, we denote
@ @u @ @uðx; tÞ d ðx; tÞ ¼ dðxÞ ; @x @x @x @x h h 2 d uðx; tÞ ¼ d x þ ; t ðuðx þ h; tÞ uðx; tÞÞ d x ; t ðuðx; tÞ uðx h; tÞÞ: d x 2 2 Luðx; tÞ ¼
Integrating both sides of (2) from t k to t kþ1 , we obtain
Z tk 1 Luðxi ; gÞ þ f ðxi ; gÞ dg CðbÞ 0 ðtkþ1 gÞ1b ðt k gÞ1b 0 Z s Z tk 1 Luðxi ; gÞ þ f ðxi ; gÞ 1 Lv ðxi ; gÞ þ gðxi ; gÞ ¼ uðxi ; tk Þ þ dg þ dg; 1b CðbÞ 0 C ðbÞ ðt kþ1 gÞ ðt k gÞ1b 0
uðxi ; tkþ1 Þ ¼ uðxi ; tk Þ þ
1 CðbÞ
Z
t kþ1
Luðxi ; gÞ þ f ðxi ; gÞ
dg
ð5Þ
where v ðx; tÞ ¼ uðx; t þ sÞ uðx; tÞ; gðx; tÞ ¼ f ðx; t þ sÞ f ðx; tÞ. We can rewrite (5) as
uðxi ; tkþ1 Þ ¼ uðxi ; tk Þ þ I1 þ I2 ; where
1 CðbÞ
I1 ¼
1 CðbÞ
I2 ¼
Z s Luðxi ; gÞ þ f ðxi ; gÞ ðtkþ1 gÞ1b
0
Z
tk
dg;
Lv ðxi ; gÞ þ gðxi ; gÞ ðtk gÞ1b
0
dg:
For convenience, we introduce the notation b
bk ¼ ðk þ 1Þb k ;
k ¼ 0; 1; 2; . . . ; K:
ð6Þ
For I1 , we can use the approximation
1 CðbÞ
Z s Luðxi ; sÞ þ f ðxi ; sÞ
sb
b
ððk þ 1Þb k Þ þ R11 ðtkþ1 gÞ sb 1 sb 1 ¼ bk 2 d d2x uðxi ; sÞ þ R12 þ f ðxi ; sÞ þ R11 ¼ bk 2 d d2x uðxi ; sÞ þ f ðxi ; sÞ þ R1 ; Cð1 þ bÞ h Cð1 þ bÞ h
I1 ¼
0
1b
dg þ R11 ¼ ½Luðxi ; sÞ þ f ðxi ; sÞ
where
Z s 1 Luðxi ; gÞ Luðxi ; sÞ þ f ðxi ; gÞ f ðxi ; sÞ dg; CðbÞ 0 ðt kþ1 gÞ1b @ @u 1 ¼ d ðxi ; sÞ 2 d d2x uðxi ; sÞ; @x @x h
R11 ¼ R12
R1 ¼ Noting that
sb Cð1 þ bÞ
bk R12 þ R11 :
Cð1 þ bÞ
1948
C. Wu, L. Lu / Applied Mathematics and Computation 216 (2010) 1945–1955
" # @ @u @ @u @2 @u @f ðxi ; g1 Þ ðg sÞ; d ðxi ; gÞ þ f ðxi ; gÞ ¼ d ðxi ; sÞ þ f ðxi ; sÞ þ ðxi ; g1 Þ þ d @x @x @x @x @x @t @x @t
Luðxi ; gÞ þ f ðxi ; gÞ ¼
where 0 6 g 6 g1 6 s, we obtain
Z s
1 CðbÞ
jR11 j 6 C 1 s
1 1b
ðt kþ1 gÞ
0
dg 6
C 1 s1þb
Cð1 þ bÞ
bk :
By Taylor’s formula it is apparent that 2
jR12 j 6 Ch : 2
Hence, we have jR1 j 6 Cbk sb ðs þ h Þ. For I2 , we can use the approximation
I2 ¼ ¼ ¼
Z
1 CðbÞ
Lv ðxi ; gÞ þ gðxi ; gÞ
tk
ðt k gÞ
0
Z
k1 X
1 CðbÞ
sb
k1 X
Cð1 þ bÞ
s¼0
dg ¼
k1 Z tsþ1 1 X Lv ðxi ; gÞ þ gðxi ; gÞ dg CðbÞ s¼0 ts ðt k gÞ1b
Lv ðxi ; t sþ1 Þ þ gðxi ; tsþ1 Þ
tsþ1
ðtk gÞ1b
ts
s¼0
1b
dg þ R21 ¼
sb
k1 X
Cð1 þ bÞ
s¼0
bks1 ½Lv ðxi ; t sþ1 Þ þ gðxi ; tsþ1 Þ þ R21
bks1 ½Lh v ðxi ; t sþ1 Þ þ gðxi ; tsþ1 Þ þ R22 þ R21 ;
where
R21 ¼ R22 ¼
k1 Z t sþ1 1 X Lv ðxi ; gÞ Lv ðxi ; tsþ1 Þ þ gðxi ; gÞ gðxi ; tsþ1 Þ dg; CðbÞ s¼0 ts ðtk gÞ1b
sb
k1 X
Cð1 þ bÞ
s¼0
Lh uðxi ; ts Þ ¼
bks1 ½Lv ðxi ; t sþ1 Þ Lh v ðxi ; t sþ1 Þ;
1
d2 uðxi ; ts Þ: 2d x
h
Noting that
Lv ðxi ; gÞ þ gðxi ; gÞ ¼
" # @ @v @2 @v @gðxi ; fÞ ðg tsþ1 Þ; d ðxi ; t sþ1 Þ þ gðxi ; t sþ1 Þ þ d ðxi ; fÞ þ @x @t @x @x @t @x
where t s 6 g 6 f 6 t sþ1 . And,
@2 @v @2 @u @2 @u @3 @u ðxi ; f þ sÞ ðxi ; fÞ ¼ ðxi ; f1 Þs; ðxi ; fÞ ¼ d d d d 2 @x @x @x @x @t @x @x @t @x @t @x @t @gðxi ; fÞ @f ðxi ; f þ sÞ @f ðxi ; fÞ @ 2 f ðxi ; f2 Þ s; ¼ ¼ @t @t @t @t 2 where f 6 f1 ; f2 6 f þ s. Consequently,
jLv ðxi ; gÞ Lv ðxi ; t sþ1 Þj 6 C 3 s2 : So we obtain
jR21 j 6 C 3 s2
1 CðbÞ
Z 0
tk
dg ðt k gÞ1b
6 C s2 :
From Taylor’s formula, it is apparent that
2 2 @ @v 1 h @3 @v h @3 @v ðxi ; tsþ1 Þ ¼ 2 d d2x uðxi ; tsþ1 Þ ðn; t ðn; t sþ1 Þ d Þ ¼ L v ðx ; t Þ d d sþ1 sþ1 i h @x @x 24 @x3 @x 24 @x3 @x h " # 2 2 h @3 @u @3 @u h s @4 @u ðn; tsþ2 Þ 3 d ðn; tsþ1 Þ ¼ Lh v ðxi ; t sþ1 Þ ðn; g2 Þ: ¼ Lh v ðxi ; tsþ1 Þ d d 3 3 @x @x @x 24 @x @x 24 @x @t
Lv ðxi ; t sþ1 Þ ¼
Thus, we have 2
jR22 j 6 C sh :
1949
C. Wu, L. Lu / Applied Mathematics and Computation 216 (2010) 1945–1955
From the above results, we obtain
uðxi ; tkþ1 Þ ¼ uðxi ; tk Þ þ r 1 bkd d2x uðxi ; sÞ þ r 1
k1 X
bks1d d2x ðuðxi ; t sþ2 Þ uðxi ; tsþ1 ÞÞ þ r 2 bk f ðxi ; sÞ þ r 2
s¼0
k1 X
bks1 ðf ðxi ; tsþ2 Þ
s¼0
; f ðxi ; t sþ1 ÞÞ þ Rkþ1 i b
b
s s where r1 ¼ Cð1þbÞh 2 ; r 2 ¼ Cð1þbÞ, and 2
jRkþ1 j 6 Cðbk sb þ sÞðs þ h Þ: i
ð7Þ
To analyze the stability, we find it worthwhile to recall here the following useful lemma associated with the property of the coefficients bk ðk ¼ 0; 1; 2; . . .Þ referred to in [21]. Lemma 2.1. The coefficients bk ðk ¼ 0; 1; 2; . . .Þ defined by (6) satisfy (1) b0 ¼ 1; bk > 0; k ¼ 0; 1; 2; . . .; (2) bk > bkþ1 ; k ¼ 0; 1; 2; . . .; (3) there exists a positive constant C > 0 such that
s 6 Cbk sb ; k ¼ 1; 2; . . . Let Rk ¼ ðRk1 ; Rk2 ; . . . ; Rkm1 ÞT . By using Lemma 2.1 and (7), we obtain the following lemma. qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Pm1 k 2 d Lemma 2.2. Suppose that kRk k2 ¼ i¼1 jRi j h, if uðx; tÞ 2 U ðKÞ is the solution of (2)–(4) , then we have 2
kRkþ1 k2 6 Cbk sb ðs þ h Þ: Let uki ; fik be the numerical approximation to uðxi ; t k Þ; f ðxi ; tk Þ; respectively, and introduce the following notation: 2 k d dx ui
¼ diþ1 ðukiþ1 uki Þ di1 ðuki uki1 Þ; 2
Dx uki ¼ ukiþ1 uki :
2
We obtain the following implicit numerical approximation scheme (INAS):
ukþ1 ¼ uki þ r 1 bkd d2x u1i þ r1 i
k1 X
bks1d d2x ðusþ2 usþ1 Þ þ r2 bk fi1 þ r 2 i i
s¼0
k1 X
bks1 ðfisþ2 fisþ1 Þ;
s¼0
which can be rewritten as
ukþ1 ¼ uki þ r 1d d2x ukþ1 þ r1 i i
k1 k1 X X ðbsþ1 bs Þd d2x uks þ r 2 fikþ1 þ r2 ðbsþ1 bs Þfiks ; i s¼0
ð8Þ
s¼0
for i ¼ 1; 2; . . . ; m 1; k ¼ 0; 1; 2; . . . ; K 1. The boundary and initial conditions are
uk0 ¼ u1 ðksÞ; ukm ¼ u0i ¼
k ¼ 1; 2; . . . ; K;
u2 ðksÞ; k ¼ 1; 2; . . . ; K; xðihÞ; i ¼ 0; 1; . . . ; m:
ð9Þ
Then we rewrite Eqs. (8) and (9) as the following matrix form:
8 1 0 1 1 > < Au ¼ u þ r 1 w þ r 2 f ; kþ1 k kþ1 Au ¼ u þ r1 v þ r 1 wkþ1 þ r2 f kþ1 þ r 2 g kþ1 ; > : 0 u ¼ ;
ð10Þ
k > 0;
where
2
uk1
3
6 k 7 6 u2 7 7 6 7 6 uk ¼ 6 ... 7; 7 6 6 k 7 4 um2 5 ukm1 and
v kþ1 ¼
Pk1
s¼0 ðbsþ1
2
f1k
3
6 k 7 6 f2 7 7 6 7 6 f k ¼ 6 ... 7; 7 6 6 k 7 4 fm2 5
2
wkþ1
6 6 6 6 ¼6 6 6 4
k fm1
bs Þd d2x uks ; g kþ1 ¼
d1 ukþ1 0 2
0 .. . 0 kþ1 dm1 um1
3 7 7 7 7 7; 7 7 5
s¼0 ðbsþ1
x1 7 6 6 x2 7
7 6 6 . 7 ¼ 6 .. 7; 7 6 7 6 4 xm2 5
2
Pk1
3
2
xm1
bs Þf ks ; xi ¼ xðihÞ.
The matrix A in (10) can be written as the following form:
1950
C. Wu, L. Lu / Applied Mathematics and Computation 216 (2010) 1945–1955
2 6 6 6 6 6 6 6 6 6 4
1 þ r1 d3 þ d1 2
2
r 1 d3
r1 d3 2 1 þ r 1 d5 þ d3
..
.
..
..
.
2
2
3
0
2
.
0
r1 dm3 2 1 þ r1 dm1 þ dm3
r1 dm3 2
2
7 7 7 7 7 7: 7 7 7 5
ð11Þ
2
From the structure of the matrix A we can see that it is a strictly diagonally dominant symmetric tridiagonal matrix with positive diagonal terms and nonpositive off-diagonal terms. Then, we obtain the following theorem. Theorem 2.1. The discretization matrix A is invertible. And the system (8) and (9) has a unique solution.
3. Stability of the INAS For u ¼ ðu1 ; u2 ; . . . ; um1 ÞT ;
ðu; v Þ ¼
m 1 X
ui v i h;
v ¼ ðv 1 ; v 2 ; . . . ; v m1 ÞT , we define the inner product
ðu; v Þd ¼
m1 X
i¼1
i¼1
diþ1 ui v i h; 2
and the norm
kuk2 ¼
m1 X
pffiffiffiffiffiffiffiffiffiffiffiffi ðu; uÞ ¼
!12 u2i h
;
kukd ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðu; uÞd ¼
i¼1
m1 X i¼1
!12 diþ1 u2i h 2
:
~ ki ði ¼ 0; 1; 2; . . . ; m; k ¼ 0; 1; 2; . . . ; KÞ is the approximation solution of (8) and (9). Then the error We suppose that u eki ¼ ue ki uki satisfies
eikþ1 ¼ eki þ r1d d2x ekþ1 þ r1 i
k1 X ðbsþ1 bs Þd d2x eiks ;
ð12Þ
s¼0
and
ek0 ¼ ekm ¼ 0; k ¼ 1; 2; . . . ; K: k
Let E ¼ ðe
k 1;
e
k 2; . . . ;
e
T k m1 Þ :
Multiplying (12) by
ð13Þ
e
kþ1 h i
kEkþ1 k22 ¼ ðEkþ1 ; Ek Þ þ r1 ðd d2x Ekþ1 ; Ekþ1 Þ þ r1
k1 X
and summing i from 1 to m 1, we have
ðbsþ1 bs Þðd d2x Eks ; Ekþ1 Þ
s¼0 k1 X ks kþ1 ðbs bsþ1 Þ d1 eks ; Dx Ekþ1 Þd : ¼ ðEkþ1 ; Ek Þ r1 d1 ðe1kþ1 Þ2 h þ ðDx Ekþ1 ; Dx Ekþ1 Þd þ r 1 1 e1 h þ ðDx E 2
2
s¼0
By applying the Cauchy–Schwarz inequality, we have
kEkþ1 k22 6
k1 r X 1 1 Þ2 h þ kDx Ekþ1 k2d þ ðbs bsþ1 Þ d1 ðe1ks Þ2 h þ d1 ðe1kþ1 Þ2 h ðkEkþ1 k22 þ kEk k22 Þ r1 d1 ðekþ1 1 2 2 2 2 2 s¼0
þ Noting that
Pk1
k1 r1 X ðbs bsþ1 ÞðkDx Eks k2d þ kDx Ekþ1 k2d Þ: 2 s¼0
s¼0 ðbs
kEkþ1 k22 6
bsþ1 Þ ¼ 1 bk and bk > 0, we have
k1 r X 1 r1 1 Þ2 h þ kDx Ekþ1 k2d þ ðbs bsþ1 Þ d1 ðeks Þ2 h þ kDx Eks k2d ; ðkEkþ1 k22 þ kEk k22 Þ ð1 þ bk Þ d1 ðekþ1 1 1 2 2 2 2 2 s¼0
i.e.,
kEkþ1 k22 þ r 1
k X
k1 X 2 ks 2 bs d1 ðe1kþ1s Þ2 h þ kDx Ekþ1s k2d 6 kEk k22 þ r1 bs d1 ðeks kd : 1 Þ h þ kDx E 2
s¼0
s¼0
Defining the energy norm
kEk k2E ¼ kEk k22 þ r1
k1 X s¼0
2 ks 2 bs d1 ðeks kd ; 1 Þ h þ kDx E 2
2
ð14Þ
C. Wu, L. Lu / Applied Mathematics and Computation 216 (2010) 1945–1955
1951
we have
kEkþ1 k22 6 kEkþ1 k2E 6 kEk k2E 6 6 kE1 k2E : As
e1i ¼ e0i þ r1 d d2x e1i ; we have 1 kE1 k22 ¼ ðE1 ; E0 Þ r 1 d1 ðe11 Þ2 h þ kDx E1 k2d 6 ðkE1 k22 þ kE0 k22 Þ r 1 d1 ðe11 Þ2 h þ kDx E1 k2d : 2 2 2
Then,
kE1 k2E ¼ kE1 k22 þ r1 d1 ðe11 Þ2 h þ kDx E1 k2d 6 kE0 k22 ; 2
ð15Þ
so that kEkþ1 k22 6 kE0 k22 : Hence, we can obtain the following theorem of stability. Theorem 3.1. The implicit numerical approximation scheme (8) is unconditionally stable.
4. Convergence of the INAS In this section, we discuss the convergence of the INAS. Suppose that u(x, t) is the solution of the FFPE (2)–(4) and uðx; tÞ 2 Ud ðKÞ. Let uðxi ; tk Þ ði ¼ 0; 1; 2; . . . ; m; k ¼ 0; 1; 2; . . . ; KÞ be the exact solution of (2)–(4) at the mesh point ðxi ; t k Þ. Define
yki ¼ uðxi ; tk Þ uki ; k
and Y ¼
i ¼ 0; 1; 2; . . . ; m; k ¼ 0; 1; 2; . . . ; K;
ðyk1 ; yk2 ; . . . ; ykm1 ÞT .
By using uki ¼ uðxi ; tk Þ yki , substitution into (8) leads to
ykþ1 ¼ yki þ r 1d d2x ykþ1 þ r1 i i
k1 X ðbsþ1 bs Þd d2x yiks þ Rkþ1 ; i
ð16Þ
s¼0
where i ¼ 1; 2; . . . ; m 1; k ¼ 0; 1; . . . ; K 1, and
y0i ¼ 0; i ¼ 0; 1; . . . ; m;
yk0 ¼ ykm ¼ 0; k ¼ 0; 1; . . . ; K:
Multiplying (16) by ykþ1 h and summing i from 1 to m 1, we obtain i
kY kþ1 k22 ¼ ðY kþ1 ; Y k Þ þ r 1 ðd d2x Y kþ1 ; Y kþ1 Þ þ r1
k1 X ðbsþ1 bs Þðd d2x Y ks ; Y kþ1 Þ þ ðRkþ1 ; Y kþ1 Þ:
ð17Þ
s¼0
For s ¼ 0; 1; . . . ; k þ 1; we have s kþ1 ðd d2x Y s ; Y kþ1 Þ ¼ d1 ys1 ykþ1 Þd ; 1 h ðDx Y ; Dx Y 2
ðY
kþ1
1 ; Y Þ 6 ðkY kþ1 k22 þ kY s k22 Þ; 2 s
and s ykþ1 1 y1 6
1 kþ1 2 ðjy j þ jys1 j2 Þ: 2 1
Taking
c ¼ min diþ1 ; i
ð18Þ
2
and using jv wj 6 rv 2 þ 41r w2 ;
r > 0; we obtain
2
jðRkþ1 ; Y kþ1 Þj 6
2r 1 h bk c kþ1 2 a2 kY k2 þ kRkþ1 k22 : 2 2 a 8r 1 h bk c
ð19Þ
Similarly to the proof of the stability, we obtain
kY kþ1 k22 6
k1 r X 1 r1 1 Þ2 h þ kDx Y kþ1 k2d þ ðbs bsþ1 Þ d1 ðyks Þ2 h þ kDx Y ks k2d ðkY kþ1 k22 þ kY k k22 Þ ð1 þ bk Þ d1 ðykþ1 1 1 2 2 2 2 2 s¼0 2
þ
2r 1 h bk c kþ1 2 a2 kRkþ1 k22 : kY k2 þ 2 a2 8r 1 h bk c
Let jyki0 j ¼ max16i6m1 jyki j, then
ð20Þ
1952
C. Wu, L. Lu / Applied Mathematics and Computation 216 (2010) 1945–1955
yki0 ¼ yk1 þ
i0 1 X
Dx yki ;
yki0 ¼
m1 X
Dx yki :
i¼i0
i¼1
Thus
2jyki0 j 6 jyk1 j þ
m1 X
jDx yki j:
i¼1
Using the Cauchy–Schwarz inequality, we have
jyk1 j2
m 1 X
4jyki0 j2
6m
k
max16i6m1 jyki j,
þ
! jDx yki j2
h
i¼1
Let kY k1 ¼
kY k k21
a
6
4h
ðjyk1 j2 h 2
a
6
2
ðjyk1 j2 h þ kDx Y k k22 Þ:
Therefore,
þ kDx Y k k22 Þ:
Then we have
kY k k22 6 akY k k21 6
a2 4h
2
ðjyk1 j2 h þ kDx Y k k22 Þ:
Noting (18), we have 2
2
2r 1 h bk c kþ1 2 2r 1 h bk c a2 r1 2 kþ1 2 kY k2 6 ðjykþ1 k2 Þ ¼ bk ðcjy1kþ1 j2 h þ ckDx Y kþ1 k22 Þ 1 j h þ kDx Y 2 a2 a2 2 4h r1 kþ1 2 kþ1 2 6 bk d1 jy1 j h þ kDx Y kd : 2 2
ð21Þ
Hence, from (20) and (21), we have
kY kþ1 k22 6
k1 r X 1 r1 1 2 2 kþ1 2 ks 2 kd þ ðbs bsþ1 Þ d1 ðyks kd ðkY kþ1 k22 þ kY k k22 Þ d1 ðykþ1 1 Þ h þ kDx Y 1 Þ h þ kDx Y 2 2 2 2 2 s¼0
a2
þ Let ck ¼ kY k k22 þ r1
2
8r1 h bk c
kRkþ1 k22 :
Pk1
ks 2 s¼0 bs ðd12 ðy1 Þ h
þ kDx Y ks k2d Þ and then
ckþ1 6 ck þ Cbk sb ðs þ h2 Þ2 ; c0 ¼ 0; Hence, we obtain
ckþ1 6 C
k X
2
bs sb ðs þ h Þ2 :
s¼0
Noting that
Pk
s¼0 bs
sb ¼ ðk þ 1Þb sb 6 T b and kY kþ1 k22 6 ckþ1 ; we have 2
kY kþ1 k22 6 CT b ðs þ h Þ2 : Consequently, the following theorem of convergence is obtained. Theorem 4.1. Let uðx; tÞ 2 Ud ðKÞ be the solution of the FFPE (2)–(4) . Then the INAS (8) is convergent, and there exists a positive constant C > 0 such that 2
kY kþ1 k2 6 Cðs þ h Þ;
k ¼ 0; 1; 2; . . . ; K 1:
5. Numerical examples In this section, two numerical examples are presented to support our theoretical analysis. Example 1. Consider the following FFPE
@uðx; tÞ @ @ ¼ D1b dðxÞ uðx; tÞ t @t @x @x with initial and boundary conditions
ð22Þ
1953
C. Wu, L. Lu / Applied Mathematics and Computation 216 (2010) 1945–1955
( uðx; 0Þ ¼
x;
x 2 ½0; 1;
3x ; 2
x 2 ½1; 3;
uð0; tÞ ¼ uð3; tÞ ¼ 0;
t 2 ½0; 1:
The INAS is used to solve Eq. (22). Firstly, We take the diffusion coefficient dðxÞ ¼ 1. b ¼ 0:7; 0 6 x 6 3; 0 6 t 6 1 and pffiffiffi t ¼ 0:5; 0 6 x 6 3; 0:1 6 b 6 0:9, then taking the diffusion coefficient dðxÞ ¼ 3 x; b ¼ 0:7; 0 6 x 6 3; 0 6 t 6 1; The numerical simulation is shown in Figs. 1–3, respectively. They show the process exhibits subdiffusion behaviors. In order to show the approximation order of the INAS, we construct an example with an analytic solution. Example 2. The following FFPE with analytic solution is considered
@uðx; tÞ @ @ ¼ D1b dðxÞ uðx; tÞ þ f ðx; tÞ ; t @t @x @x
0 < x < 1; 0 < t < 1;
ð23Þ
where dðxÞ ¼ ex ; f ðx; tÞ ¼ Cðb þ 3Þex t 2 =2 2e2x tbþ2 . with the initial and boundary conditions
uð0; tÞ ¼ t 2þb ;
0 < t < 1;
uð1; tÞ ¼ et 2þb ; uðx; 0Þ ¼ 0;
0 < t < 1;
0 6 x 6 1:
The exact solution of (23) is uðx; tÞ ¼ ex t2þb :
u(x,t)
1.5 1
0.5 0 1
0.8
0.6
t
0.4
0.2
0 0
1
2
3
x
Fig. 1. Numerical approximation for b ¼ 0:7; dðxÞ ¼ 1.
0.5
u(x,t=0.5)
0.4 0.3 0.2 0.1 0 1
0.8
0.6
β
0.4
0.2
0 0
1
x
2
Fig. 2. Numerical approximation for different b at t ¼ 0:5.
3
1954
C. Wu, L. Lu / Applied Mathematics and Computation 216 (2010) 1945–1955
1.5
u(x,t)
1 0.5 3
0 1
0.8
2
0.6 t
0.4
0.2
1
x
0 0
Fig. 3. Numerical approximation for b ¼ 0:7; dðxÞ ¼
ffiffiffi p 3 x:
Table 1 The maximum absolute error of the INAS.
s¼h
b ¼ 0:4
b ¼ 0:7
b ¼ 0:9
1 50 1 100 1 200
3.0013994E003
4.6906049E003
6.0741258E003
1.5066619E003
2.3543905E003
3.0419634E003
5.9997852E004
1.2317856E003
1.5027906E003
0.5 The exact solution The numerical solution
0.45
u(x,t=0.5)
0.4 0.35 0.3 0.25 0.2 0.15 0
0.2
0.4
x
0.6
0.8
1
Fig. 4. The exact solution and numerical solution when b ¼ 0:5; t ¼ 0:5.
The maximum absolute error between the exact solution and the numerical solutions by INAS, with spatial and temporal steps s ¼ h ¼ 1=50; 1=100; 1=200 at time t ¼ 1:0, is listed in Table 1. In Table 1 the last column is computed as the ratio of the maximum absolute error at the previous grid size to the current grid size, which shows the order of convergence of the method as all step sizes are halved. This is in good agreement with the analysis of theory. We also compare the exact solution with the numerical solution by using INAS with different time and space steps at t ¼ 0:5 and b ¼ 0:5 which is shown in Fig. 4, we can see that the numerical solution is in good agreement with the exact solution. 6. Conclusions In this paper, an implicit numerical approximation scheme for the fractional Fokker–Planck equation has been described and demonstrated. The stability, consistency and convergence of the INAS for the FFPE have been discussed. These methods and analytical techniques can also be extended to some high-dimensional fractional partial differential equations.
C. Wu, L. Lu / Applied Mathematics and Computation 216 (2010) 1945–1955
1955
Acknowledgement The authors thank the referee for many thoughtful and constructive suggestions to improve the paper. References [1] Walter Wyss, The fractional diffusion equation, J. Math. Phys. 27 (1986) 2782–5785. [2] Ralf Metzler, Joseph Klafter, The random walk’s guide to anomalous diffusion: a fractional dynamics approach, Phys. Rep. 339 (2000) 1–77. [3] Ralf Metzler, Joseph Klafter, The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics, J. Phys. A 37 (2004) R161–208. [4] I.M. Sokolov, J. Klafter, From diffusion to anomalous diffusion: a century after Einstein’s Brownian motion, Chaos 15 (2005) 026103. [5] S.B. Yuste, L. Acedo, An explicit finite difference method and a new von Neumann-type stability analysis for fractional diffusion equations, SIAM J. Numer. Anal. 42 (2004) 1862–1874. [6] T.A.M. Langlands, B.I. Henry, The accuracy and stability of an implicit solution method for the fractional diffusion equation, J. Comput. Phys. 205 (2005) 719–736. [7] M.A. Zahran, E.M. Abulwafa, S.A. Elwakil, The fractional Fokker–Planck equation on comb-like model, Physica A 323 (2003) 237–248. [8] I. Podlubny, Fractional Differential Equations, Academic Press, 1999. [9] F. Mainardi, Y. Luchko, G. Pagnini, The fundamental solution of the space-time fractional diffusion equation, Fract. Calculus Appl. Anal. 4 (2) (2001) 153–192. [10] H. Qi, M. Xu, Unsteady flow of viscoelastic fluid with fractional Maxwell model in a channel, Mech. Res. Commun. 34 (2007) 210–212. [11] S. Wang, M. Xu, Axial Couette flow of two kinds of fractional viscoelastic fluids in an annulus, Nonlinear Anal.: Real World Appl. (2008). [12] F. Liu, V. Anh, I. Turner, Numerical solution of the space fractional Fokker–Planck equation, J. Comput. Appl. Math. 166 (2004) 209–219. [13] F. Liu, V. Anh, I. Turner, P. Zhuang, Numerical solution for the solute transport in fractal porous media, ANZIAM J. 45 (E) (2004) 461–473. [14] M.M. Meerschaert, C. Tadjeran, Finite difference approximations for fractional advection–dispersion flow equations, J. Comput. Appl. Math. 172 (1) (2004) 65–77. [15] S. Shen, F. Liu, Error analysis of an explicit finite difference approximation for the space fractional diffusion, ANZIAM J. 46 (E) (2005) 871–887. [16] M.M. Meerschaert, H. Scheffler, C. Tadjeran, Finite difference methods for two-dimensional fractional dispersion equation, J. Comput. Phys. 211 (2006) 249–261. [17] J.P. Roop, Computational aspects of FEM approximation of fractional advection dispersion equation on bounded domains in R2 , J. Comput. Appl. Math. 193 (1) (2006) 243–268. [18] Q. Liu, F. Liu, I. Turner, V. Anh, Approximation of the Lévy–Feller advection–dispersion process by random walk and finite difference method, J. Comput. Phys. 222 (2007) 57–70. [19] Chang-Ming Chen, F. Liu, I. Turner, V. Anh, A Fourier method for the fractional diffusion equation describing sub-diffusion, J. Comput. Phys. 227 (2007) 886–897. [20] Chang-Ming Chen, F. Liu, I. Turner, V. Anh, A new Fourier analysis method for the Galilei invariant fractional advection diffusion equation, ANZIAM J. 48 (2007) 605–619. [21] P. Zhuang, F. Liu, V. Anh, I. Turner, New solution and analytical techniques of the implicit numerical method for the anomalous subdiffusion equation, SIAM J. Numer. Anal. 46 (2) (2008) 1079–1095.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Implicit numerical approximation scheme for the fractional Fokker–Planck equation q Chunhong Wu a,*, Linzhang Lu b a b
Department of Mathematics and Physics, Xiamen University of Technology, Xiamen 361024, PR China School of Mathematical Sciences, Xiamen University, Xiamen 361005, PR China
a r t i c l e
i n f o
a b s t r a c t In this paper, we consider an anomalous subdiffusion process, governed by fractional Fokker–Planck equation. An effective numerical method for approximating Fokker–Planck equation in a bounded domain is presented. The stability and convergence of the numerical method are analyzed. Some numerical examples are presented to show the application of the present technique. The numerical results exhibit the good performance of our theoretical analysis. Ó 2010 Elsevier Inc. All rights reserved.
Keywords: Fractional Fokker–Planck equation Subdiffusion Numerical method Stability Convergence
1. Introduction Fractional kinetic equations of the diffusion, diffusion–advection, and Fokker–Planck type have received much attention because of their practical importance [1,2]. Since Metzler and Klafter [2] firstly reviewed the fractional Fokker–Planck equation, a large number of research papers developing fractional dynamics further, or applying it to various systems have appeared, such as, Metzler and Klafter [3] discussed the theoretical modelling of sub- and superdiffusive processes and their applications, and introduced extensively the occurrence of anomalous dynamics in various fields ranging from nanoscale over biological to geophysical and environmental systems. Sokolov and Klafter [4] presented a continuous-time random walk model for approximating the subdiffusive processes, and derived two different but equivalent forms of kinetic equations, which reduce to known fractional diffusion or Fokker–Planck equations for waiting-time distributions following a power law. Yuste and Acedo [5] obtained an explicit finite difference method and a new von Neumann-type stability analysis for the fractional subdiffusion equation. Langlands and Henry [6] investigated the accuracy and stability of an implicit numerical scheme for solving the fractional diffusion equation. Zahran et al. [7] derived the counterpart of fractional Fokker–Planck equation using comb-like model as fractal background medium. The resulting equation which is the fractional Fokker–Planck equation essentially differs from the usual Fokker–Planck equation (i.e., instead of first-order derivative of time evolution operator in Fokker–Planck equation, the derivatives of fractional order b arises in fractional Fokker–Planck equation)
uðx; tÞ ¼ Db t LFP uðx; tÞ;
0 < b 6 1;
ð1Þ
where
LFP ¼
@ @ dðxÞ v ðxÞ ; @x @x
q Research supported by National Natural Science Foundation of PR China (10531080) and the Scientific Research Foundation of Xiamen University of Technology (YKJ09016R). * Corresponding author. E-mail addresses: [email protected] (C. Wu), [email protected] (L. Lu).
0096-3003/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2010.03.024
1946
C. Wu, L. Lu / Applied Mathematics and Computation 216 (2010) 1945–1955
and
Db t uðx; tÞ ¼
1 CðbÞ
Z
t
0
uðx; gÞ ðt gÞ1b
dg:
The above expression is defined as the Riemann–Liouville fractional integral operator [8], CðÞ is the gamma function, dðxÞ > 0 and v ðxÞ are the diffusion and drift coefficients related to the fractal system under consideration, respectively. In order to construct a suitable numerical method, we differentiate both sides of (1) and let v ðxÞ ¼ 0, then the following equation is obtained
@uðx; tÞ @ @ ¼ D1b dðxÞ uðx; tÞ ; t @t @x @x where the symbol D1b is the Riemann–Liouville fractional derivative operator, defined by t
Dt1b uðx; tÞ ¼
1 @ CðbÞ @t
Z 0
t
uðx; gÞ ðt gÞ1b
dg:
Without loss of generality, a source term is added, then we consider the fractional Fokker–Planck equation (FFPE) without the drift term
@uðx; tÞ @ @ ¼ D1b dðxÞ uðx; tÞ þ f ðx; tÞ ; t @t @x @x
0 6 x 6 a; 0 < t 6 T
ð2Þ
with boundary conditions
uð0; tÞ ¼ u1 ðtÞ;
0 6 t 6 T;
uða; tÞ ¼ u2 ðtÞ;
0 6 t 6 T;
ð3Þ
and initial condition
uðx; 0Þ ¼ xðxÞ;
0 6 x 6 a;
ð4Þ
where 0 < b 6 1. When the fractional differential equation is derived from various fields, it is urgent to research that how to solve it. Some researchers are interested in the analytical solutions, for example, Mainardi et al. [9] discussed the Cauchy problem for the space–time fractional diffusion-wave equation, investigated the fundamental solution by Fourier–Laplace representation and provided a general representation of the solution in terms of Mellin-Barnes integrals in the complex plane by using the Mellin transform. Qi and Xu [10] studied the unsteady flow of viscoelastic fluid with the fractional derivative Maxwell model (FDMM) in a channel, the exact solutions were obtained for an arbitrary pressure gradient by means of the finite Fourier cosine transform and the Laplace transform. Wang and Xu [11] considered the unsteady axial Couette flow of fractional second grade fluid (FSGF) and fractional Maxwell fluid (FMF) between two infinitely long concentric circular cylinders, and obtained the analytical solutions by using the theory of Laplace transform and Weber transform. As analytic solutions of most fractional differential equations are not usually expressed explicitly, many authors resort to numerical solution strategies based on convergence and stability analysis. Liu et al. [12,13] considered the space fractional partial differential equation and transformed the equation into a system of ordinary differential equations that was then solved using backward differentiation formulas. Meerschaert and Tadjeran [14] presented practical numerical methods to solve the one-dimensional fractional advection–dispersion equations with variable coefficients on a finite domain. Shen and Liu [15] proposed an explicit finite difference approximation for the space fractional diffusion equation and gave an error analysis. Meerschaert et al. [16] developed practical numerical methods to solve two-dimensional fractional dispersion equations with variable coefficients on a finite domain and obtained the first-order accuracy in space and time. And they proposed the stability analysis by using matrix analysis method. Roop [17] investigated the computational aspects for the Galerkin approximation using continuous piecewise polynomial basis functions on a regular triangulation of a bounded domain in R2 . Liu et al. [18] discussed an approximation of the Lévy–Feller advection–dispersion process by a random walk and finite difference method. Chen et al. [19,20] proposed a Fourier method for the fractional diffusion equation describing subdiffusion and the Galilei invariant fractional advection diffusion equation. Zhuang et al. [21] proposed a new implicit numerical method and two solution techniques for improving its convergence order for solving an anomalous subdiffusion equation. The rest of this paper is organized as follows. In Section 2, we construct an effective implicit numerical approximation scheme for the FFPE (2) with boundary and initial conditions (3) and (4). In Sections 3 and 4, the stability and convergence analyses are discussed. Finally, in Section 5, some numerical results are given to testify our method and theoretical analyses.
2. An implicit numerical approximation scheme for the FFPE In this section, an implicit numerical approximation scheme for the FFPE (2) with boundary and initial conditions (3), (4) is proposed.
C. Wu, L. Lu / Applied Mathematics and Computation 216 (2010) 1945–1955
1947
Let K ¼ ½0; a ½0; T, then define the function space associated with a special function d(x)
(
@
@
@4
@
)
dðxÞ uðx; tÞ 2 C 2 ðKÞ; 3 UðKÞ ¼ uðx; tÞ dðxÞ uðx; tÞ 2 CðKÞ : @x @x @x @x @t And we suppose the FFPE (2)–(4) has a solution uðx; tÞ 2 Ud ðKÞ for the special function d(x) in (2). We now discretize space and time by grid points and time instants as follows:
a ; m T k ¼ 0; 1; 2; . . . ; K; s ¼ ; K
xi ¼ ih;
i ¼ 0; 1; 2; . . . ; m; h ¼
t k ¼ k s; where h and
xiþ1 2
s are the space and time steps, respectively. And we introduce the dual subdivision points
1 h; ¼ iþ 2
i ¼ 0; 1; 2; . . . ; m 1:
Then, we can discretize the variable uðxi ; tk Þ. For convenience, we denote
@ @u @ @uðx; tÞ d ðx; tÞ ¼ dðxÞ ; @x @x @x @x h h 2 d uðx; tÞ ¼ d x þ ; t ðuðx þ h; tÞ uðx; tÞÞ d x ; t ðuðx; tÞ uðx h; tÞÞ: d x 2 2 Luðx; tÞ ¼
Integrating both sides of (2) from t k to t kþ1 , we obtain
Z tk 1 Luðxi ; gÞ þ f ðxi ; gÞ dg CðbÞ 0 ðtkþ1 gÞ1b ðt k gÞ1b 0 Z s Z tk 1 Luðxi ; gÞ þ f ðxi ; gÞ 1 Lv ðxi ; gÞ þ gðxi ; gÞ ¼ uðxi ; tk Þ þ dg þ dg; 1b CðbÞ 0 C ðbÞ ðt kþ1 gÞ ðt k gÞ1b 0
uðxi ; tkþ1 Þ ¼ uðxi ; tk Þ þ
1 CðbÞ
Z
t kþ1
Luðxi ; gÞ þ f ðxi ; gÞ
dg
ð5Þ
where v ðx; tÞ ¼ uðx; t þ sÞ uðx; tÞ; gðx; tÞ ¼ f ðx; t þ sÞ f ðx; tÞ. We can rewrite (5) as
uðxi ; tkþ1 Þ ¼ uðxi ; tk Þ þ I1 þ I2 ; where
1 CðbÞ
I1 ¼
1 CðbÞ
I2 ¼
Z s Luðxi ; gÞ þ f ðxi ; gÞ ðtkþ1 gÞ1b
0
Z
tk
dg;
Lv ðxi ; gÞ þ gðxi ; gÞ ðtk gÞ1b
0
dg:
For convenience, we introduce the notation b
bk ¼ ðk þ 1Þb k ;
k ¼ 0; 1; 2; . . . ; K:
ð6Þ
For I1 , we can use the approximation
1 CðbÞ
Z s Luðxi ; sÞ þ f ðxi ; sÞ
sb
b
ððk þ 1Þb k Þ þ R11 ðtkþ1 gÞ sb 1 sb 1 ¼ bk 2 d d2x uðxi ; sÞ þ R12 þ f ðxi ; sÞ þ R11 ¼ bk 2 d d2x uðxi ; sÞ þ f ðxi ; sÞ þ R1 ; Cð1 þ bÞ h Cð1 þ bÞ h
I1 ¼
0
1b
dg þ R11 ¼ ½Luðxi ; sÞ þ f ðxi ; sÞ
where
Z s 1 Luðxi ; gÞ Luðxi ; sÞ þ f ðxi ; gÞ f ðxi ; sÞ dg; CðbÞ 0 ðt kþ1 gÞ1b @ @u 1 ¼ d ðxi ; sÞ 2 d d2x uðxi ; sÞ; @x @x h
R11 ¼ R12
R1 ¼ Noting that
sb Cð1 þ bÞ
bk R12 þ R11 :
Cð1 þ bÞ
1948
C. Wu, L. Lu / Applied Mathematics and Computation 216 (2010) 1945–1955
" # @ @u @ @u @2 @u @f ðxi ; g1 Þ ðg sÞ; d ðxi ; gÞ þ f ðxi ; gÞ ¼ d ðxi ; sÞ þ f ðxi ; sÞ þ ðxi ; g1 Þ þ d @x @x @x @x @x @t @x @t
Luðxi ; gÞ þ f ðxi ; gÞ ¼
where 0 6 g 6 g1 6 s, we obtain
Z s
1 CðbÞ
jR11 j 6 C 1 s
1 1b
ðt kþ1 gÞ
0
dg 6
C 1 s1þb
Cð1 þ bÞ
bk :
By Taylor’s formula it is apparent that 2
jR12 j 6 Ch : 2
Hence, we have jR1 j 6 Cbk sb ðs þ h Þ. For I2 , we can use the approximation
I2 ¼ ¼ ¼
Z
1 CðbÞ
Lv ðxi ; gÞ þ gðxi ; gÞ
tk
ðt k gÞ
0
Z
k1 X
1 CðbÞ
sb
k1 X
Cð1 þ bÞ
s¼0
dg ¼
k1 Z tsþ1 1 X Lv ðxi ; gÞ þ gðxi ; gÞ dg CðbÞ s¼0 ts ðt k gÞ1b
Lv ðxi ; t sþ1 Þ þ gðxi ; tsþ1 Þ
tsþ1
ðtk gÞ1b
ts
s¼0
1b
dg þ R21 ¼
sb
k1 X
Cð1 þ bÞ
s¼0
bks1 ½Lv ðxi ; t sþ1 Þ þ gðxi ; tsþ1 Þ þ R21
bks1 ½Lh v ðxi ; t sþ1 Þ þ gðxi ; tsþ1 Þ þ R22 þ R21 ;
where
R21 ¼ R22 ¼
k1 Z t sþ1 1 X Lv ðxi ; gÞ Lv ðxi ; tsþ1 Þ þ gðxi ; gÞ gðxi ; tsþ1 Þ dg; CðbÞ s¼0 ts ðtk gÞ1b
sb
k1 X
Cð1 þ bÞ
s¼0
Lh uðxi ; ts Þ ¼
bks1 ½Lv ðxi ; t sþ1 Þ Lh v ðxi ; t sþ1 Þ;
1
d2 uðxi ; ts Þ: 2d x
h
Noting that
Lv ðxi ; gÞ þ gðxi ; gÞ ¼
" # @ @v @2 @v @gðxi ; fÞ ðg tsþ1 Þ; d ðxi ; t sþ1 Þ þ gðxi ; t sþ1 Þ þ d ðxi ; fÞ þ @x @t @x @x @t @x
where t s 6 g 6 f 6 t sþ1 . And,
@2 @v @2 @u @2 @u @3 @u ðxi ; f þ sÞ ðxi ; fÞ ¼ ðxi ; f1 Þs; ðxi ; fÞ ¼ d d d d 2 @x @x @x @x @t @x @x @t @x @t @x @t @gðxi ; fÞ @f ðxi ; f þ sÞ @f ðxi ; fÞ @ 2 f ðxi ; f2 Þ s; ¼ ¼ @t @t @t @t 2 where f 6 f1 ; f2 6 f þ s. Consequently,
jLv ðxi ; gÞ Lv ðxi ; t sþ1 Þj 6 C 3 s2 : So we obtain
jR21 j 6 C 3 s2
1 CðbÞ
Z 0
tk
dg ðt k gÞ1b
6 C s2 :
From Taylor’s formula, it is apparent that
2 2 @ @v 1 h @3 @v h @3 @v ðxi ; tsþ1 Þ ¼ 2 d d2x uðxi ; tsþ1 Þ ðn; t ðn; t sþ1 Þ d Þ ¼ L v ðx ; t Þ d d sþ1 sþ1 i h @x @x 24 @x3 @x 24 @x3 @x h " # 2 2 h @3 @u @3 @u h s @4 @u ðn; tsþ2 Þ 3 d ðn; tsþ1 Þ ¼ Lh v ðxi ; t sþ1 Þ ðn; g2 Þ: ¼ Lh v ðxi ; tsþ1 Þ d d 3 3 @x @x @x 24 @x @x 24 @x @t
Lv ðxi ; t sþ1 Þ ¼
Thus, we have 2
jR22 j 6 C sh :
1949
C. Wu, L. Lu / Applied Mathematics and Computation 216 (2010) 1945–1955
From the above results, we obtain
uðxi ; tkþ1 Þ ¼ uðxi ; tk Þ þ r 1 bkd d2x uðxi ; sÞ þ r 1
k1 X
bks1d d2x ðuðxi ; t sþ2 Þ uðxi ; tsþ1 ÞÞ þ r 2 bk f ðxi ; sÞ þ r 2
s¼0
k1 X
bks1 ðf ðxi ; tsþ2 Þ
s¼0
; f ðxi ; t sþ1 ÞÞ þ Rkþ1 i b
b
s s where r1 ¼ Cð1þbÞh 2 ; r 2 ¼ Cð1þbÞ, and 2
jRkþ1 j 6 Cðbk sb þ sÞðs þ h Þ: i
ð7Þ
To analyze the stability, we find it worthwhile to recall here the following useful lemma associated with the property of the coefficients bk ðk ¼ 0; 1; 2; . . .Þ referred to in [21]. Lemma 2.1. The coefficients bk ðk ¼ 0; 1; 2; . . .Þ defined by (6) satisfy (1) b0 ¼ 1; bk > 0; k ¼ 0; 1; 2; . . .; (2) bk > bkþ1 ; k ¼ 0; 1; 2; . . .; (3) there exists a positive constant C > 0 such that
s 6 Cbk sb ; k ¼ 1; 2; . . . Let Rk ¼ ðRk1 ; Rk2 ; . . . ; Rkm1 ÞT . By using Lemma 2.1 and (7), we obtain the following lemma. qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Pm1 k 2 d Lemma 2.2. Suppose that kRk k2 ¼ i¼1 jRi j h, if uðx; tÞ 2 U ðKÞ is the solution of (2)–(4) , then we have 2
kRkþ1 k2 6 Cbk sb ðs þ h Þ: Let uki ; fik be the numerical approximation to uðxi ; t k Þ; f ðxi ; tk Þ; respectively, and introduce the following notation: 2 k d dx ui
¼ diþ1 ðukiþ1 uki Þ di1 ðuki uki1 Þ; 2
Dx uki ¼ ukiþ1 uki :
2
We obtain the following implicit numerical approximation scheme (INAS):
ukþ1 ¼ uki þ r 1 bkd d2x u1i þ r1 i
k1 X
bks1d d2x ðusþ2 usþ1 Þ þ r2 bk fi1 þ r 2 i i
s¼0
k1 X
bks1 ðfisþ2 fisþ1 Þ;
s¼0
which can be rewritten as
ukþ1 ¼ uki þ r 1d d2x ukþ1 þ r1 i i
k1 k1 X X ðbsþ1 bs Þd d2x uks þ r 2 fikþ1 þ r2 ðbsþ1 bs Þfiks ; i s¼0
ð8Þ
s¼0
for i ¼ 1; 2; . . . ; m 1; k ¼ 0; 1; 2; . . . ; K 1. The boundary and initial conditions are
uk0 ¼ u1 ðksÞ; ukm ¼ u0i ¼
k ¼ 1; 2; . . . ; K;
u2 ðksÞ; k ¼ 1; 2; . . . ; K; xðihÞ; i ¼ 0; 1; . . . ; m:
ð9Þ
Then we rewrite Eqs. (8) and (9) as the following matrix form:
8 1 0 1 1 > < Au ¼ u þ r 1 w þ r 2 f ; kþ1 k kþ1 Au ¼ u þ r1 v þ r 1 wkþ1 þ r2 f kþ1 þ r 2 g kþ1 ; > : 0 u ¼ ;
ð10Þ
k > 0;
where
2
uk1
3
6 k 7 6 u2 7 7 6 7 6 uk ¼ 6 ... 7; 7 6 6 k 7 4 um2 5 ukm1 and
v kþ1 ¼
Pk1
s¼0 ðbsþ1
2
f1k
3
6 k 7 6 f2 7 7 6 7 6 f k ¼ 6 ... 7; 7 6 6 k 7 4 fm2 5
2
wkþ1
6 6 6 6 ¼6 6 6 4
k fm1
bs Þd d2x uks ; g kþ1 ¼
d1 ukþ1 0 2
0 .. . 0 kþ1 dm1 um1
3 7 7 7 7 7; 7 7 5
s¼0 ðbsþ1
x1 7 6 6 x2 7
7 6 6 . 7 ¼ 6 .. 7; 7 6 7 6 4 xm2 5
2
Pk1
3
2
xm1
bs Þf ks ; xi ¼ xðihÞ.
The matrix A in (10) can be written as the following form:
1950
C. Wu, L. Lu / Applied Mathematics and Computation 216 (2010) 1945–1955
2 6 6 6 6 6 6 6 6 6 4
1 þ r1 d3 þ d1 2
2
r 1 d3
r1 d3 2 1 þ r 1 d5 þ d3
..
.
..
..
.
2
2
3
0
2
.
0
r1 dm3 2 1 þ r1 dm1 þ dm3
r1 dm3 2
2
7 7 7 7 7 7: 7 7 7 5
ð11Þ
2
From the structure of the matrix A we can see that it is a strictly diagonally dominant symmetric tridiagonal matrix with positive diagonal terms and nonpositive off-diagonal terms. Then, we obtain the following theorem. Theorem 2.1. The discretization matrix A is invertible. And the system (8) and (9) has a unique solution.
3. Stability of the INAS For u ¼ ðu1 ; u2 ; . . . ; um1 ÞT ;
ðu; v Þ ¼
m 1 X
ui v i h;
v ¼ ðv 1 ; v 2 ; . . . ; v m1 ÞT , we define the inner product
ðu; v Þd ¼
m1 X
i¼1
i¼1
diþ1 ui v i h; 2
and the norm
kuk2 ¼
m1 X
pffiffiffiffiffiffiffiffiffiffiffiffi ðu; uÞ ¼
!12 u2i h
;
kukd ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðu; uÞd ¼
i¼1
m1 X i¼1
!12 diþ1 u2i h 2
:
~ ki ði ¼ 0; 1; 2; . . . ; m; k ¼ 0; 1; 2; . . . ; KÞ is the approximation solution of (8) and (9). Then the error We suppose that u eki ¼ ue ki uki satisfies
eikþ1 ¼ eki þ r1d d2x ekþ1 þ r1 i
k1 X ðbsþ1 bs Þd d2x eiks ;
ð12Þ
s¼0
and
ek0 ¼ ekm ¼ 0; k ¼ 1; 2; . . . ; K: k
Let E ¼ ðe
k 1;
e
k 2; . . . ;
e
T k m1 Þ :
Multiplying (12) by
ð13Þ
e
kþ1 h i
kEkþ1 k22 ¼ ðEkþ1 ; Ek Þ þ r1 ðd d2x Ekþ1 ; Ekþ1 Þ þ r1
k1 X
and summing i from 1 to m 1, we have
ðbsþ1 bs Þðd d2x Eks ; Ekþ1 Þ
s¼0 k1 X ks kþ1 ðbs bsþ1 Þ d1 eks ; Dx Ekþ1 Þd : ¼ ðEkþ1 ; Ek Þ r1 d1 ðe1kþ1 Þ2 h þ ðDx Ekþ1 ; Dx Ekþ1 Þd þ r 1 1 e1 h þ ðDx E 2
2
s¼0
By applying the Cauchy–Schwarz inequality, we have
kEkþ1 k22 6
k1 r X 1 1 Þ2 h þ kDx Ekþ1 k2d þ ðbs bsþ1 Þ d1 ðe1ks Þ2 h þ d1 ðe1kþ1 Þ2 h ðkEkþ1 k22 þ kEk k22 Þ r1 d1 ðekþ1 1 2 2 2 2 2 s¼0
þ Noting that
Pk1
k1 r1 X ðbs bsþ1 ÞðkDx Eks k2d þ kDx Ekþ1 k2d Þ: 2 s¼0
s¼0 ðbs
kEkþ1 k22 6
bsþ1 Þ ¼ 1 bk and bk > 0, we have
k1 r X 1 r1 1 Þ2 h þ kDx Ekþ1 k2d þ ðbs bsþ1 Þ d1 ðeks Þ2 h þ kDx Eks k2d ; ðkEkþ1 k22 þ kEk k22 Þ ð1 þ bk Þ d1 ðekþ1 1 1 2 2 2 2 2 s¼0
i.e.,
kEkþ1 k22 þ r 1
k X
k1 X 2 ks 2 bs d1 ðe1kþ1s Þ2 h þ kDx Ekþ1s k2d 6 kEk k22 þ r1 bs d1 ðeks kd : 1 Þ h þ kDx E 2
s¼0
s¼0
Defining the energy norm
kEk k2E ¼ kEk k22 þ r1
k1 X s¼0
2 ks 2 bs d1 ðeks kd ; 1 Þ h þ kDx E 2
2
ð14Þ
C. Wu, L. Lu / Applied Mathematics and Computation 216 (2010) 1945–1955
1951
we have
kEkþ1 k22 6 kEkþ1 k2E 6 kEk k2E 6 6 kE1 k2E : As
e1i ¼ e0i þ r1 d d2x e1i ; we have 1 kE1 k22 ¼ ðE1 ; E0 Þ r 1 d1 ðe11 Þ2 h þ kDx E1 k2d 6 ðkE1 k22 þ kE0 k22 Þ r 1 d1 ðe11 Þ2 h þ kDx E1 k2d : 2 2 2
Then,
kE1 k2E ¼ kE1 k22 þ r1 d1 ðe11 Þ2 h þ kDx E1 k2d 6 kE0 k22 ; 2
ð15Þ
so that kEkþ1 k22 6 kE0 k22 : Hence, we can obtain the following theorem of stability. Theorem 3.1. The implicit numerical approximation scheme (8) is unconditionally stable.
4. Convergence of the INAS In this section, we discuss the convergence of the INAS. Suppose that u(x, t) is the solution of the FFPE (2)–(4) and uðx; tÞ 2 Ud ðKÞ. Let uðxi ; tk Þ ði ¼ 0; 1; 2; . . . ; m; k ¼ 0; 1; 2; . . . ; KÞ be the exact solution of (2)–(4) at the mesh point ðxi ; t k Þ. Define
yki ¼ uðxi ; tk Þ uki ; k
and Y ¼
i ¼ 0; 1; 2; . . . ; m; k ¼ 0; 1; 2; . . . ; K;
ðyk1 ; yk2 ; . . . ; ykm1 ÞT .
By using uki ¼ uðxi ; tk Þ yki , substitution into (8) leads to
ykþ1 ¼ yki þ r 1d d2x ykþ1 þ r1 i i
k1 X ðbsþ1 bs Þd d2x yiks þ Rkþ1 ; i
ð16Þ
s¼0
where i ¼ 1; 2; . . . ; m 1; k ¼ 0; 1; . . . ; K 1, and
y0i ¼ 0; i ¼ 0; 1; . . . ; m;
yk0 ¼ ykm ¼ 0; k ¼ 0; 1; . . . ; K:
Multiplying (16) by ykþ1 h and summing i from 1 to m 1, we obtain i
kY kþ1 k22 ¼ ðY kþ1 ; Y k Þ þ r 1 ðd d2x Y kþ1 ; Y kþ1 Þ þ r1
k1 X ðbsþ1 bs Þðd d2x Y ks ; Y kþ1 Þ þ ðRkþ1 ; Y kþ1 Þ:
ð17Þ
s¼0
For s ¼ 0; 1; . . . ; k þ 1; we have s kþ1 ðd d2x Y s ; Y kþ1 Þ ¼ d1 ys1 ykþ1 Þd ; 1 h ðDx Y ; Dx Y 2
ðY
kþ1
1 ; Y Þ 6 ðkY kþ1 k22 þ kY s k22 Þ; 2 s
and s ykþ1 1 y1 6
1 kþ1 2 ðjy j þ jys1 j2 Þ: 2 1
Taking
c ¼ min diþ1 ; i
ð18Þ
2
and using jv wj 6 rv 2 þ 41r w2 ;
r > 0; we obtain
2
jðRkþ1 ; Y kþ1 Þj 6
2r 1 h bk c kþ1 2 a2 kY k2 þ kRkþ1 k22 : 2 2 a 8r 1 h bk c
ð19Þ
Similarly to the proof of the stability, we obtain
kY kþ1 k22 6
k1 r X 1 r1 1 Þ2 h þ kDx Y kþ1 k2d þ ðbs bsþ1 Þ d1 ðyks Þ2 h þ kDx Y ks k2d ðkY kþ1 k22 þ kY k k22 Þ ð1 þ bk Þ d1 ðykþ1 1 1 2 2 2 2 2 s¼0 2
þ
2r 1 h bk c kþ1 2 a2 kRkþ1 k22 : kY k2 þ 2 a2 8r 1 h bk c
Let jyki0 j ¼ max16i6m1 jyki j, then
ð20Þ
1952
C. Wu, L. Lu / Applied Mathematics and Computation 216 (2010) 1945–1955
yki0 ¼ yk1 þ
i0 1 X
Dx yki ;
yki0 ¼
m1 X
Dx yki :
i¼i0
i¼1
Thus
2jyki0 j 6 jyk1 j þ
m1 X
jDx yki j:
i¼1
Using the Cauchy–Schwarz inequality, we have
jyk1 j2
m 1 X
4jyki0 j2
6m
k
max16i6m1 jyki j,
þ
! jDx yki j2
h
i¼1
Let kY k1 ¼
kY k k21
a
6
4h
ðjyk1 j2 h 2
a
6
2
ðjyk1 j2 h þ kDx Y k k22 Þ:
Therefore,
þ kDx Y k k22 Þ:
Then we have
kY k k22 6 akY k k21 6
a2 4h
2
ðjyk1 j2 h þ kDx Y k k22 Þ:
Noting (18), we have 2
2
2r 1 h bk c kþ1 2 2r 1 h bk c a2 r1 2 kþ1 2 kY k2 6 ðjykþ1 k2 Þ ¼ bk ðcjy1kþ1 j2 h þ ckDx Y kþ1 k22 Þ 1 j h þ kDx Y 2 a2 a2 2 4h r1 kþ1 2 kþ1 2 6 bk d1 jy1 j h þ kDx Y kd : 2 2
ð21Þ
Hence, from (20) and (21), we have
kY kþ1 k22 6
k1 r X 1 r1 1 2 2 kþ1 2 ks 2 kd þ ðbs bsþ1 Þ d1 ðyks kd ðkY kþ1 k22 þ kY k k22 Þ d1 ðykþ1 1 Þ h þ kDx Y 1 Þ h þ kDx Y 2 2 2 2 2 s¼0
a2
þ Let ck ¼ kY k k22 þ r1
2
8r1 h bk c
kRkþ1 k22 :
Pk1
ks 2 s¼0 bs ðd12 ðy1 Þ h
þ kDx Y ks k2d Þ and then
ckþ1 6 ck þ Cbk sb ðs þ h2 Þ2 ; c0 ¼ 0; Hence, we obtain
ckþ1 6 C
k X
2
bs sb ðs þ h Þ2 :
s¼0
Noting that
Pk
s¼0 bs
sb ¼ ðk þ 1Þb sb 6 T b and kY kþ1 k22 6 ckþ1 ; we have 2
kY kþ1 k22 6 CT b ðs þ h Þ2 : Consequently, the following theorem of convergence is obtained. Theorem 4.1. Let uðx; tÞ 2 Ud ðKÞ be the solution of the FFPE (2)–(4) . Then the INAS (8) is convergent, and there exists a positive constant C > 0 such that 2
kY kþ1 k2 6 Cðs þ h Þ;
k ¼ 0; 1; 2; . . . ; K 1:
5. Numerical examples In this section, two numerical examples are presented to support our theoretical analysis. Example 1. Consider the following FFPE
@uðx; tÞ @ @ ¼ D1b dðxÞ uðx; tÞ t @t @x @x with initial and boundary conditions
ð22Þ
1953
C. Wu, L. Lu / Applied Mathematics and Computation 216 (2010) 1945–1955
( uðx; 0Þ ¼
x;
x 2 ½0; 1;
3x ; 2
x 2 ½1; 3;
uð0; tÞ ¼ uð3; tÞ ¼ 0;
t 2 ½0; 1:
The INAS is used to solve Eq. (22). Firstly, We take the diffusion coefficient dðxÞ ¼ 1. b ¼ 0:7; 0 6 x 6 3; 0 6 t 6 1 and pffiffiffi t ¼ 0:5; 0 6 x 6 3; 0:1 6 b 6 0:9, then taking the diffusion coefficient dðxÞ ¼ 3 x; b ¼ 0:7; 0 6 x 6 3; 0 6 t 6 1; The numerical simulation is shown in Figs. 1–3, respectively. They show the process exhibits subdiffusion behaviors. In order to show the approximation order of the INAS, we construct an example with an analytic solution. Example 2. The following FFPE with analytic solution is considered
@uðx; tÞ @ @ ¼ D1b dðxÞ uðx; tÞ þ f ðx; tÞ ; t @t @x @x
0 < x < 1; 0 < t < 1;
ð23Þ
where dðxÞ ¼ ex ; f ðx; tÞ ¼ Cðb þ 3Þex t 2 =2 2e2x tbþ2 . with the initial and boundary conditions
uð0; tÞ ¼ t 2þb ;
0 < t < 1;
uð1; tÞ ¼ et 2þb ; uðx; 0Þ ¼ 0;
0 < t < 1;
0 6 x 6 1:
The exact solution of (23) is uðx; tÞ ¼ ex t2þb :
u(x,t)
1.5 1
0.5 0 1
0.8
0.6
t
0.4
0.2
0 0
1
2
3
x
Fig. 1. Numerical approximation for b ¼ 0:7; dðxÞ ¼ 1.
0.5
u(x,t=0.5)
0.4 0.3 0.2 0.1 0 1
0.8
0.6
β
0.4
0.2
0 0
1
x
2
Fig. 2. Numerical approximation for different b at t ¼ 0:5.
3
1954
C. Wu, L. Lu / Applied Mathematics and Computation 216 (2010) 1945–1955
1.5
u(x,t)
1 0.5 3
0 1
0.8
2
0.6 t
0.4
0.2
1
x
0 0
Fig. 3. Numerical approximation for b ¼ 0:7; dðxÞ ¼
ffiffiffi p 3 x:
Table 1 The maximum absolute error of the INAS.
s¼h
b ¼ 0:4
b ¼ 0:7
b ¼ 0:9
1 50 1 100 1 200
3.0013994E003
4.6906049E003
6.0741258E003
1.5066619E003
2.3543905E003
3.0419634E003
5.9997852E004
1.2317856E003
1.5027906E003
0.5 The exact solution The numerical solution
0.45
u(x,t=0.5)
0.4 0.35 0.3 0.25 0.2 0.15 0
0.2
0.4
x
0.6
0.8
1
Fig. 4. The exact solution and numerical solution when b ¼ 0:5; t ¼ 0:5.
The maximum absolute error between the exact solution and the numerical solutions by INAS, with spatial and temporal steps s ¼ h ¼ 1=50; 1=100; 1=200 at time t ¼ 1:0, is listed in Table 1. In Table 1 the last column is computed as the ratio of the maximum absolute error at the previous grid size to the current grid size, which shows the order of convergence of the method as all step sizes are halved. This is in good agreement with the analysis of theory. We also compare the exact solution with the numerical solution by using INAS with different time and space steps at t ¼ 0:5 and b ¼ 0:5 which is shown in Fig. 4, we can see that the numerical solution is in good agreement with the exact solution. 6. Conclusions In this paper, an implicit numerical approximation scheme for the fractional Fokker–Planck equation has been described and demonstrated. The stability, consistency and convergence of the INAS for the FFPE have been discussed. These methods and analytical techniques can also be extended to some high-dimensional fractional partial differential equations.
C. Wu, L. Lu / Applied Mathematics and Computation 216 (2010) 1945–1955
1955
Acknowledgement The authors thank the referee for many thoughtful and constructive suggestions to improve the paper. References [1] Walter Wyss, The fractional diffusion equation, J. Math. Phys. 27 (1986) 2782–5785. [2] Ralf Metzler, Joseph Klafter, The random walk’s guide to anomalous diffusion: a fractional dynamics approach, Phys. Rep. 339 (2000) 1–77. [3] Ralf Metzler, Joseph Klafter, The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics, J. Phys. A 37 (2004) R161–208. [4] I.M. Sokolov, J. Klafter, From diffusion to anomalous diffusion: a century after Einstein’s Brownian motion, Chaos 15 (2005) 026103. [5] S.B. Yuste, L. Acedo, An explicit finite difference method and a new von Neumann-type stability analysis for fractional diffusion equations, SIAM J. Numer. Anal. 42 (2004) 1862–1874. [6] T.A.M. Langlands, B.I. Henry, The accuracy and stability of an implicit solution method for the fractional diffusion equation, J. Comput. Phys. 205 (2005) 719–736. [7] M.A. Zahran, E.M. Abulwafa, S.A. Elwakil, The fractional Fokker–Planck equation on comb-like model, Physica A 323 (2003) 237–248. [8] I. Podlubny, Fractional Differential Equations, Academic Press, 1999. [9] F. Mainardi, Y. Luchko, G. Pagnini, The fundamental solution of the space-time fractional diffusion equation, Fract. Calculus Appl. Anal. 4 (2) (2001) 153–192. [10] H. Qi, M. Xu, Unsteady flow of viscoelastic fluid with fractional Maxwell model in a channel, Mech. Res. Commun. 34 (2007) 210–212. [11] S. Wang, M. Xu, Axial Couette flow of two kinds of fractional viscoelastic fluids in an annulus, Nonlinear Anal.: Real World Appl. (2008). [12] F. Liu, V. Anh, I. Turner, Numerical solution of the space fractional Fokker–Planck equation, J. Comput. Appl. Math. 166 (2004) 209–219. [13] F. Liu, V. Anh, I. Turner, P. Zhuang, Numerical solution for the solute transport in fractal porous media, ANZIAM J. 45 (E) (2004) 461–473. [14] M.M. Meerschaert, C. Tadjeran, Finite difference approximations for fractional advection–dispersion flow equations, J. Comput. Appl. Math. 172 (1) (2004) 65–77. [15] S. Shen, F. Liu, Error analysis of an explicit finite difference approximation for the space fractional diffusion, ANZIAM J. 46 (E) (2005) 871–887. [16] M.M. Meerschaert, H. Scheffler, C. Tadjeran, Finite difference methods for two-dimensional fractional dispersion equation, J. Comput. Phys. 211 (2006) 249–261. [17] J.P. Roop, Computational aspects of FEM approximation of fractional advection dispersion equation on bounded domains in R2 , J. Comput. Appl. Math. 193 (1) (2006) 243–268. [18] Q. Liu, F. Liu, I. Turner, V. Anh, Approximation of the Lévy–Feller advection–dispersion process by random walk and finite difference method, J. Comput. Phys. 222 (2007) 57–70. [19] Chang-Ming Chen, F. Liu, I. Turner, V. Anh, A Fourier method for the fractional diffusion equation describing sub-diffusion, J. Comput. Phys. 227 (2007) 886–897. [20] Chang-Ming Chen, F. Liu, I. Turner, V. Anh, A new Fourier analysis method for the Galilei invariant fractional advection diffusion equation, ANZIAM J. 48 (2007) 605–619. [21] P. Zhuang, F. Liu, V. Anh, I. Turner, New solution and analytical techniques of the implicit numerical method for the anomalous subdiffusion equation, SIAM J. Numer. Anal. 46 (2) (2008) 1079–1095.