- Email: [email protected]
Numerical simulation with high order accuracy for the time fractional reaction–subdiffusion equation
Accepted Manuscript Numerical simulation with high order accuracy for the time fractional reaction-subdiffusion equation Y. Chen, Chang-Ming Chen PII: DOI: Reference:
S0378-4754(17)30084-8 http://dx.doi.org/10.1016/j.matcom.2017.03.008 MATCOM 4445
To appear in:
Mathematics and Computers in Simulation
Received date: 2 January 2016 Revised date: 28 December 2016 Accepted date: 12 March 2017 Please cite this article as: Y. Chen, C.-M. Chen, Numerical simulation with high order accuracy for the time fractional reaction-subdiffusion equation, Math. Comput. Simulation (2017), http://dx.doi.org/10.1016/j.matcom.2017.03.008 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Numerical simulation with high order accuracy for the time fractional reaction-subdiffusion equation Y. Chen a , Chang-Ming.Chen a,∗ a School
of Mathematical Sciences, Xiamen University, Xiamen 361006, China
Abstract In this paper, we propose the numerical simulation method with second order temporal accuracy and fourth order spatial accuracy for the time fractional reactionsubdiffusion equation; the stability, convergence and solvability of the numerical simulation method respectively are discussed by Fourier analysis and algebraic theory; the theoretical analysis results very consistent with the numerical experiment. Key words: The time fractional reaction-subdiffusion equation, The numerical simulation method with second order temporal accuracy and fourth order spatial accuracy, Stability, Solvability, Convergence.
1
Introduction
Fractional reaction-subdiffusion equations are noticed and researched in recent years in many areas of science and engineering [38,32,39,33,36,29,30,20,14,42,28,25,9,22,26,31,15–18,40,23,21,4,1,19]. At the same time, some authors have make out a lot of research on numerical simulation and numerical analysis of fractional reaction-subdiffusion equations, and this research continues. For one-dimensional fractional reaction-subdiffusion equation, based on operator splitting, Baeumer et al. [3] proposed a practical method for numerical solution; Chen et al. [6] developed the finite difference methods, by Fourier analysis, they also verified that these finite difference methods have first order ∗ Corresponding author. Email address: [email protected] (Chang-Ming.Chen).
Preprint submitted to Elsevier Science
4 April 2017
temporal accuracy, second order spatial accuracy and unconditionally stable; Wang et al. [35] presented a variational iteration method; Zhuang et al. [41] developed an implicit numerical method, they also demonstrated that this implicit numerical method has first order temporal accuracy, second order spatial accuracy and unconditionally stable; Rida et al. [27] studied a new application of generalized differential transform method; Das et al. [10] presented an approximate solution; Ding et al. [13] proposed the two classes of difference schemes, they also verified that one scheme has first order temporal accuracy, second order spatial accuracy and unconditionally stable, whereas another scheme has first order temporal accuracy, fourth order spatial accuracy and unconditionally stable; Sweilam et al. [34] presented the weighted average finite difference methods, but they only given stable analysis of the methods; Doha et al. [11] constructed the efficient Legendre spectral tau matrix formulation. For multi-dimensional fractional reaction-subdiffusion equation, Huang et al. [24] proposed a numerical method, they also verified that this method has (1 + γ) order temporal accuracy, second order spatial accuracy and unconditionally stable; Dehghan et al. [12] proposed a numerical technique based on a meshless method, this technique has γ order temporal accuracy and unconditionally stable; based on a new combination of alternating direction implicit (ADI) approach and interpolating element free Galerkin (EFG) method, Abbaszadeh et al. [2] researched a meshless numerical procedure, the procedure has (1 + γ) order temporal accuracy and unconditionally stable; Yu et al. [37] developed the novel compact numerical method, the method has second order temporal accuracy and fourth order spatial accuracy. For one-dimensional variable order reaction-subdiffusion equation, Chen et al. [8] discussed the numerical approximation, they also verified that this numerical approximation method has first order temporal accuracy, second order spatial accuracy and unconditionally stable, and they also developed a numerical approximation method for improving temporal accuracy. According to the literature, so far, although there are some numerical solving techniques for improving temporal accuracy of fractional reactionsubdiffusion equation, but they are only supported by the numerical experiments whereas lack of support for strict numerical analysis. Therefore, how to structure numerical simulation method with high temporal accuracy for solving fractional reaction-subdiffusion equation is very meaningful. In this paper, we research the following time fractional reaction-subdiffusion equation: ∂u(x, t) =0 Dt1−γ ∂t
!
∂ 2 u(x, t) − u(x, t) + f (x, t), 0 < t ≤ T, 0 < x < L, (1) ∂x2 2
with the following initial and boundary conditions: u(x, 0) = φ(x), 0 ≤ x ≤ L,
(2)
u(0, t) = ϕ(t), u(L, t) = ψ(t), 0 ≤ t ≤ T,
(3)
where 0 < γ < 1, and 0 Dt1−γ u(x, t) is the Riemann-Liouville fractional partial derivative of order 1 − γ defined by t
1−γ u(x, t) 0 Dt
1 ∂ Z = u(x, ν)(t − ν)γ−1 dν. Γ(γ) ∂t 0
In this article, we always assume that u(x, t) ∈ U (Ω) is the exact solution 2 ∈ C(Ω), where of the problem (1)-(3), and ∂ f∂t(x,t) 2 Ω = {(x, t)| 0 ≤ x ≤ L, 0 ≤ t ≤ T } , (
)
∂ 6 u(x, t) ∂ 4 u(x, t) , ∈ C(Ω) . U (Ω) = u(x, t) | ∂x6 ∂x2 ∂t2
2
The numerical simulation method
In this paper, we let that xj = j4x , j = 0, 1, . . . , J; tk = k4t , k = 0, 1, . . . , K, where 4x = L/J and 4t = T /K is spatial step and temporal step, respectively. Also define that δx2 skj = skj−1 − 2skj + skj+1 . Lemma 2.1. If
∂ 6 u(x,t) ∂x6
1+
Proof. See [7]. Lemma 2.2. If
∈ C(Ω), then
1 2 ∂ 2 u(xj , tm ) δx2 u(xj , tm ) + O(44x ), δx = 2 2 12 ∂x 4x
∂ 4 u(x,t) ∂x2 ∂t2
∈ C(Ω), then 3
∂ 2 u(xj , ν) tl − ν ∂ 2 u(xj , tl−1 ) ν − tl−1 ∂ 2 u(xj , tl ) = + + O(42t ), tl−1 ≤ ν ≤ tl , ∂x2 4t ∂x2 4t ∂x2 Proof. From the linear interpolation formula easy to get the conclusion. 2 2 ∈ C(Ω), then following trapezoidal formula holds Lemma 2.3. If ∂ f∂t(x,t) 2 Ztk
f (xj , ν)dν =
tk−1
4t [f (xj , tk−1 ) + f (xj , tk )] + O(43t ), 2
Proof. From the trapezoidal formula of definite integral easy to get the conclusion. 2
Integrating both sides of Eq.(1) from tk−1 to tk and take x = xj , immediately gets u(xj , tk ) − u(xj , tk−1 ) t k
1 Z = Γ(γ) 0
− =
tZk−1 0
1 Γ(γ) −
!
Ztk ∂ u(xj , ν) γ−1 − u(xj , ν) (tk−1 − ν) dν + f (xj , ν)dν ∂x2 2
k Ztl X
k−1 X
!
∂ 2 u(xj , ν) − u(xj , ν) (tk − ν)γ−1 dν ∂x2
∂ 2 u(xj , ν) − u(xj , ν) (tk − ν)γ−1 dν 2 ∂x
l=1 tl−1
Ztl
l=1 tl−1
tk−1
!
!
Ztk ∂ u(xj , ν) γ−1 f (xj , ν)dν,(4) − u(xj , ν) (tk−1 − ν) dν + ∂x2 2
tk−1
further apply Lemma 2.1, Lemma 2.2 and Lemma 2.3, we have
1+
1 2 1 δx u(xj , tk ) − 1 + δx2 u(xj , tk−1 ) 12 12
" k Ztl 1 2 1 X tl − ν = 1 + δx 12 Γ(γ) l=1 4t
tl−1
ν − tl−1 + 4t
2
!
#
∂ u(xj , tl ) − u(xj , tl ) + O(42t ) (tk − ν)γ−1 dν ∂x2
tl " X Z 1 2 1 k−1 tl − ν − 1 + δx 12 Γ(γ) l=1 4t
!
∂ 2 u(xj , tl−1 ) − u(xj , tl−1 ) ∂x2
tl−1
4
∂ 2 u(xj , tl−1 ) − u(xj , tl−1 ) ∂x2
!
!
#
ν − tl−1 ∂ 2 u(xj , tl ) + − u(xj , tl ) + O(42t ) (tk−1 − ν)γ−1 dν 4t ∂x2 1 4t 1 + δx2 [f (xj , tk−1 ) + f (xj , tk )] + O(43t ) + 2 12
" ! k Ztl tl − ν ∂ 2 u(xj , tl−1 ) 1 X 1 2 = 1 + δx − u(xj , tl−1 ) Γ(γ) l=1 4t 12 ∂x2 tl−1
ν − tl−1 1 + 1 + δx2 4t 12
#
!
∂ 2 u(xj , tl ) − u(xj , tl ) + O(42t ) (tk − ν)γ−1 dν ∂x2
! tl " X Z 1 k−1 tl − ν ∂ 2 u(xj , tl−1 ) 1 2 − 1 + δx − u(xj , tl−1 ) Γ(γ) l=1 4t 12 ∂x2 tl−1
#
!
∂ 2 u(xj , tl ) ν − tl−1 1 + 1 + δx2 − u(xj , tl ) + O(42t ) (tk−1 − ν)γ−1 dν 2 4t 12 ∂x 1 4t 1 + δx2 [f (xj , tk−1 ) + f (xj , tk )] + O(43t ) + 2 12
( " # k Ztl tl − ν δx2 u(xj , tl−1 ) 1 X 1 2 4 = + O(4x ) − 1 + δx u(xj , tl−1 ) Γ(γ) l=1 4t 42x 12 tl−1
"
#
#
)
ν − tl−1 δx2 u(xj , tl ) 1 + + O(44x ) − 1 + δx2 u(xj , tl ) + O(42t ) (tk − ν)γ−1 dν 2 4t 4x 12
" # tl ( X Z 1 2 1 k−1 tl − ν δx2 u(xj , tl−1 ) 4 + O(4x ) − 1 + δx u(xj , tl−1 ) − Γ(γ) l=1 4t 42x 12 tl−1
"
)
ν − tl−1 δx2 u(xj , tl ) 1 + + O(44x ) − 1 + δx2 u(xj , tl ) + O(42t ) (tk−1 − ν)γ−1 dν 2 τ 4x 12 4t 1 + 1 + δx2 [f (xj , tk−1 ) + f (xj , tk )] + O(43t ) 2 12 γ k k X X 1 4 4γt = 2 t λl δx2 u(xj , tk−l ) − λl 1 + δx2 u(xj , tk−l ) 4x Γ(2 + γ) l=0 Γ(2 + γ) l=0 12 1 4t 1 + δx2 [f (xj , tk−1 ) + f (xj , tk )] + Rjk , + 2 12 namely
Lx u(xj , tk ) − Lx u(xj , tk−1 ) = µ1
k X l=0
+
λl δx2 u(xj , tk−l )
− µ2
k X
λl Lx u(xj , tk−l )
l=0
4t Lx [f (xj , tk−1 ) + f (xj , tk )] + Rjk , 2 5
(5)
where
1 4γ Lx = 1 + δx2 , µ1 = 2 t > 0, 12 4x Γ(2 + γ)
21+γ − 3, λ = 1;
(1 + γ) [k γ
Rjk =
4γt > 0, Γ(2 + γ)
1, λ = 0;
λl =
µ2 =
(6) −(l − 2)1+γ + 3(l − 1)1+γ − 3l1+γ + (l + 1)1+γ , λ = 2, · · · , k − 1;
− (k − 1)γ ] − (k − 2)1+γ + 2(k − 1)1+γ − k 1+γ , λ = k,
k Ztl X 2 4 O(4t ) + O(4x )
l=1 tl−1
(tk − ν)γ−1 dν −
k−1 X
Ztl
l=1 tl−1
3 (tk−1 − ν)γ−1 dν +O(4t ).
Because k Ztl X
l=1 tl−1
k−1 X
Ztl
l=1 tl−1
γ−1
(tk − ν)
dν =
Ztk 0
(tk−1 − ν)γ−1 dν =
(tk − ν)γ−1 dν =
tZk−1 0
(tk−1 − ν)γ−1 dν =
(k4t )γ , γ ((k − 1)4t )γ , γ
and k4t ≤ T , therefore Rjk = O(42t + 44x )
(7)
According to the above analysis, we now present the following numerical simulation method for problem (1)-(3) Lx ukj − Lx uk−1 j = µ1
k X l=0
λl δx2 uk−l − µ2 j
k X
λl Lx uk−l + j
l=0
4t k−1 Lx fj + fjk , 2
(8)
k = 1, 2, . . . , K; j = 1, 2, . . . , J − 1, u0j = wj , j = 0, 1, . . . , J,
(9)
uk0 = ϕ(tk ), ukJ = ψ(tk ), k = 0, 1, . . . , K, 6
(10)
where fjl = f (xj , tl ), wj = w(xj ), whereas ulj denote the numerical approximation to the exact solution u(xj , tl ).
3
Stability of the numerical simulation method
Similar to [24], we have the following lemma: Lemma 3.1. If 0 < γ < 1, then the coefficients λl (l = 0, 1, . . . , k) possess the following properties: (1) λ0 = 1, (2)
k P
l=0
λ1 = 21+γ − 3,
l = 2, · · · , k;
λl < 0,
λl = (1 + γ) [k γ − (k − 1)γ ] > 0;
(3) −1 < −λ1 < 1, −
k P
l=2
λl < 1 + λ1 .
We now study stability of the numerical simulation method (8)-(10), considering the following difference equation
Lx ρkj − Lx ρk−1 = µ1 j
k X l=0
λl δx2 ρk−l − µ2 j
k X
λl Lx ρk−l j ,
(11)
l=0
k = 1, 2, . . . , K; j = 1, 2, . . . , J − 1,
where ρkj = ukj − uekj , whereas uekj is an approximation for ukj . Let that
h
ρk = ρk1 , ρk2 , . . . , ρkJ−1
iT
.
Similar to the derivation in [5], we have that k
kρ k2 ≡
J−1 X j=1
1
2
4x |ρkj |2
=
∞ X
1 2
l=−∞
2
|ξk (l)|
, k = 0, 1, . . . , K,
(12)
where ξk (l) is Fourier expansion coefficient of ρk (x), whereas ρk (x) is the following grid function
ρk (x) =
ρk , j
i
when x ∈ xj− 1 , xj+ 1 , j = 1, 2, . . . , J − 1;
0, when x ∈ [0,
2
4x S ] (L 2
7
2
−
4x , L]. 2
We now suppose that ρkj = ξk eiσj4x , where σ = 2πl/L. By substituting the above expression into (11) immediately !
!
σ4x σ4x 1 1 ξk − 1 − sin2 ξk−1 1 − sin2 3 2 3 2
!
k σ4x X σ4x 1 = −4µ1 sin λl ξk−l − µ2 1 − sin2 λl ξk−l , 2 l=0 3 2 2
(13)
k = 1, 2, . . . , K. For convenience, we rewrite Eq.(13) as follows "
!
#
k k σ4x X 1 σ4x X 1 Bξk−1 − 4µ1 sin2 λl ξk−l − µ2 1 − sin2 λl ξk−l , ξk = A 2 l=2 3 2 l=2 (14)
k = 1, 2, . . . , K, where
B
A = (1 + µ2 ) 1 − 13 sin2
= (1 − λ1 µ2 ) 1 − 13 sin2
σ4x 2
σ4x 2
+ 4µ1 sin2
σ4x , 2
− 4λ1 µ1 sin2
σ4x . 2
(15)
Lemma 3.2. For 0 < γ < 1, if the temporal and spatial steps satisfy that 42t = 44x , then A ≥ 1. Proof. By 0 < γ < 1 and 42t = 44x gets 4γt 42(γ−1) 1 1 x µ1 = = ≥ = , 2 Γ(2 + γ)4x Γ(2 + γ) Γ(3) 2 then
A= ≥ = ≥
!
1 σ4x σ4x (1 + µ2 ) 1 − sin2 + 4µ1 sin2 3 2 2 ! σ4x 1 σ4x 1 − sin2 + 4µ1 sin2 3 2 2 1 σ4x 4µ1 − sin2 +1 3 2 1. 8
This complete the proof of Lemma 3.2. 2 Lemma 3.3. For 0 < γ < 1, if λ1 < 0, or if 0 ≤ λ1 ≤
1 , 6µ1 +µ2
then
B ≥ 0. Proof. When λ1 < 0, the conclusion is clear. When 0 ≤ λ1 ≤
1 , 6µ1 +µ2
notice that
0 < µ2 =
4γt 1 < =1 Γ(2 + γ) Γ(2)
and −1 < −λ1 < 1,
then from this we obtain
1 − λ1 µ2 > 0,
1 σ4x B = (1 − λ1 µ2 ) 1 − sin2 3 2 2 ≥ (1 − λ1 µ2 ) − 4λ1 µ1 3 2 2 = − λ1 (6µ1 + µ2 ) 3 3 ≥ 0.
!
− 4λ1 µ1 sin2
σ4x 2
This complete the proof of Lemma 3.3. 2 Lemma 3.4. For 0 < γ < 1, if λ1 < 0, or if 0 ≤ λ1 ≤
1 , 6µ1 +µ2
then
B <1 A Proof. From Lemma 3.3 and |λ1 | < 1, we have (1 − λ1 µ2 )
!
1 σ4x σ4x B = |B| = 1 − sin2 − 4λ1 µ1 sin2 3 2 2 ! 1 2 σ4x 2 σ4x ≤ (1 − λ1 µ2 ) 1 − sin + 4 λ1 µ1 sin 3 2 2 ! σ4x 1 2 σ4x ≤ (1 + |λ1 | µ2 ) 1 − sin + 4 |λ1 | µ1 sin2 3 2 2 ! 1 σ4x σ4x < (1 + µ2 ) 1 − sin2 + 4µ1 sin2 3 2 2 = A. 9
This complete the proof of Lemma 3.4. 2 Theorem 3.5. Assume that the temporal and spatial steps satisfy that 42t = 44x , if λ1 < 0, or if 0 ≤ λ1 ≤ 6µ11+µ2 , then the numerical simulation method (8)-(10) is stable. Proof. Firstly, we verify the following conclusion by induction: |ξk | ≤ |ξ0 |, k = 1, 2, . . . , K.
(16)
where ξk (k = 1, 2, . . . , K) be the solution of (14). For k = 1, apply Lemma 3.4, from Eq. (14) gets that |ξ1 | =
1 B|ξ0 | ≤ |ξ0 |, A
Suppose that |ξn | ≤ |ξ0 |, n = 1, 2, . . . , k − 1,
according to Lemma 3.1, Lemma 3.2 and Lemma 3.3, from Eq. (14) we have that
!
k k σ4x X 1 σ4x X 1 λl ξk−l − µ2 1 − sin2 λl ξk−l |ξk | = Bξk−1 − 4µ1 sin2 A 2 l=2 3 2 l=2
"
!
#
k k 1 1 σ4x X σ4x X ≤ |λl ||ξk−l | + µ2 1 − sin2 |λl ||ξk−l | B|ξk−1 | + 4µ1 sin2 A 2 l=2 3 2 l=2
"
!
#
k k X 1 1 2 σ4x X 2 σ4x ≤ B + 4µ1 sin |λl | + µ2 1 − sin |λl | |ξ0 | A 2 l=2 3 2 l=2
"
!
#
k k 1 1 σ4x X σ4x X = B − 4µ1 sin2 λl − µ2 1 − sin2 λl |ξ0 | A 2 l=2 3 2 l=2
"
σ4x 1 σ4x 1 B + 4µ1 (1 + λ1 ) sin2 + µ2 (1 + λ1 ) 1 − sin2 ≤ A 2 3 2 " !# 1 σ4x 1 σ4x = 4µ1 sin2 + (1 + µ2 ) 1 − sin2 |ξ0 | A 2 3 2 = |ξ0 |.
!#
|ξ0 |
From this the conclusion (16) holds. Further applying (12) and (16), then the solution ρk of difference equation (11) satisfies that kρk k2 ≤ kρ0 k2 , k = 1, 2, . . . , K. This completes the proof of Theorem 3.5. 2 10
4
Solvability of the numerical simulation method
In this section, we discuss solvability of the numerical simulation method (8)-(10) . If define that Uk = [uk1 , uk2 , . . . , ukJ−1 ]T , then on the k layer, the numerical simulation method (8)-(10) can be written in the following matrix form: AUk = bk ,
(17)
where the right term bk is a (J − 1) dimension column vector witch includes the known initial and boundary values, the values of known term fjk and the before k layer known values unj (n = 0, 1, . . . , k − 1), whereas the coefficient matrix is a (J − 1) × (J − 1) known matrix:
e a a
A=
where
e a
a ae ... ... ... ae a
, e a
ae a
5 5 1 1 + 2µ1 + µ2 , ae = − µ1 − µ2 . 6 6 12 12 Theorem 4.1. The numerical simulation method (8)-(10) is uniquely solvable. Proof. It is easy to see that A is a strictly diagonally dominant matrix, then A is a nonsingular matrix, therefore the linear equation system (17) exists unique solution, i.e., the numerical simulation method (8)-(10) is uniquely solvable. 2 a=
5
Convergence of the numerical simulation method
In this section, we discuss convergence of the numerical simulation method (8)-(10). Subtracting (8) from (5), we obtain the following error equation Lx Ejk − Lx Ejk−1 11
= µ1
k X
λl δx2 Ejk−l
l=0
− µ2
k X
λl Lx Ejk−l + Rjk ,
(18)
l=0
k = 1, 2, . . . , K; j = 1, 2, . . . , J − 1, where Ejk = u(xj , tk ) − ukj . Let that h
k E k = E1k , E2k , . . . , EJ−1
iT
h
k , Rk = R1k , R2k , . . . , RJ−1
iT
.
Similar to the derivation in [5], we have that, respectively k
kE k2 ≡
J−1 X
j=1
1 2
4x |Ejk |2
=
∞ X
l=−∞
1 2
2
|αk (l)|
, k = 0, 1, . . . , K
(19)
, k = 0, 1, . . . , K,
(20)
and k
kR k2 ≡
J−1 X
j=1
1
2
4x |Rjk |2
=
∞ X
l=−∞
1 2
2
|βk (l)|
where αk (l) and βk (l) are Fourier expansion coefficients of E k (x) and Rk (x), whereas E k (x) and Rk (x) are the following grid functions, respectively
and
Ek,
when x ∈ xj− 1 , xj+ 1 , j = 1, 2, . . . , J − 1;
Rk ,
when x ∈ xj− 1 , xj+ 1 , j = 1, 2, . . . , J − 1;
E k (x) =
Rk (x) =
j
0, when x ∈ [0, j
0, when x ∈ [0,
2
4x S ] (L 2 2
4x S ] (L 2
2
i
− 2
i
−
4x , L], 2
4x , L]. 2
We now suppose that, respectively Ejk = αk eiσj4x ,
Rjk = βk eiσj4x ,
where σ = 2πl/L. By substituting the above expressions into (18), immediately !
!
σ4x 1 σ4x 1 1 − sin2 αk − 1 − sin2 αk−1 3 2 3 2
!
k k σ4x X 1 σ4x X = −4µ1 sin2 λl αk−l + βk , (21) λl αk−l − µ2 1 − sin2 2 l=0 3 2 l=0
12
k = 1, 2, . . . , K. For convenience, we rewrite equation (21) as follows "
k σ4x X 1 σ4x 1 Bαk−1 − 4µ1 sin2 λl αk−l − µ2 1 − sin2 αk = A 2 l=2 3 2 k X
#
λl αk−l + βk , k = 1, 2, . . . , K.
l=2
!
(22)
From (7), then there is a positive constant C1 , such that
|Rjk | ≤ C1 42t + 44x , k = 1, 2, . . . , K, again use the first equality of (20), result in √ kRk k2 ≤ C1 L 42t + 44x , k = 1, 2, . . . , K.
(23)
(24)
Based on the convergence of the series in RHS of (20), then there is a positive constant C2 , such that |βk | ≡ |βk (l)| ≤ C2 4t |β1 (l)| ≡ C2 4t |β1 |, k = 1, 2, . . . , K.
(25)
Theorem 5.1. Assume that u(x, t) ∈ U (Ω) is the exact solution of the 2 ∈ C(Ω), and the temporal and spatial steps satisfy problem (1)-(3), ∂ f∂t(x,t) 2 4 2 that 4t = 4x , if λ1 < 0, or if 0 ≤ λ1 ≤ 6µ11+µ2 , then the numerical simulation method (8)-(10) is convergent with the order O(42t + 44x ). Proof. At first, we verify the following conclusion by induction: |αk | ≤ C2 k4t |β1 |, k = 1, 2, . . . , K, where αk (k = 1, 2, . . . , K) be the solution of (22). For k = 1, from Eq. (22) gets 1 β1 , A applying Lemma 3.2 and (25), we obtain that α1 =
|α1 | =
1 |β1 | ≤ |β1 | ≤ C2 4t |β1 |. A 13
(26)
Suppose that |αn | ≤ C2 n4t |β1 |, n = 1, 2, . . . , k − 1, in view of Lemma 3.1, Lemma 3.2 and Lemma 3.3, from Eq. (22), we have that
!
k k X 1 2 σ4x X 1 2 σ4x λl αk−l − µ2 1 − sin λl αk−l + βk αk = Bαk−1 − 4µ1 sin A 2 l=2 3 2 l=2
"
!
#
k k X 1 1 2 σ4x X 2 σ4x ≤ B|αk−1 | + 4µ1 sin |λl | |αk−l | + µ2 1 − sin |λl | |αk−l | + |βk | A 2 l=2 3 2 l=2
"
!
k k σ4x X σ4x X 1 1 B(k − 1) + 4µ1 sin2 |λl | (k − l) + µ2 1 − sin2 |λl |( k − l) + 1 ≤ A 2 l=2 3 2 l=2 C2 4t |β1 |
(
"
!
#
k k 1 σ4x X 1 σ4x X ≤ B + 4µ1 sin2 |λl | + µ2 1 − sin2 |λl | (k − 1) + 1 A 2 l=2 3 2 l=2 C2 4t |β1 |
=
(
(
"
!
#
)
)
k k 1 1 σ4x X σ4x X B − 4µ1 sin2 λl − µ2 1 − sin2 λl (k − 1) + 1 C2 4t |β1 | A 2 l=2 3 2 l=2
"
!#
#
)
1 1 σ4x σ4x ≤ B + 4µ1 (1 + λ1 ) sin2 + µ2 (1 + λ1 ) 1 − sin2 (k − 1) + 1 C2 4t |β1 | A 2 3 2 ! # ) ( " 1 2 σ4x 1 2 σ4x (1 + µ2 ) 1 − sin + 4µ1 sin (k − 1) + 1 C2 4t |β1 | = A 3 2 2 ≤ C2 k4t |β1 |. From this the conclusion (26) holds. Further combining (19), (20) and (26), we have that kE k k2 ≤ C2 k4t kR1 k2 ≤ C(42t + 44x ), √ where C = C1 C2 LT. This completes the proof of Theorem 5.1. 2 Remark 1. In fact, because there are many numerical integration formulas with high order accuracy, to continue to structure numerical method with high temporal accuracy for solving the reaction-subdiffusion equation (1) is not difficult. However, stability and convergence analysis may be very difficult.
6
Numerical experiment
In order to certificate the theoretical analysis results, we employ the numerical simulation method (8)-(10) to solve the following initial-boundary value problem of the time fractional reaction-subdiffusion equation 14
∂u(x, t) =0 Dt1−γ ∂t
!
∂ 2 u(x, t) − u(x, t) + 2ex t, ∂x2
(27)
0 < t ≤ 1, 0 < x < 1, u(x, 0) = 0, 0 ≤ x ≤ 1.
(28)
u(0, t) = t2 , u(1, t) = et2 , 0 ≤ t ≤ 1,
(29)
The exact solution of the problem (27)-(29) is u(x, t) = ex t2 . Let that the maximum error of the numerical solutions E∞ = max
0≤k≤K
n
o
kE k k2 = O(4ζt + 4ηx ).
Suppose that in the front experiment, the maximum error, temporal step, spa(f ) (f ) ) tial step of the numerical solutions are E∞ , 4t , 4(f x , respectively, whereas in the subsequent experiment, the maximum error, temporal step, spatial step (s) (f ) , 4t , 4(s) of the numerical solutions are E∞ x , respectively, then the temporal accuracy ζ and spatial accuracy η calculated by the following formulas (s)
(s)
ζ=
ln E∞ (f ) E∞
ln
(s)
4t
,
η=
(f )
4t
ln E∞ (f ) E∞
(s)
x ln 4(f ) 4x
.
Table 1 exhibits the maximum error E∞ , temporal accuracy ζ and spatial accuracy η of the numerical simulation method (8)-(10) for solving for the problem (27)-(29) versus various γ and 42t = 44x , respectively. From Table 1, it can be seen that our theoretical analysis has been strongly supported by the numerical experiments. Fig. 1 reveals the comparison of the numerical solution using the numerical simulation method (8)-(10) and the exact solution for the problem (27)-(29) 1 . at t = 14 , 24 , 34 , 1, respectively, for γ = 0.5 and 42t = 44x = 256 Fig. 2 shows the comparison of the numerical solution using the numerical simulation method (8)-(10) and the exact solution for the problem (27)-(29) 1 at x = 14 , 42 , 34 , respectively, for γ = 0.5 and 42t = 44x = 256 . Fig. 3 illustrates the absolute error E(x, t) of the numerical solution for the problem (27)-(29) using the numerical simulation method (8)-(10) for γ = 0.5 15
and 42t = 44x =
1 , 256
where E(x, t) = u(x, t) − ue(x, t),
whereas ue(x, t) is the numerical approximation for u(x, t).
Fig. 1, 2 and 3 show that the numerical solution are very consistent with the exact solution. Table 1 The maximum error E∞ , temporal accuracy ζ and spatial accuracy η of the numerical simulation method (8)-(10) for solving for the problem (27)-(29)
42t = 44x =
γ 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
1 16
3.1544 × 10−5 3.1292 × 10−5 3.1039 × 10−5 3.0701 × 10−5 3.0364 × 10−5 3.0111 × 10−5 2.9774 × 10−5 2.9521 × 10−5 2.9184 × 10−5
42t = 44x =
1 81
6.5689 × 10−6 6.5689 × 10−6 6.4718 × 10−6 6.4718 × 10−6 6.3230 × 10−6 6.2776 × 10−6 6.2258 × 10−6 6.1287 × 10−6 6.0834 × 10−6
ζ
η
1.93 3.87 1.93 3.85 1.93 3.87 1.92 3.84 1.93 3.87 1.93 3.87 1.93 3.86 1.94 3.88 1.93 3.87
42t = 44x =
1 256
2.2508 × 10−6 2.2231 × 10−6 2.2231 × 10−6 2.1832 × 10−6 2.1832 × 10−6 2.0464 × 10−6 2.0464 × 10−6 2.0066 × 10−6 2.0066 × 10−6
ζ
1.86 3.72 1.88 3.77 1.86 3.71 1.89 3.78 1.85 3.70 1.95 3.90 1.93 3.87 1.94 3.88 1.93 3.86
Remark 2. In section 3 and 5, in order to discuss stability and convergence of the numerical simulation method, we suppose that holds the condition: λ1 < 0, or 0 ≤ λ1 ≤ 6µ11+µ2 , this condition is effectively equivalent to that some γ are restricted. When this condition is not holds, we are not able to confirm the stability and convergence of the numerical simulation method. However, the numerical experiment in this section shows that for the wide range of γ, there are no symptoms of numerical instability.
7
Conclusion
In this paper, the numerical simulation method with second order temporal accuracy and fourth order spatial accuracy has been proposed to solve the time fractional reaction-subdiffusion equation; by Fourier analysis and 16
η
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
0.6
0.8
1
x
Fig. 1. The comparison of the numerical solution (n.s.) using the numerical simulation method (8)-(10) and the exact solution (e.s.) for the problem (27)-(29) at 1 t = 14 , 24 , 34 , 1, respectively, for γ = 0.5 and 42t = 44x = 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
0.6
0.8
1
t
Fig. 2. The comparison of the numerical solution (n.s.) using the numerical simulation method (8)-(10) and the exact solution (e.s.) for the problem (27)-(29) at 1 x = 14 , 42 , 34 , respectively, for γ = 0.5 and 42t = 44x = 256 .
algebraic theory, the stability, convergence and solvability of the numerical simulation method have been discussed, respectively; the theoretical analysis results very consistent with the numerical experiment. The numerical simu17
−6
x 10
Absolute error |E(x,t)|
3 2.5 2 1.5 1 0.5 0 1 1 0.8
0.5
0.6 0.4 t
0
0.2 0
x
Fig. 3. The absolute error of the numerical solution for the problem (27)-(29) using 1 the numerical simulation method (8)-(10) for γ = 0.5 and 42t = 44x = 256 .
lation method and numerical analysis method in this paper can be extended to other fractional partial deferential equations with the Riemann-Liouville fractional partial derivative. Acknowledgments The authors would like to thank the reviewers for their careful reading of this article and their comments and suggestions which have improved the article.
References [1] E. Abad, S. B. Yuste, K. Lindenberg, Reaction-subdiffusion and reactionsuperdiffusion equations for evanescent particles performing continuous-time random walks, Phys. Rev. E 81 (2010) 031115. [2] M. Abbaszadeh, M. Dehghan, A meshless numerical procedure for solving fractional reaction subdiffusion model via a new combination of alternating direction implicit (ADI) approach and interpolating element free Galerkin (EFG) method, Comput. Math. Appl. 70 (2015) 2493-2512. [3] B. Baeumer, M. Kovcs, M. Meerschaert, Numerical solutions for fractional reaction-diffusion equations, Comput. Math. Appl. 55 (2008) 2212–2226. [4] V. K. Baranwal, R. K. Pandey, M. P. Tripathi, O. P. Singh, An analytic algorithm for time fractional nonlinear reaction-diffusion equation based on a new iterative method, Commun. Nonlinear. Sci. Numer. Simulat. 17 (2012) 3906-3921.
18
[5] C.-M. Chen, F. Liu, I. Turner, and V. Anh , Fourier method for the fractional diffusion equation describing sub-diffusion, J. Comp. Phys. 227 (2007) 886-897. [6] C.-M. Chen, F. Liu, K. Burrage, Finite difference methods and a Fourier analysis for the fractional reaction-subdiffusion equation, Appl. Math. Comp. 198 (2008) 754-769. [7] C.-M. Chen, F. Liu, V. Anh, I. Turner, Numerical schemes with high spatial accuracy for a variable-order anomalous subdiffusion equation, SIAM J. Sci. Comput. 32 (4) (2010), 1740-1760. [8] C.-M. Chen, F. Liu, I. Turner, V. Anh, Y. Chen, Numerical approximation for a variable-order nonlinear reaction-subdiffusion equation, Numer. Algor. 63 (2013) 265-290. [9] J. W. Chiu, K. -H. Chiam, Monte Carlo simulation and linear stability analysis of Turing pattern formation in reaction-subdiffusion systems. Phys. Rev. E 78 (2008) 056708. [10] S. Das, P. K. Gupta, P. Ghosh, An approximate solution of nonlinear fractional reaction-diffusion equation, Appl. Math. Model. 35 (2011) 4071-4076. [11] E. H. Doha, A. H. Bhrawy, S. S. Ezz-Eldien, An efficient Legendre spectral tau matrix formulation for solving fractional subdiffusion and reaction subdiffusion equations, J. Comput. Nonlin. Dyn. 10 (2015) DOI: 10.1115/1.4027944. [12] M. Dehghan, M. Abbaszadeh, A. Mohebbi, Error estimate for the numerical solution of fractional reaction-subdiffusion process based on a meshless method, J. Comput. Appl. Math. 280 (2015) 14-36. [13] H. F. Ding, C. P. Li, Mixed spline function method for reaction-subdiffusion equations, J. Comp. Phys. 242 (2013) 103-123. [14] V. Gafiychuk, B. Y. Datsko, Stability analysis and oscillatory structures in time-fractional reaction-diffusion systems, Phys. Rev. E 75 (2007) 055201. [15] V. Gafiychuk, B. Datsko, V. Meleshko, Mathematical modeling of time fractional reaction-diffusion systems, J. Comput. Appl. Math. 220 (2008) 215225. [16] V. Gafiychuk, B. Datsko, Inhomogeneous oscillatory solutions in fractional reaction-diffusion systems and their computer modeling, Appl. Math. Comp. 198 (2008) 251-260. [17] V. Gafiychuk, B. Datsko, V. Meleshko, D. Blackmore, Analysis of the solutions of coupled nonlinear fractional reaction-diffusion equations, Chaos Soliton. Fract. 41 (2009) 1095-1104. [18] V. Gafiychuk, B. Datsko, Mathematical modeling of different types of instabilities in time fractional reaction-diffusion systems, Comput. Math. Appl. 59 (2010) 1101-1107.
19
[19] C. Y. Gong, W. M. Bao, G. J. Tang, A parallel algorithm for the Riesz fractional reaction-diffusion equation with explicit finite difference method, Fract. Calc. Appl. Anal. 16 (2013) 654-669. [20] H. J. Haubold, A. M. Mathai, R. K. Saxena, Solutions of fractional reactiondiffusion equations in terms of the H-function, Bull. Astron. Soc. India 35 (2007) 681-689. [21] H. J. Haubold, A. M. Mathai, R. K. Saxena, Further solutions of fractional reaction-diffusion equations in terms of the H-function, J. Comput. Appl. Math. 235 (2011) 1311-1316. [22] J. M. Haugh, Analysis of reaction-diffusion systems with anomalous subdiffusion, Biophys. J. 97 (2009) 435-442. [23] B. I. Henry, S. L. Wearne, Fractional reaction-diffusion, Phys. A 276 (2000) 448-455. [24] H. Huang, X. Cao, Numerical method for two dimensional fractional reaction subdiffusion equation, Eur. Phys. J.-Spec. Top. 222 (2013) 1961-1973. [25] T. A. M. Langlands, B. I. Henry, S. L. Wearne, Anomalous subdiffusion with multispecies linear reaction dynamics, Phys. Rev. E 77(2008) 021111. [26] M. S. Mommer, D. Lebiedz, Modeling subdiffusion using reaction diffusion systems, SIAM J. Appl. Math. 70 (2009) 112-132. [27] S. Z. Rida, A. M. A. El-Sayed, A. A. M. Arafa, On the solutions of time-fractional reaction-diffusion equations, Commun. Nonlinear. Sci. Numer. Simulat. 15 (2010) 3847-3854. [28] F. Sagus, V. P. Shkilev, I. M. Sokolov, Reaction-subdiffusion equations for the A B reaction. Phys. Rev. E 77 (2008) 032102. [29] R. K. Saxena, A. M. Mathai, H. J. Haubold, Fractional reaction-diffusion equations, Astrophys. Space Sci. 305(2006) 289-296. [30] R. K. Saxena, A. M. Mathai, H. J. Haubold, Solution of generalized fractional reaction-diffusion equations, Astrophys. Space Sci. 305 (2006) 305-313. [31] H. H. Schmidt-Martens, D. Froemberg, I. M. Sokolov, Front propagation in a one-dimensional autocatalytic reaction-subdiffusion system. Phys. Rev. E 79 (2009) 041135. [32] K. Seki, M. Wojcik, M. Tachiya, Fractional reaction-diffusion equation, J. Chem. Phys. 119 (2003) 2165-2170. [33] I. M. Sokolov, M. G. W. Schmidt F. Sagus, Reaction-subdiffusion equations. Phys. Rev. E 73 (2006) 031102. [34] N. H. Sweilam, M. M. Khader, M. Adel, Weighted average finite difference methods for fractional reaction-subdiffusion equation, Walailak J. Sci. Tech. 11 (2014) 361-377.
20
[35] S. Q. Wang, J. H. He, Variational iteration method for a nonlinear reactiondiffusion process, J. Chem. Reactor Eng. 6 (2008) A37. [36] A. Yadav, Werner Horsthemke: Kinetic equations for reaction-subdiffusion systems: Derivation and stability analysis. Phys. Rev. E 74 (2006) 066118. [37] B. Yu, X. Y. Jiang, H. Y. Xu, A novel compact numerical method for solving the two-dimensional non-linear fractional reaction-subdiffusion equation, Numer. Algor. 68 (2015) 923-950. [38] S. B. Yuste, K. Lindenberg, Subdiffusion-limited A+A reactions, Phys. Rev. Lett. 87 (11) (2001) art. no.-118301. [39] S. B. Yuste, L. Acedo, Katja Lindenberg: Reaction front in an A+B →C reaction-subdiffusion process. Phys. Rev. E 69 (2004) 036126. [40] S. B. Yuste, E. Abad, K. Lindenberg, Reaction-subdiffusion model of morphogen gradient formation, Phys. Rev. E 82 (2010) 061123. [41] P. Zhuang, F. Liu, V. Anh, I. Turner, Stability and convergence of an implicit numerical method for the non-linear fractional reaction-subdiffusion process, IMA J. Appl. Math. 74 (2009) 645-667. [42] A. Zoia, Continuous-time random-walk approach to normal and anomalous reaction-diffusion processes, Phys. Rev. E 77 (2008) 041115.
21
S0378-4754(17)30084-8 http://dx.doi.org/10.1016/j.matcom.2017.03.008 MATCOM 4445
To appear in:
Mathematics and Computers in Simulation
Received date: 2 January 2016 Revised date: 28 December 2016 Accepted date: 12 March 2017 Please cite this article as: Y. Chen, C.-M. Chen, Numerical simulation with high order accuracy for the time fractional reaction-subdiffusion equation, Math. Comput. Simulation (2017), http://dx.doi.org/10.1016/j.matcom.2017.03.008 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Numerical simulation with high order accuracy for the time fractional reaction-subdiffusion equation Y. Chen a , Chang-Ming.Chen a,∗ a School
of Mathematical Sciences, Xiamen University, Xiamen 361006, China
Abstract In this paper, we propose the numerical simulation method with second order temporal accuracy and fourth order spatial accuracy for the time fractional reactionsubdiffusion equation; the stability, convergence and solvability of the numerical simulation method respectively are discussed by Fourier analysis and algebraic theory; the theoretical analysis results very consistent with the numerical experiment. Key words: The time fractional reaction-subdiffusion equation, The numerical simulation method with second order temporal accuracy and fourth order spatial accuracy, Stability, Solvability, Convergence.
1
Introduction
Fractional reaction-subdiffusion equations are noticed and researched in recent years in many areas of science and engineering [38,32,39,33,36,29,30,20,14,42,28,25,9,22,26,31,15–18,40,23,21,4,1,19]. At the same time, some authors have make out a lot of research on numerical simulation and numerical analysis of fractional reaction-subdiffusion equations, and this research continues. For one-dimensional fractional reaction-subdiffusion equation, based on operator splitting, Baeumer et al. [3] proposed a practical method for numerical solution; Chen et al. [6] developed the finite difference methods, by Fourier analysis, they also verified that these finite difference methods have first order ∗ Corresponding author. Email address: [email protected] (Chang-Ming.Chen).
Preprint submitted to Elsevier Science
4 April 2017
temporal accuracy, second order spatial accuracy and unconditionally stable; Wang et al. [35] presented a variational iteration method; Zhuang et al. [41] developed an implicit numerical method, they also demonstrated that this implicit numerical method has first order temporal accuracy, second order spatial accuracy and unconditionally stable; Rida et al. [27] studied a new application of generalized differential transform method; Das et al. [10] presented an approximate solution; Ding et al. [13] proposed the two classes of difference schemes, they also verified that one scheme has first order temporal accuracy, second order spatial accuracy and unconditionally stable, whereas another scheme has first order temporal accuracy, fourth order spatial accuracy and unconditionally stable; Sweilam et al. [34] presented the weighted average finite difference methods, but they only given stable analysis of the methods; Doha et al. [11] constructed the efficient Legendre spectral tau matrix formulation. For multi-dimensional fractional reaction-subdiffusion equation, Huang et al. [24] proposed a numerical method, they also verified that this method has (1 + γ) order temporal accuracy, second order spatial accuracy and unconditionally stable; Dehghan et al. [12] proposed a numerical technique based on a meshless method, this technique has γ order temporal accuracy and unconditionally stable; based on a new combination of alternating direction implicit (ADI) approach and interpolating element free Galerkin (EFG) method, Abbaszadeh et al. [2] researched a meshless numerical procedure, the procedure has (1 + γ) order temporal accuracy and unconditionally stable; Yu et al. [37] developed the novel compact numerical method, the method has second order temporal accuracy and fourth order spatial accuracy. For one-dimensional variable order reaction-subdiffusion equation, Chen et al. [8] discussed the numerical approximation, they also verified that this numerical approximation method has first order temporal accuracy, second order spatial accuracy and unconditionally stable, and they also developed a numerical approximation method for improving temporal accuracy. According to the literature, so far, although there are some numerical solving techniques for improving temporal accuracy of fractional reactionsubdiffusion equation, but they are only supported by the numerical experiments whereas lack of support for strict numerical analysis. Therefore, how to structure numerical simulation method with high temporal accuracy for solving fractional reaction-subdiffusion equation is very meaningful. In this paper, we research the following time fractional reaction-subdiffusion equation: ∂u(x, t) =0 Dt1−γ ∂t
!
∂ 2 u(x, t) − u(x, t) + f (x, t), 0 < t ≤ T, 0 < x < L, (1) ∂x2 2
with the following initial and boundary conditions: u(x, 0) = φ(x), 0 ≤ x ≤ L,
(2)
u(0, t) = ϕ(t), u(L, t) = ψ(t), 0 ≤ t ≤ T,
(3)
where 0 < γ < 1, and 0 Dt1−γ u(x, t) is the Riemann-Liouville fractional partial derivative of order 1 − γ defined by t
1−γ u(x, t) 0 Dt
1 ∂ Z = u(x, ν)(t − ν)γ−1 dν. Γ(γ) ∂t 0
In this article, we always assume that u(x, t) ∈ U (Ω) is the exact solution 2 ∈ C(Ω), where of the problem (1)-(3), and ∂ f∂t(x,t) 2 Ω = {(x, t)| 0 ≤ x ≤ L, 0 ≤ t ≤ T } , (
)
∂ 6 u(x, t) ∂ 4 u(x, t) , ∈ C(Ω) . U (Ω) = u(x, t) | ∂x6 ∂x2 ∂t2
2
The numerical simulation method
In this paper, we let that xj = j4x , j = 0, 1, . . . , J; tk = k4t , k = 0, 1, . . . , K, where 4x = L/J and 4t = T /K is spatial step and temporal step, respectively. Also define that δx2 skj = skj−1 − 2skj + skj+1 . Lemma 2.1. If
∂ 6 u(x,t) ∂x6
1+
Proof. See [7]. Lemma 2.2. If
∈ C(Ω), then
1 2 ∂ 2 u(xj , tm ) δx2 u(xj , tm ) + O(44x ), δx = 2 2 12 ∂x 4x
∂ 4 u(x,t) ∂x2 ∂t2
∈ C(Ω), then 3
∂ 2 u(xj , ν) tl − ν ∂ 2 u(xj , tl−1 ) ν − tl−1 ∂ 2 u(xj , tl ) = + + O(42t ), tl−1 ≤ ν ≤ tl , ∂x2 4t ∂x2 4t ∂x2 Proof. From the linear interpolation formula easy to get the conclusion. 2 2 ∈ C(Ω), then following trapezoidal formula holds Lemma 2.3. If ∂ f∂t(x,t) 2 Ztk
f (xj , ν)dν =
tk−1
4t [f (xj , tk−1 ) + f (xj , tk )] + O(43t ), 2
Proof. From the trapezoidal formula of definite integral easy to get the conclusion. 2
Integrating both sides of Eq.(1) from tk−1 to tk and take x = xj , immediately gets u(xj , tk ) − u(xj , tk−1 ) t k
1 Z = Γ(γ) 0
− =
tZk−1 0
1 Γ(γ) −
!
Ztk ∂ u(xj , ν) γ−1 − u(xj , ν) (tk−1 − ν) dν + f (xj , ν)dν ∂x2 2
k Ztl X
k−1 X
!
∂ 2 u(xj , ν) − u(xj , ν) (tk − ν)γ−1 dν ∂x2
∂ 2 u(xj , ν) − u(xj , ν) (tk − ν)γ−1 dν 2 ∂x
l=1 tl−1
Ztl
l=1 tl−1
tk−1
!
!
Ztk ∂ u(xj , ν) γ−1 f (xj , ν)dν,(4) − u(xj , ν) (tk−1 − ν) dν + ∂x2 2
tk−1
further apply Lemma 2.1, Lemma 2.2 and Lemma 2.3, we have
1+
1 2 1 δx u(xj , tk ) − 1 + δx2 u(xj , tk−1 ) 12 12
" k Ztl 1 2 1 X tl − ν = 1 + δx 12 Γ(γ) l=1 4t
tl−1
ν − tl−1 + 4t
2
!
#
∂ u(xj , tl ) − u(xj , tl ) + O(42t ) (tk − ν)γ−1 dν ∂x2
tl " X Z 1 2 1 k−1 tl − ν − 1 + δx 12 Γ(γ) l=1 4t
!
∂ 2 u(xj , tl−1 ) − u(xj , tl−1 ) ∂x2
tl−1
4
∂ 2 u(xj , tl−1 ) − u(xj , tl−1 ) ∂x2
!
!
#
ν − tl−1 ∂ 2 u(xj , tl ) + − u(xj , tl ) + O(42t ) (tk−1 − ν)γ−1 dν 4t ∂x2 1 4t 1 + δx2 [f (xj , tk−1 ) + f (xj , tk )] + O(43t ) + 2 12
" ! k Ztl tl − ν ∂ 2 u(xj , tl−1 ) 1 X 1 2 = 1 + δx − u(xj , tl−1 ) Γ(γ) l=1 4t 12 ∂x2 tl−1
ν − tl−1 1 + 1 + δx2 4t 12
#
!
∂ 2 u(xj , tl ) − u(xj , tl ) + O(42t ) (tk − ν)γ−1 dν ∂x2
! tl " X Z 1 k−1 tl − ν ∂ 2 u(xj , tl−1 ) 1 2 − 1 + δx − u(xj , tl−1 ) Γ(γ) l=1 4t 12 ∂x2 tl−1
#
!
∂ 2 u(xj , tl ) ν − tl−1 1 + 1 + δx2 − u(xj , tl ) + O(42t ) (tk−1 − ν)γ−1 dν 2 4t 12 ∂x 1 4t 1 + δx2 [f (xj , tk−1 ) + f (xj , tk )] + O(43t ) + 2 12
( " # k Ztl tl − ν δx2 u(xj , tl−1 ) 1 X 1 2 4 = + O(4x ) − 1 + δx u(xj , tl−1 ) Γ(γ) l=1 4t 42x 12 tl−1
"
#
#
)
ν − tl−1 δx2 u(xj , tl ) 1 + + O(44x ) − 1 + δx2 u(xj , tl ) + O(42t ) (tk − ν)γ−1 dν 2 4t 4x 12
" # tl ( X Z 1 2 1 k−1 tl − ν δx2 u(xj , tl−1 ) 4 + O(4x ) − 1 + δx u(xj , tl−1 ) − Γ(γ) l=1 4t 42x 12 tl−1
"
)
ν − tl−1 δx2 u(xj , tl ) 1 + + O(44x ) − 1 + δx2 u(xj , tl ) + O(42t ) (tk−1 − ν)γ−1 dν 2 τ 4x 12 4t 1 + 1 + δx2 [f (xj , tk−1 ) + f (xj , tk )] + O(43t ) 2 12 γ k k X X 1 4 4γt = 2 t λl δx2 u(xj , tk−l ) − λl 1 + δx2 u(xj , tk−l ) 4x Γ(2 + γ) l=0 Γ(2 + γ) l=0 12 1 4t 1 + δx2 [f (xj , tk−1 ) + f (xj , tk )] + Rjk , + 2 12 namely
Lx u(xj , tk ) − Lx u(xj , tk−1 ) = µ1
k X l=0
+
λl δx2 u(xj , tk−l )
− µ2
k X
λl Lx u(xj , tk−l )
l=0
4t Lx [f (xj , tk−1 ) + f (xj , tk )] + Rjk , 2 5
(5)
where
1 4γ Lx = 1 + δx2 , µ1 = 2 t > 0, 12 4x Γ(2 + γ)
21+γ − 3, λ = 1;
(1 + γ) [k γ
Rjk =
4γt > 0, Γ(2 + γ)
1, λ = 0;
λl =
µ2 =
(6) −(l − 2)1+γ + 3(l − 1)1+γ − 3l1+γ + (l + 1)1+γ , λ = 2, · · · , k − 1;
− (k − 1)γ ] − (k − 2)1+γ + 2(k − 1)1+γ − k 1+γ , λ = k,
k Ztl X 2 4 O(4t ) + O(4x )
l=1 tl−1
(tk − ν)γ−1 dν −
k−1 X
Ztl
l=1 tl−1
3 (tk−1 − ν)γ−1 dν +O(4t ).
Because k Ztl X
l=1 tl−1
k−1 X
Ztl
l=1 tl−1
γ−1
(tk − ν)
dν =
Ztk 0
(tk−1 − ν)γ−1 dν =
(tk − ν)γ−1 dν =
tZk−1 0
(tk−1 − ν)γ−1 dν =
(k4t )γ , γ ((k − 1)4t )γ , γ
and k4t ≤ T , therefore Rjk = O(42t + 44x )
(7)
According to the above analysis, we now present the following numerical simulation method for problem (1)-(3) Lx ukj − Lx uk−1 j = µ1
k X l=0
λl δx2 uk−l − µ2 j
k X
λl Lx uk−l + j
l=0
4t k−1 Lx fj + fjk , 2
(8)
k = 1, 2, . . . , K; j = 1, 2, . . . , J − 1, u0j = wj , j = 0, 1, . . . , J,
(9)
uk0 = ϕ(tk ), ukJ = ψ(tk ), k = 0, 1, . . . , K, 6
(10)
where fjl = f (xj , tl ), wj = w(xj ), whereas ulj denote the numerical approximation to the exact solution u(xj , tl ).
3
Stability of the numerical simulation method
Similar to [24], we have the following lemma: Lemma 3.1. If 0 < γ < 1, then the coefficients λl (l = 0, 1, . . . , k) possess the following properties: (1) λ0 = 1, (2)
k P
l=0
λ1 = 21+γ − 3,
l = 2, · · · , k;
λl < 0,
λl = (1 + γ) [k γ − (k − 1)γ ] > 0;
(3) −1 < −λ1 < 1, −
k P
l=2
λl < 1 + λ1 .
We now study stability of the numerical simulation method (8)-(10), considering the following difference equation
Lx ρkj − Lx ρk−1 = µ1 j
k X l=0
λl δx2 ρk−l − µ2 j
k X
λl Lx ρk−l j ,
(11)
l=0
k = 1, 2, . . . , K; j = 1, 2, . . . , J − 1,
where ρkj = ukj − uekj , whereas uekj is an approximation for ukj . Let that
h
ρk = ρk1 , ρk2 , . . . , ρkJ−1
iT
.
Similar to the derivation in [5], we have that k
kρ k2 ≡
J−1 X j=1
1
2
4x |ρkj |2
=
∞ X
1 2
l=−∞
2
|ξk (l)|
, k = 0, 1, . . . , K,
(12)
where ξk (l) is Fourier expansion coefficient of ρk (x), whereas ρk (x) is the following grid function
ρk (x) =
ρk , j
i
when x ∈ xj− 1 , xj+ 1 , j = 1, 2, . . . , J − 1;
0, when x ∈ [0,
2
4x S ] (L 2
7
2
−
4x , L]. 2
We now suppose that ρkj = ξk eiσj4x , where σ = 2πl/L. By substituting the above expression into (11) immediately !
!
σ4x σ4x 1 1 ξk − 1 − sin2 ξk−1 1 − sin2 3 2 3 2
!
k σ4x X σ4x 1 = −4µ1 sin λl ξk−l − µ2 1 − sin2 λl ξk−l , 2 l=0 3 2 2
(13)
k = 1, 2, . . . , K. For convenience, we rewrite Eq.(13) as follows "
!
#
k k σ4x X 1 σ4x X 1 Bξk−1 − 4µ1 sin2 λl ξk−l − µ2 1 − sin2 λl ξk−l , ξk = A 2 l=2 3 2 l=2 (14)
k = 1, 2, . . . , K, where
B
A = (1 + µ2 ) 1 − 13 sin2
= (1 − λ1 µ2 ) 1 − 13 sin2
σ4x 2
σ4x 2
+ 4µ1 sin2
σ4x , 2
− 4λ1 µ1 sin2
σ4x . 2
(15)
Lemma 3.2. For 0 < γ < 1, if the temporal and spatial steps satisfy that 42t = 44x , then A ≥ 1. Proof. By 0 < γ < 1 and 42t = 44x gets 4γt 42(γ−1) 1 1 x µ1 = = ≥ = , 2 Γ(2 + γ)4x Γ(2 + γ) Γ(3) 2 then
A= ≥ = ≥
!
1 σ4x σ4x (1 + µ2 ) 1 − sin2 + 4µ1 sin2 3 2 2 ! σ4x 1 σ4x 1 − sin2 + 4µ1 sin2 3 2 2 1 σ4x 4µ1 − sin2 +1 3 2 1. 8
This complete the proof of Lemma 3.2. 2 Lemma 3.3. For 0 < γ < 1, if λ1 < 0, or if 0 ≤ λ1 ≤
1 , 6µ1 +µ2
then
B ≥ 0. Proof. When λ1 < 0, the conclusion is clear. When 0 ≤ λ1 ≤
1 , 6µ1 +µ2
notice that
0 < µ2 =
4γt 1 < =1 Γ(2 + γ) Γ(2)
and −1 < −λ1 < 1,
then from this we obtain
1 − λ1 µ2 > 0,
1 σ4x B = (1 − λ1 µ2 ) 1 − sin2 3 2 2 ≥ (1 − λ1 µ2 ) − 4λ1 µ1 3 2 2 = − λ1 (6µ1 + µ2 ) 3 3 ≥ 0.
!
− 4λ1 µ1 sin2
σ4x 2
This complete the proof of Lemma 3.3. 2 Lemma 3.4. For 0 < γ < 1, if λ1 < 0, or if 0 ≤ λ1 ≤
1 , 6µ1 +µ2
then
B <1 A Proof. From Lemma 3.3 and |λ1 | < 1, we have (1 − λ1 µ2 )
!
1 σ4x σ4x B = |B| = 1 − sin2 − 4λ1 µ1 sin2 3 2 2 ! 1 2 σ4x 2 σ4x ≤ (1 − λ1 µ2 ) 1 − sin + 4 λ1 µ1 sin 3 2 2 ! σ4x 1 2 σ4x ≤ (1 + |λ1 | µ2 ) 1 − sin + 4 |λ1 | µ1 sin2 3 2 2 ! 1 σ4x σ4x < (1 + µ2 ) 1 − sin2 + 4µ1 sin2 3 2 2 = A. 9
This complete the proof of Lemma 3.4. 2 Theorem 3.5. Assume that the temporal and spatial steps satisfy that 42t = 44x , if λ1 < 0, or if 0 ≤ λ1 ≤ 6µ11+µ2 , then the numerical simulation method (8)-(10) is stable. Proof. Firstly, we verify the following conclusion by induction: |ξk | ≤ |ξ0 |, k = 1, 2, . . . , K.
(16)
where ξk (k = 1, 2, . . . , K) be the solution of (14). For k = 1, apply Lemma 3.4, from Eq. (14) gets that |ξ1 | =
1 B|ξ0 | ≤ |ξ0 |, A
Suppose that |ξn | ≤ |ξ0 |, n = 1, 2, . . . , k − 1,
according to Lemma 3.1, Lemma 3.2 and Lemma 3.3, from Eq. (14) we have that
!
k k σ4x X 1 σ4x X 1 λl ξk−l − µ2 1 − sin2 λl ξk−l |ξk | = Bξk−1 − 4µ1 sin2 A 2 l=2 3 2 l=2
"
!
#
k k 1 1 σ4x X σ4x X ≤ |λl ||ξk−l | + µ2 1 − sin2 |λl ||ξk−l | B|ξk−1 | + 4µ1 sin2 A 2 l=2 3 2 l=2
"
!
#
k k X 1 1 2 σ4x X 2 σ4x ≤ B + 4µ1 sin |λl | + µ2 1 − sin |λl | |ξ0 | A 2 l=2 3 2 l=2
"
!
#
k k 1 1 σ4x X σ4x X = B − 4µ1 sin2 λl − µ2 1 − sin2 λl |ξ0 | A 2 l=2 3 2 l=2
"
σ4x 1 σ4x 1 B + 4µ1 (1 + λ1 ) sin2 + µ2 (1 + λ1 ) 1 − sin2 ≤ A 2 3 2 " !# 1 σ4x 1 σ4x = 4µ1 sin2 + (1 + µ2 ) 1 − sin2 |ξ0 | A 2 3 2 = |ξ0 |.
!#
|ξ0 |
From this the conclusion (16) holds. Further applying (12) and (16), then the solution ρk of difference equation (11) satisfies that kρk k2 ≤ kρ0 k2 , k = 1, 2, . . . , K. This completes the proof of Theorem 3.5. 2 10
4
Solvability of the numerical simulation method
In this section, we discuss solvability of the numerical simulation method (8)-(10) . If define that Uk = [uk1 , uk2 , . . . , ukJ−1 ]T , then on the k layer, the numerical simulation method (8)-(10) can be written in the following matrix form: AUk = bk ,
(17)
where the right term bk is a (J − 1) dimension column vector witch includes the known initial and boundary values, the values of known term fjk and the before k layer known values unj (n = 0, 1, . . . , k − 1), whereas the coefficient matrix is a (J − 1) × (J − 1) known matrix:
e a a
A=
where
e a
a ae ... ... ... ae a
, e a
ae a
5 5 1 1 + 2µ1 + µ2 , ae = − µ1 − µ2 . 6 6 12 12 Theorem 4.1. The numerical simulation method (8)-(10) is uniquely solvable. Proof. It is easy to see that A is a strictly diagonally dominant matrix, then A is a nonsingular matrix, therefore the linear equation system (17) exists unique solution, i.e., the numerical simulation method (8)-(10) is uniquely solvable. 2 a=
5
Convergence of the numerical simulation method
In this section, we discuss convergence of the numerical simulation method (8)-(10). Subtracting (8) from (5), we obtain the following error equation Lx Ejk − Lx Ejk−1 11
= µ1
k X
λl δx2 Ejk−l
l=0
− µ2
k X
λl Lx Ejk−l + Rjk ,
(18)
l=0
k = 1, 2, . . . , K; j = 1, 2, . . . , J − 1, where Ejk = u(xj , tk ) − ukj . Let that h
k E k = E1k , E2k , . . . , EJ−1
iT
h
k , Rk = R1k , R2k , . . . , RJ−1
iT
.
Similar to the derivation in [5], we have that, respectively k
kE k2 ≡
J−1 X
j=1
1 2
4x |Ejk |2
=
∞ X
l=−∞
1 2
2
|αk (l)|
, k = 0, 1, . . . , K
(19)
, k = 0, 1, . . . , K,
(20)
and k
kR k2 ≡
J−1 X
j=1
1
2
4x |Rjk |2
=
∞ X
l=−∞
1 2
2
|βk (l)|
where αk (l) and βk (l) are Fourier expansion coefficients of E k (x) and Rk (x), whereas E k (x) and Rk (x) are the following grid functions, respectively
and
Ek,
when x ∈ xj− 1 , xj+ 1 , j = 1, 2, . . . , J − 1;
Rk ,
when x ∈ xj− 1 , xj+ 1 , j = 1, 2, . . . , J − 1;
E k (x) =
Rk (x) =
j
0, when x ∈ [0, j
0, when x ∈ [0,
2
4x S ] (L 2 2
4x S ] (L 2
2
i
− 2
i
−
4x , L], 2
4x , L]. 2
We now suppose that, respectively Ejk = αk eiσj4x ,
Rjk = βk eiσj4x ,
where σ = 2πl/L. By substituting the above expressions into (18), immediately !
!
σ4x 1 σ4x 1 1 − sin2 αk − 1 − sin2 αk−1 3 2 3 2
!
k k σ4x X 1 σ4x X = −4µ1 sin2 λl αk−l + βk , (21) λl αk−l − µ2 1 − sin2 2 l=0 3 2 l=0
12
k = 1, 2, . . . , K. For convenience, we rewrite equation (21) as follows "
k σ4x X 1 σ4x 1 Bαk−1 − 4µ1 sin2 λl αk−l − µ2 1 − sin2 αk = A 2 l=2 3 2 k X
#
λl αk−l + βk , k = 1, 2, . . . , K.
l=2
!
(22)
From (7), then there is a positive constant C1 , such that
|Rjk | ≤ C1 42t + 44x , k = 1, 2, . . . , K, again use the first equality of (20), result in √ kRk k2 ≤ C1 L 42t + 44x , k = 1, 2, . . . , K.
(23)
(24)
Based on the convergence of the series in RHS of (20), then there is a positive constant C2 , such that |βk | ≡ |βk (l)| ≤ C2 4t |β1 (l)| ≡ C2 4t |β1 |, k = 1, 2, . . . , K.
(25)
Theorem 5.1. Assume that u(x, t) ∈ U (Ω) is the exact solution of the 2 ∈ C(Ω), and the temporal and spatial steps satisfy problem (1)-(3), ∂ f∂t(x,t) 2 4 2 that 4t = 4x , if λ1 < 0, or if 0 ≤ λ1 ≤ 6µ11+µ2 , then the numerical simulation method (8)-(10) is convergent with the order O(42t + 44x ). Proof. At first, we verify the following conclusion by induction: |αk | ≤ C2 k4t |β1 |, k = 1, 2, . . . , K, where αk (k = 1, 2, . . . , K) be the solution of (22). For k = 1, from Eq. (22) gets 1 β1 , A applying Lemma 3.2 and (25), we obtain that α1 =
|α1 | =
1 |β1 | ≤ |β1 | ≤ C2 4t |β1 |. A 13
(26)
Suppose that |αn | ≤ C2 n4t |β1 |, n = 1, 2, . . . , k − 1, in view of Lemma 3.1, Lemma 3.2 and Lemma 3.3, from Eq. (22), we have that
!
k k X 1 2 σ4x X 1 2 σ4x λl αk−l − µ2 1 − sin λl αk−l + βk αk = Bαk−1 − 4µ1 sin A 2 l=2 3 2 l=2
"
!
#
k k X 1 1 2 σ4x X 2 σ4x ≤ B|αk−1 | + 4µ1 sin |λl | |αk−l | + µ2 1 − sin |λl | |αk−l | + |βk | A 2 l=2 3 2 l=2
"
!
k k σ4x X σ4x X 1 1 B(k − 1) + 4µ1 sin2 |λl | (k − l) + µ2 1 − sin2 |λl |( k − l) + 1 ≤ A 2 l=2 3 2 l=2 C2 4t |β1 |
(
"
!
#
k k 1 σ4x X 1 σ4x X ≤ B + 4µ1 sin2 |λl | + µ2 1 − sin2 |λl | (k − 1) + 1 A 2 l=2 3 2 l=2 C2 4t |β1 |
=
(
(
"
!
#
)
)
k k 1 1 σ4x X σ4x X B − 4µ1 sin2 λl − µ2 1 − sin2 λl (k − 1) + 1 C2 4t |β1 | A 2 l=2 3 2 l=2
"
!#
#
)
1 1 σ4x σ4x ≤ B + 4µ1 (1 + λ1 ) sin2 + µ2 (1 + λ1 ) 1 − sin2 (k − 1) + 1 C2 4t |β1 | A 2 3 2 ! # ) ( " 1 2 σ4x 1 2 σ4x (1 + µ2 ) 1 − sin + 4µ1 sin (k − 1) + 1 C2 4t |β1 | = A 3 2 2 ≤ C2 k4t |β1 |. From this the conclusion (26) holds. Further combining (19), (20) and (26), we have that kE k k2 ≤ C2 k4t kR1 k2 ≤ C(42t + 44x ), √ where C = C1 C2 LT. This completes the proof of Theorem 5.1. 2 Remark 1. In fact, because there are many numerical integration formulas with high order accuracy, to continue to structure numerical method with high temporal accuracy for solving the reaction-subdiffusion equation (1) is not difficult. However, stability and convergence analysis may be very difficult.
6
Numerical experiment
In order to certificate the theoretical analysis results, we employ the numerical simulation method (8)-(10) to solve the following initial-boundary value problem of the time fractional reaction-subdiffusion equation 14
∂u(x, t) =0 Dt1−γ ∂t
!
∂ 2 u(x, t) − u(x, t) + 2ex t, ∂x2
(27)
0 < t ≤ 1, 0 < x < 1, u(x, 0) = 0, 0 ≤ x ≤ 1.
(28)
u(0, t) = t2 , u(1, t) = et2 , 0 ≤ t ≤ 1,
(29)
The exact solution of the problem (27)-(29) is u(x, t) = ex t2 . Let that the maximum error of the numerical solutions E∞ = max
0≤k≤K
n
o
kE k k2 = O(4ζt + 4ηx ).
Suppose that in the front experiment, the maximum error, temporal step, spa(f ) (f ) ) tial step of the numerical solutions are E∞ , 4t , 4(f x , respectively, whereas in the subsequent experiment, the maximum error, temporal step, spatial step (s) (f ) , 4t , 4(s) of the numerical solutions are E∞ x , respectively, then the temporal accuracy ζ and spatial accuracy η calculated by the following formulas (s)
(s)
ζ=
ln E∞ (f ) E∞
ln
(s)
4t
,
η=
(f )
4t
ln E∞ (f ) E∞
(s)
x ln 4(f ) 4x
.
Table 1 exhibits the maximum error E∞ , temporal accuracy ζ and spatial accuracy η of the numerical simulation method (8)-(10) for solving for the problem (27)-(29) versus various γ and 42t = 44x , respectively. From Table 1, it can be seen that our theoretical analysis has been strongly supported by the numerical experiments. Fig. 1 reveals the comparison of the numerical solution using the numerical simulation method (8)-(10) and the exact solution for the problem (27)-(29) 1 . at t = 14 , 24 , 34 , 1, respectively, for γ = 0.5 and 42t = 44x = 256 Fig. 2 shows the comparison of the numerical solution using the numerical simulation method (8)-(10) and the exact solution for the problem (27)-(29) 1 at x = 14 , 42 , 34 , respectively, for γ = 0.5 and 42t = 44x = 256 . Fig. 3 illustrates the absolute error E(x, t) of the numerical solution for the problem (27)-(29) using the numerical simulation method (8)-(10) for γ = 0.5 15
and 42t = 44x =
1 , 256
where E(x, t) = u(x, t) − ue(x, t),
whereas ue(x, t) is the numerical approximation for u(x, t).
Fig. 1, 2 and 3 show that the numerical solution are very consistent with the exact solution. Table 1 The maximum error E∞ , temporal accuracy ζ and spatial accuracy η of the numerical simulation method (8)-(10) for solving for the problem (27)-(29)
42t = 44x =
γ 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
1 16
3.1544 × 10−5 3.1292 × 10−5 3.1039 × 10−5 3.0701 × 10−5 3.0364 × 10−5 3.0111 × 10−5 2.9774 × 10−5 2.9521 × 10−5 2.9184 × 10−5
42t = 44x =
1 81
6.5689 × 10−6 6.5689 × 10−6 6.4718 × 10−6 6.4718 × 10−6 6.3230 × 10−6 6.2776 × 10−6 6.2258 × 10−6 6.1287 × 10−6 6.0834 × 10−6
ζ
η
1.93 3.87 1.93 3.85 1.93 3.87 1.92 3.84 1.93 3.87 1.93 3.87 1.93 3.86 1.94 3.88 1.93 3.87
42t = 44x =
1 256
2.2508 × 10−6 2.2231 × 10−6 2.2231 × 10−6 2.1832 × 10−6 2.1832 × 10−6 2.0464 × 10−6 2.0464 × 10−6 2.0066 × 10−6 2.0066 × 10−6
ζ
1.86 3.72 1.88 3.77 1.86 3.71 1.89 3.78 1.85 3.70 1.95 3.90 1.93 3.87 1.94 3.88 1.93 3.86
Remark 2. In section 3 and 5, in order to discuss stability and convergence of the numerical simulation method, we suppose that holds the condition: λ1 < 0, or 0 ≤ λ1 ≤ 6µ11+µ2 , this condition is effectively equivalent to that some γ are restricted. When this condition is not holds, we are not able to confirm the stability and convergence of the numerical simulation method. However, the numerical experiment in this section shows that for the wide range of γ, there are no symptoms of numerical instability.
7
Conclusion
In this paper, the numerical simulation method with second order temporal accuracy and fourth order spatial accuracy has been proposed to solve the time fractional reaction-subdiffusion equation; by Fourier analysis and 16
η
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
0.6
0.8
1
x
Fig. 1. The comparison of the numerical solution (n.s.) using the numerical simulation method (8)-(10) and the exact solution (e.s.) for the problem (27)-(29) at 1 t = 14 , 24 , 34 , 1, respectively, for γ = 0.5 and 42t = 44x = 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
0.6
0.8
1
t
Fig. 2. The comparison of the numerical solution (n.s.) using the numerical simulation method (8)-(10) and the exact solution (e.s.) for the problem (27)-(29) at 1 x = 14 , 42 , 34 , respectively, for γ = 0.5 and 42t = 44x = 256 .
algebraic theory, the stability, convergence and solvability of the numerical simulation method have been discussed, respectively; the theoretical analysis results very consistent with the numerical experiment. The numerical simu17
−6
x 10
Absolute error |E(x,t)|
3 2.5 2 1.5 1 0.5 0 1 1 0.8
0.5
0.6 0.4 t
0
0.2 0
x
Fig. 3. The absolute error of the numerical solution for the problem (27)-(29) using 1 the numerical simulation method (8)-(10) for γ = 0.5 and 42t = 44x = 256 .
lation method and numerical analysis method in this paper can be extended to other fractional partial deferential equations with the Riemann-Liouville fractional partial derivative. Acknowledgments The authors would like to thank the reviewers for their careful reading of this article and their comments and suggestions which have improved the article.
References [1] E. Abad, S. B. Yuste, K. Lindenberg, Reaction-subdiffusion and reactionsuperdiffusion equations for evanescent particles performing continuous-time random walks, Phys. Rev. E 81 (2010) 031115. [2] M. Abbaszadeh, M. Dehghan, A meshless numerical procedure for solving fractional reaction subdiffusion model via a new combination of alternating direction implicit (ADI) approach and interpolating element free Galerkin (EFG) method, Comput. Math. Appl. 70 (2015) 2493-2512. [3] B. Baeumer, M. Kovcs, M. Meerschaert, Numerical solutions for fractional reaction-diffusion equations, Comput. Math. Appl. 55 (2008) 2212–2226. [4] V. K. Baranwal, R. K. Pandey, M. P. Tripathi, O. P. Singh, An analytic algorithm for time fractional nonlinear reaction-diffusion equation based on a new iterative method, Commun. Nonlinear. Sci. Numer. Simulat. 17 (2012) 3906-3921.
18
[5] C.-M. Chen, F. Liu, I. Turner, and V. Anh , Fourier method for the fractional diffusion equation describing sub-diffusion, J. Comp. Phys. 227 (2007) 886-897. [6] C.-M. Chen, F. Liu, K. Burrage, Finite difference methods and a Fourier analysis for the fractional reaction-subdiffusion equation, Appl. Math. Comp. 198 (2008) 754-769. [7] C.-M. Chen, F. Liu, V. Anh, I. Turner, Numerical schemes with high spatial accuracy for a variable-order anomalous subdiffusion equation, SIAM J. Sci. Comput. 32 (4) (2010), 1740-1760. [8] C.-M. Chen, F. Liu, I. Turner, V. Anh, Y. Chen, Numerical approximation for a variable-order nonlinear reaction-subdiffusion equation, Numer. Algor. 63 (2013) 265-290. [9] J. W. Chiu, K. -H. Chiam, Monte Carlo simulation and linear stability analysis of Turing pattern formation in reaction-subdiffusion systems. Phys. Rev. E 78 (2008) 056708. [10] S. Das, P. K. Gupta, P. Ghosh, An approximate solution of nonlinear fractional reaction-diffusion equation, Appl. Math. Model. 35 (2011) 4071-4076. [11] E. H. Doha, A. H. Bhrawy, S. S. Ezz-Eldien, An efficient Legendre spectral tau matrix formulation for solving fractional subdiffusion and reaction subdiffusion equations, J. Comput. Nonlin. Dyn. 10 (2015) DOI: 10.1115/1.4027944. [12] M. Dehghan, M. Abbaszadeh, A. Mohebbi, Error estimate for the numerical solution of fractional reaction-subdiffusion process based on a meshless method, J. Comput. Appl. Math. 280 (2015) 14-36. [13] H. F. Ding, C. P. Li, Mixed spline function method for reaction-subdiffusion equations, J. Comp. Phys. 242 (2013) 103-123. [14] V. Gafiychuk, B. Y. Datsko, Stability analysis and oscillatory structures in time-fractional reaction-diffusion systems, Phys. Rev. E 75 (2007) 055201. [15] V. Gafiychuk, B. Datsko, V. Meleshko, Mathematical modeling of time fractional reaction-diffusion systems, J. Comput. Appl. Math. 220 (2008) 215225. [16] V. Gafiychuk, B. Datsko, Inhomogeneous oscillatory solutions in fractional reaction-diffusion systems and their computer modeling, Appl. Math. Comp. 198 (2008) 251-260. [17] V. Gafiychuk, B. Datsko, V. Meleshko, D. Blackmore, Analysis of the solutions of coupled nonlinear fractional reaction-diffusion equations, Chaos Soliton. Fract. 41 (2009) 1095-1104. [18] V. Gafiychuk, B. Datsko, Mathematical modeling of different types of instabilities in time fractional reaction-diffusion systems, Comput. Math. Appl. 59 (2010) 1101-1107.
19
[19] C. Y. Gong, W. M. Bao, G. J. Tang, A parallel algorithm for the Riesz fractional reaction-diffusion equation with explicit finite difference method, Fract. Calc. Appl. Anal. 16 (2013) 654-669. [20] H. J. Haubold, A. M. Mathai, R. K. Saxena, Solutions of fractional reactiondiffusion equations in terms of the H-function, Bull. Astron. Soc. India 35 (2007) 681-689. [21] H. J. Haubold, A. M. Mathai, R. K. Saxena, Further solutions of fractional reaction-diffusion equations in terms of the H-function, J. Comput. Appl. Math. 235 (2011) 1311-1316. [22] J. M. Haugh, Analysis of reaction-diffusion systems with anomalous subdiffusion, Biophys. J. 97 (2009) 435-442. [23] B. I. Henry, S. L. Wearne, Fractional reaction-diffusion, Phys. A 276 (2000) 448-455. [24] H. Huang, X. Cao, Numerical method for two dimensional fractional reaction subdiffusion equation, Eur. Phys. J.-Spec. Top. 222 (2013) 1961-1973. [25] T. A. M. Langlands, B. I. Henry, S. L. Wearne, Anomalous subdiffusion with multispecies linear reaction dynamics, Phys. Rev. E 77(2008) 021111. [26] M. S. Mommer, D. Lebiedz, Modeling subdiffusion using reaction diffusion systems, SIAM J. Appl. Math. 70 (2009) 112-132. [27] S. Z. Rida, A. M. A. El-Sayed, A. A. M. Arafa, On the solutions of time-fractional reaction-diffusion equations, Commun. Nonlinear. Sci. Numer. Simulat. 15 (2010) 3847-3854. [28] F. Sagus, V. P. Shkilev, I. M. Sokolov, Reaction-subdiffusion equations for the A B reaction. Phys. Rev. E 77 (2008) 032102. [29] R. K. Saxena, A. M. Mathai, H. J. Haubold, Fractional reaction-diffusion equations, Astrophys. Space Sci. 305(2006) 289-296. [30] R. K. Saxena, A. M. Mathai, H. J. Haubold, Solution of generalized fractional reaction-diffusion equations, Astrophys. Space Sci. 305 (2006) 305-313. [31] H. H. Schmidt-Martens, D. Froemberg, I. M. Sokolov, Front propagation in a one-dimensional autocatalytic reaction-subdiffusion system. Phys. Rev. E 79 (2009) 041135. [32] K. Seki, M. Wojcik, M. Tachiya, Fractional reaction-diffusion equation, J. Chem. Phys. 119 (2003) 2165-2170. [33] I. M. Sokolov, M. G. W. Schmidt F. Sagus, Reaction-subdiffusion equations. Phys. Rev. E 73 (2006) 031102. [34] N. H. Sweilam, M. M. Khader, M. Adel, Weighted average finite difference methods for fractional reaction-subdiffusion equation, Walailak J. Sci. Tech. 11 (2014) 361-377.
20
[35] S. Q. Wang, J. H. He, Variational iteration method for a nonlinear reactiondiffusion process, J. Chem. Reactor Eng. 6 (2008) A37. [36] A. Yadav, Werner Horsthemke: Kinetic equations for reaction-subdiffusion systems: Derivation and stability analysis. Phys. Rev. E 74 (2006) 066118. [37] B. Yu, X. Y. Jiang, H. Y. Xu, A novel compact numerical method for solving the two-dimensional non-linear fractional reaction-subdiffusion equation, Numer. Algor. 68 (2015) 923-950. [38] S. B. Yuste, K. Lindenberg, Subdiffusion-limited A+A reactions, Phys. Rev. Lett. 87 (11) (2001) art. no.-118301. [39] S. B. Yuste, L. Acedo, Katja Lindenberg: Reaction front in an A+B →C reaction-subdiffusion process. Phys. Rev. E 69 (2004) 036126. [40] S. B. Yuste, E. Abad, K. Lindenberg, Reaction-subdiffusion model of morphogen gradient formation, Phys. Rev. E 82 (2010) 061123. [41] P. Zhuang, F. Liu, V. Anh, I. Turner, Stability and convergence of an implicit numerical method for the non-linear fractional reaction-subdiffusion process, IMA J. Appl. Math. 74 (2009) 645-667. [42] A. Zoia, Continuous-time random-walk approach to normal and anomalous reaction-diffusion processes, Phys. Rev. E 77 (2008) 041115.
21