- Email: [email protected]
Implicit compact difference schemes for the fractional cable equation
Applied Mathematical Modelling 36 (2012) 4027–4043
Contents lists available at SciVerse ScienceDirect
Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm
Implicit compact difference schemes for the fractional cable equation q Xiuling Hu a,b, Luming Zhang a,⇑ a b
College of Science, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China School of Mathematical Sciences, Xuzhou Normal University, Xuzhou 221116, China
a r t i c l e
i n f o
Article history: Received 28 April 2011 Received in revised form 29 October 2011 Accepted 2 November 2011 Available online 17 November 2011 Keywords: Fractional cable equation Compact finite difference scheme Stability Convergence
a b s t r a c t In this paper, we propose two implicit compact difference schemes for the fractional cable equation. The first scheme is proved to be stable and convergent in l1-norm with the convergence order O(s + h4) by the energy method, where new inner products defined in this paper gives great convenience for the theoretical analysis. Numerical experiments are presented to demonstrate the accuracy and effectiveness of the two compact schemes. The computational results show that the two new schemes proposed in this paper are more accurate and effective than the previous. 2011 Elsevier Inc. All rights reserved.
1. Introduction Fractional differential and integral equations have been applied to numerous fields of science, such as physics [1,2], engineering [3,4], finance [5], medicine [6], and so on. Many scholars have been attracted to study them using kinds of ways in recent years. For most of the systems that cannot be solved analytically, some numerical methods including the finite difference methods [7–13] have been considered by researchers. Chen et al. [14,15] established the finite difference schemes based on Grünwald–Letnikov approximation and analyzed the schemes by Fourier method for the fractional sub-diffusion and reaction–subdiffusion equation. Langlands and Henry [16] constructed an unconditional implicit scheme based on L1 approximation for the fractional subdiffusion but without global convergence analysis. Zhuang et al. [17] proposed a new implicit numerical method (INM) by integrating method for an anomalous subdiffusion equation and presented two techniques for improving the convergence order of the scheme and explored the stability and convergence by the energy method. Then they proposed an implicit difference scheme by using the same method for nonlinear fractional reaction–subdiffusion equation [18]. Liu et al. [19] studied a modified anomalous subdiffusion equation with a nonlinear source term and constructed a new implicit difference scheme and also analyzed the stability and convergence in l2-norm by a new energy method. Chen et al. [20] also investigated two-dimensional anomalous subdiffusion. Sun and Wu [21] explored a fractional diffusion-wave system and proposed a fully discrete scheme by the method of order reduction and proved the solvability, stability and convergence by the energy method. Because compact finite difference scheme is developed to approximate the second-order derivatives and keep the desirable trigonal nature, the compact method has been received great attention. Cui [22] considered a spacial fourth order finite difference scheme for the one-dimensional fractional diffusion equation and analyzed the stability by the Fourier method. Du et al. [23] presented a compact scheme for the time Caputo’s sense fractional diffusion-wave system on the basis of q
This work is partially supported by the National Science Foundation of China under grant 11102179.
⇑ Corresponding author.
E-mail addresses: [email protected] (X. Hu), [email protected] (L. Zhang). 0307-904X/$ - see front matter 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2011.11.027
4028
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043
[21] and prove the stability and convergence by energy method. Guanghua et al. [24] provided a compact finite difference scheme by applying L1 discretization for the fractional derivative and fourth order approximation for the second order space derivative for the fractional subdiffusion equation after a transformation of the original problem and discussed its unique solvability, stability and convergence by introducing a new inner product. Chen et al. [25] proposed a scheme with order O(s + h4) for a variable-order anomalous subdiffusion equation and analyzed the convergence, stability and solvability by the Fourier method. Then an improved scheme with the order O(s2 + h4) are explored. The cable equation is one of the most fundamental equations for modeling neuronal dynamics. Sun and Zhao [30] and Chaos and Sun [31] found a fractional cable equation involving two fractional temporal Riemann–Liouville derivatives by assuming small ionic concentration gradients and simplifying the fractional Nernst–Planck equations and investigated solutions of a fractional cable equation for finite case and infinite case, respectively. The fractional cable equation (FCE) was introduced by Liao and Sun [32] from the fractional Nernst–Planck equations for modeling the anomalous electrodiffusion of ions in spiny neuronal dendrites. Langlans et al. [26] proposed two implicit numerical methods for the fractional cable equation with l2-convergence order O(s + h2) and O(s2 + h2). Because the l2-norm depends on the spatial mesh step, it cannot portray the degree of the closeness of the numerical solution to the exact solution better than the l1-norm. Recently, the maximum problem has also been paid to much attention [27–29]. This paper studies the fractional cable equation and proposes two implicit compact schemes and prove the l1-stability and convergency of the first scheme by the energy method. The following of this paper is organized as follows. In Section 2, an implicit compact finite difference scheme is derived. In Section 3, the l1-stability of the compact scheme is proved by the energy method where new inner products are introduced. Then the l1-convergency of the compact scheme is proved in Section 4. In Section 5, an improved implicit compact finite difference scheme is derived and two illustrative numerical examples are carried out to demonstrate the accuracy and effectiveness of the two schemes in Section 6. In the end, some conclusions are made in Section 7. 2. An implicit compact finite difference scheme Consider the following fractional cable equation:
! @u @ 2 uðx; tÞ 1c1 1c ¼ 0 Dt a b0 Dt 2 uðx; tÞ þ f ðx; tÞ; @t @x2 uð0; tÞ ¼ uðtÞ;
uðL; tÞ ¼ wðtÞ;
uðx; 0Þ ¼ gðxÞ;
x 2 ½0; L;
1c uðx; tÞ 0 Dt
¼
1
@
CðcÞ @t
ð1Þ
t 2 ð0; T;
1c uðx; tÞ 0 Dt
t
0
t 2 ð0; T;
ð2Þ ð3Þ
where 0 < c1, c2 < 1, a, b > 0, and
Z
x 2 ð0; LÞ;
uðx; nÞ ðt nÞ1c
is the Rieman–Liouville fractional partial derivative of order 1 c defined by
dn;
ð4Þ
Firstly, we denote xj = jh, tn = ns, Xh = {xjj0 6 j 6 M}, Xs = {tnj0 6 n 6 N}, and Xsh ¼ Xh Xs ; where h = L/M, s = T/N are the uniform spatial and temporal mesh sizes respectively, and M, N are two positive integers. In the rest of the paper, C denotes a general positive constant, which may have different values in different occurrences. Suppose U ¼ fU nj j0 6 j 6 M; 0 6 n 6 Ng be a grid function on Xsh : Introduce the following notations:
unj ¼ uðxj ; tn Þ; U nj uðxj ; t n Þ; 1 n 1 U U nj ; d2x U nj ¼ 2 U njþ1 2U nj þ U nj1 U nj ¼ x h jþ1 h Let Uh = {UjU = (U0, U1, . . . , UM), U0 = UM = 0}. For V 2 Uh, U 2 Uh, we firstly introduce the following usual discrete l2 inner product (, ), discrete l2-norm kk and discrete l1-norm kk1:
ðU n ; V n Þ ¼ h
M 1 X
U nj V nj ;
kU n k2 ¼ ðU n ; U n Þ;
kU nx k2 ¼ h
j¼1
M 1 X j¼0
U nj
2 x
;
kU n k1 ¼ max jU nj j: 16j6M1
Lemma 2.1 [17]. Let ðcÞ
c
bk ¼ ðk þ 1Þc k ; then
ðcÞ bk
k ¼ 0; 1; . . . ;
ð5Þ
satisfies
ðcÞ
ðcÞ
(1) b0 ¼ 1; bk > 0; k ¼ 0; 1; . . . ðcÞ
ðcÞ
(2) bk > bkþ1 ; k ¼ 0; 1; . . . (3) There exists a positive constant C > 0, such that Pk ðcÞ c c c c (4) j¼0 bj s ¼ ðk þ 1Þ s 6 T :
s 6 CbkðcÞ sc ; k ¼ 1; 2; 3 . . .
4029
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043
Assume that problem (1)–(3) has a smooth solution uðx; tÞ 2 C 7;3 x;t ð½0; L ½0; TÞ; and denote
Lu ¼
@ 2 uðx; tÞ : @x2
ð6Þ
From Taylor expansion, we can easily know [22] 4
Luðxi ; tn Þ ¼
@ 2 uðxi ; tn Þ d2x h @ 4 uðxi ; t n Þ 6 ¼ uðxi ; tn Þ þ þ Oðh Þ: 2 2 2 1 @x @x4 240 1 þ h dx
ð7Þ
12
Now, we consider the compact scheme for the problem (1)–(3). Integrating both sides of Eq. (1) from tn to tn+1, we have
uðxi ; tnþ1 Þ uðxi ; t n Þ ¼
a Cðc1 Þ
Z
t nþ1
Luðxi ; nÞ 1c1
ðt nþ1 nÞ
0
dn
a Cðc1 Þ
Z
tn
0
Luðxi ; nÞ ðtn nÞ
1c1
dn
b Cðc2 Þ
Z 0
tnþ1
uðxi ; nÞ ðt nþ1 nÞ1c2
dn
Z tnþ1 Z tn b uðxi ; nÞ dn þ f ðxi ; nÞdn Cðc2 Þ 0 ðt n nÞ1c2 tn Z s Z tn Z s a Luðxi ; nÞ a Luðxi ; n þ sÞ Luðxi ; nÞ b uðxi ; nÞ ¼ dn þ dn dn 1c1 1c2 Cðc1 Þ 0 ðtnþ1 nÞ1c1 Cðc1 Þ 0 C ð c Þ ðt n nÞ 0 ðt nþ1 nÞ 2 Z tnþ1 Z tn b uðxi ; n þ sÞ uðxi ; nÞ dn þ f ðxi ; nÞdn ¼ I1 I2 þ I3 ; ð8Þ Cðc2 Þ 0 ðt n nÞ1c2 tn þ
where
I1 ¼ I2 ¼ I3 ¼
a Cðc1 Þ b Cðc2 Þ Z tnþ1
Z s 0
Z s 0
Luðxi ; nÞ 1c1
ðt nþ1 nÞ
uðxi ; nÞ 1c2
ðt nþ1 nÞ
dn þ dn þ
a Cðc1 Þ b Cðc2 Þ
Z
tn
Luðxi ; n þ sÞ Luðxi ; nÞ ðt n nÞ1c1
0
Z
tn
uðxi ; n þ sÞ uðxi ; nÞ
0
ðt n nÞ1c2
dn;
dn;
ð9Þ
f ðxi ; nÞdn:
tn
According to Taylor expansion and (7), for I1 ; we have
I1 ¼
a Cðc1 Þ
Z s 0
Luðxi ; sÞ ðtnþ1 nÞ
1c1
dn þ R11 þ
n1 a X Cðc1 Þ j¼0
Z
t jþ1
Luðxi ; n þ sÞ Luðxi ; nÞ
tj
ðtn nÞ1c1
dn
¼
n1 asc1 asc1 X ðc Þ ðc Þ b 1 Luðxi ; t jþ2 Þ Luðxi ; tjþ1 Þ þ R12 bn 1 Luðxi ; t 1 Þ þ R11 þ Cð1 þ c1 Þ Cð1 þ c1 Þ j¼0 nj1
¼
asc1 d2x asc1 ðc Þ uðxi ; t 1 Þ þ R13 þ R11 þ bn 1 2 2 h Cð1 þ c1 Þ Cð1 þ c1 Þ 1 þ 12 dx
"
n1 X
ðc1 Þ bnj1
"
#
d2x
d2x
#
uðxi ; t jþ2 Þ uðxi ; tjþ1 Þ þ R14 þ R12 h2 2 h2 2 1 þ 12 dx 1 þ 12 dx " # " # n1 asc1 d2x asc1 X d2x d2x ðc Þ ðc1 Þ ¼ uðx ; t Þ þ b uðx ; t Þ uðx ; t Þ þ R1 ; bn 1 i 1 i jþ2 i jþ1 h2 2 h2 2 Cð1 þ c1 Þ Cð1 þ c1 Þ j¼0 nj1 1 þ h2 d2x 1 þ 12 dx 1 þ 12 dx 12
j¼0
ðcÞ
where bj
ð10Þ
satisfies (5) and
R1 ¼ R11 þ R12 þ
n1 asc1 asc1 X ðc Þ ðc Þ b 1 R : bn 1 R13 þ Cð1 þ c1 Þ Cð1 þ c1 Þ j¼0 nj1 14
ð11Þ
For R11, from Taylor expansion, we have
" # a Z s a Z s 1 1 @ 3 uðxi ; n1 Þ jR11 j ¼ ½Luðxi ; nÞ Luðxi ; sÞdn ¼ ðn sÞ dn 1 c 1 c 2 Cðc1 Þ 0 ðt nþ1 nÞ 1 Cðc1 Þ 0 ðt nþ1 nÞ 1 @x @t Z 3 as @ uðxi ; n1 Þ s 1 ðc Þ 6 dn 6 C s1þc1 bn 1 ; 0 < n < n1 < s: Cðc1 Þ @x2 @t 0 ðt nþ1 nÞ1c1
ð12Þ
4030
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043
For R12, from Taylor expansion, we have
R12 ¼
n1 a X Cðc1 Þ j¼0
Z
t jþ1
½Luðxi ; n þ sÞ Luðxi ; nÞ ½Luðxi ; t jþ2 Þ Luðxi ; t jþ1 Þ ðt n nÞ1c1
tj
dn
" # @ 3 uðxi ; n2 Þ @ 3 uðxi ; n3 Þ ðt jþ2 n sÞ ðt jþ1 nÞ ðt n nÞc1 1 dn @x2 @t @x2 @t tj " # n1 Z t jþ1 a X @ 3 uðxi ; n2 Þ @ 3 uðxi ; n3 Þ ðtjþ1 nÞðtn nÞc1 1 dn ¼ @x2 @t @x2 @t Cðc1 Þ j¼0 tj
n1 a X ¼ Cðc1 Þ j¼0
¼
n1 a X Cðc1 Þ j¼0
Z
t jþ1
Z
t jþ1 tj
@ 4 uðxi ; n4 Þ ðn2 n3 Þðtjþ1 nÞðtn nÞc1 1 dn; @x2 @t 2
where
n þ s < n2 < tjþ2 ;
tj < n < t jþ1 ;
n < n3 < t jþ1 ;
n3 < n4 < n2 :
Thus,
tjþ1 < n < tj ;
tjþ1 < n3 < n:
So we have
0 ¼ tj þ s t jþ1 < n2 n3 < tjþ2 n < tjþ2 tj < 2s;
0 < t jþ1 n < t jþ1 tj ¼ s:
Then we can get
a X X n1 Z t jþ1 n1 Z tjþ1 @ 4 uðxi ; n4 Þ c1 1 c1 1 2 6 C jR12 j ¼ ðn n Þðt nÞðt nÞ dn s ðt nÞ dn n n jþ1 2 3 2 2 Cðc1 Þ j¼0 tj j¼0 tj @x @t Z t n ¼ C s2 ðt n nÞc1 1 dn 6 C s2 : 0
And from Lemma 2.1 (3), we have ðc Þ
jR12 j 6 C s1þc1 bn 1 :
ð13Þ
From (7), we have 4
jR13 j 6 Ch ;
ð14Þ
and
d2x d2x jR14 j ¼ Luðxi ; tjþ2 Þ Luðxi ; t jþ1 Þ uðx ; t Þ þ uðx ; t Þ i jþ2 i jþ1 2 2 2 2 h h 1 þ 12 dx 1 þ 12 dx 4 h4 @ 4 uðx ; t Þ 4 6 6 h @ uðxi ; t jþ1 Þ 6 @ uðxi ; t jþ2 Þ 6 @ uðxi ; t jþ1 Þ i jþ2 ¼ þ Ch Ch þ 4 4 6 6 240 @x @x @x @x 240 4 sh @ 5 uðx ; gÞ 7 6 @ uðxi ; g1 Þ 4 i ¼ þ C sh þ 6 C sh ; tjþ1 < g; g1 < tjþ2 : 240 @x4 @t @x6 @t Again from Lemma 2.1 (3), we have ðc Þ 4
jR14 j 6 C sc1 bn 1 h :
ð15Þ
So from (12)–(15) and Lemma 2.1 (1), (2) and (4), we can get
n1 asc1 asc1 X ðc1 Þ ðc1 Þ jR1 j ¼ R11 þ R12 þ bnj1 R14 bn R13 þ Cð1 þ c1 Þ Cð1 þ c1 Þ j¼0 ðc Þ
ðc Þ
ðc Þ 4
ðc Þ 4
ðc Þ
4
6 C s1þc1 bn 1 þ C s1þc1 bn 1 þ C sc1 bn 1 h þ C sc1 bn 1 h 6 C sc1 bn 1 ðs þ h Þ:
ð16Þ
Similarly, for I2 ; we have
I2 ¼
n1 X bsc2 bsc2 ðc Þ ðc Þ b 2 ½uðxi ; t jþ2 Þ uðxi ; t jþ1 Þ þ R2 : bn 2 uðxi ; t 1 Þ þ Cð1 þ c2 Þ Cð1 þ c2 Þ j¼0 nj1
ð17Þ
4031
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043
where ðc Þ 1þc2
s
jR2 j 6 Cbn 2
ð18Þ
:
For I3 , we have
I3 ¼
Z
t nþ1
f ðxi ; nÞdn ¼
tn
s 2
½f ðxi ; tn Þ þ f ðxi ; t nþ1 Þ þ R3 :
ð19Þ
From Taylor expansion and Lemma 2.13, we have ðc Þ 1þc2
jR3 j 6 C s2 6 Cbn 2
s
ð20Þ
:
Thus, we can obtain ðc Þ
uðxi ; tnþ1 Þ uðxi ; t n Þ ¼ r 1 bn 1
(
d2x
uðxi ; t 1 Þ þ r 1 h2 2 1 þ 12 dx ðc Þ
r 2 bn 2 uðxi ; t1 Þ þ
n1 X
n1 X
ðc Þ
"
1 bnj1
j¼0
d2x
d2x uðxi ; t jþ2 Þ uðxi ; tjþ1 Þ 2 2 h h2 2 1 þ 12 dx 1 þ 12 dx )
ðc Þ
2 bnj1 ½uðxi ; t jþ2 Þ uðxi ; t jþ1 Þ
j¼0
#
s
þ ½f ðxi ; t nþ1 Þ þ f ðxi ; t n Þ þ R1 þ R2 þ R3 ; 2 ð21Þ
or
1þ
! ! 2 2 n1 X h 2 h 2 ðc Þ ðc Þ bjþ11 bj 1 d2x uðxi ; tnj Þ r2 1 þ dx ½uðxi ; t nþ1 Þ uðxi ; t n Þ ¼ r1 d2x uðxi ; tnþ1 Þ þ r1 dx uðxi ; tnþ1 Þ 12 12 j¼0 ! ! 2 2 n1 h 2 X ðc2 Þ s h 2 ðc2 Þ b bj d 1þ d ½f ðxi ; t nþ1 Þ uðxi ; tnj Þ þ r2 1 þ 12 x j¼0 jþ1 2 12 x þ f ðxi ; t n Þ þ R;
ð22Þ
where
r1 ¼
asc1 Cð1 þ c1 Þ 2
R¼
1þ
bsc2 ; r2 ¼ ; Cð1 þ c2 Þ !
h 2 d ½R1 þ R2 þ R3 ¼ 12 x
ð23Þ 2
1þ
!
h 2 d ðR1 þ R2 Þ; 12 x
ð24Þ
ðcÞ
and bj satisfies (5). Then omitting the small terms in (22) and discretizing the initial and boundary conditions, we can get the following implicit compact finite difference scheme (ICFDS) for the problem (1)–(3):
! ! 2 2 n1 X h 2 nþ1 h 2 ðc1 Þ ðc1 Þ n 2 nþ1 2 nj 1þ bjþ1 bj d d U nþ1 dx U i r2 1 þ U i U i ¼ r1 dx U i þ r1 i 12 x 12 x j¼0 ! ! 2 2 n1 h 2 X s h 2 ðc Þ ðc Þ r2 1 þ bjþ12 bj 2 U nj þ dx 1þ d ½f ðxi ; t nþ1 Þ þ f ðxi ; tn Þ; i 12 2 12 x j¼0 1 6 i 6 M 1; U n0 ¼ uðnsÞ; U 0i
0 6 n 6 N 1; U nM ¼ wðnsÞ;
¼ gðihÞ;
0 6 n 6 N;
0 6 i 6 M:
ð25Þ ð26Þ ð27Þ
3. l‘-stability of the implicit compact finite difference scheme In this section, we will show the stability of the ICFDS (25)–(27). Firstly, for the convenience of the subsequent proof, we here define the following new inner products: For V 2 Uh, U 2 Uh,
hU; Vi ¼ h
M 1 X
2
Ui V i
i¼1
hU x ; V x i ¼ h
M1 X
1 X h M ðU i Þx ðV i Þx ; h 12 i¼0
ðU i Þx ðV i Þx
i¼0
ð28Þ
2
X 2 h M1 ðd U i Þðd2x V i Þ: h 12 i¼1 x
For the new and usual inner products, we have the following lemma:
ð29Þ
4032
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043
Lemma 3.1. For "U, V 2 Uh, we have
h
M 1 X i¼1
h
! ! 2 h 2 dx U i V i ¼ hU; Vi; 12 ! ! 2 h 2 1þ dx U i d2x V i ¼ hU x ; V x i; 12
1þ
M 1 X i¼1
ð30Þ ð31Þ
2 kUk2 6 hU; Ui 6 kUk2 ; 3 2 kU x k2 6 hU x ; U x i 6 kU x k2 : 3
ð32Þ ð33Þ
For (30) and (31), one can prove them easily using summation formula by parts and (28). For (32) and (33), one can use the inequality (a ± b)2 6 2a2 + 2b2 and refer to [24] to obtain the proof. Lemma 3.2 [33]. If a P 0, b P 0, there’s the following Young’s inequality q
ab 6
ap b þ ; p q
ab 6
ðeaÞp þ p
or
1 q eb ; q
where e > 0, p > 1, q > 1 and 1p þ 1q ¼ 1: Lemma 3.3 [34] Discrete Sobolev’s inequality. Suppose that dependent on e such that
n o unj are mesh functions. Given e > 0, there exists a constant C
kuk1 6 ekun k þ CðeÞkunx k: Theorem 3.1. The implicit compact finite difference scheme (ICFDS) (25)–(27) is unconditionally stable to the initial data in l1norm. e n is the approximation solution of the scheme (25)–(27) arised from the initial error. Then the error Proof. Suppose U i
e n U n satisfies eni ¼ U i i
1þ
! ! 2 2 n1 X h 2 nþ1 h 2 nþ1 ðc1 Þ ðc1 Þ 2 nj þ r b b e r 1 þ e dx ei eni ¼ r1 d2x enþ1 d d 1 2 x i i jþ1 j 12 12 x i j¼0 ! 2 n1 h 2 X ðc Þ ðc Þ r2 1 þ bjþ12 bj 2 enj 1 6 i 6 M 1; 0 6 n 6 N 1; dx i ; 12 j¼0
en0 ¼ 0; enM ¼ 0; 0 6 n 6 N;
ð35Þ
e0i ¼ qi ; 0 6 i 6 M:
ð36Þ
where qi is the initial error. Multiplying (34) by he
h
M1 X i¼1
2
1þ
ð34Þ
!
h 2 nþ1 ei eni enþ1 ¼h d i 12 x
M1 X
and summing up for i from 1 to M 1, we have
M1 n1 X X nþ1 ðc Þ ðc Þ r 1 d2x enþ1 ei þ h r1 bjþ11 bj 1 enþ1 d2x enj i i i
i¼1
h
nþ1 i
i¼1
M 1 X
j¼0
! ! 2 2 M 1 n1 X h 2 nþ1 2 h 2 X ðc Þ ðc Þ nþ1 1þ h r2 1 þ bjþ12 bj 2 enj : ð37Þ dx ei dx i ei 12 12 i¼1 j¼0
r2
i¼1
According to (30) and (35), we can get
henþ1 ; enþ1 i hen ; enþ1 i ¼ r 1 kenþ1 k2 þ r 1 h x
M1 n1 XX i¼0
r2
n1 X j¼0
ðc Þ
ðc1 Þ
bj
ðc2 Þ
bjþ12 bj
ðc Þ
bjþ11
j¼0
nþ1 i: henj i ; ei
enj ðenþ1 Þx r2 henþ1 ; enþ1 i i i x
ð38Þ
4033
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043
P ðcÞ ðcÞ ðcÞ ðcÞ Then applying Lemma 3.2, noticing that n1 j¼0 bj bjþ1 ¼ 1 bn and bn > 0; from the above equality, it follows
henþ1 ; enþ1 i 6
n1 n1 1 n n 1 r1 X r1 X ðc Þ ðc Þ ðc Þ ðc Þ 2 k2 þ bj 1 bjþ11 kenj bj 1 bjþ11 kenþ1 k2 he ; e i þ henþ1 ; enþ1 i r 1 kenþ1 x x Þk þ x 2 2 2 j¼0 2 j¼0
r 2 henþ1 ; enþ1 i þ
6
n1 1 n n 1 r1 X r1 ðc Þ ðc Þ ðc Þ 2 k2 þ bj 1 bjþ11 kenj k2 he ; e i þ henþ1 ; enþ1 i r 1 kenþ1 1 bn 1 kenþ1 x x Þk þ x 2 2 2 j¼0 2
r 2 henþ1 ; enþ1 i þ
6
n1 n1 r2 X r2 X ðc Þ ðc Þ ðc Þ ðc Þ bj 2 bjþ12 henj ; enj i þ bj 2 bjþ12 henþ1 ; enþ1 i 2 j¼0 2 j¼0
n1 r2 X r2 ðc Þ ðc Þ ðc Þ bj 2 bjþ12 henj ; enj i þ 1 bn 2 henþ1 ; enþ1 i 2 j¼0 2
ðc Þ n1 n 1 n n 1 r1 X r1 X r 1 bn 1 nþ1 2 ðc Þ ðc Þ 2 bj 1 kenj bj 1 kexnþ1j k2 he ; e i þ henþ1 ; enþ1 i þ kex k x k 2 2 2 j¼0 2 j¼0 2
þ
ðc Þ n1 n r2 X r2 X r2 bn 2 nþ1 nþ1 ðc Þ ðc Þ bj 2 henj ; enj i bj 2 henþ1j ; enþ1j i he ; e i; 2 j¼0 2 j¼0 2
ð39Þ
which means
henþ1 ; enþ1 i þ r 1
n X
ðc1 Þ
bj
kexnþ1j k2 þ r 2
n X
j¼0
6 hen ; en i þ r 1
ðc2 Þ
bj
ðc Þ
ðc Þ
henþ1j ; enþ1j i þ r1 bn 1 kenþ1 k2 þ r2 bn 2 henþ1 ; enþ1 i x
j¼0
n1 X
ðc1 Þ
bj
2 kenj x k þ r2
j¼0
n1 X
ðc2 Þ
bj
henj ; enj i:
ð40Þ
j¼0
Defining the new energy norm
ken k2E ¼ hen ; en i þ r 1
n1 X
ðc1 Þ
bj
2 kenj x k þ r2
j¼0
n1 X
ðc2 Þ
bj
henj ; enj i;
ð41Þ
j¼0
and applying Lemma 3.1 (32), we have
2 nþ1 2 ke k 6 henþ1 ; enþ1 i 6 kenþ1 k2E 6 ken k2E 6 6 ke0 k2E 6 ke0 k2 ; 3
ð42Þ
that is
kenþ1 k 6 Cke0 k:
ð43Þ hd2x nþ1 i
e
Then multiplying (34) by
and summing up for i from 1 to M 1, noticing (31) and (35), we can get
henþ1 ; enþ1 i ¼ henx ; enþ1 i r 1 kd2x enþ1 k2 r1 h x x x
M1 n1 XX
r2
ðc Þ
ðc1 Þ
ðbjþ11 bj
i¼1
2 nþ1 Þd2x enj r2 henþ1 ; enþ1 i x x i dx ei
j¼0
n1 X ðc Þ ðc Þ nþ1 bjþ12 bj 2 henj i: x ; ex
ð44Þ
j¼0
Similarly as the derivation as (40), using Young’s inequality (Lemma 3.2) and noticing that ðcÞ bn > 0; we get
henþ1 ; enþ1 i þ r1 x x
n X
ðc1 Þ
bj
kd2x enþ1j k2 þ r 2
j¼0
6 henx ; enx i þ r 1
n1 X
n X
ðc2 Þ
bj
ðc Þ
Pn1 j¼0
ðcÞ ðcÞ ðcÞ bj bjþ1 ¼ 1 bn and
ðc Þ
henþ1j ; exnþ1j i þ r1 bn 1 kd2x enþ1 k2 þ r2 bn 2 henþ1 ; enþ1 i x x x
j¼0 ðc1 Þ
bj
kd2x enj k2 þ r 2
j¼0
n1 X
ðc2 Þ
bj
nj henj x ; ex i:
ð45Þ
j¼0
Defining the new energy norm
kenþ1 k2E ¼ henþ1 ; enþ1 i þ r1 x x x
n X
ðc1 Þ
bj
j¼0
and applying Lemma 3.1 (33), we have
kd2x enþ1j k2 þ r 2
n X j¼0
ðc2 Þ
bj
hexnþ1j ; enþ1j i; x
ð46Þ
4034
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043
2 nþ1 2 ; enþ1 i 6 kenþ1 k2E 6 kenx k2E 6 6 ke0x k2E 6 ke0x k2 ; ke k 6 henþ1 x x x 3 x
ð47Þ
which means
kenþ1 k 6 Cke0x k: x
ð48Þ
Then combining (43) and (48) and Lemma 3.3, we can obtain
kek1 6 Cðke0 k þ ke0x kÞ: This completes the proof. h
4. l‘ convergence of the implicit compact finite difference scheme In this section, we shall discuss the convergence of the ICFDS (25)–(27). Lemma 4.1 [35]. For "U 2 Uh, we have
L kUk 6 pffiffiffi kU x k: 6
n Theorem 4.1. Let uðx; tÞ 2 C 7;3 x;t be the solution of the problem (1)–(3), then for ns 6 T, the solution U i j0 6 i 6 M; 0 6 n 6 N of 4 the compact scheme (25)–(27) is convergent to the solution u(x, t) with the convergence order O(s + h ) in the l1-norm. Proof. Let eni ¼ uni U ni ; 0 6 i 6 M; 0 6 n 6 N. From (22) and (25), we can obtain the following error equations:
1þ
! ! 2 2 n1 X h 2 nþ1 h 2 nþ1 ðc1 Þ ðc1 Þ 2 nj e þ r b b e r 1 þ dx ei eni ¼ r 1 d2x enþ1 d d 1 2 x i i jþ1 j 12 12 x i j¼0 ! ! 2 2 n1 h 2 X h 2 ðc2 Þ ðc2 Þ nj r2 1 þ b bj d d ðR1 þ R2 Þ; ei þ 1 þ 12 x j¼0 jþ1 12 x
1 6 i 6 M 1;
0 6 n 6 N 1;
ð49Þ
0 6 n 6 N;
ð50Þ
en0 ¼ 0;
enM ¼ 0;
e0i ¼ 0;
0 6 i 6 M:
ð51Þ
Similarly as the derivation as (40), multiplying (49) by (30) and Lemma 3.2, we have
henþ1 ; enþ1 i þ r1
n X
ðc1 Þ
bj
kexnþ1j k2 þ r2
j¼0
6 hen ; en i þ r1
n1 X
n X
ðc2 Þ
bj
nþ1 hei
and summing up for i from 1 to M 1, then applying Lemma 3.1
henþ1j ; enþ1j i
j¼0 ðc1 Þ
bj
2 kenj x k þ r2
j¼0
n1 X
ðc2 Þ
bj
ðc Þ
ðc Þ
henj ; enj i r 1 bn 1 kenþ1 k2 r2 bn 2 henþ1 ; enþ1 i þ 2hR1 þ R2 ; enþ1 i: x
ð52Þ
j¼0
According to Lemma 4.1 and the left side of Lemma 3.1 (32), we have ðc Þ
ðc Þ
r 1 bn 1 kenþ1 k2 6 x
6r 1 bn 1
kenþ1 k2 ; L2 ðc Þ 2r2 bn 2 nþ1 2 ðc Þ r 2 bn 2 henþ1 ; enþ1 i 6 ke k : 3
ð53Þ ð54Þ
Substituting the above two inequalities into (52) and defining the new energy form
ken k2E ¼ hen ; en i þ r 1
n1 X
ðc1 Þ
bj
2 kenj x k þ r2
j¼0
n1 X
ðc2 Þ
henj ; enj i;
ð55Þ
2r2 bn 2 kenþ1 k2 þ 2hR1 þ R2 ; enþ1 i: 3
ð56Þ
bj
j¼0
we have ðc Þ
kenþ1 k2E 6 ken k2E
6r 1 bn 1 2
L
ðc Þ
kenþ1 k2
4035
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043
According to Lemma 3.2 and the right side of (32), we obtain
hR1 þ R2 ; enþ1 i ¼ hR1 ; enþ1 i þ hR2 ; enþ1 i 6
ðc Þ
12r1 bn 1
hR1 ; R1 i þ
3r 1 bn 1 L2
he
nþ1
nþ1
;e
ðc Þ
L2
6
!
ðc Þ
L2
ðc Þ 12r 1 bn 1
kR1 k2 þ
3r1 bn 1
kenþ1 k2 þ
2
L
ðc Þ
3
r 2 bn 2 nþ1 nþ1 i þ hR2 ; R2 i þ he ; e i ðc2 Þ 3 4r 2 bn
!
ðc Þ
3 ðc Þ 4r2 bn 2
kR2 k2 þ
r2 bn 2 kenþ1 k2 : 3
ð57Þ
Substituting (57) into (56), then using (32), (16), (18), (23) and (51), we can get
2 nþ1 2 L2 3 ðc Þ 4 ðc Þ kR1 k2 þ kR2 k2 6 ken k2E þ Cbn 1 sc1 ðs þ h Þ2 þ Cbn 2 s2þc2 6 ke k 6 ken k2E þ ðc Þ ðc Þ 3 2r 2 bn 2 6r 1 bn 1 n X
6 ke0 k2E þ C
n X
ðc1 Þ c1
s ðs þ h4 Þ2 þ C
bj
j¼0
6C
n X
s
6 ke0 k2 þ C
j¼0
ðc1 Þ c1
s ðs þ h4 Þ2 þ C
bj
ðc2 Þ 2þc2
bj
j¼0
n X
ðc1 Þ c1
s ðs þ h4 Þ2 þ C
bj
n X
j¼0
ðc2 Þ 2þc2
s
bj
n X
ðc2 Þ 2þc2
bj
s
j¼0
ð58Þ
:
j¼0
From Lemma 2.1 (4), we obtain 4
kenþ1 k2 6 Cðs þ h Þ2 :
ð59Þ
hd2x enþ1 i
Then multiplying (49) by we can get
and summing up for i from 1 to M 1, then applying Lemma 3.1 (31) and noticing (50),
henþ1 ; enþ1 i ¼ henx ; enþ1 i r1 kd2x enþ1 k2 r 1 h x x x
M 1 X n1 X i¼1
r2
n1 X
ðc Þ
ðc1 Þ
ðbjþ11 bj
2 nþ1 Þd2x enj r 2 henþ1 ; enþ1 i x x i dx ei
j¼0
ðc Þ ðc Þ nþ1 bjþ12 bj 2 henj i hR1 þ R2 ; d2x enþ1 i: x ; ex
ð60Þ
j¼0
Similarly as the derivation as (40), using Lemma 3.2, and noting that
henþ1 ; enþ1 i þ r1 x x
n X
ðc1 Þ
bj
kd2x enjþ1 k2 þ r 2
j¼0
6 henx ; enx i þ r 1
n X
ðc2 Þ
bj
Pn1 ðcÞ ðcÞ ðcÞ ðcÞ j¼0 bj bjþ1 ¼ 1 bn and bn > 0; we get
ðc Þ
ðc Þ
kenjþ1 k2 þ r 1 bn 1 kd2x enþ1 k2 þ r2 bn 2 henþ1 ; enþ1 i x x x
j¼0
n1 X
ðc1 Þ
bj
kd2x enj k2 þ r2
j¼0
n1 X
ðc2 Þ
bj
2 2 nþ1 kenj i: x k 2hR1 þ R2 ; dx e
ð61Þ
j¼0
Now we estimate the last term of (61). From (28) and (18), summation formula by parts, (29) and the smoothness of the u(x, t), we have
hR2 ; d2x enþ1 i ¼ h
M1 X
2
ðR2 Þi d2x enþ1 i
i¼1
2
M1 X X X 2 h M1 h M1 ððR2 Þi Þx ðd2x enþ1 Þx ¼ h ððR2 Þi Þx enþ1 þ dx ððR2 Þi Þx d2x enþ1 h h i i i x 12 i¼0 12 i¼0 i¼1
¼ hðR2 Þx ; enþ1 i: x So according to Lemma 3.2, (16), (18), (23), we obtain
hR1 þ R2 ; d2x enþ1 i ¼ hR1 ; d2x enþ1 i þ hðR2 Þx ; enþ1 i 6 hR1 ; R1 ihd2x enþ1 ; d2x enþ1 i þ henþ1 ; enþ1 ihðR2 Þx ; ðR2 Þx i x x x ðc Þ
1
6
ðc
2r 1 bn 1
kR1 k2 þ Þ
ðc Þ c1
6 Cbn 1
ðc Þ
r1 bn 1 1 r 2 bn 2 nþ1 nþ1 kðR2 Þx k2 þ kd2x enþ1 k2 þ hex ; ex i ðc2 Þ 2 2 2r 2 bn ðc Þ
s ðs þ h4 Þ2 þ
ðc Þ
r 1 bn 1 r2 bn 2 nþ1 nþ1 ðc Þ kd2x enþ1 k2 þ Cbn 2 sc2 s2 þ hex ; ex i: 2 2
ð62Þ
Then substituting the above inequality into (61), we have
henþ1 ; enþ1 i þ r1 x x
n X j¼0
ðc1 Þ
bj
kd2x enjþ1 k2 þ r 2
n X
ðc2 Þ
bj
kenjþ1 k2 6 henx ; enx i þ r 1 x
j¼0
n1 X
ðc1 Þ
kd2x enj k2 þ r 2
4 2
ðc Þ Cbn 2 c2
bj
j¼0
þ
ðc Þ Cbn 1 c1 ð
n1 X
ðc2 Þ
bj
2 kenj x k
j¼0
s sþh Þ þ
2
s s:
ð63Þ
4036
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043
Defining the new energy norm
kenþ1 k2E ¼ henþ1 ; enþ1 i þ r1 x x x
n X
ðc1 Þ
bj
kd2x enþ1j k2 þ r2
n X
j¼0
ðc2 Þ
bj
kexnþ1j k2 ;
ð64Þ
j¼0
then according to Lemma 3.1 (33) and (51), for ns 6 T, we get
2 nþ1 2 ðc Þ 4 ðc Þ k2E 6 kenx k2E þ Cbn 1 sc1 ðs þ h Þ2 þ Cbn 1 sc1 s2 6 ke k 6 kenþ1 x 3 x n n n n X X X X ðc Þ 4 ðc Þ ðc Þ 4 ðc Þ bj 1 sc1 ðs þ h Þ2 þ C b0 1 sc1 s2 6 ke0x k2 þ C bj 1 sc1 ðs þ h Þ2 þ C b0 1 sc1 s2 6 ke0x k2E þ C j¼0
6C
n X
j¼0
ðc1 Þ c1
s ðs þ h4 Þ2 þ C
bj
j¼0
n X
j¼0
j¼0
ðc Þ c1
b0 1
s s2 :
ð65Þ
j¼0
From Lemma 2.1 (4), we have 4
kenþ1 k2 6 cðs þ h Þ2 : x
ð66Þ
Finally, according to (59) and (66) and Lemma 3.3, we get 4
kenþ1 k1 6 Cðs þ h Þ: This completes the proof. h
5. Another improved implicit compact scheme
Lemma 5.1 [29]. Suppose u(t) is sufficiently smooth, then
uðnÞ ¼
ðt jþ1 nÞuðtj Þ þ ðn t j Þuðtjþ1 Þ
s
þ Oðs2 Þ:
ð67Þ
From the first equality of (8), we have
a Cðc1 Þ
uðxi ; t nþ1 Þ uðxi ; t n Þ ¼
Z
Luðxi ; nÞ ðtnþ1 nÞ
0
b Cðc2 Þ
þ
t nþ1
Z
tn 0
1c1
uðxi ; nÞ ðt n nÞ
1c2
a Cðc1 Þ
dn
dn þ
Z
Z
tn
0
Luðxi ; nÞ 1c1
ðt n nÞ
dn
b Cðc2 Þ
Z
t nþ1
f ðxi ; nÞdn ¼ I1 I2 I3 þ I4 þ
tn
t nþ1
Z
uðxi ; nÞ ðt nþ1 nÞ1c2
0
dn
t nþ1
f ðxi ; nÞdn;
ð68Þ
tn
where
I1 ¼ I2 ¼ I3 ¼ I4 ¼
a Cðc1 Þ a Cðc1 Þ b Cðc2 Þ b Cðc2 Þ
Z
t nþ1
Luðxi ; nÞ ðtnþ1 nÞ1c1
0
Z
tn
Z
Luðxi ; nÞ ðt n nÞ1c1
0
t nþ1
Z
dn;
uðxi ; nÞ ðtnþ1 nÞ1c2
0 tn
0
uðxi ; nÞ ðt n nÞ1c2
dn;
dn;
dn:
ð69Þ
From Lemma 5.1 and (7), we can get
I1 ¼
n a X Cðc1 Þ j¼0
¼ r1
n X j¼0
"
Z
t jþ1
ðt jþ1 nÞLuðxi ; tj Þ þ ðn t j ÞLuðxi ; tjþ1 Þ
s
tj
xðj c1 Þ
d2x 2
ðc1 Þ
uðxi ; tnj Þ þ tj 2
h 1 þ 12 dx
d2x 2
h 1 þ 12 d2x
1 ðtnþ1 nÞ1c1 #
uðxi ; tnjþ1 Þ þ R1 :
dn
ð70Þ
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043
4037
Similarly, we have
I2 ¼ r 1
n1 X
"
xðj c1 Þ
d2x
ðc1 Þ
uðxi ; tnj1 Þ þ tj h2 2 1 þ 12 dx
#
d2x
uðxi ; t nj Þ þ R2 ; h2 2 1 þ 12 dx j¼0 n h i X xðj c2 Þ uðxi ; tnj Þ þ tðj c2 Þ uðxi ; tnjþ1 Þ þ R3 ; I3 ¼ r 2
ð71Þ ð72Þ
j¼0
I4 ¼ r 2
n1 h X
i
xðj c2 Þ uðxi ; tnj1 Þ þ tjðc2 Þ uðxi ; tnj Þ þ R4 ;
ð73Þ
i 1 h ðj þ 1Þ1þc ðjÞ1þc ; 1þc h i 1þc c ðj þ 1Þ1þc j j ;
ð74Þ
j¼0
where
xjðcÞ ¼ ðj þ 1Þc
tjðcÞ ¼
1 1þc
4
2
R1 ; R2 ¼ Oðs þ h Þ; R3 ; R4 ¼ Oðs2 Þ; and r1, r2 satisfy (23). Substituting (70)–(73) and (19) into (68), omitting the small terms and rearranging the equality, we can obtain the following improved implicit compact finite difference scheme (IICFDS)
! ( ) 2 n1 h i X h 2 nþ1 ðc1 Þ 2 n ðc1 Þ 2 nþ1 ðc1 Þ ðc1 Þ ðc1 Þ ðc1 Þ n 2 nj1 2 nj U i U i ¼ r1 x0 dx U i þ t0 dx U i þ 1þ xjþ1 xj dx U i þ tjþ1 tj dx U i d 12 x j¼0 ) !( 2 n1 h i X h 2 ðc2 Þ ðc2 Þ ðc2 Þ ðc2 Þ nj1 nj xð0c2 Þ U ni þ tð0c2 Þ U nþ1 þ ð x x ÞU þ t t dx U r2 1 þ i jþ1 j jþ1 j i i 12 j¼0 ! 2 s h 2 nþ1 fi þ fin ; 1 6 i 6 M 1; 0 6 n 6 N 1; þ 1þ d ð75Þ 2 12 x U n0 ¼ uðnsÞ; U 0i ¼ gðihÞ;
U nM ¼ wðnsÞ;
0 6 n 6 N;
ð76Þ
0 6 i 6 M:
ð77Þ
6. Numerical experiments In this section, numerical results are presented to demonstrate the accuracy and effectiveness of the two compact schemes (25)–(27) and (75)–(77). Example 1 [26]. Consider the following initial and boundary problem of the fractional cable equation: 2 @uðx; tÞ 1c @ uðx; tÞ 1c ¼ 0 Dt 1 0 Dt 2 uðx; tÞ þ f ðx; tÞ; @t @x2 uð0; tÞ ¼ 0; uðL; tÞ ¼ 0; t 2 ð0; T; uðx; 0Þ ¼ 0; x 2 ½0; 1;
x 2 ð0; 1Þ;
t 2 ð0; T;
ð78Þ ð79Þ ð80Þ
where f ðx; tÞ ¼ 2ðt þ p2 t1þc1 =Cð2 þ c1 Þ þ t 1þc2 =Cð2 þ c2 ÞÞ sin px: The exact solution of the system (78)–(80) is u(x, t) = t2 sinpx, which can be obtained by evaluating the fractional derivative formula. In [26], the authors gave the following implicit numerical method (INM) for the problem (1)–(3)
U nþ1 ¼ U ni þ r 1 d2x U nþ1 þ r1 i i
n1 X
ðc Þ
ðc1 Þ
bjþ11 bj
r 2 U nþ1 d2x U nj i i
j¼0 n1 X s ðc Þ ðc Þ r2 bjþ12 bj 2 U nj þ ½f ðxi ; t nþ1 Þ þ f ðxi ; t n Þ 1 6 i 6 M 1; i 2 j¼0
U n0 ¼ uðnsÞ; U 0i ¼ gðihÞ;
U nM ¼ wðnsÞ; 0 6 i 6 M:
0 6 n 6 N;
0 6 n 6 N 1;
ð81Þ
ð82Þ ð83Þ
4038
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043 ðcÞ
where r1, r2 satisfy (23) and bj satisfies (5), and proved the stability and convergency of the scheme in l2-norm both theoretically and numerically. For the l2-norm depends on the mesh step, the l1-norm can describe the approximation degree of the numerical solutions better. By using the similar method and Lemma 3.3 of this paper, the scheme INM (81)–(83) can be proved to be stable and convergent in l1-norm theoretically. The details are not discussed any more. Here we only give the l1-convergency from numerical experiments, as well as the scheme IINM of [26]. Firstly, we prove the n l1-convergence order of the o four schemes INM, IINM of [26] and ICFDS (25)–(27) and IICFDS (75)– (77) of this paper. Let U nj j1 6 j 6 M 1; 0 6 n 6 N be the solution of the numerical scheme and denote
e1 ðs; hÞ ¼ max juNj U Nj j: 16j6M1
Suppose that p
e1 ðs; hÞ ¼ Oðsq þ h Þ: Then we use
log 2
ðe1 ð2p=q s; 2hÞÞ p: e1 ðs; hÞ
to compute the convergence order. For the case c1 = c2, we take c1 = c2 = 0.5. For the case c1 – c2, we take c1 = 0.9, c2 = 0.3 and c1 = 0.2, c2 = 0.8 respectively. Tables 1–3 provide some computational results with the schemes INM, IINM [26] and our two schemes ICFDS amd IICFDS at T = 1, which shows that the convergence order of our schemes ICFDS amd IICFDS is O(s + h4) and O(s2 + h4), respectively, while the order of INM and IINM [26] is only O(s + h2) and O(s2 + h2). Then we give some accuracy comparisons. Tables 4 and 5 provides some computational results of the maximum errors of INM [26] versus ICFDS, and IINM [26] versus IICFDS at T = 1 with the same h and s, respectively. From these data we can see that our two schemes are actually more accurate than the two schemes proposed in [26]. To test and verify the effectiveness of our two new schemes for short time, Tables 6 and 7 present some numerical maximum error comparisons for 0 6 t 6 1 with h = 1/10, s = 1/20,000 (Table 6) and h = 1/ 100, s = 1/1000 (Table 7) respectively. Similarly to test and verify the effectiveness of our two new schemes for long time, Tables 8 and 9 present some numerical maximum error comparisons for 1 6 t 6 10 with h = 1/10, s = 1/1000 (Table 8) and h = 1/100, s = 1/1000 (Table 9) respectively. From these tables, we see that our two compact schemes can reproduce the exact analytical solutions and also work better than the two schemes proposed in [26] for either short time or long time, especially IICFDS.
Table 1 Convergence order of scheme INM [26] and ICFDS for c1 = c2 = 0.5 at T = 1. h, s
e1(s, h)
e1 ð2p=q s;2hÞ e1 ðs;hÞ
s;2hÞÞ log 2 ðe1eð21 ðs;hÞ
Scheme INM [26] s = h = 1/8 h = 1/16, s = 1/32 h = 1/32, s = 1/128 h = 1/64, s = 1/516
– 7.846988e02 2.280452e02 6.074373e03 1.563170e03
– – 3.441 3.754 3.886
– – 1.783 1.909 1.958
Scheme ICFDS
– 8.786068e02 6.705252e03 4.358651e04 2.749260e05
– – 13.103 15.384 15.854
– – 3.720 3.943 3.987
s = h = 1/8 h = 1/16, s = 1/128 h = 1/32, s = 1/2048 h = 1/64, s = 1/32,768
p=q
Table 2 Convergence order of scheme IINM [26] and IICFDS for c1 = c2 = 0.5 at T = 1. h, s
e1(s, h)
e1 ð2p=q s;2hÞ e1 ðs;hÞ
s;2hÞ log 2 e1eð21 ðs;hÞ
Scheme IINM [26] s = h = 1/8 h = s = 1/16 h = s = 1/32 h = s = 1/64
– 1.01789e02 2.532596e03 6.324977e04 1.583664e04
– – 4.002 4.001 3.994
– – 2.008 2.000 1.998
Scheme IICFDS
– 7.939284e05 4.938659e05 3.083094e07 1.926444e08
– – 16.076 16.019 16.004
– – 4.007 4.002 4.000
s = h = 1/8 h = 1/16, s = 1/32 h = 1/32, s = 1/128 h = 1/64, s = 1/512
p=q
4039
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043 Table 3 Convergence order of scheme ICFDS for c1 = 0.9, c2 = 0.3 and c1 = 0.2, c2 = 0.8 at T = 1. h, s
e1(s, h)
e1 ð16s;2hÞ e1 ðs;hÞ
s;2hÞ log 2 e1eð16 1 ðs;hÞ
c1 = 0.9, c2 = 0.3 s = h = 1/8 h = 1/16, s = 1/128 h = 1/32, s = 1/2048 h = 1/64, s = 1/32,768
– 9.98366e02 6.890772e03 6.074373e04 2.743784e05
– – 14.488 15.777 15.929
– – 3.857 3.979 3.994
c1 = 0.2, c2 = 0.8 s = h = 1/7 h = 1/14, s = 1/112 h = 1/28, s = 1/1792 h = 1/56, s = 1/28,672
– 5.893423e02 5.676444e03 4.196787e04 2.853184e05
– – 10.382 13.526 14.709
– – 3.376 3.758 3.879
Table 4 Comparison of some maximum errors for INM [26] and ICFDS with c1 = c2 = 0.5 at T = 1.
s = 1/100,000
INM [26]
ICFDS
h = 1/4 h = 1/8 h = 1/16 h = 1/32
4.1759e2 1.0285e2 2.3523e3 6.3069e4
1.2845e3 7.0353e5 4.0881e6 8.7184e6
Table 5 Comparison of some maximum errors for IINM [26] and IICFDS with c1 = c2 = 0.5 at T = 1.
s = 1/5000
IINM [26]
IICFDS
h = 1/10 h = 1/20 h = 1/40 h = 1/80
4.5612e3 1.2907e3 3.6420e4 9.6578e5
3.2443e5 2.0217e5 1.2627e7 7.8954e9
Table 6 Comparison of some computational maximum errors with c1 = c2 = 0.5. t
INM [26]
ICFDS
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
4.7796e5 2.1914e4 5.2286e4 9.6227e4 1.5392e3 2.2552e3 3.1110e3 4.1015e3 5.2452e3 6.5246e3
3.4436e6 6.8604e6 9.8036e6 1.2163e5 1.3893e5 1.4974e5 1.5394e5 1.5141e5 1.4211e5 1.2596e5
Table 7 Comparison of some computational maximum errors with c1 = c2 = 0.5. t
IINM [26]
IICFDS
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
5.1147e7 2.2494e6 5.2883e6 9.6540e6 1.5360e5 2.2415e5 3.0824e5 4.0590e5 5.1716e5 6.4203e5
2.5351e11 1.1179e10 2.6388e10 4.8168e10 7.6760e10 1.1218e9 1.5446e9 2.0368e9 2.5978e9 3.2287e9
4040
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043
Table 8 Comparison of some computational maximum errors for h = 1/10, s = 1/1000. t
1 2 3 4 5 6 7 8 9 10
c1 = 0.25, c2 = 0.75
c1 = 0.75, c2 = 0.25
INM [26]
ICFDS
INM [26]
ICFDS
6.1580e3 2.6030e2 5.9273e2 1.0056e1 1.6460e1 2.3615e1 3.1198e1 4.1588e1 5.2367e1 6.4318e1
7.5722e4 1.5180e3 2.3338e3 2.8977e3 3.5075e3 4.0622e3 4.5615e3 5.0057e3 5.3949e3 5.7294e3
5.2246e3 2.5152e2 6.5865e2 1.1176e1 1.7880e1 2.6179e1 3.6079e1 4.7584e1 6.0698e1 7.5424e1
8.6978e4 1.7457e3 2.5562e3 3.2968e3 3.9960e3 4.5630e3 5.0874e3 5.5390e3 5.9174e3 6.2227e3
Table 9 Comparison of some computational maximum errors for h = 1/100, s = 1/1000. t
c1 = c2 = 0.25
c1 = c2 = 0.75
IINM [26]
IICFDS
IINM [26]
IICFDS
1 2 3 4 5 6 7 8 9 10
6.3581e5 2.4896e4 5.4610e4 9.4570e4 1.4396e3 2.0194e3 2.6784e3 3.4091e3 4.2055e3 5.0607e3
3.2287e9 1.3397e8 3.0644e8 5.5161e8 8.6717e8 1.2551e7 1.7159e7 2.2501e7 2.8580e7 3.5195e7
6.3596e5 2.7024e4 6.2248e4 1.2083e3 1.7655e3 2.5567e3 3.4943e3 4.5785e3 5.8093e3 7.1866e3
3.1523e9 1.3412e8 3.0920e8 5.5710e8 8.7801e8 1.2720e7 1.7392e7 2.2797e7 2.8935e7 3.5805e7
Finally we draw some error curves. Take M = 32, N = 32, 64, 128, 256, 512, respectively. Fig. 1 plots the curves of the max errors of the ICFNS solutions for c1 = c2 = 0.25, c1 = c2 = 0.75, c1 = 0.8, c2 = 0.3 and c1 = 0.2, c2 = 0.7, respectively, at T = 1. Then
0.02
0.03
N=32 N=64 N=128 N=256 N=512
0.018 0.016 0.014
0.02
0.012 0.01
0.015
0.008
0.01
0.006 0.004
0.005
0.002 0
0
0.2
0.4
0.6
0.8
0.03
1
0
0
0.2
0.4
0.6
0.8
0.018
N=32 N=64 N=128 N=256 N=512
0.025
1
N=32 N=64 N=128 N=256 N=512
0.016 0.014
0.02
0.012 0.01
0.015
0.008
0.01
0.006 0.004
0.005 0
N=32 N=64 N=128 N=256 N=512
0.025
0.002 0
0.2
0.4
0.6
0.8
1
0
0
0.2
0.4
0.6
0.8
1
Fig. 1. The curves of the max errors of the ICFDS at T = 1. Top left: c1 = c2 = 0.25. Top right: c1 = c2 = 0.75. Bottom left: c1 = 0.8, c2 = 0.3. Bottom right: c1 = 0.2, c2 = 0.7.
4041
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043
x 10 −8
x 10 −7 4
2 M=M=32
3 2
1
1
0.5
0 1
0.8
0.6
0.4 x
0.2
0 0
0.4
0.2
0.6
0.8
M=64,N=128
1.5
0 1
1
0.8
0.6
0.4 x
t
0.2
0 0
0.6
0.4
0.2
0.8
1
t
x 10 −8
x 10 −7 2
3 M=N=32 2
M=64,N=128
1.5 1
1 0 1
0.5
0.8
0.6
0.4 x
0.2
0 0
0.4
0.2
0.6
0.8
1
0 1
0.8
0.6
0.4 x
t
0.2
0 0
0.2
0.4
0.6
0.8
1
t
Fig. 2. The max error-surface of the IICFDS for c1 = c2 = 0.3 (up) and c1 = c2 = 0.8 (down) at T = 1.
take M = N = 32, and M = 64, N = 128. Fig. 2 plots max error curved surface of IICFNS solutions for c1 = c2 = 0.3 (up) and c1 = c2 = 0.8 (down), respectively. From these figures, we can see that the max errors become smaller as the time step s becomes smaller, which implies the two methods are actually efficient for the computation. Example 2. Consider the following fractional cable equation [26,32]: 2 @uðx; tÞ 1c @ uðx; tÞ 1c ¼ 0 Dt 1 0:50 Dt 2 uðx; tÞ; @t @x2 uð0; tÞ ¼ 0; uðL; tÞ ¼ 0; t 2 ð0; 0:1;
uðx; 0Þ ¼ 10dðx 5Þ;
x 2 ð0; 10Þ;
t 2 ð0; 0:1;
ð84Þ ð85Þ
x 2 ½0; 10;
ð86Þ
approximation lg(u(x,t))
where d(x) is the Dirac delta function. This problem does not exists an analytical solution. In [32], Liao and Sun did not consider the case of c1 – c2. Here take M = 32, N = 160. Fig. 3 plots the changes of u(x, t) for fixed c1 = 0.5 and different c2 = 0.2, 0.4, 0.6, 0.8, 1.0 at t = 0.1, respectively. And Fig. 4 plots the changes of u(x, t) for fixed c2 = 0.5 and different c1 = 0.2, 0.4, 0.6, 0.8, 1.0 at t = 0.1. These figures show the impacts of c1 and c2 on the anomalous behavior. We can see that there is a crossover between more and less anomalous behaviors, which is coincidence with the results of [26].
0.4 0.2 0 −0.2 −0.4 −0.6 4
4.5
5 x
5.5
6
Fig. 3. u(x, t) for fixed c1 = 0.5 and different c2 = 0.2, 0.4, 0.6, 0.8, 1.0 at t = 0.1.
4042
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043
approximation log2(u(x,t))
5
0
−5
−10
−15
−20
3
3.5
4
4.5
5 x
5.5
6
6.5
7
7.5
Fig. 4. u(x, t) for fixed c2 = 0.5 and different c1 = 0.2, 0.4, 0.6, 0.8, 1.0 at t = 0.1.
7. Conclusion In this paper, we build two implicit compact finite difference schemes for the fractional cable equation. For the l1-norm can describe the approximation degree of the numerical solutions, we prove the l1-stability and convergence of the first scheme by defining the new inner product. The proof techniques can also be applied to the other linear or nonlinear problems. Numerical experiments confirm the l1 convergence order of the two schemes and the numerical comparisons show that our two new schemes are more accurate than the two previous schemes. And the two compact schemes are effective and also work better for either short time or long time than the previous. Among the four schemes, the second compact scheme IICFDS can enhance the accuracy greatly and work best. Acknowledgment The authors are grateful to the referees for their constructive and valuable suggestions on this manuscript. References [1] R. Metzler, J. Klafter, The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics, J. Phys. A 37 (2004) R161–R208. [2] I.M. Sokolov, J. Klafter, From diffusion to anomalous diffusion: A century after Einstein’s Brownian motion, Chaos 15 (2005) 026103–026103-7. [3] J.A. Tenreiro Machado, Analysis and design of fractional-order digital control systems, Syst. Anal. Model. Simul. 27 (1997) 107–122. [4] D. Baleanu, O. Defterli, O.P. Agrawal, A central difference numerical scheme for fractional optimal control problems, J. Vip. Contr. 15 (2009) 583–597. [5] F. Mainardi, M. Raberto, R. Gorenflo, E. Scalas, Fractional calculus and continuous-time finance I: the waiting-time distribution, Physica A. 287 (2000) 468–481. [6] R. Magin, X. Feng, D. Baleanu, Solving the fractional order Bloch equation, Conc. Mag. Reson. Part A 34A (1) (2009) 16–23. [7] T.A.M. Langlands, B.I. Henry, The accuracy and stability of an implicit solution method for the fractional diffusion equation, J. Comput. Phys. 205 (2005) 719–736. [8] S.B. Yuste, Weighted average finite difference methods for fractional diffusion equations, J. Comput. Phys. 216 (2006) 264–274. [9] D.A. Murio, Implicit finite difference approximation for time fractional diffusion equation, Comput. Math. Appl. 56 (2008) 1138–1145. [10] Y. Zhang, A finite difference method for fractional partial differential equation, Appl. Math. Comput. 215 (2009) 524–529. [11] M.M. Meerschaert, H.P. Scheffler, C. Tadjeran, Finite difference methods for two-dimensional fractional dispersion equation, J. Comput. Phys. 211 (2006) 249–261. [12] C. Tadjeran, M.M. Meerschaert, H.P. Scheffler, A second-order accurate numerical approximation for the fractional diffusion equation, J. Comput. Phys. 213 (2006) 205–213. [13] C. Tadjeran, M.M. Meerschaert, A second-order accurate numerical method for the two-dimensional fractional diffusion equation, J. Comput. Phys. 220 (2007) 813–823. [14] C. Chen, F. Liu, I. Turner, V. Anh, A Fourier method for the fractional diffusion equation describing sub-diffusion, J. Comput. Phys. 227 (2007) 886–897. [15] C. Chen, F. Liu, K. Burrage, Finite difference methods and a Fourier analysis for the fractional reaction–subdiffusion equation, Appl. Math. Comput. 198 (2008) 754–769. [16] T. Langlands, B. Henry, The accuracy and stability of an implicit solution method for the fractional diffusion equation, J. Comput. Phys. 205 (2005) 719– 736. [17] P. Zhuang, F. Liu, V. Anh, I. Turner, New solution and analytical techniques of the implicit numerical method for the anomalous subdiffusion equation, SIAM J. Numer. Anal. 46 (2008) 1079–1095. [18] 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. [19] 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. [20] C. Chen, F. Liu, I. Turner, V. Anh, Numerical schemes and multivariate extrapolation of a two-dimensional anomalous sub-diffusion equation, Numer. Algor. 54 (2010) 1–21. [21] Z.Z. Sun, X. Wu, A fully difference scheme for a diffusion-wave system, Appl. Numer. Math. 56 (2006) 193–209. [22] M. Cui, Compact finite difference method for the fractional diffusion equation, J. Comput. Phys. 228 (2009) 7792–7804.
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043
4043
[23] R. Du, W.R. Cao, Z.Z. Sun, A compact difference scheme for the fractional diffusion-wave equation, Appl. Math. Model. 34 (2010) 2998–3007. [24] Guanghua, Gao, Z.Z. Sun, A compact difference scheme for the fractional sub-diffusion equation, J. Comput. Phys. 230 (2011) 586–595. [25] C. 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 (2010) 1740–1760. [26] T.A.M. Langlans, B. Henry,S. Wearne, Solution of a fractional cable equation: finite case,. [27] T.A.M. Langlans, B. Henry,S. Wearne, Solution of a fractional cable equation: finite case, .B.I. [28] Henry, T.A.M. Langlands, S.L. Wearne, Fractional cable models for spiny neuronal dendrites, Phys. Rev. Lett. 100 (2008) 128103–128106. [29] F. Liu, Q. Yang, I. Turner, Two new implicit numerical methods for the fractional cable equation, J. Comput. Nonlinear Dyn. 6 (2011) 0110091–0110097. [30] Z.Z. Sun, D. Zhao, On the l1 convergence of a difference scheme for the coupled schrödinger equations, Comput. Math. Appl. 59 (2010) 3286–3330. [31] W. Cao, Z.Z. Sun, Maximum norm error estimates of the Cranck–Niclson scheme for solving a linear moving boundary problem, J. Comput. Appl. Math. 234 (2010) 2578–2586. [32] H. Liao, Z.Z. Sun, Maximum norm error estimates of efficient difference schemes for second-order wave equations, J. Comput. Appl. Math. 235 (2011) 2217–2233. [33] T. Wang, B. Guo, A robust semi-explicit difference scheme for the Kuramoto–Tsuzuki equation, J. Comput. Appl. Math. 233 (2009) 878–888. [34] Y. Zhou, Application of Discrete Functional Analysis of the Finite Difference Method, International Academic Publishers, Beijing, 1990. [35] Z.Z. Sun, Numerical Methods of the Partial Differential Equations, Science Press, Beijing, 2005.
Contents lists available at SciVerse ScienceDirect
Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm
Implicit compact difference schemes for the fractional cable equation q Xiuling Hu a,b, Luming Zhang a,⇑ a b
College of Science, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China School of Mathematical Sciences, Xuzhou Normal University, Xuzhou 221116, China
a r t i c l e
i n f o
Article history: Received 28 April 2011 Received in revised form 29 October 2011 Accepted 2 November 2011 Available online 17 November 2011 Keywords: Fractional cable equation Compact finite difference scheme Stability Convergence
a b s t r a c t In this paper, we propose two implicit compact difference schemes for the fractional cable equation. The first scheme is proved to be stable and convergent in l1-norm with the convergence order O(s + h4) by the energy method, where new inner products defined in this paper gives great convenience for the theoretical analysis. Numerical experiments are presented to demonstrate the accuracy and effectiveness of the two compact schemes. The computational results show that the two new schemes proposed in this paper are more accurate and effective than the previous. 2011 Elsevier Inc. All rights reserved.
1. Introduction Fractional differential and integral equations have been applied to numerous fields of science, such as physics [1,2], engineering [3,4], finance [5], medicine [6], and so on. Many scholars have been attracted to study them using kinds of ways in recent years. For most of the systems that cannot be solved analytically, some numerical methods including the finite difference methods [7–13] have been considered by researchers. Chen et al. [14,15] established the finite difference schemes based on Grünwald–Letnikov approximation and analyzed the schemes by Fourier method for the fractional sub-diffusion and reaction–subdiffusion equation. Langlands and Henry [16] constructed an unconditional implicit scheme based on L1 approximation for the fractional subdiffusion but without global convergence analysis. Zhuang et al. [17] proposed a new implicit numerical method (INM) by integrating method for an anomalous subdiffusion equation and presented two techniques for improving the convergence order of the scheme and explored the stability and convergence by the energy method. Then they proposed an implicit difference scheme by using the same method for nonlinear fractional reaction–subdiffusion equation [18]. Liu et al. [19] studied a modified anomalous subdiffusion equation with a nonlinear source term and constructed a new implicit difference scheme and also analyzed the stability and convergence in l2-norm by a new energy method. Chen et al. [20] also investigated two-dimensional anomalous subdiffusion. Sun and Wu [21] explored a fractional diffusion-wave system and proposed a fully discrete scheme by the method of order reduction and proved the solvability, stability and convergence by the energy method. Because compact finite difference scheme is developed to approximate the second-order derivatives and keep the desirable trigonal nature, the compact method has been received great attention. Cui [22] considered a spacial fourth order finite difference scheme for the one-dimensional fractional diffusion equation and analyzed the stability by the Fourier method. Du et al. [23] presented a compact scheme for the time Caputo’s sense fractional diffusion-wave system on the basis of q
This work is partially supported by the National Science Foundation of China under grant 11102179.
⇑ Corresponding author.
E-mail addresses: [email protected] (X. Hu), [email protected] (L. Zhang). 0307-904X/$ - see front matter 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2011.11.027
4028
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043
[21] and prove the stability and convergence by energy method. Guanghua et al. [24] provided a compact finite difference scheme by applying L1 discretization for the fractional derivative and fourth order approximation for the second order space derivative for the fractional subdiffusion equation after a transformation of the original problem and discussed its unique solvability, stability and convergence by introducing a new inner product. Chen et al. [25] proposed a scheme with order O(s + h4) for a variable-order anomalous subdiffusion equation and analyzed the convergence, stability and solvability by the Fourier method. Then an improved scheme with the order O(s2 + h4) are explored. The cable equation is one of the most fundamental equations for modeling neuronal dynamics. Sun and Zhao [30] and Chaos and Sun [31] found a fractional cable equation involving two fractional temporal Riemann–Liouville derivatives by assuming small ionic concentration gradients and simplifying the fractional Nernst–Planck equations and investigated solutions of a fractional cable equation for finite case and infinite case, respectively. The fractional cable equation (FCE) was introduced by Liao and Sun [32] from the fractional Nernst–Planck equations for modeling the anomalous electrodiffusion of ions in spiny neuronal dendrites. Langlans et al. [26] proposed two implicit numerical methods for the fractional cable equation with l2-convergence order O(s + h2) and O(s2 + h2). Because the l2-norm depends on the spatial mesh step, it cannot portray the degree of the closeness of the numerical solution to the exact solution better than the l1-norm. Recently, the maximum problem has also been paid to much attention [27–29]. This paper studies the fractional cable equation and proposes two implicit compact schemes and prove the l1-stability and convergency of the first scheme by the energy method. The following of this paper is organized as follows. In Section 2, an implicit compact finite difference scheme is derived. In Section 3, the l1-stability of the compact scheme is proved by the energy method where new inner products are introduced. Then the l1-convergency of the compact scheme is proved in Section 4. In Section 5, an improved implicit compact finite difference scheme is derived and two illustrative numerical examples are carried out to demonstrate the accuracy and effectiveness of the two schemes in Section 6. In the end, some conclusions are made in Section 7. 2. An implicit compact finite difference scheme Consider the following fractional cable equation:
! @u @ 2 uðx; tÞ 1c1 1c ¼ 0 Dt a b0 Dt 2 uðx; tÞ þ f ðx; tÞ; @t @x2 uð0; tÞ ¼ uðtÞ;
uðL; tÞ ¼ wðtÞ;
uðx; 0Þ ¼ gðxÞ;
x 2 ½0; L;
1c uðx; tÞ 0 Dt
¼
1
@
CðcÞ @t
ð1Þ
t 2 ð0; T;
1c uðx; tÞ 0 Dt
t
0
t 2 ð0; T;
ð2Þ ð3Þ
where 0 < c1, c2 < 1, a, b > 0, and
Z
x 2 ð0; LÞ;
uðx; nÞ ðt nÞ1c
is the Rieman–Liouville fractional partial derivative of order 1 c defined by
dn;
ð4Þ
Firstly, we denote xj = jh, tn = ns, Xh = {xjj0 6 j 6 M}, Xs = {tnj0 6 n 6 N}, and Xsh ¼ Xh Xs ; where h = L/M, s = T/N are the uniform spatial and temporal mesh sizes respectively, and M, N are two positive integers. In the rest of the paper, C denotes a general positive constant, which may have different values in different occurrences. Suppose U ¼ fU nj j0 6 j 6 M; 0 6 n 6 Ng be a grid function on Xsh : Introduce the following notations:
unj ¼ uðxj ; tn Þ; U nj uðxj ; t n Þ; 1 n 1 U U nj ; d2x U nj ¼ 2 U njþ1 2U nj þ U nj1 U nj ¼ x h jþ1 h Let Uh = {UjU = (U0, U1, . . . , UM), U0 = UM = 0}. For V 2 Uh, U 2 Uh, we firstly introduce the following usual discrete l2 inner product (, ), discrete l2-norm kk and discrete l1-norm kk1:
ðU n ; V n Þ ¼ h
M 1 X
U nj V nj ;
kU n k2 ¼ ðU n ; U n Þ;
kU nx k2 ¼ h
j¼1
M 1 X j¼0
U nj
2 x
;
kU n k1 ¼ max jU nj j: 16j6M1
Lemma 2.1 [17]. Let ðcÞ
c
bk ¼ ðk þ 1Þc k ; then
ðcÞ bk
k ¼ 0; 1; . . . ;
ð5Þ
satisfies
ðcÞ
ðcÞ
(1) b0 ¼ 1; bk > 0; k ¼ 0; 1; . . . ðcÞ
ðcÞ
(2) bk > bkþ1 ; k ¼ 0; 1; . . . (3) There exists a positive constant C > 0, such that Pk ðcÞ c c c c (4) j¼0 bj s ¼ ðk þ 1Þ s 6 T :
s 6 CbkðcÞ sc ; k ¼ 1; 2; 3 . . .
4029
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043
Assume that problem (1)–(3) has a smooth solution uðx; tÞ 2 C 7;3 x;t ð½0; L ½0; TÞ; and denote
Lu ¼
@ 2 uðx; tÞ : @x2
ð6Þ
From Taylor expansion, we can easily know [22] 4
Luðxi ; tn Þ ¼
@ 2 uðxi ; tn Þ d2x h @ 4 uðxi ; t n Þ 6 ¼ uðxi ; tn Þ þ þ Oðh Þ: 2 2 2 1 @x @x4 240 1 þ h dx
ð7Þ
12
Now, we consider the compact scheme for the problem (1)–(3). Integrating both sides of Eq. (1) from tn to tn+1, we have
uðxi ; tnþ1 Þ uðxi ; t n Þ ¼
a Cðc1 Þ
Z
t nþ1
Luðxi ; nÞ 1c1
ðt nþ1 nÞ
0
dn
a Cðc1 Þ
Z
tn
0
Luðxi ; nÞ ðtn nÞ
1c1
dn
b Cðc2 Þ
Z 0
tnþ1
uðxi ; nÞ ðt nþ1 nÞ1c2
dn
Z tnþ1 Z tn b uðxi ; nÞ dn þ f ðxi ; nÞdn Cðc2 Þ 0 ðt n nÞ1c2 tn Z s Z tn Z s a Luðxi ; nÞ a Luðxi ; n þ sÞ Luðxi ; nÞ b uðxi ; nÞ ¼ dn þ dn dn 1c1 1c2 Cðc1 Þ 0 ðtnþ1 nÞ1c1 Cðc1 Þ 0 C ð c Þ ðt n nÞ 0 ðt nþ1 nÞ 2 Z tnþ1 Z tn b uðxi ; n þ sÞ uðxi ; nÞ dn þ f ðxi ; nÞdn ¼ I1 I2 þ I3 ; ð8Þ Cðc2 Þ 0 ðt n nÞ1c2 tn þ
where
I1 ¼ I2 ¼ I3 ¼
a Cðc1 Þ b Cðc2 Þ Z tnþ1
Z s 0
Z s 0
Luðxi ; nÞ 1c1
ðt nþ1 nÞ
uðxi ; nÞ 1c2
ðt nþ1 nÞ
dn þ dn þ
a Cðc1 Þ b Cðc2 Þ
Z
tn
Luðxi ; n þ sÞ Luðxi ; nÞ ðt n nÞ1c1
0
Z
tn
uðxi ; n þ sÞ uðxi ; nÞ
0
ðt n nÞ1c2
dn;
dn;
ð9Þ
f ðxi ; nÞdn:
tn
According to Taylor expansion and (7), for I1 ; we have
I1 ¼
a Cðc1 Þ
Z s 0
Luðxi ; sÞ ðtnþ1 nÞ
1c1
dn þ R11 þ
n1 a X Cðc1 Þ j¼0
Z
t jþ1
Luðxi ; n þ sÞ Luðxi ; nÞ
tj
ðtn nÞ1c1
dn
¼
n1 asc1 asc1 X ðc Þ ðc Þ b 1 Luðxi ; t jþ2 Þ Luðxi ; tjþ1 Þ þ R12 bn 1 Luðxi ; t 1 Þ þ R11 þ Cð1 þ c1 Þ Cð1 þ c1 Þ j¼0 nj1
¼
asc1 d2x asc1 ðc Þ uðxi ; t 1 Þ þ R13 þ R11 þ bn 1 2 2 h Cð1 þ c1 Þ Cð1 þ c1 Þ 1 þ 12 dx
"
n1 X
ðc1 Þ bnj1
"
#
d2x
d2x
#
uðxi ; t jþ2 Þ uðxi ; tjþ1 Þ þ R14 þ R12 h2 2 h2 2 1 þ 12 dx 1 þ 12 dx " # " # n1 asc1 d2x asc1 X d2x d2x ðc Þ ðc1 Þ ¼ uðx ; t Þ þ b uðx ; t Þ uðx ; t Þ þ R1 ; bn 1 i 1 i jþ2 i jþ1 h2 2 h2 2 Cð1 þ c1 Þ Cð1 þ c1 Þ j¼0 nj1 1 þ h2 d2x 1 þ 12 dx 1 þ 12 dx 12
j¼0
ðcÞ
where bj
ð10Þ
satisfies (5) and
R1 ¼ R11 þ R12 þ
n1 asc1 asc1 X ðc Þ ðc Þ b 1 R : bn 1 R13 þ Cð1 þ c1 Þ Cð1 þ c1 Þ j¼0 nj1 14
ð11Þ
For R11, from Taylor expansion, we have
" # a Z s a Z s 1 1 @ 3 uðxi ; n1 Þ jR11 j ¼ ½Luðxi ; nÞ Luðxi ; sÞdn ¼ ðn sÞ dn 1 c 1 c 2 Cðc1 Þ 0 ðt nþ1 nÞ 1 Cðc1 Þ 0 ðt nþ1 nÞ 1 @x @t Z 3 as @ uðxi ; n1 Þ s 1 ðc Þ 6 dn 6 C s1þc1 bn 1 ; 0 < n < n1 < s: Cðc1 Þ @x2 @t 0 ðt nþ1 nÞ1c1
ð12Þ
4030
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043
For R12, from Taylor expansion, we have
R12 ¼
n1 a X Cðc1 Þ j¼0
Z
t jþ1
½Luðxi ; n þ sÞ Luðxi ; nÞ ½Luðxi ; t jþ2 Þ Luðxi ; t jþ1 Þ ðt n nÞ1c1
tj
dn
" # @ 3 uðxi ; n2 Þ @ 3 uðxi ; n3 Þ ðt jþ2 n sÞ ðt jþ1 nÞ ðt n nÞc1 1 dn @x2 @t @x2 @t tj " # n1 Z t jþ1 a X @ 3 uðxi ; n2 Þ @ 3 uðxi ; n3 Þ ðtjþ1 nÞðtn nÞc1 1 dn ¼ @x2 @t @x2 @t Cðc1 Þ j¼0 tj
n1 a X ¼ Cðc1 Þ j¼0
¼
n1 a X Cðc1 Þ j¼0
Z
t jþ1
Z
t jþ1 tj
@ 4 uðxi ; n4 Þ ðn2 n3 Þðtjþ1 nÞðtn nÞc1 1 dn; @x2 @t 2
where
n þ s < n2 < tjþ2 ;
tj < n < t jþ1 ;
n < n3 < t jþ1 ;
n3 < n4 < n2 :
Thus,
tjþ1 < n < tj ;
tjþ1 < n3 < n:
So we have
0 ¼ tj þ s t jþ1 < n2 n3 < tjþ2 n < tjþ2 tj < 2s;
0 < t jþ1 n < t jþ1 tj ¼ s:
Then we can get
a X X n1 Z t jþ1 n1 Z tjþ1 @ 4 uðxi ; n4 Þ c1 1 c1 1 2 6 C jR12 j ¼ ðn n Þðt nÞðt nÞ dn s ðt nÞ dn n n jþ1 2 3 2 2 Cðc1 Þ j¼0 tj j¼0 tj @x @t Z t n ¼ C s2 ðt n nÞc1 1 dn 6 C s2 : 0
And from Lemma 2.1 (3), we have ðc Þ
jR12 j 6 C s1þc1 bn 1 :
ð13Þ
From (7), we have 4
jR13 j 6 Ch ;
ð14Þ
and
d2x d2x jR14 j ¼ Luðxi ; tjþ2 Þ Luðxi ; t jþ1 Þ uðx ; t Þ þ uðx ; t Þ i jþ2 i jþ1 2 2 2 2 h h 1 þ 12 dx 1 þ 12 dx 4 h4 @ 4 uðx ; t Þ 4 6 6 h @ uðxi ; t jþ1 Þ 6 @ uðxi ; t jþ2 Þ 6 @ uðxi ; t jþ1 Þ i jþ2 ¼ þ Ch Ch þ 4 4 6 6 240 @x @x @x @x 240 4 sh @ 5 uðx ; gÞ 7 6 @ uðxi ; g1 Þ 4 i ¼ þ C sh þ 6 C sh ; tjþ1 < g; g1 < tjþ2 : 240 @x4 @t @x6 @t Again from Lemma 2.1 (3), we have ðc Þ 4
jR14 j 6 C sc1 bn 1 h :
ð15Þ
So from (12)–(15) and Lemma 2.1 (1), (2) and (4), we can get
n1 asc1 asc1 X ðc1 Þ ðc1 Þ jR1 j ¼ R11 þ R12 þ bnj1 R14 bn R13 þ Cð1 þ c1 Þ Cð1 þ c1 Þ j¼0 ðc Þ
ðc Þ
ðc Þ 4
ðc Þ 4
ðc Þ
4
6 C s1þc1 bn 1 þ C s1þc1 bn 1 þ C sc1 bn 1 h þ C sc1 bn 1 h 6 C sc1 bn 1 ðs þ h Þ:
ð16Þ
Similarly, for I2 ; we have
I2 ¼
n1 X bsc2 bsc2 ðc Þ ðc Þ b 2 ½uðxi ; t jþ2 Þ uðxi ; t jþ1 Þ þ R2 : bn 2 uðxi ; t 1 Þ þ Cð1 þ c2 Þ Cð1 þ c2 Þ j¼0 nj1
ð17Þ
4031
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043
where ðc Þ 1þc2
s
jR2 j 6 Cbn 2
ð18Þ
:
For I3 , we have
I3 ¼
Z
t nþ1
f ðxi ; nÞdn ¼
tn
s 2
½f ðxi ; tn Þ þ f ðxi ; t nþ1 Þ þ R3 :
ð19Þ
From Taylor expansion and Lemma 2.13, we have ðc Þ 1þc2
jR3 j 6 C s2 6 Cbn 2
s
ð20Þ
:
Thus, we can obtain ðc Þ
uðxi ; tnþ1 Þ uðxi ; t n Þ ¼ r 1 bn 1
(
d2x
uðxi ; t 1 Þ þ r 1 h2 2 1 þ 12 dx ðc Þ
r 2 bn 2 uðxi ; t1 Þ þ
n1 X
n1 X
ðc Þ
"
1 bnj1
j¼0
d2x
d2x uðxi ; t jþ2 Þ uðxi ; tjþ1 Þ 2 2 h h2 2 1 þ 12 dx 1 þ 12 dx )
ðc Þ
2 bnj1 ½uðxi ; t jþ2 Þ uðxi ; t jþ1 Þ
j¼0
#
s
þ ½f ðxi ; t nþ1 Þ þ f ðxi ; t n Þ þ R1 þ R2 þ R3 ; 2 ð21Þ
or
1þ
! ! 2 2 n1 X h 2 h 2 ðc Þ ðc Þ bjþ11 bj 1 d2x uðxi ; tnj Þ r2 1 þ dx ½uðxi ; t nþ1 Þ uðxi ; t n Þ ¼ r1 d2x uðxi ; tnþ1 Þ þ r1 dx uðxi ; tnþ1 Þ 12 12 j¼0 ! ! 2 2 n1 h 2 X ðc2 Þ s h 2 ðc2 Þ b bj d 1þ d ½f ðxi ; t nþ1 Þ uðxi ; tnj Þ þ r2 1 þ 12 x j¼0 jþ1 2 12 x þ f ðxi ; t n Þ þ R;
ð22Þ
where
r1 ¼
asc1 Cð1 þ c1 Þ 2
R¼
1þ
bsc2 ; r2 ¼ ; Cð1 þ c2 Þ !
h 2 d ½R1 þ R2 þ R3 ¼ 12 x
ð23Þ 2
1þ
!
h 2 d ðR1 þ R2 Þ; 12 x
ð24Þ
ðcÞ
and bj satisfies (5). Then omitting the small terms in (22) and discretizing the initial and boundary conditions, we can get the following implicit compact finite difference scheme (ICFDS) for the problem (1)–(3):
! ! 2 2 n1 X h 2 nþ1 h 2 ðc1 Þ ðc1 Þ n 2 nþ1 2 nj 1þ bjþ1 bj d d U nþ1 dx U i r2 1 þ U i U i ¼ r1 dx U i þ r1 i 12 x 12 x j¼0 ! ! 2 2 n1 h 2 X s h 2 ðc Þ ðc Þ r2 1 þ bjþ12 bj 2 U nj þ dx 1þ d ½f ðxi ; t nþ1 Þ þ f ðxi ; tn Þ; i 12 2 12 x j¼0 1 6 i 6 M 1; U n0 ¼ uðnsÞ; U 0i
0 6 n 6 N 1; U nM ¼ wðnsÞ;
¼ gðihÞ;
0 6 n 6 N;
0 6 i 6 M:
ð25Þ ð26Þ ð27Þ
3. l‘-stability of the implicit compact finite difference scheme In this section, we will show the stability of the ICFDS (25)–(27). Firstly, for the convenience of the subsequent proof, we here define the following new inner products: For V 2 Uh, U 2 Uh,
hU; Vi ¼ h
M 1 X
2
Ui V i
i¼1
hU x ; V x i ¼ h
M1 X
1 X h M ðU i Þx ðV i Þx ; h 12 i¼0
ðU i Þx ðV i Þx
i¼0
ð28Þ
2
X 2 h M1 ðd U i Þðd2x V i Þ: h 12 i¼1 x
For the new and usual inner products, we have the following lemma:
ð29Þ
4032
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043
Lemma 3.1. For "U, V 2 Uh, we have
h
M 1 X i¼1
h
! ! 2 h 2 dx U i V i ¼ hU; Vi; 12 ! ! 2 h 2 1þ dx U i d2x V i ¼ hU x ; V x i; 12
1þ
M 1 X i¼1
ð30Þ ð31Þ
2 kUk2 6 hU; Ui 6 kUk2 ; 3 2 kU x k2 6 hU x ; U x i 6 kU x k2 : 3
ð32Þ ð33Þ
For (30) and (31), one can prove them easily using summation formula by parts and (28). For (32) and (33), one can use the inequality (a ± b)2 6 2a2 + 2b2 and refer to [24] to obtain the proof. Lemma 3.2 [33]. If a P 0, b P 0, there’s the following Young’s inequality q
ab 6
ap b þ ; p q
ab 6
ðeaÞp þ p
or
1 q eb ; q
where e > 0, p > 1, q > 1 and 1p þ 1q ¼ 1: Lemma 3.3 [34] Discrete Sobolev’s inequality. Suppose that dependent on e such that
n o unj are mesh functions. Given e > 0, there exists a constant C
kuk1 6 ekun k þ CðeÞkunx k: Theorem 3.1. The implicit compact finite difference scheme (ICFDS) (25)–(27) is unconditionally stable to the initial data in l1norm. e n is the approximation solution of the scheme (25)–(27) arised from the initial error. Then the error Proof. Suppose U i
e n U n satisfies eni ¼ U i i
1þ
! ! 2 2 n1 X h 2 nþ1 h 2 nþ1 ðc1 Þ ðc1 Þ 2 nj þ r b b e r 1 þ e dx ei eni ¼ r1 d2x enþ1 d d 1 2 x i i jþ1 j 12 12 x i j¼0 ! 2 n1 h 2 X ðc Þ ðc Þ r2 1 þ bjþ12 bj 2 enj 1 6 i 6 M 1; 0 6 n 6 N 1; dx i ; 12 j¼0
en0 ¼ 0; enM ¼ 0; 0 6 n 6 N;
ð35Þ
e0i ¼ qi ; 0 6 i 6 M:
ð36Þ
where qi is the initial error. Multiplying (34) by he
h
M1 X i¼1
2
1þ
ð34Þ
!
h 2 nþ1 ei eni enþ1 ¼h d i 12 x
M1 X
and summing up for i from 1 to M 1, we have
M1 n1 X X nþ1 ðc Þ ðc Þ r 1 d2x enþ1 ei þ h r1 bjþ11 bj 1 enþ1 d2x enj i i i
i¼1
h
nþ1 i
i¼1
M 1 X
j¼0
! ! 2 2 M 1 n1 X h 2 nþ1 2 h 2 X ðc Þ ðc Þ nþ1 1þ h r2 1 þ bjþ12 bj 2 enj : ð37Þ dx ei dx i ei 12 12 i¼1 j¼0
r2
i¼1
According to (30) and (35), we can get
henþ1 ; enþ1 i hen ; enþ1 i ¼ r 1 kenþ1 k2 þ r 1 h x
M1 n1 XX i¼0
r2
n1 X j¼0
ðc Þ
ðc1 Þ
bj
ðc2 Þ
bjþ12 bj
ðc Þ
bjþ11
j¼0
nþ1 i: henj i ; ei
enj ðenþ1 Þx r2 henþ1 ; enþ1 i i i x
ð38Þ
4033
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043
P ðcÞ ðcÞ ðcÞ ðcÞ Then applying Lemma 3.2, noticing that n1 j¼0 bj bjþ1 ¼ 1 bn and bn > 0; from the above equality, it follows
henþ1 ; enþ1 i 6
n1 n1 1 n n 1 r1 X r1 X ðc Þ ðc Þ ðc Þ ðc Þ 2 k2 þ bj 1 bjþ11 kenj bj 1 bjþ11 kenþ1 k2 he ; e i þ henþ1 ; enþ1 i r 1 kenþ1 x x Þk þ x 2 2 2 j¼0 2 j¼0
r 2 henþ1 ; enþ1 i þ
6
n1 1 n n 1 r1 X r1 ðc Þ ðc Þ ðc Þ 2 k2 þ bj 1 bjþ11 kenj k2 he ; e i þ henþ1 ; enþ1 i r 1 kenþ1 1 bn 1 kenþ1 x x Þk þ x 2 2 2 j¼0 2
r 2 henþ1 ; enþ1 i þ
6
n1 n1 r2 X r2 X ðc Þ ðc Þ ðc Þ ðc Þ bj 2 bjþ12 henj ; enj i þ bj 2 bjþ12 henþ1 ; enþ1 i 2 j¼0 2 j¼0
n1 r2 X r2 ðc Þ ðc Þ ðc Þ bj 2 bjþ12 henj ; enj i þ 1 bn 2 henþ1 ; enþ1 i 2 j¼0 2
ðc Þ n1 n 1 n n 1 r1 X r1 X r 1 bn 1 nþ1 2 ðc Þ ðc Þ 2 bj 1 kenj bj 1 kexnþ1j k2 he ; e i þ henþ1 ; enþ1 i þ kex k x k 2 2 2 j¼0 2 j¼0 2
þ
ðc Þ n1 n r2 X r2 X r2 bn 2 nþ1 nþ1 ðc Þ ðc Þ bj 2 henj ; enj i bj 2 henþ1j ; enþ1j i he ; e i; 2 j¼0 2 j¼0 2
ð39Þ
which means
henþ1 ; enþ1 i þ r 1
n X
ðc1 Þ
bj
kexnþ1j k2 þ r 2
n X
j¼0
6 hen ; en i þ r 1
ðc2 Þ
bj
ðc Þ
ðc Þ
henþ1j ; enþ1j i þ r1 bn 1 kenþ1 k2 þ r2 bn 2 henþ1 ; enþ1 i x
j¼0
n1 X
ðc1 Þ
bj
2 kenj x k þ r2
j¼0
n1 X
ðc2 Þ
bj
henj ; enj i:
ð40Þ
j¼0
Defining the new energy norm
ken k2E ¼ hen ; en i þ r 1
n1 X
ðc1 Þ
bj
2 kenj x k þ r2
j¼0
n1 X
ðc2 Þ
bj
henj ; enj i;
ð41Þ
j¼0
and applying Lemma 3.1 (32), we have
2 nþ1 2 ke k 6 henþ1 ; enþ1 i 6 kenþ1 k2E 6 ken k2E 6 6 ke0 k2E 6 ke0 k2 ; 3
ð42Þ
that is
kenþ1 k 6 Cke0 k:
ð43Þ hd2x nþ1 i
e
Then multiplying (34) by
and summing up for i from 1 to M 1, noticing (31) and (35), we can get
henþ1 ; enþ1 i ¼ henx ; enþ1 i r 1 kd2x enþ1 k2 r1 h x x x
M1 n1 XX
r2
ðc Þ
ðc1 Þ
ðbjþ11 bj
i¼1
2 nþ1 Þd2x enj r2 henþ1 ; enþ1 i x x i dx ei
j¼0
n1 X ðc Þ ðc Þ nþ1 bjþ12 bj 2 henj i: x ; ex
ð44Þ
j¼0
Similarly as the derivation as (40), using Young’s inequality (Lemma 3.2) and noticing that ðcÞ bn > 0; we get
henþ1 ; enþ1 i þ r1 x x
n X
ðc1 Þ
bj
kd2x enþ1j k2 þ r 2
j¼0
6 henx ; enx i þ r 1
n1 X
n X
ðc2 Þ
bj
ðc Þ
Pn1 j¼0
ðcÞ ðcÞ ðcÞ bj bjþ1 ¼ 1 bn and
ðc Þ
henþ1j ; exnþ1j i þ r1 bn 1 kd2x enþ1 k2 þ r2 bn 2 henþ1 ; enþ1 i x x x
j¼0 ðc1 Þ
bj
kd2x enj k2 þ r 2
j¼0
n1 X
ðc2 Þ
bj
nj henj x ; ex i:
ð45Þ
j¼0
Defining the new energy norm
kenþ1 k2E ¼ henþ1 ; enþ1 i þ r1 x x x
n X
ðc1 Þ
bj
j¼0
and applying Lemma 3.1 (33), we have
kd2x enþ1j k2 þ r 2
n X j¼0
ðc2 Þ
bj
hexnþ1j ; enþ1j i; x
ð46Þ
4034
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043
2 nþ1 2 ; enþ1 i 6 kenþ1 k2E 6 kenx k2E 6 6 ke0x k2E 6 ke0x k2 ; ke k 6 henþ1 x x x 3 x
ð47Þ
which means
kenþ1 k 6 Cke0x k: x
ð48Þ
Then combining (43) and (48) and Lemma 3.3, we can obtain
kek1 6 Cðke0 k þ ke0x kÞ: This completes the proof. h
4. l‘ convergence of the implicit compact finite difference scheme In this section, we shall discuss the convergence of the ICFDS (25)–(27). Lemma 4.1 [35]. For "U 2 Uh, we have
L kUk 6 pffiffiffi kU x k: 6
n Theorem 4.1. Let uðx; tÞ 2 C 7;3 x;t be the solution of the problem (1)–(3), then for ns 6 T, the solution U i j0 6 i 6 M; 0 6 n 6 N of 4 the compact scheme (25)–(27) is convergent to the solution u(x, t) with the convergence order O(s + h ) in the l1-norm. Proof. Let eni ¼ uni U ni ; 0 6 i 6 M; 0 6 n 6 N. From (22) and (25), we can obtain the following error equations:
1þ
! ! 2 2 n1 X h 2 nþ1 h 2 nþ1 ðc1 Þ ðc1 Þ 2 nj e þ r b b e r 1 þ dx ei eni ¼ r 1 d2x enþ1 d d 1 2 x i i jþ1 j 12 12 x i j¼0 ! ! 2 2 n1 h 2 X h 2 ðc2 Þ ðc2 Þ nj r2 1 þ b bj d d ðR1 þ R2 Þ; ei þ 1 þ 12 x j¼0 jþ1 12 x
1 6 i 6 M 1;
0 6 n 6 N 1;
ð49Þ
0 6 n 6 N;
ð50Þ
en0 ¼ 0;
enM ¼ 0;
e0i ¼ 0;
0 6 i 6 M:
ð51Þ
Similarly as the derivation as (40), multiplying (49) by (30) and Lemma 3.2, we have
henþ1 ; enþ1 i þ r1
n X
ðc1 Þ
bj
kexnþ1j k2 þ r2
j¼0
6 hen ; en i þ r1
n1 X
n X
ðc2 Þ
bj
nþ1 hei
and summing up for i from 1 to M 1, then applying Lemma 3.1
henþ1j ; enþ1j i
j¼0 ðc1 Þ
bj
2 kenj x k þ r2
j¼0
n1 X
ðc2 Þ
bj
ðc Þ
ðc Þ
henj ; enj i r 1 bn 1 kenþ1 k2 r2 bn 2 henþ1 ; enþ1 i þ 2hR1 þ R2 ; enþ1 i: x
ð52Þ
j¼0
According to Lemma 4.1 and the left side of Lemma 3.1 (32), we have ðc Þ
ðc Þ
r 1 bn 1 kenþ1 k2 6 x
6r 1 bn 1
kenþ1 k2 ; L2 ðc Þ 2r2 bn 2 nþ1 2 ðc Þ r 2 bn 2 henþ1 ; enþ1 i 6 ke k : 3
ð53Þ ð54Þ
Substituting the above two inequalities into (52) and defining the new energy form
ken k2E ¼ hen ; en i þ r 1
n1 X
ðc1 Þ
bj
2 kenj x k þ r2
j¼0
n1 X
ðc2 Þ
henj ; enj i;
ð55Þ
2r2 bn 2 kenþ1 k2 þ 2hR1 þ R2 ; enþ1 i: 3
ð56Þ
bj
j¼0
we have ðc Þ
kenþ1 k2E 6 ken k2E
6r 1 bn 1 2
L
ðc Þ
kenþ1 k2
4035
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043
According to Lemma 3.2 and the right side of (32), we obtain
hR1 þ R2 ; enþ1 i ¼ hR1 ; enþ1 i þ hR2 ; enþ1 i 6
ðc Þ
12r1 bn 1
hR1 ; R1 i þ
3r 1 bn 1 L2
he
nþ1
nþ1
;e
ðc Þ
L2
6
!
ðc Þ
L2
ðc Þ 12r 1 bn 1
kR1 k2 þ
3r1 bn 1
kenþ1 k2 þ
2
L
ðc Þ
3
r 2 bn 2 nþ1 nþ1 i þ hR2 ; R2 i þ he ; e i ðc2 Þ 3 4r 2 bn
!
ðc Þ
3 ðc Þ 4r2 bn 2
kR2 k2 þ
r2 bn 2 kenþ1 k2 : 3
ð57Þ
Substituting (57) into (56), then using (32), (16), (18), (23) and (51), we can get
2 nþ1 2 L2 3 ðc Þ 4 ðc Þ kR1 k2 þ kR2 k2 6 ken k2E þ Cbn 1 sc1 ðs þ h Þ2 þ Cbn 2 s2þc2 6 ke k 6 ken k2E þ ðc Þ ðc Þ 3 2r 2 bn 2 6r 1 bn 1 n X
6 ke0 k2E þ C
n X
ðc1 Þ c1
s ðs þ h4 Þ2 þ C
bj
j¼0
6C
n X
s
6 ke0 k2 þ C
j¼0
ðc1 Þ c1
s ðs þ h4 Þ2 þ C
bj
ðc2 Þ 2þc2
bj
j¼0
n X
ðc1 Þ c1
s ðs þ h4 Þ2 þ C
bj
n X
j¼0
ðc2 Þ 2þc2
s
bj
n X
ðc2 Þ 2þc2
bj
s
j¼0
ð58Þ
:
j¼0
From Lemma 2.1 (4), we obtain 4
kenþ1 k2 6 Cðs þ h Þ2 :
ð59Þ
hd2x enþ1 i
Then multiplying (49) by we can get
and summing up for i from 1 to M 1, then applying Lemma 3.1 (31) and noticing (50),
henþ1 ; enþ1 i ¼ henx ; enþ1 i r1 kd2x enþ1 k2 r 1 h x x x
M 1 X n1 X i¼1
r2
n1 X
ðc Þ
ðc1 Þ
ðbjþ11 bj
2 nþ1 Þd2x enj r 2 henþ1 ; enþ1 i x x i dx ei
j¼0
ðc Þ ðc Þ nþ1 bjþ12 bj 2 henj i hR1 þ R2 ; d2x enþ1 i: x ; ex
ð60Þ
j¼0
Similarly as the derivation as (40), using Lemma 3.2, and noting that
henþ1 ; enþ1 i þ r1 x x
n X
ðc1 Þ
bj
kd2x enjþ1 k2 þ r 2
j¼0
6 henx ; enx i þ r 1
n X
ðc2 Þ
bj
Pn1 ðcÞ ðcÞ ðcÞ ðcÞ j¼0 bj bjþ1 ¼ 1 bn and bn > 0; we get
ðc Þ
ðc Þ
kenjþ1 k2 þ r 1 bn 1 kd2x enþ1 k2 þ r2 bn 2 henþ1 ; enþ1 i x x x
j¼0
n1 X
ðc1 Þ
bj
kd2x enj k2 þ r2
j¼0
n1 X
ðc2 Þ
bj
2 2 nþ1 kenj i: x k 2hR1 þ R2 ; dx e
ð61Þ
j¼0
Now we estimate the last term of (61). From (28) and (18), summation formula by parts, (29) and the smoothness of the u(x, t), we have
hR2 ; d2x enþ1 i ¼ h
M1 X
2
ðR2 Þi d2x enþ1 i
i¼1
2
M1 X X X 2 h M1 h M1 ððR2 Þi Þx ðd2x enþ1 Þx ¼ h ððR2 Þi Þx enþ1 þ dx ððR2 Þi Þx d2x enþ1 h h i i i x 12 i¼0 12 i¼0 i¼1
¼ hðR2 Þx ; enþ1 i: x So according to Lemma 3.2, (16), (18), (23), we obtain
hR1 þ R2 ; d2x enþ1 i ¼ hR1 ; d2x enþ1 i þ hðR2 Þx ; enþ1 i 6 hR1 ; R1 ihd2x enþ1 ; d2x enþ1 i þ henþ1 ; enþ1 ihðR2 Þx ; ðR2 Þx i x x x ðc Þ
1
6
ðc
2r 1 bn 1
kR1 k2 þ Þ
ðc Þ c1
6 Cbn 1
ðc Þ
r1 bn 1 1 r 2 bn 2 nþ1 nþ1 kðR2 Þx k2 þ kd2x enþ1 k2 þ hex ; ex i ðc2 Þ 2 2 2r 2 bn ðc Þ
s ðs þ h4 Þ2 þ
ðc Þ
r 1 bn 1 r2 bn 2 nþ1 nþ1 ðc Þ kd2x enþ1 k2 þ Cbn 2 sc2 s2 þ hex ; ex i: 2 2
ð62Þ
Then substituting the above inequality into (61), we have
henþ1 ; enþ1 i þ r1 x x
n X j¼0
ðc1 Þ
bj
kd2x enjþ1 k2 þ r 2
n X
ðc2 Þ
bj
kenjþ1 k2 6 henx ; enx i þ r 1 x
j¼0
n1 X
ðc1 Þ
kd2x enj k2 þ r 2
4 2
ðc Þ Cbn 2 c2
bj
j¼0
þ
ðc Þ Cbn 1 c1 ð
n1 X
ðc2 Þ
bj
2 kenj x k
j¼0
s sþh Þ þ
2
s s:
ð63Þ
4036
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043
Defining the new energy norm
kenþ1 k2E ¼ henþ1 ; enþ1 i þ r1 x x x
n X
ðc1 Þ
bj
kd2x enþ1j k2 þ r2
n X
j¼0
ðc2 Þ
bj
kexnþ1j k2 ;
ð64Þ
j¼0
then according to Lemma 3.1 (33) and (51), for ns 6 T, we get
2 nþ1 2 ðc Þ 4 ðc Þ k2E 6 kenx k2E þ Cbn 1 sc1 ðs þ h Þ2 þ Cbn 1 sc1 s2 6 ke k 6 kenþ1 x 3 x n n n n X X X X ðc Þ 4 ðc Þ ðc Þ 4 ðc Þ bj 1 sc1 ðs þ h Þ2 þ C b0 1 sc1 s2 6 ke0x k2 þ C bj 1 sc1 ðs þ h Þ2 þ C b0 1 sc1 s2 6 ke0x k2E þ C j¼0
6C
n X
j¼0
ðc1 Þ c1
s ðs þ h4 Þ2 þ C
bj
j¼0
n X
j¼0
j¼0
ðc Þ c1
b0 1
s s2 :
ð65Þ
j¼0
From Lemma 2.1 (4), we have 4
kenþ1 k2 6 cðs þ h Þ2 : x
ð66Þ
Finally, according to (59) and (66) and Lemma 3.3, we get 4
kenþ1 k1 6 Cðs þ h Þ: This completes the proof. h
5. Another improved implicit compact scheme
Lemma 5.1 [29]. Suppose u(t) is sufficiently smooth, then
uðnÞ ¼
ðt jþ1 nÞuðtj Þ þ ðn t j Þuðtjþ1 Þ
s
þ Oðs2 Þ:
ð67Þ
From the first equality of (8), we have
a Cðc1 Þ
uðxi ; t nþ1 Þ uðxi ; t n Þ ¼
Z
Luðxi ; nÞ ðtnþ1 nÞ
0
b Cðc2 Þ
þ
t nþ1
Z
tn 0
1c1
uðxi ; nÞ ðt n nÞ
1c2
a Cðc1 Þ
dn
dn þ
Z
Z
tn
0
Luðxi ; nÞ 1c1
ðt n nÞ
dn
b Cðc2 Þ
Z
t nþ1
f ðxi ; nÞdn ¼ I1 I2 I3 þ I4 þ
tn
t nþ1
Z
uðxi ; nÞ ðt nþ1 nÞ1c2
0
dn
t nþ1
f ðxi ; nÞdn;
ð68Þ
tn
where
I1 ¼ I2 ¼ I3 ¼ I4 ¼
a Cðc1 Þ a Cðc1 Þ b Cðc2 Þ b Cðc2 Þ
Z
t nþ1
Luðxi ; nÞ ðtnþ1 nÞ1c1
0
Z
tn
Z
Luðxi ; nÞ ðt n nÞ1c1
0
t nþ1
Z
dn;
uðxi ; nÞ ðtnþ1 nÞ1c2
0 tn
0
uðxi ; nÞ ðt n nÞ1c2
dn;
dn;
dn:
ð69Þ
From Lemma 5.1 and (7), we can get
I1 ¼
n a X Cðc1 Þ j¼0
¼ r1
n X j¼0
"
Z
t jþ1
ðt jþ1 nÞLuðxi ; tj Þ þ ðn t j ÞLuðxi ; tjþ1 Þ
s
tj
xðj c1 Þ
d2x 2
ðc1 Þ
uðxi ; tnj Þ þ tj 2
h 1 þ 12 dx
d2x 2
h 1 þ 12 d2x
1 ðtnþ1 nÞ1c1 #
uðxi ; tnjþ1 Þ þ R1 :
dn
ð70Þ
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043
4037
Similarly, we have
I2 ¼ r 1
n1 X
"
xðj c1 Þ
d2x
ðc1 Þ
uðxi ; tnj1 Þ þ tj h2 2 1 þ 12 dx
#
d2x
uðxi ; t nj Þ þ R2 ; h2 2 1 þ 12 dx j¼0 n h i X xðj c2 Þ uðxi ; tnj Þ þ tðj c2 Þ uðxi ; tnjþ1 Þ þ R3 ; I3 ¼ r 2
ð71Þ ð72Þ
j¼0
I4 ¼ r 2
n1 h X
i
xðj c2 Þ uðxi ; tnj1 Þ þ tjðc2 Þ uðxi ; tnj Þ þ R4 ;
ð73Þ
i 1 h ðj þ 1Þ1þc ðjÞ1þc ; 1þc h i 1þc c ðj þ 1Þ1þc j j ;
ð74Þ
j¼0
where
xjðcÞ ¼ ðj þ 1Þc
tjðcÞ ¼
1 1þc
4
2
R1 ; R2 ¼ Oðs þ h Þ; R3 ; R4 ¼ Oðs2 Þ; and r1, r2 satisfy (23). Substituting (70)–(73) and (19) into (68), omitting the small terms and rearranging the equality, we can obtain the following improved implicit compact finite difference scheme (IICFDS)
! ( ) 2 n1 h i X h 2 nþ1 ðc1 Þ 2 n ðc1 Þ 2 nþ1 ðc1 Þ ðc1 Þ ðc1 Þ ðc1 Þ n 2 nj1 2 nj U i U i ¼ r1 x0 dx U i þ t0 dx U i þ 1þ xjþ1 xj dx U i þ tjþ1 tj dx U i d 12 x j¼0 ) !( 2 n1 h i X h 2 ðc2 Þ ðc2 Þ ðc2 Þ ðc2 Þ nj1 nj xð0c2 Þ U ni þ tð0c2 Þ U nþ1 þ ð x x ÞU þ t t dx U r2 1 þ i jþ1 j jþ1 j i i 12 j¼0 ! 2 s h 2 nþ1 fi þ fin ; 1 6 i 6 M 1; 0 6 n 6 N 1; þ 1þ d ð75Þ 2 12 x U n0 ¼ uðnsÞ; U 0i ¼ gðihÞ;
U nM ¼ wðnsÞ;
0 6 n 6 N;
ð76Þ
0 6 i 6 M:
ð77Þ
6. Numerical experiments In this section, numerical results are presented to demonstrate the accuracy and effectiveness of the two compact schemes (25)–(27) and (75)–(77). Example 1 [26]. Consider the following initial and boundary problem of the fractional cable equation: 2 @uðx; tÞ 1c @ uðx; tÞ 1c ¼ 0 Dt 1 0 Dt 2 uðx; tÞ þ f ðx; tÞ; @t @x2 uð0; tÞ ¼ 0; uðL; tÞ ¼ 0; t 2 ð0; T; uðx; 0Þ ¼ 0; x 2 ½0; 1;
x 2 ð0; 1Þ;
t 2 ð0; T;
ð78Þ ð79Þ ð80Þ
where f ðx; tÞ ¼ 2ðt þ p2 t1þc1 =Cð2 þ c1 Þ þ t 1þc2 =Cð2 þ c2 ÞÞ sin px: The exact solution of the system (78)–(80) is u(x, t) = t2 sinpx, which can be obtained by evaluating the fractional derivative formula. In [26], the authors gave the following implicit numerical method (INM) for the problem (1)–(3)
U nþ1 ¼ U ni þ r 1 d2x U nþ1 þ r1 i i
n1 X
ðc Þ
ðc1 Þ
bjþ11 bj
r 2 U nþ1 d2x U nj i i
j¼0 n1 X s ðc Þ ðc Þ r2 bjþ12 bj 2 U nj þ ½f ðxi ; t nþ1 Þ þ f ðxi ; t n Þ 1 6 i 6 M 1; i 2 j¼0
U n0 ¼ uðnsÞ; U 0i ¼ gðihÞ;
U nM ¼ wðnsÞ; 0 6 i 6 M:
0 6 n 6 N;
0 6 n 6 N 1;
ð81Þ
ð82Þ ð83Þ
4038
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043 ðcÞ
where r1, r2 satisfy (23) and bj satisfies (5), and proved the stability and convergency of the scheme in l2-norm both theoretically and numerically. For the l2-norm depends on the mesh step, the l1-norm can describe the approximation degree of the numerical solutions better. By using the similar method and Lemma 3.3 of this paper, the scheme INM (81)–(83) can be proved to be stable and convergent in l1-norm theoretically. The details are not discussed any more. Here we only give the l1-convergency from numerical experiments, as well as the scheme IINM of [26]. Firstly, we prove the n l1-convergence order of the o four schemes INM, IINM of [26] and ICFDS (25)–(27) and IICFDS (75)– (77) of this paper. Let U nj j1 6 j 6 M 1; 0 6 n 6 N be the solution of the numerical scheme and denote
e1 ðs; hÞ ¼ max juNj U Nj j: 16j6M1
Suppose that p
e1 ðs; hÞ ¼ Oðsq þ h Þ: Then we use
log 2
ðe1 ð2p=q s; 2hÞÞ p: e1 ðs; hÞ
to compute the convergence order. For the case c1 = c2, we take c1 = c2 = 0.5. For the case c1 – c2, we take c1 = 0.9, c2 = 0.3 and c1 = 0.2, c2 = 0.8 respectively. Tables 1–3 provide some computational results with the schemes INM, IINM [26] and our two schemes ICFDS amd IICFDS at T = 1, which shows that the convergence order of our schemes ICFDS amd IICFDS is O(s + h4) and O(s2 + h4), respectively, while the order of INM and IINM [26] is only O(s + h2) and O(s2 + h2). Then we give some accuracy comparisons. Tables 4 and 5 provides some computational results of the maximum errors of INM [26] versus ICFDS, and IINM [26] versus IICFDS at T = 1 with the same h and s, respectively. From these data we can see that our two schemes are actually more accurate than the two schemes proposed in [26]. To test and verify the effectiveness of our two new schemes for short time, Tables 6 and 7 present some numerical maximum error comparisons for 0 6 t 6 1 with h = 1/10, s = 1/20,000 (Table 6) and h = 1/ 100, s = 1/1000 (Table 7) respectively. Similarly to test and verify the effectiveness of our two new schemes for long time, Tables 8 and 9 present some numerical maximum error comparisons for 1 6 t 6 10 with h = 1/10, s = 1/1000 (Table 8) and h = 1/100, s = 1/1000 (Table 9) respectively. From these tables, we see that our two compact schemes can reproduce the exact analytical solutions and also work better than the two schemes proposed in [26] for either short time or long time, especially IICFDS.
Table 1 Convergence order of scheme INM [26] and ICFDS for c1 = c2 = 0.5 at T = 1. h, s
e1(s, h)
e1 ð2p=q s;2hÞ e1 ðs;hÞ
s;2hÞÞ log 2 ðe1eð21 ðs;hÞ
Scheme INM [26] s = h = 1/8 h = 1/16, s = 1/32 h = 1/32, s = 1/128 h = 1/64, s = 1/516
– 7.846988e02 2.280452e02 6.074373e03 1.563170e03
– – 3.441 3.754 3.886
– – 1.783 1.909 1.958
Scheme ICFDS
– 8.786068e02 6.705252e03 4.358651e04 2.749260e05
– – 13.103 15.384 15.854
– – 3.720 3.943 3.987
s = h = 1/8 h = 1/16, s = 1/128 h = 1/32, s = 1/2048 h = 1/64, s = 1/32,768
p=q
Table 2 Convergence order of scheme IINM [26] and IICFDS for c1 = c2 = 0.5 at T = 1. h, s
e1(s, h)
e1 ð2p=q s;2hÞ e1 ðs;hÞ
s;2hÞ log 2 e1eð21 ðs;hÞ
Scheme IINM [26] s = h = 1/8 h = s = 1/16 h = s = 1/32 h = s = 1/64
– 1.01789e02 2.532596e03 6.324977e04 1.583664e04
– – 4.002 4.001 3.994
– – 2.008 2.000 1.998
Scheme IICFDS
– 7.939284e05 4.938659e05 3.083094e07 1.926444e08
– – 16.076 16.019 16.004
– – 4.007 4.002 4.000
s = h = 1/8 h = 1/16, s = 1/32 h = 1/32, s = 1/128 h = 1/64, s = 1/512
p=q
4039
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043 Table 3 Convergence order of scheme ICFDS for c1 = 0.9, c2 = 0.3 and c1 = 0.2, c2 = 0.8 at T = 1. h, s
e1(s, h)
e1 ð16s;2hÞ e1 ðs;hÞ
s;2hÞ log 2 e1eð16 1 ðs;hÞ
c1 = 0.9, c2 = 0.3 s = h = 1/8 h = 1/16, s = 1/128 h = 1/32, s = 1/2048 h = 1/64, s = 1/32,768
– 9.98366e02 6.890772e03 6.074373e04 2.743784e05
– – 14.488 15.777 15.929
– – 3.857 3.979 3.994
c1 = 0.2, c2 = 0.8 s = h = 1/7 h = 1/14, s = 1/112 h = 1/28, s = 1/1792 h = 1/56, s = 1/28,672
– 5.893423e02 5.676444e03 4.196787e04 2.853184e05
– – 10.382 13.526 14.709
– – 3.376 3.758 3.879
Table 4 Comparison of some maximum errors for INM [26] and ICFDS with c1 = c2 = 0.5 at T = 1.
s = 1/100,000
INM [26]
ICFDS
h = 1/4 h = 1/8 h = 1/16 h = 1/32
4.1759e2 1.0285e2 2.3523e3 6.3069e4
1.2845e3 7.0353e5 4.0881e6 8.7184e6
Table 5 Comparison of some maximum errors for IINM [26] and IICFDS with c1 = c2 = 0.5 at T = 1.
s = 1/5000
IINM [26]
IICFDS
h = 1/10 h = 1/20 h = 1/40 h = 1/80
4.5612e3 1.2907e3 3.6420e4 9.6578e5
3.2443e5 2.0217e5 1.2627e7 7.8954e9
Table 6 Comparison of some computational maximum errors with c1 = c2 = 0.5. t
INM [26]
ICFDS
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
4.7796e5 2.1914e4 5.2286e4 9.6227e4 1.5392e3 2.2552e3 3.1110e3 4.1015e3 5.2452e3 6.5246e3
3.4436e6 6.8604e6 9.8036e6 1.2163e5 1.3893e5 1.4974e5 1.5394e5 1.5141e5 1.4211e5 1.2596e5
Table 7 Comparison of some computational maximum errors with c1 = c2 = 0.5. t
IINM [26]
IICFDS
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
5.1147e7 2.2494e6 5.2883e6 9.6540e6 1.5360e5 2.2415e5 3.0824e5 4.0590e5 5.1716e5 6.4203e5
2.5351e11 1.1179e10 2.6388e10 4.8168e10 7.6760e10 1.1218e9 1.5446e9 2.0368e9 2.5978e9 3.2287e9
4040
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043
Table 8 Comparison of some computational maximum errors for h = 1/10, s = 1/1000. t
1 2 3 4 5 6 7 8 9 10
c1 = 0.25, c2 = 0.75
c1 = 0.75, c2 = 0.25
INM [26]
ICFDS
INM [26]
ICFDS
6.1580e3 2.6030e2 5.9273e2 1.0056e1 1.6460e1 2.3615e1 3.1198e1 4.1588e1 5.2367e1 6.4318e1
7.5722e4 1.5180e3 2.3338e3 2.8977e3 3.5075e3 4.0622e3 4.5615e3 5.0057e3 5.3949e3 5.7294e3
5.2246e3 2.5152e2 6.5865e2 1.1176e1 1.7880e1 2.6179e1 3.6079e1 4.7584e1 6.0698e1 7.5424e1
8.6978e4 1.7457e3 2.5562e3 3.2968e3 3.9960e3 4.5630e3 5.0874e3 5.5390e3 5.9174e3 6.2227e3
Table 9 Comparison of some computational maximum errors for h = 1/100, s = 1/1000. t
c1 = c2 = 0.25
c1 = c2 = 0.75
IINM [26]
IICFDS
IINM [26]
IICFDS
1 2 3 4 5 6 7 8 9 10
6.3581e5 2.4896e4 5.4610e4 9.4570e4 1.4396e3 2.0194e3 2.6784e3 3.4091e3 4.2055e3 5.0607e3
3.2287e9 1.3397e8 3.0644e8 5.5161e8 8.6717e8 1.2551e7 1.7159e7 2.2501e7 2.8580e7 3.5195e7
6.3596e5 2.7024e4 6.2248e4 1.2083e3 1.7655e3 2.5567e3 3.4943e3 4.5785e3 5.8093e3 7.1866e3
3.1523e9 1.3412e8 3.0920e8 5.5710e8 8.7801e8 1.2720e7 1.7392e7 2.2797e7 2.8935e7 3.5805e7
Finally we draw some error curves. Take M = 32, N = 32, 64, 128, 256, 512, respectively. Fig. 1 plots the curves of the max errors of the ICFNS solutions for c1 = c2 = 0.25, c1 = c2 = 0.75, c1 = 0.8, c2 = 0.3 and c1 = 0.2, c2 = 0.7, respectively, at T = 1. Then
0.02
0.03
N=32 N=64 N=128 N=256 N=512
0.018 0.016 0.014
0.02
0.012 0.01
0.015
0.008
0.01
0.006 0.004
0.005
0.002 0
0
0.2
0.4
0.6
0.8
0.03
1
0
0
0.2
0.4
0.6
0.8
0.018
N=32 N=64 N=128 N=256 N=512
0.025
1
N=32 N=64 N=128 N=256 N=512
0.016 0.014
0.02
0.012 0.01
0.015
0.008
0.01
0.006 0.004
0.005 0
N=32 N=64 N=128 N=256 N=512
0.025
0.002 0
0.2
0.4
0.6
0.8
1
0
0
0.2
0.4
0.6
0.8
1
Fig. 1. The curves of the max errors of the ICFDS at T = 1. Top left: c1 = c2 = 0.25. Top right: c1 = c2 = 0.75. Bottom left: c1 = 0.8, c2 = 0.3. Bottom right: c1 = 0.2, c2 = 0.7.
4041
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043
x 10 −8
x 10 −7 4
2 M=M=32
3 2
1
1
0.5
0 1
0.8
0.6
0.4 x
0.2
0 0
0.4
0.2
0.6
0.8
M=64,N=128
1.5
0 1
1
0.8
0.6
0.4 x
t
0.2
0 0
0.6
0.4
0.2
0.8
1
t
x 10 −8
x 10 −7 2
3 M=N=32 2
M=64,N=128
1.5 1
1 0 1
0.5
0.8
0.6
0.4 x
0.2
0 0
0.4
0.2
0.6
0.8
1
0 1
0.8
0.6
0.4 x
t
0.2
0 0
0.2
0.4
0.6
0.8
1
t
Fig. 2. The max error-surface of the IICFDS for c1 = c2 = 0.3 (up) and c1 = c2 = 0.8 (down) at T = 1.
take M = N = 32, and M = 64, N = 128. Fig. 2 plots max error curved surface of IICFNS solutions for c1 = c2 = 0.3 (up) and c1 = c2 = 0.8 (down), respectively. From these figures, we can see that the max errors become smaller as the time step s becomes smaller, which implies the two methods are actually efficient for the computation. Example 2. Consider the following fractional cable equation [26,32]: 2 @uðx; tÞ 1c @ uðx; tÞ 1c ¼ 0 Dt 1 0:50 Dt 2 uðx; tÞ; @t @x2 uð0; tÞ ¼ 0; uðL; tÞ ¼ 0; t 2 ð0; 0:1;
uðx; 0Þ ¼ 10dðx 5Þ;
x 2 ð0; 10Þ;
t 2 ð0; 0:1;
ð84Þ ð85Þ
x 2 ½0; 10;
ð86Þ
approximation lg(u(x,t))
where d(x) is the Dirac delta function. This problem does not exists an analytical solution. In [32], Liao and Sun did not consider the case of c1 – c2. Here take M = 32, N = 160. Fig. 3 plots the changes of u(x, t) for fixed c1 = 0.5 and different c2 = 0.2, 0.4, 0.6, 0.8, 1.0 at t = 0.1, respectively. And Fig. 4 plots the changes of u(x, t) for fixed c2 = 0.5 and different c1 = 0.2, 0.4, 0.6, 0.8, 1.0 at t = 0.1. These figures show the impacts of c1 and c2 on the anomalous behavior. We can see that there is a crossover between more and less anomalous behaviors, which is coincidence with the results of [26].
0.4 0.2 0 −0.2 −0.4 −0.6 4
4.5
5 x
5.5
6
Fig. 3. u(x, t) for fixed c1 = 0.5 and different c2 = 0.2, 0.4, 0.6, 0.8, 1.0 at t = 0.1.
4042
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043
approximation log2(u(x,t))
5
0
−5
−10
−15
−20
3
3.5
4
4.5
5 x
5.5
6
6.5
7
7.5
Fig. 4. u(x, t) for fixed c2 = 0.5 and different c1 = 0.2, 0.4, 0.6, 0.8, 1.0 at t = 0.1.
7. Conclusion In this paper, we build two implicit compact finite difference schemes for the fractional cable equation. For the l1-norm can describe the approximation degree of the numerical solutions, we prove the l1-stability and convergence of the first scheme by defining the new inner product. The proof techniques can also be applied to the other linear or nonlinear problems. Numerical experiments confirm the l1 convergence order of the two schemes and the numerical comparisons show that our two new schemes are more accurate than the two previous schemes. And the two compact schemes are effective and also work better for either short time or long time than the previous. Among the four schemes, the second compact scheme IICFDS can enhance the accuracy greatly and work best. Acknowledgment The authors are grateful to the referees for their constructive and valuable suggestions on this manuscript. References [1] R. Metzler, J. Klafter, The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics, J. Phys. A 37 (2004) R161–R208. [2] I.M. Sokolov, J. Klafter, From diffusion to anomalous diffusion: A century after Einstein’s Brownian motion, Chaos 15 (2005) 026103–026103-7. [3] J.A. Tenreiro Machado, Analysis and design of fractional-order digital control systems, Syst. Anal. Model. Simul. 27 (1997) 107–122. [4] D. Baleanu, O. Defterli, O.P. Agrawal, A central difference numerical scheme for fractional optimal control problems, J. Vip. Contr. 15 (2009) 583–597. [5] F. Mainardi, M. Raberto, R. Gorenflo, E. Scalas, Fractional calculus and continuous-time finance I: the waiting-time distribution, Physica A. 287 (2000) 468–481. [6] R. Magin, X. Feng, D. Baleanu, Solving the fractional order Bloch equation, Conc. Mag. Reson. Part A 34A (1) (2009) 16–23. [7] T.A.M. Langlands, B.I. Henry, The accuracy and stability of an implicit solution method for the fractional diffusion equation, J. Comput. Phys. 205 (2005) 719–736. [8] S.B. Yuste, Weighted average finite difference methods for fractional diffusion equations, J. Comput. Phys. 216 (2006) 264–274. [9] D.A. Murio, Implicit finite difference approximation for time fractional diffusion equation, Comput. Math. Appl. 56 (2008) 1138–1145. [10] Y. Zhang, A finite difference method for fractional partial differential equation, Appl. Math. Comput. 215 (2009) 524–529. [11] M.M. Meerschaert, H.P. Scheffler, C. Tadjeran, Finite difference methods for two-dimensional fractional dispersion equation, J. Comput. Phys. 211 (2006) 249–261. [12] C. Tadjeran, M.M. Meerschaert, H.P. Scheffler, A second-order accurate numerical approximation for the fractional diffusion equation, J. Comput. Phys. 213 (2006) 205–213. [13] C. Tadjeran, M.M. Meerschaert, A second-order accurate numerical method for the two-dimensional fractional diffusion equation, J. Comput. Phys. 220 (2007) 813–823. [14] C. Chen, F. Liu, I. Turner, V. Anh, A Fourier method for the fractional diffusion equation describing sub-diffusion, J. Comput. Phys. 227 (2007) 886–897. [15] C. Chen, F. Liu, K. Burrage, Finite difference methods and a Fourier analysis for the fractional reaction–subdiffusion equation, Appl. Math. Comput. 198 (2008) 754–769. [16] T. Langlands, B. Henry, The accuracy and stability of an implicit solution method for the fractional diffusion equation, J. Comput. Phys. 205 (2005) 719– 736. [17] P. Zhuang, F. Liu, V. Anh, I. Turner, New solution and analytical techniques of the implicit numerical method for the anomalous subdiffusion equation, SIAM J. Numer. Anal. 46 (2008) 1079–1095. [18] 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. [19] 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. [20] C. Chen, F. Liu, I. Turner, V. Anh, Numerical schemes and multivariate extrapolation of a two-dimensional anomalous sub-diffusion equation, Numer. Algor. 54 (2010) 1–21. [21] Z.Z. Sun, X. Wu, A fully difference scheme for a diffusion-wave system, Appl. Numer. Math. 56 (2006) 193–209. [22] M. Cui, Compact finite difference method for the fractional diffusion equation, J. Comput. Phys. 228 (2009) 7792–7804.
X. Hu, L. Zhang / Applied Mathematical Modelling 36 (2012) 4027–4043
4043
[23] R. Du, W.R. Cao, Z.Z. Sun, A compact difference scheme for the fractional diffusion-wave equation, Appl. Math. Model. 34 (2010) 2998–3007. [24] Guanghua, Gao, Z.Z. Sun, A compact difference scheme for the fractional sub-diffusion equation, J. Comput. Phys. 230 (2011) 586–595. [25] C. 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 (2010) 1740–1760. [26] T.A.M. Langlans, B. Henry,S. Wearne, Solution of a fractional cable equation: finite case,