- Email: [email protected]
Numerical scheme with high order accuracy for solving a modified fractional diffusion equation
Applied Mathematics and Computation 244 (2014) 772–782
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Numerical scheme with high order accuracy for solving a modified fractional diffusion equation Y. Chen, Chang-Ming Chen ⇑ School of Mathematical Sciences, Xiamen University, Xiamen 361006, China
a r t i c l e
i n f o
Keywords: Modified fractional diffusion equation Numerical scheme with second order temporal accuracy and fourth order spatial accuracy Convergence Stability Solvability Fourier analysis
a b s t r a c t In recent years, some researchers have developed various numerical schemes to solve the modified fractional diffusion equation. For the numerical solutions of the modified fractional diffusion equation, there are already some important progresses. In this paper, a numerical scheme with second order temporal accuracy and fourth order spatial accuracy is developed to solve a modified fractional diffusion equation; the convergence, stability and solvability of the numerical scheme are analyzed by Fourier analysis; the theoretical results extremely consistent with the numerical experiment. Ó 2014 Elsevier Inc. All rights reserved.
1. Introduction In recent years, the time or space or time–space fractional diffusion equations are widely used to describe anomalous diffusion processes in numerous physical and biological systems, many authors have developed various numerical schemes to solve these fractional diffusion equations. By the inclusion of a secondary fractional time derivative acting on a diffusion operator, a model for describing processes that become less anomalous as time progresses has been proposed [3,10,11]
@pðx; tÞ ¼ @t
A
! @ 1a @ 1b @ 2 pðx; tÞ þ B ; @x2 @t 1a @t 1b
ð1Þ
for this modified fractional diffusion equation, Langlands et al. proposed the solution with the form of an infinite series of Fox functions on an infinite domain [6]; Merdan et al. presented the numerical solution of time-fraction modified equal width wave equation, and show how an application of fractional two dimensional differential transformation method obtained approximate analytical solution of time-fraction modified equal width wave equation [9]; Liu et al. discussed the numerical method and analytical technique of the modified fractional diffusion equation with a nonlinear source term, they proved that this method has first order temporal accuracy and second order spatial accuracy using a new energy method [7]; Liu et al. researched the finite element approximation for a modified fractional diffusion equation, they proved that this approximation has first order temporal accuracy and m (here m is the degree of the piecewise polynomials) order spatial accuracy [8]; Zhang et al. studied the finite difference/element method for a two-dimensional modified fractional diffusion equation, they proved that this method has ð1 þ minfa; bgÞ order temporal accuracy and m (here m is the degree of the piecewise polynomials) order spatial accuracy [12]; Chen researched the numerical method ⇑ Corresponding author. E-mail address: [email protected] (C.-M. Chen). http://dx.doi.org/10.1016/j.amc.2014.07.044 0096-3003/Ó 2014 Elsevier Inc. All rights reserved.
773
Y. Chen, C.-M. Chen / Applied Mathematics and Computation 244 (2014) 772–782
for solving a two-dimensional variable-order modified fractional diffusion equation, this method has first order temporal accuracy and fourth order spatial accuracy, further, a numerical method for improving temporal accuracy have also been developed [4]; apply the fourth-order compact scheme of the second-order space partial derivative and the Grünwald–Letnikov discretization of the Rieman–Liouille fractional time derivative, Abbaszadeh et al. proposed a high-order and unconditionally stable scheme for the modified anomalous fractional sub-diffusion equation with a nonlinear source term, they proved that this scheme has first order temporal accuracy and fourth order spatial accuracy using the Fourier method [1], and they also presented a fourth-order compact solution of the two-dimensional modified anomalous fractional sub-diffusion equation with a nonlinear source term, they proved that this solution has first order temporal accuracy and fourth order spatial accuracy using the Fourier method [2]. On the basis of the literature, for the cases of one-dimensional and multidimensional of the modified fractional diffusion equation Eq. (1), although the numerical schemes with ideal spatial accuracy can be constructed by the fourth-order compact scheme of the second-order space partial derivative or the finite element approximation, but unfortunately that numerical scheme with second-order temporal accuracy is still unable to achieve. Chen constructed a numerical method for improving temporal accuracy for solving a two-dimensional variable-order modified fractional diffusion equation, but it should be point out that it only is supported by the numerical experiment, whereas it is not supported by rigorous numerical analysis [4]. Therefore to date, numerical scheme with second order temporal accuracy for solving modified fractional diffusions supported by rigorous numerical analysis does not exist. So, how to improve temporal accuracy of numerical scheme for solving modified fractional diffusions is a exciting and meaningful work. Easy to understand that the key problem is that convergence and stability analysis of numerical scheme are very difficult. We expect that this article will change this situation. In this paper, we study numerical scheme and perform related numerical analysis for solving the following modified fractional diffusion equation (MFDE)
! @ 1a @ 1b @ 2 pðx; tÞ þ þ f ðx; tÞ; @x2 @t 1a @t 1b
@pðx; tÞ ¼ @t
ð2Þ
with the following initial and boundary conditions:
pðx; 0Þ ¼ wðxÞ;
0 6 x 6 L;
pð0; tÞ ¼ uðtÞ;
pðL; tÞ ¼ wðtÞ;
where 0 < a < b < 1, and 1c pðx; ;tÞ 0 Dt
¼
1
ð3Þ 0 6 t 6 T;
1c pðx; tÞ 0 Dt
@
Z
CðcÞ @t
t 0
ð4Þ
is the Rieman–Liouille fractional partial derivative of order 1 c defined by
pðx; sÞ ðt sÞ1c
ds:
In this paper, we always suppose that pðx; tÞ 2 PðXÞ is the exact solution of the problem (1)–(3) and
@ 2 f ðx;tÞ @t 2
X ¼ fðx; tÞj0 6 x 6 L; 0 6 t 6 T g; ( PðXÞ ¼
) @ 6 pðx; tÞ @ 4 pðx; y; tÞ pðx; tÞj ; 2 CðXÞ : @x6 @x2 @t2
2. A numerical scheme for solving MFDE In this paper, we let
xj ¼ jh;
j ¼ 0; 1; . . . ; J;
where h ¼ L=J and
d2x
v
i j
¼v
i j1
tk ¼ ks;
k ¼ 0; 1; . . . ; K;
s ¼ T=K is space step and time step, respectively. And define
2v ij þ v ijþ1 :
Integrating both sides of Eq. (2) from t k1 to t k and take x ¼ xj , get that
1 pðxj ; t k Þ pðxj ; tk1 Þ ¼ CðaÞ
"Z
1 þ CðbÞ
tk
0
"Z 0
@ 2 pðxj ; sÞ ds @x2 ðtk sÞ1a tk
Z
t k1
0
@ 2 pðxj ; sÞ ds @x2 ðt k sÞ1b
Z 0
@ 2 pðxj ; sÞ ds @x2 ðt k1 sÞ1a
t k1
#
# Z tk @ 2 pðxj ; sÞ ds þ f ðxj ; sÞds 1b @x2 ðt k1 sÞ t k1
2 CðXÞ, where
774
Y. Chen, C.-M. Chen / Applied Mathematics and Computation 244 (2014) 772–782
from this we have that
"Z # Z tk1 2 tk 2 1 1 1 1 @ pðxj ; sÞ ds @ pðxj ; sÞ ds 1 þ d2x pðxj ;t k Þ 1 þ d2x pðxj ;t k1 Þ ¼ 1 þ d2x 12 12 12 @x2 ðt k sÞ1a @x2 ðt k1 sÞ1a CðaÞ 0 0 " # Z tk 2 Z tk1 2 1 2 1 @ pðxj ;sÞ ds @ pðxj ; sÞ ds þ 1 þ dx 12 @x2 ðt k sÞ1b @x2 ðt k1 sÞ1b CðbÞ 0 0 Z tk 1 f ðxj ;sÞds: þ 1 þ d2x 12 t k1 Because pðx; tÞ 2 PðXÞ gives
1 2 @ 2 pðxj ; t i Þ d2x pðxj ; t i Þ 4 1þ dx ¼ þ Oðh Þ; 2 12 @x2 h and
@ 2 pðxj ; sÞ t l s @ 2 pðxj ; t l1 Þ s tl1 @ 2 pðxj ; t l Þ ¼ þ þ Oðs2 Þ; @x2 @x2 @x2 s s
t l1 6 s 6 tl ;
then
"Z # Z tk1 2 tk 1 2 1 @ 2 pðxj ; sÞ ds @ pðxj ; sÞ ds 1þ dx 12 @x2 @x2 CðaÞ 0 ðt k sÞ1a ðt k1 sÞ1a 0 ! k Z tl 1 2 1 X t l s @ 2 pðxj ; t l1 Þ s t l1 @ 2 pðxj ; tl Þ ds 2 d þ þ Oðs Þ ¼ 1þ 12 x CðaÞ l¼1 tl1 @x2 @x2 s s ðt k sÞ1a ! k1 Z t l 1 2 1 X t l s @ 2 pðxj ; t l1 Þ s t l1 @ 2 pðxj ; tl Þ ds 2 1þ d þ þ Oðs Þ 12 x CðaÞ l¼1 tl1 @x2 @x2 s s ðt k1 sÞ1a " # Z k tl 1 X tl s 1 2 @ 2 pðxj ; t l1 Þ s t l1 1 2 @ 2 pðxj ; t l Þ ds 2 ¼ 1þ dx þ 1 þ d þ Oð s Þ 12 12 x CðaÞ l¼1 tl1 s s @x2 @x2 ðt k sÞ1a " # k1 Z t l 1 X tl s 1 2 @ 2 pðxj ; t l1 Þ s t l1 1 2 @ 2 pðxj ; t l Þ ds 2 1þ d þ 1þ d þ Oðs Þ 12 x @x2 12 x @x2 CðaÞ l¼1 tl1 s s ðt k1 sÞ1a " ! ! # k Z tl 1 X t l s d2 pðxj ; tl1 Þ s tl1 d2x pðxj ; t l Þ ds 4 4 2 ¼ þ Oðh Þ þ þ Oðh Þ þ Oðs Þ 2 2 CðaÞ l¼1 tl1 s s ðt k sÞ1a h h " ! ! # k1 Z t l 1 X t l s d2 pðxj ; tl1 Þ s tl1 d2x pðxj ; t l Þ ds 4 4 2 þ þ Oð þ Oðh Þ þ Oðh Þ s Þ 2 2 CðaÞ l¼1 tl1 s s ðt k1 sÞ1a h h ¼ la
k X 2 k kðlÞ a dx pðxj ; t kl Þ þ Rj ; l¼0
where
la ¼ Rkj
sa Cð2 þ aÞh2
;
" k Z X 4 2 ¼ Oðs Þ þ Oðh Þ l¼1
tl
ds
t l1
ðtk sÞ1a
k1 Z X l¼1
tl
ds
t l1
ðt k1 sÞ1a
# ;
and
kðlÞ c ¼
8 1; > > > > 1þc < 2 3; 1þc
k ¼ 0; k ¼ 1; 1þc
1þc
1þc
ðl 2Þ þ 3ðl 1Þ 3l þ ðl þ 1Þ ; k ¼ 2; . . . ; k 1; > > > > : ð1 þ cÞkc ðk 1Þc ðk 2Þ1þc þ 2ðk 1Þ1þc k1þc ; k ¼ k;
Similarly, we also have that
"Z # Z tk1 2 k tk X 1 2 1 @ 2 pðxj ; sÞ ds @ pðxj ; sÞ ds ðlÞ ek; ¼ 1þ dx l kb d2x pðxj ; tkl Þ þ R j b 12 @x2 @x2 CðbÞ 0 ðt k sÞ1b ðt k1 sÞ1b 0 l¼0
ð5Þ
Y. Chen, C.-M. Chen / Applied Mathematics and Computation 244 (2014) 772–782
775
where
sb
lb ¼
Cð2 þ bÞh2
;
and
" k Z X 4 k 2 e R j ¼ Oðs Þ þ Oðh Þ l¼1
By
@ 2 f ðx;tÞ @t 2
Z
tl
ds
t l1
ðtk sÞ1b
k1 Z X
l¼1
tl
ds
t l1
ðtk1 sÞ1b
# :
2 CðXÞ, the following Trapezaoidal formula holds
tk
s
f ðxj ; sÞds ¼
2
t k1
f ðxj ; tk1 Þ þ f ðxj ; tk Þ þ Oðs3 Þ:
According the above analysis, we have that
k k X X 1 2 1 2 ðlÞ 2 1þ dx pðxj ; t k Þ 1 þ dx pðxj ; t k1 Þ ¼ la kðlÞ kb d2x pðxj ; tkl Þ a dx pðxj ; t kl Þ þ lb 12 12 l¼0 l¼0 s 1 2 dx f ðxj ; tk1 Þ þ f ðxj ; tk Þ þ Rkj ; 1þ þ 12 2
ð6Þ
where
e k þ Oðs3 Þ Rkj ¼ Rkj þ R j
" k Z X 4 ¼ Oðs2 Þ þ Oðh Þ l¼1
4 þ Oðs2 Þ þ Oðh Þ
tl
ds
t l1
" k Z X l¼1
ðtk sÞ
1a
k1 Z X l¼1
tl
ds
t l1
ðtk sÞ1b
tl
ds
t l1
ðt k1 sÞ1a
k1 Z X l¼1
#
tl
ds
t l1
ðtk1 sÞ1b
# þ Oðs3 Þ:
Notice that k Z X l¼1 k1 Z X l¼1 k Z X l¼1
tl tl1
tl tl1
tl tl1
ds 1a
ðt k sÞ
¼
1a
ds 1b
ðt k sÞ
tk
¼
¼
ds ðtk sÞ
0
ds ðt k1 sÞ
Z
Z
tk1
0
Z
tk
a
a
ds
¼
;
ððk 1ÞsÞ
ðt k1 sÞ
1a
ds
ðksÞ ; b
a
a
;
b
1b
ðt k sÞ
0
ðksÞ
¼
1a
¼
and k1 Z X l¼1
tl tl1
ds ðt k1 sÞ
1b
¼
Z
t k1
0
b
ds 1b
ðtk1 sÞ
¼
ððk 1ÞsÞ ; b
whereas ks 6 T, then 4
Rkj ¼ Oðs2 þ h Þ
ð7Þ
We now present the following numerical scheme for solving problem (2)–(4)
k k X X 1 2 k 1 2 k1 s 1 2 k1 ðlÞ 2 kl 2 kl fj þ fjk ; 1þ dx pj 1 þ dx pj ¼ la kðlÞ d p þ l k d p þ d 1 þ x b x j a x j b 12 12 12 2 l¼0 l¼0 k ¼ 1; 2; . . . ; K; p0j ¼ wj ;
ð8Þ
j ¼ 1; 2; . . . ; J 1;
j ¼ 0; 1; . . . ; J;
pk0 ¼ uðtk Þ; pkJ ¼ wðt k Þ;
ð9Þ k ¼ 0; 1; . . . ; K;
where fjl ¼ f ðxj ; tl Þ; wj ¼ wðxj Þ, whereas plj denote the numerical approximation to the exact solution pðxj ; tl Þ.
ð10Þ
776
Y. Chen, C.-M. Chen / Applied Mathematics and Computation 244 (2014) 772–782
3. Convergence of the numerical scheme for solving MFDE In this section, we study convergence of the numerical scheme (8)–(10). Subtracting (8) from (6), we obtain the following error equation
k k X X 1 2 k 1 2 k1 ðlÞ 2 kl 1þ dx Ej 1 þ dx Ej ¼ la kðlÞ kb d2x Ekl þ Rkj ; a d x E j þ lb j 12 12 l¼0 l¼0 k ¼ 1; 2; . . . ; K;
ð11Þ
j ¼ 1; 2; . . . ; J 1;
where Ekj ¼ uðxj ; t k Þ ukj . For k ¼ 0; 1; . . . ; K, define the following grid functions:
8 i < Ek ; when x 2 x 1 ;x 1 ; j ¼ 1; 2; . . . ; J 1; j2 jþ2 j E ðxÞ ¼ : 0; when x 2 0; h S L h ; L; 2 2 k
and
8 i < Rk ; when x 2 x 1 ;x 1 ; j ¼ 1; 2; . . . ; J 1; j2 jþ2 j R ðxÞ ¼ : 0; when x 2 0; h S L h ; L; 2 2 k
respectively. Then, Ek ðxÞ and Rk ðxÞ can be developed Fourier series, respectively 1 X
Ek ðxÞ ¼
nk ðlÞei2plx=L ;
k ¼ 0; 1; . . . ; K;
l¼1
and
Rk ðxÞ ¼
1 X
gk ðlÞei2plx=L ; k ¼ 0; 1; . . . ; K;
l¼1
where
nk ðlÞ ¼
1 L
Z
L
Ek ðxÞei2plx=L dx;
gk ðlÞ ¼
0
Z
1 L
L
Rk ðxÞei2plx=L dx:
0
Let that
h iT Ek ¼ Ek1 ; Ek2 ; . . . ; EkJ1 ;
h iT Rk ¼ Rk1 ; Rk2 ; . . . ; RkJ1 :
Applying Parseval equalities:
Z
L
jEk ðxÞj2 dx ¼
0
1 X
jnk ðlÞj2 ;
Z
L
jRk ðxÞj2 dx ¼
0
l¼1
1 X
jgk ðlÞj2 ;
l¼1
we have, respectively k
kE k2
J1 X hjEkj j2
!12 ¼
j¼1
kR k2
!12 jnk ðlÞj
2
;
k ¼ 0; 1; . . . ; K
ð12Þ
k ¼ 0; 1; . . . ; K:
ð13Þ
l¼1
and k
1 X
J1 X hjRkj j2
!12 ¼
j¼1
1 X
!12 jgk ðlÞj
2
;
l¼1
We now assume that, respectively
Ekj ¼ nk eirjh ; where
Rkj ¼ gk eirjh ;
r ¼ 2pl=L. By substituting the above expressions into (11), we get k k X X 1 1 2 rh 2 rh 2 rh 2 rh ðlÞ nk 1 sin nk1 ¼ 4la sin 1 sin kðlÞ k n þ gk ; a nkl 4lb sin 3 2 3 2 2 l¼0 2 l¼0 b kl k ¼ 1; 2; . . . ; K;
ð14Þ
Y. Chen, C.-M. Chen / Applied Mathematics and Computation 244 (2014) 772–782
777
We now arrange Eq. (14) as follows
nk ¼
"
1
U
Wnk1 4la sin2
k rh X
2
kðlÞ a nkl 4lb sin
2
l¼2
k rh X
2
# ðlÞ
kb nkl þ gk ;
ð15Þ
l¼2
k ¼ 1; 2; . . . ; K; where
U¼1
1 2 rh 2 rh sin þ 4 la þ lb sin ; 3 2 2
W¼1
1 2 rh 2 rh ð1Þ sin 4 la kð1Þ : sin a þ lb kb 3 2 2
Similar to [5], we have the following lemma: ðlÞ
Lemma 1. If 0 < c < 1, the coefficients kc ðl ¼ 0; 1; . . . ; kÞ possess the following properties: 1þc (1) kð0Þ kð1Þ 3; kðlÞ l ¼ 2; . . . ; k; c ¼ 1; c ¼ 2 c < 0; P c c k ðlÞ (2) k ¼ ð1 þ c Þ k ðk 1Þ > 0; l¼0 c P k ðlÞ ð1Þ (3) 1 < kð1Þ c < 1, l¼2 kc < 1 þ kc .
Lemma 2. For 0 < a < b < 1, if the temporal and spatial steps satisfy
s2 ¼ h4 , then
U P 1: Proof. By 0 < a < b < 1 and
la þ lb ¼
sa Cð2 þ aÞh
2
s2 ¼ h4 , we have that
þ
2ða1Þ
sb 2
Cð2 þ bÞh
¼
2ðb1Þ
h h 1 1 P þ þ ¼1 Cð2 þ aÞ Cð2 þ bÞ Cð3Þ Cð3Þ
from this we complete the proof of Lemma 2. h Lemma 3. For 0 < a < b < 1, if
ð1Þ ð1Þ ð1Þ 1 la kð1Þ a þ lb kb 0, or if 0 < la ka þ lb kb 6 6, then
WP0 From (7), then there is a positive constant C 1 , such that 4
jRkj j 6 C 1 ðs2 þ h Þ;
k ¼ 1; 2; . . . ; K;
ð16Þ
again noticing the first equality of (13), we have
pffiffiffi 4 kRk k2 6 C 1 Lðs2 þ h Þ;
k ¼ 1; 2; . . . ; K:
ð17Þ
In terms of the convergence of the series in RHS of (13), then there is a positive constant C 2 , such that
jgk j jgk ðlÞj 6 C 2 sjg1 ðlÞj C 2 sjg1 j;
k ¼ 1; 2; . . . ; K:
ð18Þ 2
f ðx;tÞ Theorem 1. Suppose that pðx; tÞ 2 PðXÞ is the exact solution of the problem (2)–(4), @ @t 2 CðXÞ, and the temporal and spatial 2 4 ð1Þ ð1Þ ð1Þ ð1Þ 2 1 steps satisfy s ¼ h , if la ka þ lb kb 0, or if 0 < la ka þ lb kb 6 6, then the numerical scheme (8)–(10) is convergent with the 4 order Oðs2 þ h Þ.
Proof. Firstly, we certify the following conclusion: there is a positive constant C, such that
jnk j 6 Ckjg1 j;
k ¼ 1; 2; . . . ; K;
where nk ðk ¼ 1; 2; . . . ; KÞ be the solution of (15). For k ¼ 1, from E0 ¼ 0 lead to n0 ¼ 0, further applying Lemma 2 and (18), from Eq. (15) immediately gets
jn1 j ¼
1
U
jg1 j 6 jg1 j 6 C 2 sjg1 j:
ð19Þ
778
Y. Chen, C.-M. Chen / Applied Mathematics and Computation 244 (2014) 772–782
Suppose that
jnn j 6 C 2 nsjg1 j;
n ¼ 1; 2; . . . ; k 1;
in term of Lemmas 1–3, from Eq. (15), we have that
nk ¼
6
6
6
¼ 6
k k X X 1
2 rh 2 rh ðlÞ kðlÞ n 4 l sin k n þ g
Wnk1 4la sin kl k
b kl a b
2 l¼2 2 l¼2 U
1
Wjnk1 j þ 4la sin
U 1
2
Wðk 1Þ þ 4la sin ("
W þ 4la sin2
U 1
(" 2
W 4la sin
U 1
2
"
U 1
k
rh X
ðlÞ
2
U
l¼2
2
ka jnkl j þ 4lb sin
l¼2
2
ka ðk lÞ þ 4lb sin
k
rh X
ðlÞ
2
rh
2
2
k
rh X
ðlÞ
2
2
ka þ 4lb sin
l¼2 k X
2
ðlÞ
ka 4lb sin
l¼2
2 W þ 4la 1 þ kð1Þ sin a
rh 2
kb jnkl j þ jgk j
l¼2
l¼2
k rh X
2
#
k
rh X
ðlÞ
2
l¼2
#
k
rh X
ðlÞ
2
!
k
rh X
ðlÞ
kb ðk lÞ þ 1 C 2 sjg1 j )
kb ðk 1Þ þ 1 C 2 sjg1 j #
ðlÞ kb
)
ðk 1Þ þ 1 C 2 sjg1 j
l¼2
2 ð1Þ þ 4lb 1 þ kb sin
rh 2
1 ðk 1Þ þ 1 C 2 sjg1 j ¼ k 1 þ C 2 sjg1 j 6 C 2 ksjg1 j
U
This proves the conclusion (19) by induction. Applying the conclusion (19), (12), (13) and (17), we have that 4
kEk k2 6 C 2 kskR1 k2 6 Cðs2 þ h Þ:
pffiffiffi where C ¼ C 1 C 2 LT. This completes the proof of Theorem 1.
h
Remark 1. It is well known that since there are many numerical integration formulas with high order accuracy, thus to continue to construct numerical scheme with high temporal accuracy for solving the modified fractional diffusion Eq. (2) is not difficult. However, the problem is that convergence and stability analysis may be very difficult.
4. Stability and solvability of the numerical scheme for solving MFDE 4.1. Stability of the numerical scheme for solving MFDE In this subsection, we research stability of the numerical scheme (8)–(10), to this end consider the following difference equation
k k X X 1 2 k 1 2 k1 ðlÞ 2 kl 1þ dx qj 1 þ dx qj ¼ la kðlÞ kb d2x qkl j ; a dx qj þ lb 12 12 l¼0 l¼0 k ¼ 1; 2; . . . ; K;
j ¼ 1; 2; . . . ; J 1;
p kj is an approximation for pkj . where q ¼ whereas e For k ¼ 0; 1; . . . ; K, define the following grid function: k j
e p kj ,
pkj
8 i < qk ; when x 2 x 1 ;x 1 ; j ¼ 1; 2; . . . ; J 1; j2 jþ2 j q ðxÞ ¼ : 0; when x 2 0; h S L h ; L; 2 2 k
then qk ðxÞ can be developed Fourier series 1 X
qk ðxÞ ¼
fk ðlÞei2plx=L ;
k ¼ 0; 1; . . . ; K;
l¼1
where
fk ðlÞ ¼
1 L
Z
L
uk ðxÞei2plx=L dx:
0
Let that
h
iT
qk ¼ qk1 ; qk2 ; . . . ; qkJ1 ;
ð20Þ
779
Y. Chen, C.-M. Chen / Applied Mathematics and Computation 244 (2014) 772–782
applying Parseval equality
Z
L
jqk ðxÞj2 dx ¼
0
1 X
jfk ðlÞj2 ;
l¼1
we have that J1 X hjqkj j2
k
kq k2
!12
1 X
¼
j¼1
!12 2
jfk ðlÞj
;
k ¼ 0; 1; . . . ; K:
ð21Þ
l¼1
We now assume that
qkj ¼ fk eirjh ; where
r ¼ 2pl=L. By substituting the above expression into (20) gives k k X X 1 1 2 rh 2 rh 2 rh 2 rh ðlÞ fk 1 sin fk1 ¼ 4la sin 1 sin kðlÞ k f ; a fkl 4lb sin 3 2 3 2 2 l¼0 2 l¼0 b kl
ð22Þ
k ¼ 1; 2; . . . ; K: We rewrite Eq. (22) as follows
fk ¼
"
1
U
Wfk1 4la sin2
k rh X
2
2
kðlÞ a fkl 4lb sin
l¼2
k rh X
2
# ðlÞ
ð23Þ
kb fkl ;
l¼2
k ¼ 1; 2; . . . ; K;
4
ð1Þ
ð1Þ
ð1Þ
ð1Þ
Theorem 2. Suppose that the temporal and spatial steps satisfy s2 ¼ h , if la ka þ lb kb 0, or if 0 < la ka þ lb kb 6 16, then the numerical scheme (8)–(10) is stable.
Proof. At first, we demonstrate the following conclusion:
jfk j 6 jf0 j;
k ¼ 1; 2; . . . ; K;
ð24Þ
where fk ðk ¼ 1; 2; . . . ; KÞ be the solution of (23). For k ¼ 1, by Lemmas 1–3, from Eq. (23) immediately gives
jf1 j ¼
1
U
Wjf0 j 6 jf0 j;
Suppose that
jfn j 6 jf0 j;
n ¼ 1; 2; . . . ; k 1:
Combining Lemmas 1–3, from Eq. (23) gets that
k k X X 1
2 rh 2 rh ðlÞ ðlÞ jfk j ¼ Wfk1 4la sin ka fkl 4lb sin kb fkl
2 l¼2 2 l¼2 U
6
6
1
Wjfk1 j þ 4la sin2
U 1
2
W þ 4la sin
U
rh 2
k
rh X
ðlÞ
2
l¼2
k
X l¼2
2
ka jfkl j þ 4lb sin
ðlÞ
2
ka þ 4lb sin
rh 2
k
rh X
ðlÞ
2 l¼2 !
1
ðlÞ
kb jf0 j ¼
k
X l¼2
!
kb jfkl j
U
2
W 4la sin
1 2 rh 2 rh ð1Þ jf0 j ¼ jf0 j: 6 W þ 4la 1 þ kð1Þ þ 4 l 1 þ k sin sin b a b 2 2 U
k rh X
2
ðlÞ
2
ka 4lb sin
l¼2
k rh X
2
ðlÞ kb
! jf0 j
l¼2
So, the conclusion (24) is proved by induction. Applying the conclusion (24) and (21), then the solution qk of difference Eq. (20) satisfies that
kqk k2 6 kq0 k2 ;
k ¼ 1; 2; . . . ; K:
This prove the Theorem 2.
h
780
Y. Chen, C.-M. Chen / Applied Mathematics and Computation 244 (2014) 772–782
4.2. Solvability of the numerical scheme for solving MFDE In this subsection, we study solvability of the numerical scheme (8)–(10). It can be seen that the corresponding homogeneous linear algebraic equations for the numerical scheme (8)–(10) is
1 2 k 1 2 k1 1þ dx pj 1 þ dx p j 12 12 k k X X ðlÞ 2 kl kb d2x pkl ¼ la kðlÞ j ; a dx pj þ lb l¼0
ð25Þ ð26Þ
l¼0
k ¼ 1; 2; . . . ; K;
j ¼ 1; 2; . . . ; J 1;
p0j ¼ 0;
j ¼ 0; 1; . . . ; J:
pk0 ¼ 0;
pkJ ¼ 0;
ð27Þ
k ¼ 0; 1; . . . ; K:
ð28Þ k
Similar to the proof of Theorem 2, we can also verify that the solution p of the homogeneous linear algebraic Eqs. (26)– (28) satisfies that
kpk k2 6 kp0 k2 ;
k ¼ 1; 2; . . . ; K:
0
By p ¼ 0 immediately gets
pk ¼ 0;
k ¼ 1; 2 . . . ; K:
This reveals that the corresponding homogeneous equations for the numerical scheme (8)–(10) have only zero solution. So, we have the following result: 4
ð1Þ
ð1Þ
ð1Þ
ð1Þ
Theorem 3. Suppose that the temporal and spatial steps satisfy s2 ¼ h , if la ka þ lb kb 0, or if 0 < la ka þ lb kb 6 16, then the numerical scheme (8)–(10) is uniquely solvable.
5. Numerical experiment In order to verify the theoretical analysis results, we now apply the numerical scheme (8)–(10) to solve the following initial-boundary value problem of a modified fractional diffusion equation
@pðx; tÞ ¼ @t 0 < t 6 1;
! @ 1a @ 1b @ 2 pðx; tÞ þ þ f ðx; tÞ; @x2 @t 1a @t 1b 0 < x < 1;
pðx; 0Þ ¼ 0; pð0; tÞ ¼ t 2 ;
ð29Þ
0 6 x 6 1:
ð30Þ
pð1; tÞ ¼ et2 ;
0 6 t 6 1;
ð31Þ
where f ðx; tÞ ¼ 2ex t Cð2þaÞ Cð2þbÞ . The exact solution of the problem (29)–(31) is t 1þa
t 1þb
pðx; tÞ ¼ ex t 2 :
Table 1 The maximum error E1 , temporal accuracy a and spatial accuracy b of the numerical scheme (8)–(10) for solving for the problem (29)–(31).
a
b
0.1
0.7
0.2
0.8
0.3
a
b
3:5759 10
6
a
b
3.87
6
7:4432 10
1.94
3:5591 105
7:4432 106
1.93
2:4953 10
1.90
3.80
3.86
2:4953 106
1.90
0.9
3:5422 105
7:3461 106
3.80
1.94
3.88
2:2508 106
2.06
0.4
0.5
3:5675 105
4.11
7:4432 106
1.93
3.87
2:4272 106
1.95
0.5
0.6
3.90
3:5506 105
7:3461 106
1.94
3.89
2:4272 106
1.92
0.6
3.85
0.7
3:5253 105
7:2489 106
1.95
3.90
2:4272 106
1.90
3.80
0.7
0.71
3:5169 105
7:2489 106
1.95
3.90
2:2508 106
2.03
4.07
0.8
0.82
3:4916 105
7:3461 106
1.93
3.84
2:4272 106
1.92
3.85
0.9
0.93
3:4748 105
7:1571 106
1.95
3.90
2:0066 106
2.21
4.42
1 s2 ¼ h4 ¼ 16
1 s2 ¼ h4 ¼ 81 5
1 s2 ¼ h4 ¼ 256
Y. Chen, C.-M. Chen / Applied Mathematics and Computation 244 (2014) 772–782
781
3 E.S.(t=1/4) N.S.(t=1/4) E.S.(t=2/4) N.S.(t=2/4) E.S.(t=3/4) N.S.(t=3/4) E.S.(t=1) N.S.(t=1)
2.5
2
1.5
1
0.5
0
0
0.2
0.4
x
0.6
0.8
1
Fig. 1. The comparison of the numerical solution (N.S.) using the numerical scheme (8)–(10) and the exact solution (E.S.) for the problem (29)–(31) at 4 1 t ¼ 14 ; 24 ; 34 ; 1, respectively, for a ¼ 0:5; b ¼ 0:6 and s2 ¼ h ¼ 256 .
2.5 E.S.(x=1/4) N.S.(x=1/4) E.S.(x=2/4) N.S.(x=2/4) E.S.(x=3/4) N.S.(x=3/4)
2
1.5
1
0.5
0
0
0.2
0.4
t
0.6
0.8
1
Fig. 2. The comparison of the numerical solution (N.S.) using the numerical scheme (8)–(10) and the exact solution (E.S.) for the problem (29)–(31) at 4 1 x ¼ 14 ; 24 ; 34, respectively, for a ¼ 0:5; b ¼ 0:6 and s2 ¼ h ¼ 256 .
Let that the maximum error of the numerical solutions
n o b E1 ¼ max kEk k2 ¼ Oðsa þ h Þ: 06k6K
ðf Þ Assume that in the front experiment, the spatial step, temporal step, maximum error of the numerical solutions are hf ; sf ; E1 , respectively, whereas in the subsequent experiment, the spatial step, temporal step, maximum error of the numerical solu-
tions are hs ; ss ; EðsÞ 1 , respectively, then the temporal accuracy a and spatial accuracy b calculated by the following formulas
a ¼ log sss f
ðsÞ E1
Eðf1Þ
;
b ¼ log hs hf
ðsÞ E1
Eðf1Þ
:
782
Y. Chen, C.-M. Chen / Applied Mathematics and Computation 244 (2014) 772–782
−6
x 10
Absolute error |E(x,t)|
4 3 2 1 0 1 0.5
t
0
0
0.2
0.4
0.6
0.8
1
x
Fig. 3. The absolute error of the numerical solution for the problem (29)–(31) using the numerical scheme (8)–(10) for a ¼ 0:5; b ¼ 0:6 and
1 s2 ¼ h4 ¼ 256 .
Table 1 explains the maximum error of the numerical solutions for the problem (29)–(31) using the numerical scheme 4 (8)–(10) versus various a; b and s2 ¼ h , respectively. From Table 1, it can be seen that the numerical experiment extremely consistent with our theoretical analysis results. Fig. 1 illustrates the comparison of the numerical solution using the numerical scheme (8)–(10) and the exact solution for 4
1 . Fig. 2 reveals the comparison of the the problem (29)–(31) at t ¼ 14 ; 24 ; 34 ; 1, respectively, for a ¼ 0:5; b ¼ 0:6 and s2 ¼ h ¼ 256 numerical solution using the numerical scheme (8)–(10) and the exact solution for the problem (29)–(31) at x ¼ 14 ; 24 ; 34, 1 s2 ¼ h4 ¼ 256 . Fig. 3 illustrates the absolute error Eðx; tÞ of the numerical solution for 4 1 , where the problem (29)–(31) using the numerical scheme (8)–(10) for a ¼ 0:5; b ¼ 0:6 and s2 ¼ h ¼ 256
respectively, for a ¼ 0:5; b ¼ 0:6 and
Eðx; tÞ ¼ pðx; tÞ e p ðx; tÞ; whereas e p ðx; tÞ is the numerical approximation for pðx; tÞ. Fig. 1–3 show that the numerical solution are very consistent with the exact solution. 6. Conclusion In this paper, a numerical scheme with second order temporal accuracy and fourth order spatial accuracy for solving a modified fractional diffusion equation has been developed; the convergence, stability and solvability of the numerical scheme have been analyzed by Fourier analysis; the numerical experiment extremely consistent with our theoretical results. References [1] A. Mohebbi, M. Abbaszadeh, M. Dehghan, A high-order and unconditionally stable scheme for the modified anomalous fractional sub-diffusion equation with a nonlinear source term, J. Comput. Phys. 240 (2013) 36–48. [2] M. Abbaszadeh, A. Mohebbi, A fourth-order compact solution of the two-dimensional modified anomalous fractional sub-diffusion equation with a nonlinear source term, Comput. Math. Appl. 66 (2013) 1345–1359. [3] A. Chechkin, R. Gorenflo, I. Sokolov, V. Gonchar, Distributed order time fractional diffusion equation, Fract. Calc. Appl. Anal. 6 (3) (2003) 259–279. [4] Chang-Ming Chen, Numerical methods for solving a two-dimensional variable-order modified diffusion equation, Appl. Math. Comput. 225 (2013) 62– 78. [5] H. Huang, X. Cao, Numerical method for two dimensional fractional reaction subdiffusion equation, Eur. Phys. J. Spec. Top. 222 (2013) 1961–1973. [6] T.A.M. Langlands, Solution of a modified fractional diffusion equation, Physica A 367 (2006) 136–144. [7] F. Liu, C. Yang, K. Burrage, Numerical method and analytical technique of the modified anomalous subdiffusion equation with a nonlinear source term, J. Comput. Appl. Math. 231 (2009) 160–176. [8] Q. Liu, F. Liu, I. Turner, V. Anh, Finite element approximation for a modified anomalous subdiffusion equation, Appl. Math. Model. 35 (2011) 4103– 4116. [9] M. Merdan, A. Yildirim, A. Gokdogan, Numerical solution of time-fraction modified equal width wave equation, Eng. Comput. 29 (7–8) (2012) 766–777. [10] I. Sokolov, A. Chechkin, J. Klafter, Distributed-order fractional kinetics, Acta. Phys. Pol. B 35 (2004) 1323–1341. [11] I. Sokolov, J. Klafter, From diffusion to anomalous diffusion: a century after Einstein’s Brownian motion, Chaos 15 (2) (2005) 026103. [12] N. Zhang, W.H. Deng, Y.J. Wu, Finite difference/element method for a two-dimensional modified fractional diffusion equation, Adv. Appl. Math. Mech. 4 (2012) 496–518.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Numerical scheme with high order accuracy for solving a modified fractional diffusion equation Y. Chen, Chang-Ming Chen ⇑ School of Mathematical Sciences, Xiamen University, Xiamen 361006, China
a r t i c l e
i n f o
Keywords: Modified fractional diffusion equation Numerical scheme with second order temporal accuracy and fourth order spatial accuracy Convergence Stability Solvability Fourier analysis
a b s t r a c t In recent years, some researchers have developed various numerical schemes to solve the modified fractional diffusion equation. For the numerical solutions of the modified fractional diffusion equation, there are already some important progresses. In this paper, a numerical scheme with second order temporal accuracy and fourth order spatial accuracy is developed to solve a modified fractional diffusion equation; the convergence, stability and solvability of the numerical scheme are analyzed by Fourier analysis; the theoretical results extremely consistent with the numerical experiment. Ó 2014 Elsevier Inc. All rights reserved.
1. Introduction In recent years, the time or space or time–space fractional diffusion equations are widely used to describe anomalous diffusion processes in numerous physical and biological systems, many authors have developed various numerical schemes to solve these fractional diffusion equations. By the inclusion of a secondary fractional time derivative acting on a diffusion operator, a model for describing processes that become less anomalous as time progresses has been proposed [3,10,11]
@pðx; tÞ ¼ @t
A
! @ 1a @ 1b @ 2 pðx; tÞ þ B ; @x2 @t 1a @t 1b
ð1Þ
for this modified fractional diffusion equation, Langlands et al. proposed the solution with the form of an infinite series of Fox functions on an infinite domain [6]; Merdan et al. presented the numerical solution of time-fraction modified equal width wave equation, and show how an application of fractional two dimensional differential transformation method obtained approximate analytical solution of time-fraction modified equal width wave equation [9]; Liu et al. discussed the numerical method and analytical technique of the modified fractional diffusion equation with a nonlinear source term, they proved that this method has first order temporal accuracy and second order spatial accuracy using a new energy method [7]; Liu et al. researched the finite element approximation for a modified fractional diffusion equation, they proved that this approximation has first order temporal accuracy and m (here m is the degree of the piecewise polynomials) order spatial accuracy [8]; Zhang et al. studied the finite difference/element method for a two-dimensional modified fractional diffusion equation, they proved that this method has ð1 þ minfa; bgÞ order temporal accuracy and m (here m is the degree of the piecewise polynomials) order spatial accuracy [12]; Chen researched the numerical method ⇑ Corresponding author. E-mail address: [email protected] (C.-M. Chen). http://dx.doi.org/10.1016/j.amc.2014.07.044 0096-3003/Ó 2014 Elsevier Inc. All rights reserved.
773
Y. Chen, C.-M. Chen / Applied Mathematics and Computation 244 (2014) 772–782
for solving a two-dimensional variable-order modified fractional diffusion equation, this method has first order temporal accuracy and fourth order spatial accuracy, further, a numerical method for improving temporal accuracy have also been developed [4]; apply the fourth-order compact scheme of the second-order space partial derivative and the Grünwald–Letnikov discretization of the Rieman–Liouille fractional time derivative, Abbaszadeh et al. proposed a high-order and unconditionally stable scheme for the modified anomalous fractional sub-diffusion equation with a nonlinear source term, they proved that this scheme has first order temporal accuracy and fourth order spatial accuracy using the Fourier method [1], and they also presented a fourth-order compact solution of the two-dimensional modified anomalous fractional sub-diffusion equation with a nonlinear source term, they proved that this solution has first order temporal accuracy and fourth order spatial accuracy using the Fourier method [2]. On the basis of the literature, for the cases of one-dimensional and multidimensional of the modified fractional diffusion equation Eq. (1), although the numerical schemes with ideal spatial accuracy can be constructed by the fourth-order compact scheme of the second-order space partial derivative or the finite element approximation, but unfortunately that numerical scheme with second-order temporal accuracy is still unable to achieve. Chen constructed a numerical method for improving temporal accuracy for solving a two-dimensional variable-order modified fractional diffusion equation, but it should be point out that it only is supported by the numerical experiment, whereas it is not supported by rigorous numerical analysis [4]. Therefore to date, numerical scheme with second order temporal accuracy for solving modified fractional diffusions supported by rigorous numerical analysis does not exist. So, how to improve temporal accuracy of numerical scheme for solving modified fractional diffusions is a exciting and meaningful work. Easy to understand that the key problem is that convergence and stability analysis of numerical scheme are very difficult. We expect that this article will change this situation. In this paper, we study numerical scheme and perform related numerical analysis for solving the following modified fractional diffusion equation (MFDE)
! @ 1a @ 1b @ 2 pðx; tÞ þ þ f ðx; tÞ; @x2 @t 1a @t 1b
@pðx; tÞ ¼ @t
ð2Þ
with the following initial and boundary conditions:
pðx; 0Þ ¼ wðxÞ;
0 6 x 6 L;
pð0; tÞ ¼ uðtÞ;
pðL; tÞ ¼ wðtÞ;
where 0 < a < b < 1, and 1c pðx; ;tÞ 0 Dt
¼
1
ð3Þ 0 6 t 6 T;
1c pðx; tÞ 0 Dt
@
Z
CðcÞ @t
t 0
ð4Þ
is the Rieman–Liouille fractional partial derivative of order 1 c defined by
pðx; sÞ ðt sÞ1c
ds:
In this paper, we always suppose that pðx; tÞ 2 PðXÞ is the exact solution of the problem (1)–(3) and
@ 2 f ðx;tÞ @t 2
X ¼ fðx; tÞj0 6 x 6 L; 0 6 t 6 T g; ( PðXÞ ¼
) @ 6 pðx; tÞ @ 4 pðx; y; tÞ pðx; tÞj ; 2 CðXÞ : @x6 @x2 @t2
2. A numerical scheme for solving MFDE In this paper, we let
xj ¼ jh;
j ¼ 0; 1; . . . ; J;
where h ¼ L=J and
d2x
v
i j
¼v
i j1
tk ¼ ks;
k ¼ 0; 1; . . . ; K;
s ¼ T=K is space step and time step, respectively. And define
2v ij þ v ijþ1 :
Integrating both sides of Eq. (2) from t k1 to t k and take x ¼ xj , get that
1 pðxj ; t k Þ pðxj ; tk1 Þ ¼ CðaÞ
"Z
1 þ CðbÞ
tk
0
"Z 0
@ 2 pðxj ; sÞ ds @x2 ðtk sÞ1a tk
Z
t k1
0
@ 2 pðxj ; sÞ ds @x2 ðt k sÞ1b
Z 0
@ 2 pðxj ; sÞ ds @x2 ðt k1 sÞ1a
t k1
#
# Z tk @ 2 pðxj ; sÞ ds þ f ðxj ; sÞds 1b @x2 ðt k1 sÞ t k1
2 CðXÞ, where
774
Y. Chen, C.-M. Chen / Applied Mathematics and Computation 244 (2014) 772–782
from this we have that
"Z # Z tk1 2 tk 2 1 1 1 1 @ pðxj ; sÞ ds @ pðxj ; sÞ ds 1 þ d2x pðxj ;t k Þ 1 þ d2x pðxj ;t k1 Þ ¼ 1 þ d2x 12 12 12 @x2 ðt k sÞ1a @x2 ðt k1 sÞ1a CðaÞ 0 0 " # Z tk 2 Z tk1 2 1 2 1 @ pðxj ;sÞ ds @ pðxj ; sÞ ds þ 1 þ dx 12 @x2 ðt k sÞ1b @x2 ðt k1 sÞ1b CðbÞ 0 0 Z tk 1 f ðxj ;sÞds: þ 1 þ d2x 12 t k1 Because pðx; tÞ 2 PðXÞ gives
1 2 @ 2 pðxj ; t i Þ d2x pðxj ; t i Þ 4 1þ dx ¼ þ Oðh Þ; 2 12 @x2 h and
@ 2 pðxj ; sÞ t l s @ 2 pðxj ; t l1 Þ s tl1 @ 2 pðxj ; t l Þ ¼ þ þ Oðs2 Þ; @x2 @x2 @x2 s s
t l1 6 s 6 tl ;
then
"Z # Z tk1 2 tk 1 2 1 @ 2 pðxj ; sÞ ds @ pðxj ; sÞ ds 1þ dx 12 @x2 @x2 CðaÞ 0 ðt k sÞ1a ðt k1 sÞ1a 0 ! k Z tl 1 2 1 X t l s @ 2 pðxj ; t l1 Þ s t l1 @ 2 pðxj ; tl Þ ds 2 d þ þ Oðs Þ ¼ 1þ 12 x CðaÞ l¼1 tl1 @x2 @x2 s s ðt k sÞ1a ! k1 Z t l 1 2 1 X t l s @ 2 pðxj ; t l1 Þ s t l1 @ 2 pðxj ; tl Þ ds 2 1þ d þ þ Oðs Þ 12 x CðaÞ l¼1 tl1 @x2 @x2 s s ðt k1 sÞ1a " # Z k tl 1 X tl s 1 2 @ 2 pðxj ; t l1 Þ s t l1 1 2 @ 2 pðxj ; t l Þ ds 2 ¼ 1þ dx þ 1 þ d þ Oð s Þ 12 12 x CðaÞ l¼1 tl1 s s @x2 @x2 ðt k sÞ1a " # k1 Z t l 1 X tl s 1 2 @ 2 pðxj ; t l1 Þ s t l1 1 2 @ 2 pðxj ; t l Þ ds 2 1þ d þ 1þ d þ Oðs Þ 12 x @x2 12 x @x2 CðaÞ l¼1 tl1 s s ðt k1 sÞ1a " ! ! # k Z tl 1 X t l s d2 pðxj ; tl1 Þ s tl1 d2x pðxj ; t l Þ ds 4 4 2 ¼ þ Oðh Þ þ þ Oðh Þ þ Oðs Þ 2 2 CðaÞ l¼1 tl1 s s ðt k sÞ1a h h " ! ! # k1 Z t l 1 X t l s d2 pðxj ; tl1 Þ s tl1 d2x pðxj ; t l Þ ds 4 4 2 þ þ Oð þ Oðh Þ þ Oðh Þ s Þ 2 2 CðaÞ l¼1 tl1 s s ðt k1 sÞ1a h h ¼ la
k X 2 k kðlÞ a dx pðxj ; t kl Þ þ Rj ; l¼0
where
la ¼ Rkj
sa Cð2 þ aÞh2
;
" k Z X 4 2 ¼ Oðs Þ þ Oðh Þ l¼1
tl
ds
t l1
ðtk sÞ1a
k1 Z X l¼1
tl
ds
t l1
ðt k1 sÞ1a
# ;
and
kðlÞ c ¼
8 1; > > > > 1þc < 2 3; 1þc
k ¼ 0; k ¼ 1; 1þc
1þc
1þc
ðl 2Þ þ 3ðl 1Þ 3l þ ðl þ 1Þ ; k ¼ 2; . . . ; k 1; > > > > : ð1 þ cÞkc ðk 1Þc ðk 2Þ1þc þ 2ðk 1Þ1þc k1þc ; k ¼ k;
Similarly, we also have that
"Z # Z tk1 2 k tk X 1 2 1 @ 2 pðxj ; sÞ ds @ pðxj ; sÞ ds ðlÞ ek; ¼ 1þ dx l kb d2x pðxj ; tkl Þ þ R j b 12 @x2 @x2 CðbÞ 0 ðt k sÞ1b ðt k1 sÞ1b 0 l¼0
ð5Þ
Y. Chen, C.-M. Chen / Applied Mathematics and Computation 244 (2014) 772–782
775
where
sb
lb ¼
Cð2 þ bÞh2
;
and
" k Z X 4 k 2 e R j ¼ Oðs Þ þ Oðh Þ l¼1
By
@ 2 f ðx;tÞ @t 2
Z
tl
ds
t l1
ðtk sÞ1b
k1 Z X
l¼1
tl
ds
t l1
ðtk1 sÞ1b
# :
2 CðXÞ, the following Trapezaoidal formula holds
tk
s
f ðxj ; sÞds ¼
2
t k1
f ðxj ; tk1 Þ þ f ðxj ; tk Þ þ Oðs3 Þ:
According the above analysis, we have that
k k X X 1 2 1 2 ðlÞ 2 1þ dx pðxj ; t k Þ 1 þ dx pðxj ; t k1 Þ ¼ la kðlÞ kb d2x pðxj ; tkl Þ a dx pðxj ; t kl Þ þ lb 12 12 l¼0 l¼0 s 1 2 dx f ðxj ; tk1 Þ þ f ðxj ; tk Þ þ Rkj ; 1þ þ 12 2
ð6Þ
where
e k þ Oðs3 Þ Rkj ¼ Rkj þ R j
" k Z X 4 ¼ Oðs2 Þ þ Oðh Þ l¼1
4 þ Oðs2 Þ þ Oðh Þ
tl
ds
t l1
" k Z X l¼1
ðtk sÞ
1a
k1 Z X l¼1
tl
ds
t l1
ðtk sÞ1b
tl
ds
t l1
ðt k1 sÞ1a
k1 Z X l¼1
#
tl
ds
t l1
ðtk1 sÞ1b
# þ Oðs3 Þ:
Notice that k Z X l¼1 k1 Z X l¼1 k Z X l¼1
tl tl1
tl tl1
tl tl1
ds 1a
ðt k sÞ
¼
1a
ds 1b
ðt k sÞ
tk
¼
¼
ds ðtk sÞ
0
ds ðt k1 sÞ
Z
Z
tk1
0
Z
tk
a
a
ds
¼
;
ððk 1ÞsÞ
ðt k1 sÞ
1a
ds
ðksÞ ; b
a
a
;
b
1b
ðt k sÞ
0
ðksÞ
¼
1a
¼
and k1 Z X l¼1
tl tl1
ds ðt k1 sÞ
1b
¼
Z
t k1
0
b
ds 1b
ðtk1 sÞ
¼
ððk 1ÞsÞ ; b
whereas ks 6 T, then 4
Rkj ¼ Oðs2 þ h Þ
ð7Þ
We now present the following numerical scheme for solving problem (2)–(4)
k k X X 1 2 k 1 2 k1 s 1 2 k1 ðlÞ 2 kl 2 kl fj þ fjk ; 1þ dx pj 1 þ dx pj ¼ la kðlÞ d p þ l k d p þ d 1 þ x b x j a x j b 12 12 12 2 l¼0 l¼0 k ¼ 1; 2; . . . ; K; p0j ¼ wj ;
ð8Þ
j ¼ 1; 2; . . . ; J 1;
j ¼ 0; 1; . . . ; J;
pk0 ¼ uðtk Þ; pkJ ¼ wðt k Þ;
ð9Þ k ¼ 0; 1; . . . ; K;
where fjl ¼ f ðxj ; tl Þ; wj ¼ wðxj Þ, whereas plj denote the numerical approximation to the exact solution pðxj ; tl Þ.
ð10Þ
776
Y. Chen, C.-M. Chen / Applied Mathematics and Computation 244 (2014) 772–782
3. Convergence of the numerical scheme for solving MFDE In this section, we study convergence of the numerical scheme (8)–(10). Subtracting (8) from (6), we obtain the following error equation
k k X X 1 2 k 1 2 k1 ðlÞ 2 kl 1þ dx Ej 1 þ dx Ej ¼ la kðlÞ kb d2x Ekl þ Rkj ; a d x E j þ lb j 12 12 l¼0 l¼0 k ¼ 1; 2; . . . ; K;
ð11Þ
j ¼ 1; 2; . . . ; J 1;
where Ekj ¼ uðxj ; t k Þ ukj . For k ¼ 0; 1; . . . ; K, define the following grid functions:
8 i < Ek ; when x 2 x 1 ;x 1 ; j ¼ 1; 2; . . . ; J 1; j2 jþ2 j E ðxÞ ¼ : 0; when x 2 0; h S L h ; L; 2 2 k
and
8 i < Rk ; when x 2 x 1 ;x 1 ; j ¼ 1; 2; . . . ; J 1; j2 jþ2 j R ðxÞ ¼ : 0; when x 2 0; h S L h ; L; 2 2 k
respectively. Then, Ek ðxÞ and Rk ðxÞ can be developed Fourier series, respectively 1 X
Ek ðxÞ ¼
nk ðlÞei2plx=L ;
k ¼ 0; 1; . . . ; K;
l¼1
and
Rk ðxÞ ¼
1 X
gk ðlÞei2plx=L ; k ¼ 0; 1; . . . ; K;
l¼1
where
nk ðlÞ ¼
1 L
Z
L
Ek ðxÞei2plx=L dx;
gk ðlÞ ¼
0
Z
1 L
L
Rk ðxÞei2plx=L dx:
0
Let that
h iT Ek ¼ Ek1 ; Ek2 ; . . . ; EkJ1 ;
h iT Rk ¼ Rk1 ; Rk2 ; . . . ; RkJ1 :
Applying Parseval equalities:
Z
L
jEk ðxÞj2 dx ¼
0
1 X
jnk ðlÞj2 ;
Z
L
jRk ðxÞj2 dx ¼
0
l¼1
1 X
jgk ðlÞj2 ;
l¼1
we have, respectively k
kE k2
J1 X hjEkj j2
!12 ¼
j¼1
kR k2
!12 jnk ðlÞj
2
;
k ¼ 0; 1; . . . ; K
ð12Þ
k ¼ 0; 1; . . . ; K:
ð13Þ
l¼1
and k
1 X
J1 X hjRkj j2
!12 ¼
j¼1
1 X
!12 jgk ðlÞj
2
;
l¼1
We now assume that, respectively
Ekj ¼ nk eirjh ; where
Rkj ¼ gk eirjh ;
r ¼ 2pl=L. By substituting the above expressions into (11), we get k k X X 1 1 2 rh 2 rh 2 rh 2 rh ðlÞ nk 1 sin nk1 ¼ 4la sin 1 sin kðlÞ k n þ gk ; a nkl 4lb sin 3 2 3 2 2 l¼0 2 l¼0 b kl k ¼ 1; 2; . . . ; K;
ð14Þ
Y. Chen, C.-M. Chen / Applied Mathematics and Computation 244 (2014) 772–782
777
We now arrange Eq. (14) as follows
nk ¼
"
1
U
Wnk1 4la sin2
k rh X
2
kðlÞ a nkl 4lb sin
2
l¼2
k rh X
2
# ðlÞ
kb nkl þ gk ;
ð15Þ
l¼2
k ¼ 1; 2; . . . ; K; where
U¼1
1 2 rh 2 rh sin þ 4 la þ lb sin ; 3 2 2
W¼1
1 2 rh 2 rh ð1Þ sin 4 la kð1Þ : sin a þ lb kb 3 2 2
Similar to [5], we have the following lemma: ðlÞ
Lemma 1. If 0 < c < 1, the coefficients kc ðl ¼ 0; 1; . . . ; kÞ possess the following properties: 1þc (1) kð0Þ kð1Þ 3; kðlÞ l ¼ 2; . . . ; k; c ¼ 1; c ¼ 2 c < 0; P c c k ðlÞ (2) k ¼ ð1 þ c Þ k ðk 1Þ > 0; l¼0 c P k ðlÞ ð1Þ (3) 1 < kð1Þ c < 1, l¼2 kc < 1 þ kc .
Lemma 2. For 0 < a < b < 1, if the temporal and spatial steps satisfy
s2 ¼ h4 , then
U P 1: Proof. By 0 < a < b < 1 and
la þ lb ¼
sa Cð2 þ aÞh
2
s2 ¼ h4 , we have that
þ
2ða1Þ
sb 2
Cð2 þ bÞh
¼
2ðb1Þ
h h 1 1 P þ þ ¼1 Cð2 þ aÞ Cð2 þ bÞ Cð3Þ Cð3Þ
from this we complete the proof of Lemma 2. h Lemma 3. For 0 < a < b < 1, if
ð1Þ ð1Þ ð1Þ 1 la kð1Þ a þ lb kb 0, or if 0 < la ka þ lb kb 6 6, then
WP0 From (7), then there is a positive constant C 1 , such that 4
jRkj j 6 C 1 ðs2 þ h Þ;
k ¼ 1; 2; . . . ; K;
ð16Þ
again noticing the first equality of (13), we have
pffiffiffi 4 kRk k2 6 C 1 Lðs2 þ h Þ;
k ¼ 1; 2; . . . ; K:
ð17Þ
In terms of the convergence of the series in RHS of (13), then there is a positive constant C 2 , such that
jgk j jgk ðlÞj 6 C 2 sjg1 ðlÞj C 2 sjg1 j;
k ¼ 1; 2; . . . ; K:
ð18Þ 2
f ðx;tÞ Theorem 1. Suppose that pðx; tÞ 2 PðXÞ is the exact solution of the problem (2)–(4), @ @t 2 CðXÞ, and the temporal and spatial 2 4 ð1Þ ð1Þ ð1Þ ð1Þ 2 1 steps satisfy s ¼ h , if la ka þ lb kb 0, or if 0 < la ka þ lb kb 6 6, then the numerical scheme (8)–(10) is convergent with the 4 order Oðs2 þ h Þ.
Proof. Firstly, we certify the following conclusion: there is a positive constant C, such that
jnk j 6 Ckjg1 j;
k ¼ 1; 2; . . . ; K;
where nk ðk ¼ 1; 2; . . . ; KÞ be the solution of (15). For k ¼ 1, from E0 ¼ 0 lead to n0 ¼ 0, further applying Lemma 2 and (18), from Eq. (15) immediately gets
jn1 j ¼
1
U
jg1 j 6 jg1 j 6 C 2 sjg1 j:
ð19Þ
778
Y. Chen, C.-M. Chen / Applied Mathematics and Computation 244 (2014) 772–782
Suppose that
jnn j 6 C 2 nsjg1 j;
n ¼ 1; 2; . . . ; k 1;
in term of Lemmas 1–3, from Eq. (15), we have that
nk ¼
6
6
6
¼ 6
k k X X 1
2 rh 2 rh ðlÞ kðlÞ n 4 l sin k n þ g
Wnk1 4la sin kl k
b kl a b
2 l¼2 2 l¼2 U
1
Wjnk1 j þ 4la sin
U 1
2
Wðk 1Þ þ 4la sin ("
W þ 4la sin2
U 1
(" 2
W 4la sin
U 1
2
"
U 1
k
rh X
ðlÞ
2
U
l¼2
2
ka jnkl j þ 4lb sin
l¼2
2
ka ðk lÞ þ 4lb sin
k
rh X
ðlÞ
2
rh
2
2
k
rh X
ðlÞ
2
2
ka þ 4lb sin
l¼2 k X
2
ðlÞ
ka 4lb sin
l¼2
2 W þ 4la 1 þ kð1Þ sin a
rh 2
kb jnkl j þ jgk j
l¼2
l¼2
k rh X
2
#
k
rh X
ðlÞ
2
l¼2
#
k
rh X
ðlÞ
2
!
k
rh X
ðlÞ
kb ðk lÞ þ 1 C 2 sjg1 j )
kb ðk 1Þ þ 1 C 2 sjg1 j #
ðlÞ kb
)
ðk 1Þ þ 1 C 2 sjg1 j
l¼2
2 ð1Þ þ 4lb 1 þ kb sin
rh 2
1 ðk 1Þ þ 1 C 2 sjg1 j ¼ k 1 þ C 2 sjg1 j 6 C 2 ksjg1 j
U
This proves the conclusion (19) by induction. Applying the conclusion (19), (12), (13) and (17), we have that 4
kEk k2 6 C 2 kskR1 k2 6 Cðs2 þ h Þ:
pffiffiffi where C ¼ C 1 C 2 LT. This completes the proof of Theorem 1.
h
Remark 1. It is well known that since there are many numerical integration formulas with high order accuracy, thus to continue to construct numerical scheme with high temporal accuracy for solving the modified fractional diffusion Eq. (2) is not difficult. However, the problem is that convergence and stability analysis may be very difficult.
4. Stability and solvability of the numerical scheme for solving MFDE 4.1. Stability of the numerical scheme for solving MFDE In this subsection, we research stability of the numerical scheme (8)–(10), to this end consider the following difference equation
k k X X 1 2 k 1 2 k1 ðlÞ 2 kl 1þ dx qj 1 þ dx qj ¼ la kðlÞ kb d2x qkl j ; a dx qj þ lb 12 12 l¼0 l¼0 k ¼ 1; 2; . . . ; K;
j ¼ 1; 2; . . . ; J 1;
p kj is an approximation for pkj . where q ¼ whereas e For k ¼ 0; 1; . . . ; K, define the following grid function: k j
e p kj ,
pkj
8 i < qk ; when x 2 x 1 ;x 1 ; j ¼ 1; 2; . . . ; J 1; j2 jþ2 j q ðxÞ ¼ : 0; when x 2 0; h S L h ; L; 2 2 k
then qk ðxÞ can be developed Fourier series 1 X
qk ðxÞ ¼
fk ðlÞei2plx=L ;
k ¼ 0; 1; . . . ; K;
l¼1
where
fk ðlÞ ¼
1 L
Z
L
uk ðxÞei2plx=L dx:
0
Let that
h
iT
qk ¼ qk1 ; qk2 ; . . . ; qkJ1 ;
ð20Þ
779
Y. Chen, C.-M. Chen / Applied Mathematics and Computation 244 (2014) 772–782
applying Parseval equality
Z
L
jqk ðxÞj2 dx ¼
0
1 X
jfk ðlÞj2 ;
l¼1
we have that J1 X hjqkj j2
k
kq k2
!12
1 X
¼
j¼1
!12 2
jfk ðlÞj
;
k ¼ 0; 1; . . . ; K:
ð21Þ
l¼1
We now assume that
qkj ¼ fk eirjh ; where
r ¼ 2pl=L. By substituting the above expression into (20) gives k k X X 1 1 2 rh 2 rh 2 rh 2 rh ðlÞ fk 1 sin fk1 ¼ 4la sin 1 sin kðlÞ k f ; a fkl 4lb sin 3 2 3 2 2 l¼0 2 l¼0 b kl
ð22Þ
k ¼ 1; 2; . . . ; K: We rewrite Eq. (22) as follows
fk ¼
"
1
U
Wfk1 4la sin2
k rh X
2
2
kðlÞ a fkl 4lb sin
l¼2
k rh X
2
# ðlÞ
ð23Þ
kb fkl ;
l¼2
k ¼ 1; 2; . . . ; K;
4
ð1Þ
ð1Þ
ð1Þ
ð1Þ
Theorem 2. Suppose that the temporal and spatial steps satisfy s2 ¼ h , if la ka þ lb kb 0, or if 0 < la ka þ lb kb 6 16, then the numerical scheme (8)–(10) is stable.
Proof. At first, we demonstrate the following conclusion:
jfk j 6 jf0 j;
k ¼ 1; 2; . . . ; K;
ð24Þ
where fk ðk ¼ 1; 2; . . . ; KÞ be the solution of (23). For k ¼ 1, by Lemmas 1–3, from Eq. (23) immediately gives
jf1 j ¼
1
U
Wjf0 j 6 jf0 j;
Suppose that
jfn j 6 jf0 j;
n ¼ 1; 2; . . . ; k 1:
Combining Lemmas 1–3, from Eq. (23) gets that
k k X X 1
2 rh 2 rh ðlÞ ðlÞ jfk j ¼ Wfk1 4la sin ka fkl 4lb sin kb fkl
2 l¼2 2 l¼2 U
6
6
1
Wjfk1 j þ 4la sin2
U 1
2
W þ 4la sin
U
rh 2
k
rh X
ðlÞ
2
l¼2
k
X l¼2
2
ka jfkl j þ 4lb sin
ðlÞ
2
ka þ 4lb sin
rh 2
k
rh X
ðlÞ
2 l¼2 !
1
ðlÞ
kb jf0 j ¼
k
X l¼2
!
kb jfkl j
U
2
W 4la sin
1 2 rh 2 rh ð1Þ jf0 j ¼ jf0 j: 6 W þ 4la 1 þ kð1Þ þ 4 l 1 þ k sin sin b a b 2 2 U
k rh X
2
ðlÞ
2
ka 4lb sin
l¼2
k rh X
2
ðlÞ kb
! jf0 j
l¼2
So, the conclusion (24) is proved by induction. Applying the conclusion (24) and (21), then the solution qk of difference Eq. (20) satisfies that
kqk k2 6 kq0 k2 ;
k ¼ 1; 2; . . . ; K:
This prove the Theorem 2.
h
780
Y. Chen, C.-M. Chen / Applied Mathematics and Computation 244 (2014) 772–782
4.2. Solvability of the numerical scheme for solving MFDE In this subsection, we study solvability of the numerical scheme (8)–(10). It can be seen that the corresponding homogeneous linear algebraic equations for the numerical scheme (8)–(10) is
1 2 k 1 2 k1 1þ dx pj 1 þ dx p j 12 12 k k X X ðlÞ 2 kl kb d2x pkl ¼ la kðlÞ j ; a dx pj þ lb l¼0
ð25Þ ð26Þ
l¼0
k ¼ 1; 2; . . . ; K;
j ¼ 1; 2; . . . ; J 1;
p0j ¼ 0;
j ¼ 0; 1; . . . ; J:
pk0 ¼ 0;
pkJ ¼ 0;
ð27Þ
k ¼ 0; 1; . . . ; K:
ð28Þ k
Similar to the proof of Theorem 2, we can also verify that the solution p of the homogeneous linear algebraic Eqs. (26)– (28) satisfies that
kpk k2 6 kp0 k2 ;
k ¼ 1; 2; . . . ; K:
0
By p ¼ 0 immediately gets
pk ¼ 0;
k ¼ 1; 2 . . . ; K:
This reveals that the corresponding homogeneous equations for the numerical scheme (8)–(10) have only zero solution. So, we have the following result: 4
ð1Þ
ð1Þ
ð1Þ
ð1Þ
Theorem 3. Suppose that the temporal and spatial steps satisfy s2 ¼ h , if la ka þ lb kb 0, or if 0 < la ka þ lb kb 6 16, then the numerical scheme (8)–(10) is uniquely solvable.
5. Numerical experiment In order to verify the theoretical analysis results, we now apply the numerical scheme (8)–(10) to solve the following initial-boundary value problem of a modified fractional diffusion equation
@pðx; tÞ ¼ @t 0 < t 6 1;
! @ 1a @ 1b @ 2 pðx; tÞ þ þ f ðx; tÞ; @x2 @t 1a @t 1b 0 < x < 1;
pðx; 0Þ ¼ 0; pð0; tÞ ¼ t 2 ;
ð29Þ
0 6 x 6 1:
ð30Þ
pð1; tÞ ¼ et2 ;
0 6 t 6 1;
ð31Þ
where f ðx; tÞ ¼ 2ex t Cð2þaÞ Cð2þbÞ . The exact solution of the problem (29)–(31) is t 1þa
t 1þb
pðx; tÞ ¼ ex t 2 :
Table 1 The maximum error E1 , temporal accuracy a and spatial accuracy b of the numerical scheme (8)–(10) for solving for the problem (29)–(31).
a
b
0.1
0.7
0.2
0.8
0.3
a
b
3:5759 10
6
a
b
3.87
6
7:4432 10
1.94
3:5591 105
7:4432 106
1.93
2:4953 10
1.90
3.80
3.86
2:4953 106
1.90
0.9
3:5422 105
7:3461 106
3.80
1.94
3.88
2:2508 106
2.06
0.4
0.5
3:5675 105
4.11
7:4432 106
1.93
3.87
2:4272 106
1.95
0.5
0.6
3.90
3:5506 105
7:3461 106
1.94
3.89
2:4272 106
1.92
0.6
3.85
0.7
3:5253 105
7:2489 106
1.95
3.90
2:4272 106
1.90
3.80
0.7
0.71
3:5169 105
7:2489 106
1.95
3.90
2:2508 106
2.03
4.07
0.8
0.82
3:4916 105
7:3461 106
1.93
3.84
2:4272 106
1.92
3.85
0.9
0.93
3:4748 105
7:1571 106
1.95
3.90
2:0066 106
2.21
4.42
1 s2 ¼ h4 ¼ 16
1 s2 ¼ h4 ¼ 81 5
1 s2 ¼ h4 ¼ 256
Y. Chen, C.-M. Chen / Applied Mathematics and Computation 244 (2014) 772–782
781
3 E.S.(t=1/4) N.S.(t=1/4) E.S.(t=2/4) N.S.(t=2/4) E.S.(t=3/4) N.S.(t=3/4) E.S.(t=1) N.S.(t=1)
2.5
2
1.5
1
0.5
0
0
0.2
0.4
x
0.6
0.8
1
Fig. 1. The comparison of the numerical solution (N.S.) using the numerical scheme (8)–(10) and the exact solution (E.S.) for the problem (29)–(31) at 4 1 t ¼ 14 ; 24 ; 34 ; 1, respectively, for a ¼ 0:5; b ¼ 0:6 and s2 ¼ h ¼ 256 .
2.5 E.S.(x=1/4) N.S.(x=1/4) E.S.(x=2/4) N.S.(x=2/4) E.S.(x=3/4) N.S.(x=3/4)
2
1.5
1
0.5
0
0
0.2
0.4
t
0.6
0.8
1
Fig. 2. The comparison of the numerical solution (N.S.) using the numerical scheme (8)–(10) and the exact solution (E.S.) for the problem (29)–(31) at 4 1 x ¼ 14 ; 24 ; 34, respectively, for a ¼ 0:5; b ¼ 0:6 and s2 ¼ h ¼ 256 .
Let that the maximum error of the numerical solutions
n o b E1 ¼ max kEk k2 ¼ Oðsa þ h Þ: 06k6K
ðf Þ Assume that in the front experiment, the spatial step, temporal step, maximum error of the numerical solutions are hf ; sf ; E1 , respectively, whereas in the subsequent experiment, the spatial step, temporal step, maximum error of the numerical solu-
tions are hs ; ss ; EðsÞ 1 , respectively, then the temporal accuracy a and spatial accuracy b calculated by the following formulas
a ¼ log sss f
ðsÞ E1
Eðf1Þ
;
b ¼ log hs hf
ðsÞ E1
Eðf1Þ
:
782
Y. Chen, C.-M. Chen / Applied Mathematics and Computation 244 (2014) 772–782
−6
x 10
Absolute error |E(x,t)|
4 3 2 1 0 1 0.5
t
0
0
0.2
0.4
0.6
0.8
1
x
Fig. 3. The absolute error of the numerical solution for the problem (29)–(31) using the numerical scheme (8)–(10) for a ¼ 0:5; b ¼ 0:6 and
1 s2 ¼ h4 ¼ 256 .
Table 1 explains the maximum error of the numerical solutions for the problem (29)–(31) using the numerical scheme 4 (8)–(10) versus various a; b and s2 ¼ h , respectively. From Table 1, it can be seen that the numerical experiment extremely consistent with our theoretical analysis results. Fig. 1 illustrates the comparison of the numerical solution using the numerical scheme (8)–(10) and the exact solution for 4
1 . Fig. 2 reveals the comparison of the the problem (29)–(31) at t ¼ 14 ; 24 ; 34 ; 1, respectively, for a ¼ 0:5; b ¼ 0:6 and s2 ¼ h ¼ 256 numerical solution using the numerical scheme (8)–(10) and the exact solution for the problem (29)–(31) at x ¼ 14 ; 24 ; 34, 1 s2 ¼ h4 ¼ 256 . Fig. 3 illustrates the absolute error Eðx; tÞ of the numerical solution for 4 1 , where the problem (29)–(31) using the numerical scheme (8)–(10) for a ¼ 0:5; b ¼ 0:6 and s2 ¼ h ¼ 256
respectively, for a ¼ 0:5; b ¼ 0:6 and
Eðx; tÞ ¼ pðx; tÞ e p ðx; tÞ; whereas e p ðx; tÞ is the numerical approximation for pðx; tÞ. Fig. 1–3 show that the numerical solution are very consistent with the exact solution. 6. Conclusion In this paper, a numerical scheme with second order temporal accuracy and fourth order spatial accuracy for solving a modified fractional diffusion equation has been developed; the convergence, stability and solvability of the numerical scheme have been analyzed by Fourier analysis; the numerical experiment extremely consistent with our theoretical results. References [1] A. Mohebbi, M. Abbaszadeh, M. Dehghan, A high-order and unconditionally stable scheme for the modified anomalous fractional sub-diffusion equation with a nonlinear source term, J. Comput. Phys. 240 (2013) 36–48. [2] M. Abbaszadeh, A. Mohebbi, A fourth-order compact solution of the two-dimensional modified anomalous fractional sub-diffusion equation with a nonlinear source term, Comput. Math. Appl. 66 (2013) 1345–1359. [3] A. Chechkin, R. Gorenflo, I. Sokolov, V. Gonchar, Distributed order time fractional diffusion equation, Fract. Calc. Appl. Anal. 6 (3) (2003) 259–279. [4] Chang-Ming Chen, Numerical methods for solving a two-dimensional variable-order modified diffusion equation, Appl. Math. Comput. 225 (2013) 62– 78. [5] H. Huang, X. Cao, Numerical method for two dimensional fractional reaction subdiffusion equation, Eur. Phys. J. Spec. Top. 222 (2013) 1961–1973. [6] T.A.M. Langlands, Solution of a modified fractional diffusion equation, Physica A 367 (2006) 136–144. [7] F. Liu, C. Yang, K. Burrage, Numerical method and analytical technique of the modified anomalous subdiffusion equation with a nonlinear source term, J. Comput. Appl. Math. 231 (2009) 160–176. [8] Q. Liu, F. Liu, I. Turner, V. Anh, Finite element approximation for a modified anomalous subdiffusion equation, Appl. Math. Model. 35 (2011) 4103– 4116. [9] M. Merdan, A. Yildirim, A. Gokdogan, Numerical solution of time-fraction modified equal width wave equation, Eng. Comput. 29 (7–8) (2012) 766–777. [10] I. Sokolov, A. Chechkin, J. Klafter, Distributed-order fractional kinetics, Acta. Phys. Pol. B 35 (2004) 1323–1341. [11] I. Sokolov, J. Klafter, From diffusion to anomalous diffusion: a century after Einstein’s Brownian motion, Chaos 15 (2) (2005) 026103. [12] N. Zhang, W.H. Deng, Y.J. Wu, Finite difference/element method for a two-dimensional modified fractional diffusion equation, Adv. Appl. Math. Mech. 4 (2012) 496–518.