Numerical approaches to fractional calculus and fractional ordinary ...
Journal of Computational Physics 230 (2011) 3352–3368
Contents lists available at ScienceDirect
Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp
Numerical approaches to fractional calculus and fractional ordinary differential equation q Changpin Li ⇑, An Chen, Junjie Ye Department of Mathematics, Shanghai University, Shanghai 200444, China
a r t i c l e
i n f o
Article history: Received 19 July 2010 Received in revised form 10 November 2010 Accepted 22 January 2011 Available online 31 January 2011 Keywords: Numerical approach Fractional calculus Fractional differential equations Piecewise interpolation Simpson method
a b s t r a c t Nowadays, fractional calculus are used to model various different phenomena in nature, but due to the non-local property of the fractional derivative, it still remains a lot of improvements in the present numerical approaches. In this paper, some new numerical approaches based on piecewise interpolation for fractional calculus, and some new improved approaches based on the Simpson method for the fractional differential equations are proposed. We use higher order piecewise interpolation polynomial to approximate the fractional integral and fractional derivatives, and use the Simpson method to design a higher order algorithm for the fractional differential equations. Error analyses and stability analyses are also given, and the numerical results show that these constructed numerical approaches are efficient. Ó 2011 Elsevier Inc. All rights reserved.
1. Introduction The history of fractional calculus [1–4] is more than 300 years old, but only in the recent decades did the applied scientists and the engineers realize that such fractional differential equations provided a better approach to describe the complex phenomena in nature, such as non-Brownian motion, signal processing, systems identification, control, viscoelastic materials and polymers. Because of the non-local property of the fractional derivative, it can be used to describe those complex systems which involve long-memory in time in a better way. Based on these requirements, the numerical approach has become a very desired method to analysis the experimental data which are described in a fractional way. Some papers have discussed the numerical approaches for fractional calculus and fractional differential equations. Here let us review the previous and present numerical methods. In [1], Oldham and Spanier introduced some numerical approaches to deal with the fractional calculus. They used the definitions of the Grünwald–Letnikov and Riemann–Liouville derivatives respectively. The main ideas of those numerical algorithms are that they used a Riemann sum to approximate the fractional derivative. Later on more researches are based on the above idea and also do certain theoretical analysis, such as [5,6]. In [2], Podlubny introduced the numerical approach for the arbitrary order derivative by using the definition of Riemann–Liouville based on the relationship between Grünwald–Letnikov derivative and the Riemann–Liouville derivative. The convergence order he got is O(h). Also he suggested that the short memory principle can be used to reduce the complexity at the expense of lacking certain precision. And he also pointed out that the convergence order can be improved to a more higher order by using the method of Lubich [7]. In [8], Deng described how to use the short memory principle combined with predictor–corrector approach in fractional differential equations in details. q The work was partially supported by the Natural Science Foundation of China under Grant No. 10872119 and Shanghai Leading Academic Discipline under Grant No. S30104. ⇑ Corresponding author. E-mail address: [email protected] (C. Li).
0021-9991/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.jcp.2011.01.030
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368
3353
In [9], Yuan and Agrawal designed a novel numerical approach which can avoid the harmful effect of non-local property in fractional derivative. This method was improved by Diethelm [10]. The main idea is that it transforms the original problem into an integral and then evaluates the numerical solution by Gauss integral formula. In [11], Murio designed a stable numerical algorithm to the fractional derivative in Caputo sense based on mollification techniques. Using a series of Chebyshev polynomials to evaluate the fractional calculus was posed in [12] and it seems that the method is highly effective. However that paper only dealt with the case 0 < a < 1. For the fractional differential equations, there are also some papers dealing with the numerical solutions. Some numerical algorithms such as predictor–corrector method [13,14,8,15], the Adomian decomposition [16], and the finite element [17] have been studied. Early in 1993, Tang studied the numerical scheme of the partial integro-differential equations in which the time derivative term was approximated by a Crank–Nicolson time-stepping and the integral term was treated by the product trapezoidal method [18]. The numerical scheme by the finite difference approach was further studied and an optimal convergence rate was also obtained, see [19]. Podlubny in [20,21] used the matrix expression to unify the formula where the fractional derivatives are derived from a finite difference. Liu also did some interesting works on the numerical approach for the fractional differential equations, see [22,23]. And Li and Xu in [24] constructed an efficient spectral method for the numerical approximation of the space–time fractional diffusion equation based on their proposed weak solution formulation, which is a significant difference to the direct discrete way. Based on the above review, we can see that there is not enough efficient numerical approach for the application of technology. Therefore in this paper, we pose the novel algorithms for a better convergence. The paper is organized as follows: In the present section we design some efficient approaches for the fractional calculus. In Section 2, we introduce the treatment for the fractional differential equations, where the error analyses and stability analyses are also given. Finally, we will give some numerical examples. First of all, let us introduce the definitions of fractional integral and derivatives. Definition 1. The definition of fractional integral of a function y reads as
J aa;t yðtÞ ¼
1 CðaÞ
Z
t
ðt sÞa1 yðsÞds;
a > 0:
ð1Þ
a
Definition 2. The Riemann–Liouville fractional derivative of y is defined as n
a
RL Da;t yðtÞ
¼
1 d Cðn aÞ dxn
Z
t
ðt sÞna1 yðsÞds;
t > a; n 1 < a < n 2 Z þ :
ð2Þ
a
Definition 3. The Grünwald–Letnikov fractional derivative of y is defined as
( ) N1 ðta Þa X Cðj aÞ ta N ; y tj GL Da;t yðtÞ ¼ lim N!1 N CðaÞ j¼0 Cðj þ 1Þ a
t > a; n 1 < a < n 2 Z þ :
ð3Þ
Definition 4. The Caputo fractional derivative of y is defined as a C Da;t yðtÞ ¼
1 Cðn aÞ
Z
t
ðt sÞna1 yðnÞ ðsÞds;
t > a; n 1 < a < n 2 Z þ :
ð4Þ
a
2. Numerical approach for the fractional calculus In this section, we divide the given interval [0, b] into:
D : 0 ¼ t0 < t1 < < t n ¼ t ¼ b: And ti = ih, in which h ¼
t ; n
ð5Þ
i ¼ 0; 1; . . . ; n.
2.1. Fractional integral Firstly, we use the piecewise linear interpolation function to approximate y,
yn ðtÞ ¼
t t iþ1 t ti yðt i Þ þ yðt iþ1 Þ; ti tiþ1 tiþ1 t i
where ti 6 t 6 ti+1, i = 0, 1, 2, . . . , n 1.
ð6Þ
3354
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368
So substitution in (1) with (6) gives, a
J 0;t yðtÞ
n1 X i¼0
¼
n1 X
a
J 0;tn yn ðtÞ ¼
¼
t iþ1
ti
) ðt n sÞa1 yn ðsÞds CðaÞ ! ! ) Z tiþ1 s tiþ1 ðt n sÞa1 s t i dt yðt i Þ þ dt yðtiþ1 Þ t i t iþ1 CðaÞ tiþ1 t i ti
t iþ1 ti
i¼0
( Z
i¼0
(Z n1 X
a1
ðt n sÞ CðaÞ
a n X h W yðt Þ; Cða þ 2Þ i¼0 1;i i
ð7Þ
where
W 1;0 ¼ na ða þ 1 nÞ þ ðn 1Þaþ1 ; W 1;i ¼ ðn i 1Þ
aþ1
þ ðn i þ 1Þ
aþ1
ð8Þ aþ1
2ðn iÞ
;
i ¼ 1; 2; . . . ; n 1;
ð9Þ
W 1;n ¼ 1:
ð10Þ
Thus we get the following numerical approach for fractional integral. Theorem 5. Suppose that y(t) 2 C2[0, b], the fractional integral has the following numerical approach:
J a0;tn yðtÞ ¼
a n X h 2 W yðt Þ þ Oðh Þ; Cða þ 2Þ i¼0 1;i i
ð11Þ
where W1,i are defined as (8)–(10). Proof. Let
EðtÞ ¼ yðtÞ yn ðtÞ ¼
yð2Þ ðtÞ ðt t i Þðt t iþ1 Þ; 2!
where t 2 (ti,ti+1). Thus
kJ a0;tn yn ðtÞ J a0;tn yðtÞk1 ¼ kJ a0;tn ðyn ðtÞ yðtÞÞk1 6 kyn ðtÞ yðtÞk1
1 CðaÞ
Z
tn
ðt n sÞa1 ds ¼
0
t an ky ðtÞ yðtÞk1 Cða þ 1Þ n
a
tan b kyð2Þ k1 2 2 ¼ kyn ðtÞ yðtÞk1 ¼ Oðh Þ ¼ Oðh Þ: Cða þ 1Þ Cða þ 1Þ
Next, we apply some much higher convergence order for y. So the piecewise cubic interpolation function will be our choice. We have the piecewise cubic interpolation function as follows:
yn ðtÞ ¼
2 2 2 t ti t tiþ1 t t iþ1 t ti t tiþ1 12 yðti Þ þ 1 2 yðtiþ1 Þ þ ðt ti Þ y0 ðti Þ t i t iþ1 t i t iþ1 t iþ1 ti t iþ1 ti ti tiþ1 þ ðt t iþ1 Þ
t ti tiþ1 t i
2
y0 ðt iþ1 Þ:
ð12Þ
Substituting (12) into (1) yields a
J 0;t yðtÞ
n1 X
a
J 0;tn yn ðtÞ ¼
i¼0
þ
i¼0
Z
tiþ1
ti
þ
n1 X
Z ti
t iþ1
( Z ti
t iþ1
! ðt n sÞa1 s ti s tiþ1 2 ds yðt i Þ 12 CðaÞ t i t iþ1 t i t iþ1
! ! Z tiþ1 ðtn sÞa1 s tiþ1 s ti 2 ðt n sÞa1 s tiþ1 2 ds yðtiþ1 Þ þ ds y0 ðti Þ 12 ðs ti Þ CðaÞ tiþ1 t i tiþ1 t i CðaÞ ti tiþ1 ti ! ) a aþ1 n n X X ðtn sÞa1 s ti 2 h h 0 f y0 ðt Þ; W ds y ðt iþ1 Þ ¼ W 2;i yðt i Þ þ ðs tiþ1 Þ CðaÞ t iþ1 ti Cða þ 4Þ i¼0 Cða þ 4Þ i¼0 2;i i ð13Þ
3355
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368
where
W 2;0 ¼ 6ðn 1Þ2þa ð1 þ 2n þ aÞ þ na 12n3 6n2 ð3 þ aÞ þ ð1 þ aÞð2 þ aÞð3 þ aÞ ; W 2;i ¼ 6 4ðn iÞ3þa þ ðn i 1Þ2þa ð1 þ 2i 2n aÞ þ ð1 i þ nÞ2þa ð1 þ 2i 2n þ aÞ ;
ð14Þ i ¼ 1; 2; . . . ; n 1;
W 2;n ¼ 6ð1 þ aÞ; f 2;0 ¼ 2ðn 1Þ2þa ð3n þ aÞ þ n1þa 6n2 4nð3 þ aÞ þ ð2 þ aÞð3 þ aÞ ; W
ð15Þ ð16Þ ð17Þ
f 2;i ¼ 2ðn i 1Þaþ2 ð3i 3n aÞ 2ðn i þ 1Þaþ2 ð3i 3n þ aÞ ðn iÞaþ2 ð24 þ 8aÞ; W f 2;n ¼ 2a: W
i ¼ 1; 2; . . . ; n 1;
ð18Þ ð19Þ
So we can also get the other following numerical approach for fractional integral. Theorem 6. Suppose y(t) 2 C4[0, b], the fractional integral has the following numerical approach:
J a0;tn yðtÞ ¼
a aþ1 n n X X h h f y0 ðt Þ þ Oðh4 Þ; W W 2;i yðt i Þ þ Cða þ 4Þ i¼0 Cða þ 4Þ i¼0 2;i i
ð20Þ
f 2;i are defined as (14)–(19). where W2,i and W Proof. The numerical formula (20) is obvious. In effect, set
EðtÞ ¼ yðtÞ yn ðtÞ ¼
yð4Þ ðtÞ ðt ti Þ2 ðt t iþ1 Þ2 ; 4!
ð21Þ
where t 2 (ti,ti+1). Thus
kJ a0;tn yn ðtÞ J a0;tn yðtÞk1 ¼ kJ a0;tn ðyn ðtÞ yðtÞÞk1 6 kyn ðtÞ yðtÞk1
1 CðaÞ
Z
tn
ðt n sÞa1 ds ¼
0
t an ky ðtÞ yðtÞk1 Cða þ 1Þ n
a
tan b kyð4Þ k1 4 4 ¼ kyn ðtÞ yðtÞk1 ¼ Oðh Þ ¼ Oðh Þ: Cða þ 1Þ Cða þ 1Þ
In (20), the first-order derivative can be approximated by the back forward difference or the central difference. Since the spline function has good properties to fit the experimental data, in the following we use the spline function to approximate y. The cubic interpolation spline function reads as follows:
2 2 2 t ti t tiþ1 t tiþ1 t ti t t iþ1 yn ðtÞ ¼ 1 2 yðt i Þ þ 1 2 yðtiþ1 Þ þ ðt ti Þ mi ti tiþ1 ti tiþ1 t iþ1 ti t iþ1 ti ti tiþ1 2 t ti þ ðt tiþ1 Þ miþ1 ; t iþ1 t i
ð22Þ
where mi (i = 0, 1, 2, . . . , n 1) are the undetermined real numbers. Then one has
J a0;tn yðtÞ
Z
tn
0
ðt n sÞa1 yn ðsÞds: CðaÞ
ð23Þ
That is, a
J 0;tn yðtÞ
n1 Z X i¼0
þ þ
Z Z
t iþ1 ti tiþ1
ti
ti
t iþ1
! ( Z n1 t iþ1 X ðt n sÞa1 ðt n sÞa1 s ti s tiþ1 2 ds yðti Þ yn ðsÞds ¼ 12 CðaÞ CðaÞ ti t iþ1 t i t iþ1 ti i¼0 ! ! Z tiþ1 ðtn sÞa1 s tiþ1 s ti 2 ðt n sÞa1 s tiþ1 2 ds yðt iþ1 Þ þ ds mi 12 ðs ti Þ CðaÞ tiþ1 t i tiþ1 t i CðaÞ ti tiþ1 ti ) 2 ! a aþ1 a1 n n X X ðtn sÞ s ti h h f m; W ds miþ1 ¼ W 3;i yðt i Þ þ ð24Þ ðs t iþ1 Þ CðaÞ t iþ1 ti Cða þ 4Þ i¼0 Cða þ 4Þ i¼0 3;i i
f 3;i ¼ W f 2;i ; i ¼ 0; 1; 2; . . . ; n. where W 3;i ¼ W 2;i ; W For the undetermined real number mi, i = 1, 2, . . . , n 1, the expressions are as follows:
ki mi1 þ 2mi þ li miþ1 ¼ C i ;
i ¼ 1; 2; . . . ; n 1;
in which
t iþ1 t i ; t iþ1 ti1 t t li ¼ 1 ki ¼ i i1 ; t iþ1 ti1 C i ¼ 3½ki yðti1 ; ti Þ þ li yðti ; t iþ1 Þ: ki ¼
ð25Þ
3356
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368
The above formulas have only n 1 linear equations, but the unknown numbers are n + 1. In order to determine the unknown mi, here we present two ways: 1. Taking first-order derivative of spline interpolation function at boundary as m0 ¼ fn0 ðx0 Þ; mn ¼ fn0 ðxn Þ, then we get the following m1,i:
m1;i
8 0 i ¼ 0; > < f ðx0 Þ; ¼ f 0 ðxn Þ; i ¼ n; > : the solution mi of ð25Þ; i ¼ 1; 2; . . . ; n 1:
ð26Þ
In (25) mi is substituted by m1,i, C1 and Cn1 are updated by C1 k1m1,0 and Cn1 ln1m1,n respectively. Under these conditions, we get a system of Eq. (25) whose coefficient matrix is a triangular one. So the method of forward elimination and backward substitution can be used. 2. Taking fn00 ðx0 Þ ¼ f 00 ðx0 Þ; fn00 ðxn Þ ¼ f 00 ðxn Þ, then we get the following equations with respect to m1,i:
8 x1 x0 00 > < 2m1;0 þ m1;1 ¼ 3f ðx0 ; x1 Þ 2 f ðx0 Þ; ki m1;i1 þ 2m1;i þ li m1;iþ1 ¼ C i ; i ¼ 1; 2; . . . ; n 1; > : m1;n1 þ 2m1;n ¼ 3f ðxn1 ; xn Þ þ xn x2 n1 f 00 ðxn Þ:
ð27Þ
Thus we also get a system of equations whose coefficient matrix is also a triangular one. Thus, we have another numerical approach for fractional integral. Theorem 7. Suppose y(t) 2 C2[0, b]. Then the fractional integral has the following numerical approach:
J a0;tn yðtÞ ¼
a aþ1 n n X X h h f m þ Oðh4 Þ; W W 3;i yðt i Þ þ Cða þ 4Þ i¼0 Cða þ 4Þ i¼0 3;i 1;i
ð28Þ
f 3;i are defined as (14)–(19), and m1,i can be determined by the above suggestion (26) or (27). where W3,i and W
Proof. According to [25], we have
jyn ðtÞ yðtÞj 6
7 4 Lh ; 8
where L is a constant. Thus
kJ a0;tn yn ðtÞ J a0;tn yðtÞk1 ¼ kJ a0;tn ðyn ðtÞ yðtÞÞk1 6 kyn ðtÞ yðtÞk1
1 CðaÞ
Z
tn
ðt n sÞa1 ds ¼
0
t an ky ðtÞ yðtÞk1 Cða þ 1Þ n
a
tan 7b Lkyð4Þ k1 4 4 ¼ kyn ðtÞ yðtÞk1 ¼ Oðh Þ ¼ Oðh Þ: Cða þ 1Þ 8Cða þ 1Þ
2.2. Fractional derivative Due to a wide application of the Caputo derivative, we put our attention to the numerical approach of it. Oldham and Spanier [1] have done some researches on the numerical approach to Grünwald–Letnikov derivative and Riemann–Liouville derivative. Besides, Liu et al. used the finite difference to approximate Riemann–Liouville derivative based on the relationship between Grünwald–Letnikov derivative and Riemann–Liouville derivative. It shows that the approaches are very efficient. Here we approximate the Caputo derivative which is similar to the above discussions for the fractional integral. In the following, in order to compute the numerical solution y(m), we always assume that it makes sense in the small neighborhood of the two endpoints. Theorem 8. When y(t) 2 C2+m[0, b], m = [a] + 1, then the Caputo derivative has the following numerical approach: a
¼
ma n X h 2 W yðmÞ ðt i Þ þ Oðh Þ; Cðm a þ 2Þ i¼0 5;i
ð29Þ
a
¼
ma maþ1 n n X X h h f m þ Oðh4 Þ: W W 6;i yðmÞ ðt i Þ þ Cðm a þ 4Þ i¼0 Cðm a þ 4Þ i¼0 6;i 2;i
ð30Þ
C D0;t n yðtÞ
C D0;t n yðtÞ
And when y(t) 2 C4+m[0, b], m = [a] + 1, then the Caputo derivative has the following numerical approach: a
C D0;t n yðtÞ
¼
ma maþ1 n n X X h h f yðmþ1Þ ðt i Þ þ Oðh4 Þ; W W 7;i yðmÞ ðt i Þ þ Cðm a þ 4Þ i¼0 Cðm a þ 4Þ i¼0 7;i
ð31Þ
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368
3357
f 6;i are defined as follows: where W5,i, W6,i and W
8 > nma ðm a þ 1 nÞ þ ðn 1Þmaþ1 ; i ¼ 0; > > > < ðn i 1Þmaþ1 þ ðn i þ 1Þmaþ1 ; W 5;i ¼ > > 2ðn iÞmaþ1 ; i ¼ 1; . . . ; n 1; > > : 1; i ¼ n; 8 6ðn 1Þ2þma ð1 þ 2n þ m aÞ þ nma ð12n3 6n2 ð3 þ m aÞ > > > > > > < þð1 þ m aÞð2 þ m aÞð3 þ m aÞÞ; i ¼ 0; W 6;i ¼ 6ð4ðn iÞ3þma þ ðn i 1Þ2þma ð1 þ 2i 2n m þ aÞ > > > > þð1 i þ nÞ2þma ð1 þ 2i 2n þ m aÞÞ; i ¼ 1; . . . ; n 1; > > : 6ð1 þ m aÞ; i ¼ n; 8 2ðn 1Þ2þma ð3n þ m aÞ > > > > > þn1þma 6n2 4nð3 þ m aÞ þ ð2 þ m aÞð3 þ m aÞ ; i ¼ 0; > > > > maþ2 < ð3i 3n m þ aÞ f 6;i ¼ 2ðn i 1Þ W maþ2 > > ð3i 3n þ m aÞ 2ðn i þ 1Þ > > > > maþ2 > ð24 þ 8ðm aÞÞ; i ¼ 1; 2; . . . ; n 1; ðn iÞ > > : 2a; i ¼ n:
As to (29), y(m)(ti) have been given by [26,27]
8 yðt þhÞyðt hÞ 2 i i þ Oðh Þ; > 2h > > > 2 > yðt i þhÞ2yðt i Þþyðt i hÞ > þ Oðh Þ; > > h2 > > > 2n > P < 2 ð2nÞ! 1 ð1Þk k!ð2nkÞ! yðt i þ ðn kÞhÞ þ Oðh Þ; ym ðti Þ :¼ um ðhÞ ¼ h2n k¼0 > > > 2nþ1 > P > ð2nþ1Þ! 1 > > 2nþ1 ð1Þk k!ð2nþ1kÞ! > h > > k¼0 > > : 2 12 ðyðt i þ ðn kÞhÞ þ yðt i þ ðn þ 1 kÞhÞÞ þ Oðh Þ; And combining the above formula with Richardson extrapolation, y
m ¼ 1; m ¼ 2; m ¼ 2j;
m ¼ 2j þ 1; j ¼ 1; 2; . . .
(m)
(ti) to (30) and (31) are obtained as follows:
4 h 1 4 um ðhÞ þ Oðh Þ: yðmÞ ðt i Þ ¼ um 3 2 3 f 7;i ¼ W f 6;i , and m2,i are those values m1,i given in Theorem 7 in which y(t) is substituted by y(m)(t). And W7,i = W6,i, W Proof. The proof is simple. Once combining the definition of Caputo derivative and the already designed numerical approaches for the fractional integral, the remain thing is only the substitution. h
3. Numerical approach for the fractional differential equations In this section, we discuss a new numerical scheme for solving the following system: a
C D0;t yðtÞ
¼ f ðt; yðtÞÞ;
ðkÞ
yðkÞ ð0Þ ¼ y0 ;
k ¼ 0; . . . ; ½a:
ð32Þ
It is well known that the initial value problem (32) is equivalent to a Volterra integral equation
yðtÞ ¼
½a X
m¼0
ðmÞ
y0
tm 1 þ 2 tt ðt uÞa1 f ðu; yðuÞÞdu m! CðaÞ 0
ð33Þ
in the sense that a continuous function is a solution of (32) if and only if it is also a solution of (33). Here, we indicate the approach briefly. The main problem is to solve the integral on the right-hand side of (33) by numerical method. For the consequence of the non-local structure of the fractional-order differential operators, we use the compound quadrature formula with respect to the weight function (tk+1 )a1 to replace the integral naturally. The product trapezoidal quadrature formula is a good choice. The scheme constructed by using the trapezoidal quadrature formula to replace the integral is proved to work [28,6]. We may also construct the higher order of scheme by using the compound Simpson formula. We need not only nodes tj (j = 0,1,. . .,k + 1) but also t jþ1 ðj ¼ 0; 1; . . . ; kÞ. In other words, we apply the approximation 2
3358
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368
Z
t kþ1
ðtkþ1 zÞa1 gðzÞdz ¼
0
Z
tkþ1
ðt kþ1 zÞa1 g~kþ1 ðzÞdz;
ð34Þ
0
where g~kþ1 is the piecewise quadratic interpolation for g with nodes and knots chosen at the tj ; t jþ1 ðj ¼ 0; 1; 2; . . . ; k þ 1Þ. 2 Using the standard techniques of the Simpson theory, we write the integral on the right-hand side of (34) as
Z
t kþ1
ðtkþ1 zÞa1 g~kþ1 ðzÞdz ¼
0
kþ1 X
aj;kþ1 gðt j Þ þ
j¼0
k X
bj;kþ1 gðt jþ1 Þ;
ð35Þ
2
j¼0
where
aj;kþ1
h i 8 aþ2 aþ2 4ha > ðk þ 2Þ ðk þ 1Þ > a ð a þ1Þð a þ2Þ > > > > > h i > a a > > aðahþ1Þ ðk þ 1Þaþ1 3ðk þ 2Þaþ1 þ ha ðk þ 2Þa ; j ¼ 0; > > > > > < h i a aþ2 4ha ¼ ðk jÞaþ2 aðahþ1Þ aðaþ1Þðaþ2Þ ðk þ 2 jÞ > > > > > h i > > > ðk þ 2 jÞaþ1 þ 6ðk þ 1 jÞaþ1 þ ðk jÞaþ1 ; 1 6 j 6 k; > > > > > > > : 2ha a aðaþ1Þðaþ2Þ ð1 2Þ; j ¼ k þ 1; h i h i 8h 4h ðk þ 2 jÞaþ2 ðk þ 1 jÞaþ2 ðk þ 2 jÞaþ1 þ ðk þ 1 jÞaþ1 ; aða þ 1Þða þ 2Þ aða þ 1Þ a
bj;kþ1 ¼
ð36Þ
a
0 6 j 6 k: ð37Þ
In order to construct the explicit scheme avoiding iterations. We replace the integral on the right-hand side of Eq. (33) by the product rectangle rule to calculate the predictor of yk+1, then
ykþ1
" # ½a k k X X t jkþ1 ðjÞ 1 X P ¼ a f ðt ; y Þ þ akþ1;kþ1 f ðt kþ1 ; ykþ1 Þ þ bj;kþ1 f ðt jþ1 ; ykþ1 Þ ; y þ 2 2 j! 0 CðaÞ j¼0 j;kþ1 j j j¼0 j¼0
yPkþ1 ¼
ð38Þ
½a k X t jkþ1 ðjÞ 1 X c f ðt ; y Þ; y0 þ j! CðaÞ j¼0 j;kþ1 j j j¼0
ð39Þ
where aj,k+1, bj,k+1 are shown in (36) and (37) respectively, and
h
cj;kþ1 ¼
a
a
½ðk þ 1 jÞa ðk jÞa :
ð40Þ
The remaining problem is to determine the values ykþ1 . We use the product rectangle rule to calculate the ykþ1 . Thus 2
ykþ1 ¼ 2
j ½a t X kþ1 2
j¼0
j!
ðjÞ
y0 þ
2
k 1 X d f ðt ; y Þ; CðaÞ j¼0 j;kþ1 j j
ð41Þ
where
dj;kþ1 ¼
h
a
kþ
a
1 j 2
a
a 1 ; k j 2
0 6 j 6 k:
ð42Þ
Thus a new scheme is described by (38), (39) and (41) with the weights aj,k+1, bj,k+1, cj,k+1, dj,k+1 defined as above, respectively. 3.1. Improved algorithms In this section, we want to improve the above designed algorithm. Consider the predictor–corrector format to the intermediate nodes, we have
ykþ1 ¼ 2
yPkþ1 ¼ 2
j ½a t X kþ1 2
j¼0
j!
2
j!
þ
" k X
# P ej;kþ1 f ðt j ; yj Þ þ ekþ1;kþ1 f t kþ1 ; ykþ1 ;
j¼0
j ½a t 1 X kþ j¼0
ðjÞ y0
ðjÞ
y0 þ
k 1 X d f ðt ; y Þ; CðaÞ j¼0 j;kþ1 j j
2
2
ð43Þ
ð44Þ
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368
3359
where dj,k+1 are defined as (42) and ej,k+1 are given by
8 1a a h aþ1 i ðkþ2Þ h ha 1 aþ1 > > k þ 12 ; j ¼ 0; > Cðaþ1Þ þ Cðaþ2Þ k 2 > > h i > > a aþ1 aþ1 > 1 aþ1 a þ1 h 3 > 3ð12 Þ j ¼ k; þ Cðahþ1Þ 12 ; > > > Cðaþ2Þ 2 > > a > 1 aþ1 : 2h ; j ¼ k þ 1:
ej;kþ1
ð45Þ
Cðaþ2Þ 2
At the beginning, we apply the product rectangle rule to calculate it
a aha 1 ; Cða þ 2Þ 2
e0;1 ¼ e1;1
a a h 1 ¼ : Cða þ 2Þ 2
ð46Þ ð47Þ
We call the above algorithm Improved Algorithm I (descried by (38), (39), (43) and (44)), namely, using the intermediate nodes with the predictor–corrector method (similar to the classical Adams methods), which can be upgraded to ensure an overall accuracy. Improved accuracy of this algorithm found in the middle node is still inadequate, so we can use the formula of Simpson to calculate the intermediate nodes. Thus
ykþ1 ¼ 2
yPkþ1 ¼ 2
½a k k X X tjkþ1 ðjÞ X ej;kþ1 f ðt j1 ; yj1 Þ þ ekþ1;kþ1 f tkþ1 ; yPkþ1 þ fj;kþ1 f ðt j ; yj Þ þ fkþ1;kþ1 f t kþ1 ; yPkþ1 ; y0 þ 2 2 2 2 j! j¼1 j¼0 j¼0 j ½a t 1 X kþ 2
j¼0
j!
ðjÞ
y0 þ
k 1 X d f ðt ; y Þ; CðaÞ j¼0 j;kþ1 j j
ð48Þ
ð49Þ
where
dj;kþ1
8 h i < ha k þ 1 j a k 1 j a ; 0 6 j 6 k 1; 2 2 a ¼ : ha 1a ; j ¼ k; a 2
ð50Þ
ej;kþ1
8
a aþ1 aþ1 a aþ2 aþ2 > < 4 ðkþ32jÞ þðkjþ12Þ h ðkþ32jÞ ðkjþ12Þ 8h ; 1 6 j 6 k; Cðaþ2Þ Cðaþ3Þ ¼ > : 2ha 1 aþ1 ; j ¼ k þ 1; Cðaþ2Þ 2
ð51Þ
fj;kþ1
8 aþ1 aþ1 a aþ2 aþ2 a 3ðkþ12Þ þðk12Þ h 4 ðk12Þ ðkþ12Þ h > ha 1 a > k þ ; j ¼ 0; > > Cðaþ1Þ 2 Cðaþ2Þ Cðaþ3Þ > <
aþ2 aþ2 a aþ1 aþ1 aþ1 a 3 1 3 1 1 ¼ 4 ðkþ2jÞ ðkj2Þ h ðkþ2jÞ þ6ðkjþ2Þ þðkj2Þ h ; 1 6 j 6 k; > Cðaþ3Þ Cðaþ2Þ > > h h > 1aþ1 i aþ2 i > 3 aþ1 4ha 3 aþ2 : ha þ 5 12 þ ; j ¼ k þ 1: 2 2 2 Cðaþ2Þ Cðaþ3Þ
ð52Þ
The above algorithm is called Improved Algorithm II (described by (38), (39), (48) and (49)), namely, using the Simpson formats to calculate the intermediate nodes in order to enhance accuracy. In order to further improve the algorithm’s overall accuracy, there is still room for improvement in the predictor. Observing the above algorithm, one can find at each step whatever format we used in the predictor, its main format is not changed. If the value after correction continues to serve as a predictor, one can improve the accuracy of the whole algorithm, so as to enhance the overall accuracy of the format. However the computational complexity is almost not increased. Thus
yPkþ1 ¼
" # ½a k k X X t jkþ1 ðjÞ 1 X aj;kþ1 f ðtj ; yj Þ þ akþ1;kþ1 f ðt kþ1 ; e y Pkþ1 Þ þ bj;kþ1 f tjþ1 ; ykþ1 ; y0 þ 2 2 j! CðaÞ j¼0 j¼0 j¼0
ð53Þ
ykþ1 ¼
½a k k X X tjkþ1 ðjÞ X ~Pkþ1 Þ; ej;kþ1 f t j1 ; yk1 þ ekþ1;kþ1 f t kþ1 ; ypkþ1 þ fj;kþ1 f ðtj ; yj Þ þ fkþ1;kþ1 f ðt kþ1 ; y y0 þ 2 2 2 2 j! j¼1 j¼0 j¼0
ð54Þ
ypkþ1 ¼
½a k k X X tjkþ1 ðjÞ X ej;kþ1 f t j1 ; yk1 þ ekþ1;kþ1 f t kþ1 ; yPkþ1 þ fj;kþ1 f ðtj ; yj Þ þ fkþ1;kþ1 f ðt kþ1 ; e y Pkþ1 Þ; y0 þ 2 2 2 2 j! j¼1 j¼0 j¼0
ð55Þ
2
2
3360
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368
where e y Pkþ1 is updated by (39). The remaining format is the same as that of Improved Algorithm II. The above algorithm is called Improved Algorithm III (described by (38), (53), (54) and (55)), namely, based on the Improved Algorithm II, the value after correction continues to serve as the predictor value by using the Simpson formula to calculate the intermediate nodes. 3.2. Error analysis First of all, we focus on the errors of the compound Simpson formulas and a corresponding result for compound rectangle rule and compound trapezoidal formula. They are easy to be verified by using the Simpson theory, we omit the details. Lemma 1. Let z 2 C3[0,T], then
Z kzð3Þ k tkþ1 kþ1 k X X a1 1 a1 4 ðt kþ1 tÞ zðtÞdt aj;kþ1 zðt j Þ bj;kþ1 zðt jþ1 Þ 6 t kþ1 h : 2 0 4! j¼0 j¼0
ð56Þ
Lemma 2. Let z 2 C1[0,T], then
Z kz0 k tkþ1 k X a1 1 a ðt kþ1 tÞ zðtÞdt cj;kþ1 zðt j Þ 6 t h: 0 a kþ1 j¼0
ð57Þ
Lemma 3. Let z 2 C1[0,T], then
Z kz0 k tkþ1 k X a1 1 a ðt kþ1 tÞ zðtÞdt dj;kþ1 zðt j Þ 6 t h: 0 a kþ1 j¼0
ð58Þ
Then we give the error analysis for the (39). Lemma 4. Assume that the initial value problem solution y satisfies
Z tkþ1 k X d ðt kþ1 tÞa1 C Da yðtÞdt cj;kþ1C Da yðtj Þ 6 Ctmkþ1 h 0 j¼0
ð59Þ
with m P 0 and d > 0. Then we have
yðt kþ1 Þ yP 6 C 1 hd þ C 2 max yðt j Þ y ; j kþ1 06j6k
where C1, C2 > 0. Proof. By the construction of yPkþ1 , we find that
yðt kþ1 Þ yP j ¼ kþ1
Z k X 1 tkþ1 ðt kþ1 tÞa1 f ðt; yðtÞÞdt cj;kþ1 f ðt j ; yj Þj CðaÞ 0 j¼0 Z k k X 1 tkþ1 1 X a1 a a 6 ðt kþ1 tÞ CD yðtÞdt cj;kþ1CD yðtj Þ þ c jf ðtj ; yðt j ÞÞ f ðt j ; yj Þj CðaÞ j¼0 j;kþ1 CðaÞ 0 j¼0 6
k k X C 1 tmkþ1 d 1 X d d cj;kþ1 jyðt j Þ yj j 6 C 1 h þ C 2 max jyðt j Þ yj j cj;kþ1 6 C 1 h þ C 2 max jyðt j Þ yj j: h þ 06j6k 06j6k CðaÞ CðaÞ j¼0 j¼0
Based on the error estimate of the preceding section, we now present a general convergence result for this new scheme. Theorem 9. Assume that the initial value problem solution y satisfies
Z tkþ1 k X d a1 a a 1 ðtkþ1 tÞ C D yðtÞdt cj;kþ1C D yðt j Þ 6 C 1 tmkþ1 h 1; 0 j¼0 Z tkþ1 k X d a1 2 a a 2 ðtkþ1 tÞ C D yðtÞdt dj;kþ1C D yðtj Þ 6 C 2 tmkþ1 h 2; 0 j¼0 Z tkþ1 k k X X d 3 ðtkþ1 tÞa1 C Da yðtÞdt aj;kþ1C Da yðt j Þ akþ1;kþ1C Da yðt j Þ bj;kþ1C Da yðt jþ1 Þ 6 C 3 t mkþ1 h 3; 2 0 j¼0 j¼0
3361
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368
with m1, m2, m3 P 0 and d1, d2, d3 > 0. Then we have q
max jyðtj Þ yj j ¼ Oðh Þ;
06j6N
where q = min{d1 + a, d2, d3} and N = [T/h]. Proof. Firstly by Lemma 4, we have d
jyðtkþ1 Þ yPkþ1 j 6 C 1 h 1 þ C 2 max jyðtj Þ yj j: 06j6k
Secondly, from formula of ykþ1 and Lemma 4, we easily have 2
d
jyðtkþ1 Þ ykþ1 j 6 C 3 h 2 þ C 4 2
2
So
k X
d
dj;kþ1 jyðtj Þ yj j 6 C 3 h 2 þ C 4 max jyðt j Þ yj j: 06j6k
j¼0
Z k k X X 1 tkþ1 ðt kþ1 tÞa1 f ðt; yðtÞÞdt aj;kþ1 f ðt j ; yj Þ akþ1;kþ1 f ðtkþ1 ; yPkþ1 Þ bj;kþ1 f t jþ1 ; yjþ1 2 2 CðaÞ 0 j¼0 j¼0 Z kþ1 k k X X 1 tkþ1 1 X 6 ðt kþ1 tÞa1 C Da yðtÞdt aj;kþ1 C Da yðtj Þ bj;kþ1CD a y t jþ1 þ a ðf ðt j ; yðt j ÞÞ 2 CðaÞ 0 CðaÞ j¼0 j;kþ1 j¼0 j¼0
jyðtkþ1 Þ ykþ1 j ¼
k 1 1 X bj;kþ1 f tjþ1 ; yðt jþ1 Þ f t jþ1 ; yjþ1 akþ1;kþ1 ðf ðtkþ1 ; yðtkþ1 ÞÞ f ðtkþ1 ; yPkþ1 ÞÞ þ 2 2 2 2 CðaÞ CðaÞ j¼0 k X d d 6 C5h 3 þ C6 aj;kþ1 jyðtj Þ yj j þ C 7 akþ1;kþ1 C 1 h 1 þ C 2 max jyðtj Þ yj j f ðt j ; yj ÞÞ þ
06j6k
j¼0
þ C8
k X
d
bj;kþ1 C 3 h 2 þ C 4 max jyðt j Þ yj j
06j6k
j¼0 d1 þa
d
6 C5h 3 þ C7h
minfd3 ;d1 þa;d2 g
d
þ C 8 h 2 þ C 9 max jyðt j Þ yj j 6 Ch 06j6k
:
In the following, we mainly deal with the error analysis for the improved algorithms I–III respectively. Theorem 10. Assume that the initial value problem solution y satisfies
Z tkþ1 k X d 2 ðt kþ1 tÞa1 C Da yðtÞdt dj;kþ1C Da yðtj Þ 6 C 1 t mkþ1 1 h 1 ; 2 0 2 j¼0 Z tkþ1 k X d a1 a 2 ðt kþ1 tÞ C D yðtÞdt ej;kþ1C Da yðt j Þ ekþ1;kþ1C Da yðt jþ1 Þ 6 C 2 t mkþ2 1 h 2 ; 2 2 0 2 j¼0 with q
max jyðtj Þ yj j ¼ Oðh Þ;
06j6N
where m1, m2 P 0 and d1, d2 > 0 for some n, T > 0. Then we have p
max jyðt jþ1 Þ yjþ1 j ¼ Oðh Þ;
06j6N1
2
2
where p = min(q, d1 + a, d2) and N = [T/h]. Proof. Firstly, we find that
Z k X 1 tkþ12 a1 ðt kþ1 tÞ f ðt; yðtÞÞdt dj;kþ1 f ðtj ; yj Þ 2 CðaÞ 0 j¼0 Z k k X X 1 tkþ12 d 6 ðtkþ1 tÞa1 Da yðtÞdt dj;kþ1 Da yðtj Þ 6 C 1 h 1 þ C 2 max jyðtj Þ yj j dj;kþ1 2 06j6k CðaÞ 0 j¼0 j¼0
yðt kþ1 Þ yPkþ1 ¼ 2
2
d
6 C 1 h 1 þ C 2 max jyðt j Þ yj j: 06j6k
3362
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368
Finally, one gets that
Z k k a1 X 1 tkþ12 1 X tkþ1 t f ðt; yðtÞÞdt ej;kþ1 f ðt j ; yj Þ ekþ1;kþ1 f t kþ1 ; ykþ1 þ e jf ðt ; yðt j ÞÞ 2 2 2 CðaÞ 0 CðaÞ j¼0 j;kþ1 j j¼0 1 f ðt j ; yj Þj þ ekþ1;kþ1 f ðt kþ1 ; y t kþ1 Þ f t kþ1 ; yPkþ1 2 2 2 2 CðaÞ d a 6 C 3 h 2 þ C 4 max jyðtj Þ yj j þ C 5 h jy tkþ1 ykþ1 j
jyðtkþ1 Þ ykþ1 j ¼ 2
2
2
06j6k
d2
6 C 3 h þ C 4 max jyðtj Þ yj j þ C 5 h
aþd1
06j6k
p
6 Ch :
2
minfd1 þa;d2 g
þ C 1 max jyðt j Þ yj j 6 C 1 h 06j6k
þ C 2 max jyðt j Þ yj j 06j6k
Theorem 11. Assume that the initial value problem solution y satisfies
Z tkþ1 k X d a1 a a 2 ðtkþ1 tÞ C D yðtÞdt dj;kþ1C D yðtj Þ 6 C 1 t mkþ1 1 h 1 ; 2 0 2 j¼0 Z tkþ1 k X 2 ðtkþ1 tÞa1 C Da yðtÞdt ej;kþ1C Da yðt j1 Þ ekþ1;kþ1C Da y t kþ1 2 2 2 0 j¼1 k X d fj;kþ1C Da yðt j Þ fkþ1;kþ1 C Da yðt kþ1 Þ 6 C 2 tmkþ2 1 h 2 ; 2 j¼0 with q
max jyðtj Þ yj j ¼ Oðh 1 Þ;
06j6N
and q
max jyðtj Þ yPj j ¼ Oðh 2 Þ;
06j6N
where m1, m2 P 0 and d1, d2 > 0. Then for n, T > 0, we have
p max yðtjþ1 Þ yjþ1 ¼ Oðh Þ; 2 2
06j6N1
where p = min(q1, q1 + a, q2 + a, d1 + a, d2) and N = [T/h]. Proof. Based on the above discussion, we find that
d yðt kþ1 Þ yPkþ1 6 C 1 h 1 þ C 2 max jyðt j Þ yj j: 2
06j6k
2
In the following, one gets
yðt kþ1 Þ ykþ1 ¼ 2
2
Z k X 1 tkþ12 ðtkþ1 tÞa1 f ðt; yðtÞÞdt fj;kþ1 f ðt j ; yj Þ fkþ1;kþ1 f ðtkþ1 ; yPkþ1 Þ 2 CðaÞ 0 j¼0 Z k X 1 tkþ12 P ej;kþ1 f ðtj1 ; yj1 Þ ekþ1;kþ1 f t kþ1 ; ykþ1 6 ðtkþ1 tÞa1 f ðt; yðtÞÞdt 2 2 2 2 2 CðaÞ 0 j¼1
k X j¼0
fj;kþ1 f ðtj ; yðt j ÞÞ fkþ1;kþ1 f ðt kþ1 ; yðt kþ1 ÞÞ
k X j¼1
1 ej;kþ1 f ðt j1 ; yðt j1 ÞÞ ekþ1;kþ1 f ðtkþ1 ; yðt kþ1 ÞÞ þ 2 2 2 2 CðaÞ
k 1 X 1 ej;kþ1 f ðtj1 ; yj1 Þ f ðt j1 ; yðtj1 ÞÞ þ fj;kþ1 jf ðt j ; yj Þ f ðtj ; yðt j ÞÞj þ fkþ1;kþ1 f tkþ1 ; yPkþ1 2 2 2 2 C ð a Þ C ð a Þ j¼1 j¼0 1 d f ðt kþ1 ; yðt kþ1 ÞÞj þ ekþ1;kþ1 f t jþ1 ; yPjþ1 f ðt jþ1 ; yðt jþ1 ÞÞ 6 C 1 h 2 þ C 2 max yðt jþ1 Þ yjþ1 2 2 2 2 2 2 06j6k CðaÞ a a d þ C 3 max jyðt j Þ yj j þ C 4 yðtkþ1 Þ yPkþ1 h þ C 5 yðt kþ1 Þ yPkþ1 h 6 C 1 h 2 þ C 2 max yðt jþ1 Þ yjþ1 2 2 2 2 06j6k 06j6k q aþq aþd aþq p p þ C 3 h 1 þ C 4 h 2 þ C 5 h 1 þ C 6 h 1 6 Ch þ C 2 max yðt jþ1 Þ yjþ1 6 Ch : k X
06j6k
2
2
3363
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368
Theorem 12. Assume that the initial value problem solution y satisfies
Z tkþ1 k X d a1 a a 1 ðt kþ1 tÞ C D yðtÞdt cj;kþ1C D yðt j Þ 6 C 1 tmkþ1 h 1; 0 j¼0 Z tkþ1 k k X X d a1 a a a 2 ðt kþ1 tÞ C D yðtÞdt aj;kþ1C D yðtj Þ bj;kþ1 D yðt jþ1 Þ 6 C 3 tmkþ1 h 2; 2 0 j¼0 j¼0 q max yðt jþ1 Þ yjþ1 ¼ Oðh Þ; 06j6k
2
2
where 0 6 j 6 N, m1, m2, m3 P 0 and d1, d2, d3 > 0 for n, T > 0. Then we have
yðt k Þ yP ¼ Oðhp Þ; k
where p = min{d1 + a, q, d2}. Proof. From Lemma 4, we have
d yðt kþ1 Þ e y Pkþ1 6 C 1 h 1 þ C 2 max jyðtj Þ yj j: 06j6k
Next we focus on the error for yPkþ1 . It follows
Z k k X X 1 tkþ1 a1 P e ðtkþ1 tÞ f ðt; yðtÞÞdt aj;kþ1 f ðt j ; yj Þ akþ1;kþ1 f ðt kþ1 ; y kþ1 Þ bj;kþ1 f ðtjþ1 ; yjþ1 Þ 2 2 CðaÞ 0 j¼0 j¼0 Z kþ1 k k t kþ1 X X 1 1 X 6 ðt kþ1 tÞa1 C Da yðtÞdt aj;kþ1 C Da yðt j Þ bj;kþ1C Da yðtjþ1 Þ þ a jf ðt j ; yðt j ÞÞ 2 CðaÞ j¼0 j;kþ1 CðaÞ 0 j¼0 j¼0
yðt kþ1 Þ yP ¼ kþ1
k 1 1 X y Pkþ1 þ bj;kþ1 f ðt jþ1 ; yðt jþ1 ÞÞ f ðtjþ1 ; yjþ1 Þ akþ1;kþ1 f ðtkþ1 ; yðtkþ1 ÞÞ f t kþ1 ; e 2 2 2 2 CðaÞ CðaÞ j¼0 a d 6 C 1 h 2 þ C 2 max jyðt j Þ yj j þ C 3 yðt kþ1 Þ e y Pkþ1 h þ C 4 max yðt jþ1 Þ yjþ1 f ðt j ; yj Þj þ
06j6k
2
06j6k
aþd1
d2
6 C 1 h þ C 2 max jyðt j Þ yj j þ C 3 h 06j6k
q
þ C 4 max jyðt j Þ yj j þ C 5 h 6 Ch 06j6k
2
minfd1 þa;q;d2 g
:
3.3. Stability analysis An interpretation of stability of a numerical scheme is that for small errors in the initial conditions cause small errors in the solution [29,23]. In the following, we firstly deal with the stability analysis of the new scheme (38), (39), (41). Theorem 13. Suppose yk+1 and y0kþ1 are numerical solutions in (38), respectively. Then
ðjÞ
0ðjÞ
and the initial conditions are given by y0 and y0
jykþ1 y0kþ1 j 6 Kky0 y00 k1 for any k, i.e. the new scheme (38), (39) and (41) are numerically stable. Sketch of the proof. According to the expression of (38), we easily get
jykþ1
y0kþ1 j
" # ½a X k k X tjkþ1 ðjÞ 1 X P ¼ a f ðt ; y Þ þ akþ1;kþ1 f ðtkþ1 ; ykþ1 Þ þ bj;kþ1 f ðtjþ1 ; ykþ1 Þ y þ 2 2 j¼0 j! 0 CðaÞ j¼0 j;kþ1 j j j¼0 " #! j ½ a k k Xt X 1 X kþ1 0ðjÞ 0 0P 0 a f ðt ; y Þ þ akþ1;kþ1 f ðtkþ1 ; ykþ1 Þ þ bj;kþ1 f t jþ1 ; ykþ1 y þ 2 2 j! 0 CðaÞ j¼0 j;kþ1 j j j¼0 j¼0 6
½a k X X tjkþ1 1 aj;kþ1 jf ðtj ; yj Þ f ðt j ; y0j Þj þ akþ1;kþ1 f ðt kþ1 ; yPkþ1 Þ f ðt kþ1 ; y0P ky0 y00 k1 þ kþ1 Þ j! C ð a Þ j¼0 j¼0 ! k X bj;kþ1 f ðtjþ1 ; ykþ1 Þ f ðt jþ1 ; y0kþ1 Þ 6 Kky0 y00 k1 þ C 1 max jyj y0j j þ C 2 jyPkþ1 y0P þ kþ1 j þ C 3 max jykþ1
j¼0 y0kþ1 j 2
2
2
2
2
6 K 1 ky0 y00 k1 þ K 2 max jyj y0j j 6 Kky0 y00 k1 :
06j6k
06j6k
Next, we deal with the stability analysis for the Improved Algorithms I–III respectively.
06j6k
2
3364
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368
Theorem 14. Suppose ykþ1 and y0kþ1 are numerical solutions in (43), 2 2 respectively. Then
ðjÞ
0ðjÞ
ðjÞ
0ðjÞ
and the initial conditions are given by y0 and y0
ykþ1 y0kþ1 6 Kky0 y00 k1 2
2
for any k, i.e. the Improved Algorithm I is numerically stable. Proof. The proof is similar to that of Theorem 13, so is omitted here. h By almost the same reasoning, we can get Theorem 15. Theorem 15. Suppose ykþ1 and y0kþ1 are numerical solutions in (48), 2 2 respectively. Then
and the initial conditions are given by y0 and y0
ykþ1 y0kþ1 6 Kky0 y00 k1 2
2
for any k, i.e. the Improved Algorithm II is numerically stable. ðjÞ
0ðjÞ
Theorem 16. Suppose yPkþ1 and y0P kþ1 are numerical solutions in (53), and the initial conditions are given by y0 and y0 respectively. Then
P y y0P 6 Kky y0 k 0 kþ1 kþ1 0 1 for any k, i.e. the Improved Algorithm III is numerically stable. Sketch of the proof. According to the expression of (53), we get " # ½a j k k X P X tkþ1 ðjÞ 1 X P y y0P ¼ e a f ðt ; y Þ þ akþ1;kþ1 f ðtkþ1 ; y kþ1 Þ þ bj;kþ1 f ðtjþ1 ; ykþ1 Þ y þ kþ1 kþ1 2 2 j¼0 j! 0 CðaÞ j¼0 j;kþ1 j j j¼0
( ½a j Xt
kþ1
j¼0
6
½a j X t
kþ1
j¼0
þ
k X j¼0
j!
j!
0ðjÞ
y0 þ
" #) k k X 1 X 0 aj;kþ1 f ðtj ; y0j Þ þ akþ1;kþ1 f ðtkþ1 ; e y 0P Þ þ b f ðt ; y Þ 1 1 j;kþ1 jþ2 kþ kþ1 2 CðaÞ j¼0 j¼0
ky0 y00 k1 þ
k X 1 a jf ðt ; y Þ f ðtj ; y0j Þj þ akþ1;kþ1 f tkþ1 ; e y Pkþ1 f tkþ1 ; e y 0P kþ1 CðaÞ j¼0 j;kþ1 j j
! 0 bj;kþ1 f ðtjþ1 ; ykþ1 Þ f ðtjþ1 ; y0kþ1 Þ 6 Kky0 y00 k1 þ C 1 max jyj y0j j þ C 2 j e y Pkþ1 e y 0P kþ1 j þ C 3 max ykþ1 ykþ1 : 2
2
2
06j6k
2
06j6k
2
2
Because e y Pkþ1 and e y 0P kþ1 are updated by (39), so
P e y
kþ1
0 0 e y 0P kþ1 6 Kky0 y0 k1 þ C max jyj yj j: 06j6k
Besides, according to Theorem 15, we have
jykþ1 y0kþ1 j 6 Kky0 y00 k1 : 2
2
Thus
P y y0P 6 Kky y0 k þ C 1 max jy y0 j þ C 2 Kky y0 k þ C max jy y0 j þ C 3 ky y0 k 0 j 0 j 0 kþ1 kþ1 0 1 j 0 1 j 0 1 06j6k
6 Kky0 y00 k1 þ C max jyj y0j j 6 Kky0 y00 k1 :
06j6k
06j6k
The sum of coefficients in our designed algorithm are bounded. Besides, in the process of above proof, we do not distinguish the constants C1–C9 and K1–K2 for convenience. 4. Numerical examples In order to verify our constructed numerical schemes for the fractional calculus and fractional differential equations, we give the following numerical examples. And the numerical results show that the numerical computational algorithms are efficient.
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368
3365
Example 1. Consider the function f(x) = x4, x 2 [0, 1]. Then kf ð2Þ k1 ¼ 12; kf ð4Þ k1 ¼ 24. The numerical results are shown in Tables 1–4. Example 2. Consider the function f(x) = x4 and f(x) = x6, x 2 [0, 1]. The numerical results are shown in Tables 5–8. Example 3.
( a
C D0;t yðtÞ
¼
2
Cð3aÞ t
2a
2a 2 Cð3aÞ t
yðtÞ þ t 2 t;
1a 1 Cð2aÞ t
a > 1; yðtÞ þ t t; a 6 1:
ð60Þ
2
The numerical results are shown in Tables 9–14. Example 4. 0:5 C D0;t yðtÞ
¼
2
Cð2:5Þ
t 1:5 yðtÞ þ t2 :
ð61Þ
The initial value we chosen as y(0) = 0. Its true solution is y(t) = t2. The numerical results are shown in Table 15.
Table 1 The absolute error of J a0;1 x4 for verifying Theorem 5. h
a = 0.8
a=1
a = 1.6
0.1 0.05 0.025 0.0125 0.00625
4.161861E003 1.051866E003 2.644498E004 6.631336E005 1.660625E005
3.330000E003 8.331250E004 2.083203E004 5.208252E005 1.302078E005
1.505989E003 3.745245E004 9.348082E005 2.335869E005 5.838784E006
Table 2 The absolute error of J a0;1 x4 for verifying Theorem 6. h
a = 0.8
a=1
a = 1.6
0.1 0.05 0.025 0.0125 0.00625
3.542874E006 2.223854E007 1.393338E008 8.711287E010 5.008000E011
3.333333E006 2.083333E007 1.302083E008 8.138021E010 5.086262E011
2.334466E006 1.457874E007 9.110534E009 5.961440E010 5.476738E010
Table 3 The absolute error of J a0;1 x4 for verifying Theorem 7 (on the first suggestion). h
a = 0.8
a=1
a = 1.6
0.1 0.05 0.025 0.0125 0.00625
3.552355E006 2.225207E007 1.393541E008 8.711598E010 5.008050E011
3.333333E006 2.083333E007 1.302083E008 8.138021E010 5.086262E011
2.301434E006 1.452862E007 9.102809E009 5.960241E010 5.476719E010
Table 4 The absolute error of J a0;1 x4 for verifying Theorem 7 (on the second suggestion). h
a = 0.8
a=1
a = 1.6
0.1 0.05 0.025 0.0125 0.00625
9.953847E005 2.881589E005 8.580664E006 2.563213E006 7.636424E007
4.399626E005 1.031004E005 2.536443E006 6.316066E007 1.577471E007
1.048313E005 3.435521E006 9.687469E007 2.600681E007 6.763113E008
3366
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368 Table 5 The absolute error of C Da0;1 x4 for verifying (29). h
a = 0.8
a=1
a = 1.6
0.1 0.05 0.025 0.0125 0.00625
1.055717E002 2.892040E003 7.770508E004 2.059158E004 5.400317E005
0 0 0 0 0
1.994496E002 5.143291E003 1.315600E003 3.345430E004 8.470497E005
Table 6 The absolute error of C Da0;1 x4 for verifying (30) (on the first suggestion). h
a = 0.8
a=1
a = 1.6
0.1 0.05 0.025 0.0125 0.00625
1.441004E007 3.973507E009 1.180234E010 2.718403E011 9.706680E011
0 0 0 0 0
2.055812E006 1.159297E007 6.887261E009 4.459828E010 2.237222E009
Table 7 The absolute error of C Da0;1 x4 for verifying (30) (on the second suggestion). h
a = 0.8
a=1
a = 1.6
0.1 0.05 0.025 0.0125 0.00625
8.270386E004 3.710286E004 1.647685E004 7.272723E005 3.195694E005
0 0 0 0 0
1.189139E003 5.834701E004 2.408542E004 9.464361E005 3.653170E005
Table 8 The absolute error of C Da0;1 x6 for verifying (31). h
a = 0.8
a=1
a = 1.6
0.1 0.05 0.025 0.0125 0.00625
4.945389E005 3.431437E006 2.327757E007 1.554474E008 1.084208E009
0 0 0 0 0
9.759985E005 6.328401E006 4.063441E007 2.582014E008 2.267807E009
Table 9 Errors for Eq. (60) with t = 1, a < 1. h
a = 0.3
a = 0.5
a = 0.7
0.1 0.05 0.025 0.0125 0.00625 0.003125
0.0225 0.0124 0.0061 0.0029 0.0014 6.6099E004
0.0041 0.0031 0.0018 0.0010 5.2963E004 2.7256E004
0.0059 0.0018 4.9076E004 8.2064E005 2.2522E005 3.6356E005
Table 10 Errors for Eq. (60) with t = 1, a > 1. h
a = 1.25
a = 1.5
a = 1.85
0.1 0.05 0.025 0.0125 0.00625 0.003125
0.0255 0.0127 0.0063 0.0031 0.0016 7.8007E004
0.0200 0.0101 0.005 0.0025 0.0012 6.1623E004
0.0042 0.0031 0.0018 0.001 5.4928E004 2.8978E004
3367
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368 Table 11 (Improved Algorithm I) Errors for Eq. (60) with t = 1, a < 1. h
a = 0.3
a = 0.5
a = 0.7
0.1 0.05 0.025 0.0125 0.00625 0.003125
0.0021 4.5357E004 8.6128E005 1.1965E005 3.33131E007 1.2418E006
0.0018 7.1792E004 2.7651E004 1.0326E004 3.7853E005 1.3714E005
0.0066 0.0031 0.0014 5.7564E004 2.4134E004 1.002E004
Table 12 (Improved Algorithm I) Errors for Eq. (60) with t = 1, a > 1. h
a = 1.25
a = 1.5
a = 1.85
0.1 0.05 0.025 0.0125 0.00625 0.003125
0.0014 4.4287E004 1.3188E004 4.3188E005 1.3340E005 4.968E006
0.0018 7.1792E004 2.7651E004 1.0326E004 3.7853E005 1.3714E005
0.0073 0.0031 0.0013 5.9123E004 2.6345E004 1.18E004
Table 13 (Improved Algorithm II) Errors for Eq. (60) with t = 1, a < 1. h
a = 0.3
a = 0.5
a = 0.7
0.1 0.05 0.025 0.0125 0.00625 0.003125
0.0161 0.0059 0.0021 7.5212E004 2.7028E004 9.8093E005
0.0036 0.0011 3.402E004 1.0651E004 3.4101E005 1.1161E005
9.2062E005 1.1114E004 6.5809E005 3.2349E005 1.4721E005 6.4434E006
Table 14 (Improved Algorithm II) Errors for Eq. (60) with t = 1, a > 1. h
a = 1.25
a = 1.5
a = 1.85
0.1 0.05 0.025 0.0125 0.00625 0.003125
3.7015E004 9.3442E005 2.4233E005 6.4536E006 1.7602E006 4.8997E007
0.0015 5.2098E004 1.8259E004 6.4261E005 2.22667E005 8.0045E006
0.0130 0.0058 0.0026 0.0012 5.3348E004 2.4039E004
Table 15 Errors for Eq. (61) with t = 10. h
Improved Algorithm II
Improved Algorithm III
0.1 0.05 0.025 0.0125
0.0361 0.0117 0.0039 0.0013
0.0032 6.7576E004 1.5012E004 3.4362E005
References [1] K.B. Oldham, J. Spanier, The Fractional Calculus: Theory and Applications of Differentiation and Integration to Arbitrary Order, Academic Press, New York, 1974. [2] I. Podlubny, Fractional Differential Equations, Academic Press, 1999. [3] S.G. Samko, A.A. Kilbas, O.I. Marichev, Fractional Integrals and Derivatives Theory and Applications, Gordon and Breach, New York, 1993. [4] A.A. Kilbas, H.M. Srivastava, J.J. Trujillo, Theory and Application of Fractional Differential Equations, Elsevier, 2006. [5] Z. Odibat, Approximations of fractional integral and Caputo fractional derivatives, Appl. Math. Comput. 178 (2006) 527–533. [6] K. Diethelm, N.J. Ford, A.D. Freed, Y. Luchko, Algorithms for the fractional calculus: a selection of numerical methods, Comput. Methods Appl. Mech. Eng. 194 (6–8) (2005) 743–773. [7] C. Lubich, Discretized fractional calculus, SIAM J. Math. Anal. 17 (3) (1986) 704–719. [8] W. Deng, Short memory principle and a predictor–corrector approach for fractional differential equations, J. Comput. Appl. Math. 206 (2007) 174–188. [9] L. Yuan, O.P. Agrawal, A numerical scheme for dynamic systems containing fractional derivatives, ASME J. Vibr. Acoust. 124 (2002) 321–324. [10] K. Diethelm, An improvement of a nonclassical numerical method for the computation of fractional derivatives, Numer. Algorithms 47 (2008) 190– 361.
3368
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368
[11] D.A. Murio, On the stable numerical evaluation of Caputo fractional derivatives, Comput. Math. Appl. 51 (2006) 1539–1550. [12] T. Miyakoda, Discretized fractional calculus with a series of chebyshev polynomial, Electron. Notes Theor. Comput. Sci. 225 (2009) 239–244. [13] K. Dithelm, N.J. Ford, A.D. Freed, A predictor–corrector approach for the numerical solution of fractional differential equations, Nonlinear Dyn. 29 (2002) 3–22. [14] K. Dithelm, N.J. Ford, A.D. Freed, Detailed error analysis for a fractional adams method, Numer. Algorithms 36 (2004) 31–52. [15] C.P. Li, C.X. Tao, On the fractional adams method, Comput. Math. Appl. 58 (8) (2009) 1573–1588. [16] C.P. Li, Y.H. Wang, Numerical algorithm based on adomian decomposition for fractional differential equations, Comput. Math. Appl. 57 (2009) 1672– 1681. [17] O.P. Agrawal, A general finite element formulation for fractional variational problems, J. Math. Anal. Appl. 337 (2008) 1–12. [18] T. Tang, A finite difference scheme for partial integro-differential equations with weakly singular kernel, Appl. Numer. Math. 11 (1993) 309–319. [19] Q. Sheng, T. Tang, Optimal convergence of an euler and finite difference method for nonlinear partial integro-differential equations, Math. Comput. Model. 21 (1995) 1–11. [20] I. Podlubny, Matrix approach to discrete fractional calculus, Int. J. Theor. Appl. 3 (2000) 359–386. [21] I. Podlubny, A. Chechkin, T. Skovranek, Y. Chen, B.M.V. Jara, Matrix approach to discrete fractional calculus. II: Partial fractional differential equations, J. Comput. Phys. 228 (2009) 3137–3153. [22] Q. Yang, F. Liu, I. Turner, Numerical methods for fractional partial differential equations with riesz space fractional derivatives, Appl. Math. Model. 34 (2010) 200–218. [23] R. Lin, F. Liu, Fractional high order methods for the nonlinear fractional ordinary differential equation, Nonlinear Anal. 66 (2007) 856–869. [24] X. Li, C. Xu, Existence and uniqueness of the weak solution of the space–time fractional diffusion equation and a spectral method approximation, Commun. Comput. Phys. 8 (2010) 1016–1051. [25] M. Xi, Numerical Analysis Method, University of Science and Technology of China Press, 2004. [26] Z. Odibat, Computational algorithms for computing the fractional derivatives of functions, Math. Comput. Simul. 79 (7) (2009) 2013–2020. [27] C. Gerald, P. Wheatley, Applied Numerical Analysis, Addison-Wesley, USA, 2004. [28] K. Diethelm, An algorithm for the numerical solution of differential equations of fractional order, Electron. Trans. Numer. Anal. 5 (1997) 1–6. [29] J.W. Thomas, Numerical Partial Differential Equations: Finite Difference Methods, Springer, 1995.
Contents lists available at ScienceDirect
Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp
Numerical approaches to fractional calculus and fractional ordinary differential equation q Changpin Li ⇑, An Chen, Junjie Ye Department of Mathematics, Shanghai University, Shanghai 200444, China
a r t i c l e
i n f o
Article history: Received 19 July 2010 Received in revised form 10 November 2010 Accepted 22 January 2011 Available online 31 January 2011 Keywords: Numerical approach Fractional calculus Fractional differential equations Piecewise interpolation Simpson method
a b s t r a c t Nowadays, fractional calculus are used to model various different phenomena in nature, but due to the non-local property of the fractional derivative, it still remains a lot of improvements in the present numerical approaches. In this paper, some new numerical approaches based on piecewise interpolation for fractional calculus, and some new improved approaches based on the Simpson method for the fractional differential equations are proposed. We use higher order piecewise interpolation polynomial to approximate the fractional integral and fractional derivatives, and use the Simpson method to design a higher order algorithm for the fractional differential equations. Error analyses and stability analyses are also given, and the numerical results show that these constructed numerical approaches are efficient. Ó 2011 Elsevier Inc. All rights reserved.
1. Introduction The history of fractional calculus [1–4] is more than 300 years old, but only in the recent decades did the applied scientists and the engineers realize that such fractional differential equations provided a better approach to describe the complex phenomena in nature, such as non-Brownian motion, signal processing, systems identification, control, viscoelastic materials and polymers. Because of the non-local property of the fractional derivative, it can be used to describe those complex systems which involve long-memory in time in a better way. Based on these requirements, the numerical approach has become a very desired method to analysis the experimental data which are described in a fractional way. Some papers have discussed the numerical approaches for fractional calculus and fractional differential equations. Here let us review the previous and present numerical methods. In [1], Oldham and Spanier introduced some numerical approaches to deal with the fractional calculus. They used the definitions of the Grünwald–Letnikov and Riemann–Liouville derivatives respectively. The main ideas of those numerical algorithms are that they used a Riemann sum to approximate the fractional derivative. Later on more researches are based on the above idea and also do certain theoretical analysis, such as [5,6]. In [2], Podlubny introduced the numerical approach for the arbitrary order derivative by using the definition of Riemann–Liouville based on the relationship between Grünwald–Letnikov derivative and the Riemann–Liouville derivative. The convergence order he got is O(h). Also he suggested that the short memory principle can be used to reduce the complexity at the expense of lacking certain precision. And he also pointed out that the convergence order can be improved to a more higher order by using the method of Lubich [7]. In [8], Deng described how to use the short memory principle combined with predictor–corrector approach in fractional differential equations in details. q The work was partially supported by the Natural Science Foundation of China under Grant No. 10872119 and Shanghai Leading Academic Discipline under Grant No. S30104. ⇑ Corresponding author. E-mail address: [email protected] (C. Li).
0021-9991/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.jcp.2011.01.030
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368
3353
In [9], Yuan and Agrawal designed a novel numerical approach which can avoid the harmful effect of non-local property in fractional derivative. This method was improved by Diethelm [10]. The main idea is that it transforms the original problem into an integral and then evaluates the numerical solution by Gauss integral formula. In [11], Murio designed a stable numerical algorithm to the fractional derivative in Caputo sense based on mollification techniques. Using a series of Chebyshev polynomials to evaluate the fractional calculus was posed in [12] and it seems that the method is highly effective. However that paper only dealt with the case 0 < a < 1. For the fractional differential equations, there are also some papers dealing with the numerical solutions. Some numerical algorithms such as predictor–corrector method [13,14,8,15], the Adomian decomposition [16], and the finite element [17] have been studied. Early in 1993, Tang studied the numerical scheme of the partial integro-differential equations in which the time derivative term was approximated by a Crank–Nicolson time-stepping and the integral term was treated by the product trapezoidal method [18]. The numerical scheme by the finite difference approach was further studied and an optimal convergence rate was also obtained, see [19]. Podlubny in [20,21] used the matrix expression to unify the formula where the fractional derivatives are derived from a finite difference. Liu also did some interesting works on the numerical approach for the fractional differential equations, see [22,23]. And Li and Xu in [24] constructed an efficient spectral method for the numerical approximation of the space–time fractional diffusion equation based on their proposed weak solution formulation, which is a significant difference to the direct discrete way. Based on the above review, we can see that there is not enough efficient numerical approach for the application of technology. Therefore in this paper, we pose the novel algorithms for a better convergence. The paper is organized as follows: In the present section we design some efficient approaches for the fractional calculus. In Section 2, we introduce the treatment for the fractional differential equations, where the error analyses and stability analyses are also given. Finally, we will give some numerical examples. First of all, let us introduce the definitions of fractional integral and derivatives. Definition 1. The definition of fractional integral of a function y reads as
J aa;t yðtÞ ¼
1 CðaÞ
Z
t
ðt sÞa1 yðsÞds;
a > 0:
ð1Þ
a
Definition 2. The Riemann–Liouville fractional derivative of y is defined as n
a
RL Da;t yðtÞ
¼
1 d Cðn aÞ dxn
Z
t
ðt sÞna1 yðsÞds;
t > a; n 1 < a < n 2 Z þ :
ð2Þ
a
Definition 3. The Grünwald–Letnikov fractional derivative of y is defined as
( ) N1 ðta Þa X Cðj aÞ ta N ; y tj GL Da;t yðtÞ ¼ lim N!1 N CðaÞ j¼0 Cðj þ 1Þ a
t > a; n 1 < a < n 2 Z þ :
ð3Þ
Definition 4. The Caputo fractional derivative of y is defined as a C Da;t yðtÞ ¼
1 Cðn aÞ
Z
t
ðt sÞna1 yðnÞ ðsÞds;
t > a; n 1 < a < n 2 Z þ :
ð4Þ
a
2. Numerical approach for the fractional calculus In this section, we divide the given interval [0, b] into:
D : 0 ¼ t0 < t1 < < t n ¼ t ¼ b: And ti = ih, in which h ¼
t ; n
ð5Þ
i ¼ 0; 1; . . . ; n.
2.1. Fractional integral Firstly, we use the piecewise linear interpolation function to approximate y,
yn ðtÞ ¼
t t iþ1 t ti yðt i Þ þ yðt iþ1 Þ; ti tiþ1 tiþ1 t i
where ti 6 t 6 ti+1, i = 0, 1, 2, . . . , n 1.
ð6Þ
3354
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368
So substitution in (1) with (6) gives, a
J 0;t yðtÞ
n1 X i¼0
¼
n1 X
a
J 0;tn yn ðtÞ ¼
¼
t iþ1
ti
) ðt n sÞa1 yn ðsÞds CðaÞ ! ! ) Z tiþ1 s tiþ1 ðt n sÞa1 s t i dt yðt i Þ þ dt yðtiþ1 Þ t i t iþ1 CðaÞ tiþ1 t i ti
t iþ1 ti
i¼0
( Z
i¼0
(Z n1 X
a1
ðt n sÞ CðaÞ
a n X h W yðt Þ; Cða þ 2Þ i¼0 1;i i
ð7Þ
where
W 1;0 ¼ na ða þ 1 nÞ þ ðn 1Þaþ1 ; W 1;i ¼ ðn i 1Þ
aþ1
þ ðn i þ 1Þ
aþ1
ð8Þ aþ1
2ðn iÞ
;
i ¼ 1; 2; . . . ; n 1;
ð9Þ
W 1;n ¼ 1:
ð10Þ
Thus we get the following numerical approach for fractional integral. Theorem 5. Suppose that y(t) 2 C2[0, b], the fractional integral has the following numerical approach:
J a0;tn yðtÞ ¼
a n X h 2 W yðt Þ þ Oðh Þ; Cða þ 2Þ i¼0 1;i i
ð11Þ
where W1,i are defined as (8)–(10). Proof. Let
EðtÞ ¼ yðtÞ yn ðtÞ ¼
yð2Þ ðtÞ ðt t i Þðt t iþ1 Þ; 2!
where t 2 (ti,ti+1). Thus
kJ a0;tn yn ðtÞ J a0;tn yðtÞk1 ¼ kJ a0;tn ðyn ðtÞ yðtÞÞk1 6 kyn ðtÞ yðtÞk1
1 CðaÞ
Z
tn
ðt n sÞa1 ds ¼
0
t an ky ðtÞ yðtÞk1 Cða þ 1Þ n
a
tan b kyð2Þ k1 2 2 ¼ kyn ðtÞ yðtÞk1 ¼ Oðh Þ ¼ Oðh Þ: Cða þ 1Þ Cða þ 1Þ
Next, we apply some much higher convergence order for y. So the piecewise cubic interpolation function will be our choice. We have the piecewise cubic interpolation function as follows:
yn ðtÞ ¼
2 2 2 t ti t tiþ1 t t iþ1 t ti t tiþ1 12 yðti Þ þ 1 2 yðtiþ1 Þ þ ðt ti Þ y0 ðti Þ t i t iþ1 t i t iþ1 t iþ1 ti t iþ1 ti ti tiþ1 þ ðt t iþ1 Þ
t ti tiþ1 t i
2
y0 ðt iþ1 Þ:
ð12Þ
Substituting (12) into (1) yields a
J 0;t yðtÞ
n1 X
a
J 0;tn yn ðtÞ ¼
i¼0
þ
i¼0
Z
tiþ1
ti
þ
n1 X
Z ti
t iþ1
( Z ti
t iþ1
! ðt n sÞa1 s ti s tiþ1 2 ds yðt i Þ 12 CðaÞ t i t iþ1 t i t iþ1
! ! Z tiþ1 ðtn sÞa1 s tiþ1 s ti 2 ðt n sÞa1 s tiþ1 2 ds yðtiþ1 Þ þ ds y0 ðti Þ 12 ðs ti Þ CðaÞ tiþ1 t i tiþ1 t i CðaÞ ti tiþ1 ti ! ) a aþ1 n n X X ðtn sÞa1 s ti 2 h h 0 f y0 ðt Þ; W ds y ðt iþ1 Þ ¼ W 2;i yðt i Þ þ ðs tiþ1 Þ CðaÞ t iþ1 ti Cða þ 4Þ i¼0 Cða þ 4Þ i¼0 2;i i ð13Þ
3355
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368
where
W 2;0 ¼ 6ðn 1Þ2þa ð1 þ 2n þ aÞ þ na 12n3 6n2 ð3 þ aÞ þ ð1 þ aÞð2 þ aÞð3 þ aÞ ; W 2;i ¼ 6 4ðn iÞ3þa þ ðn i 1Þ2þa ð1 þ 2i 2n aÞ þ ð1 i þ nÞ2þa ð1 þ 2i 2n þ aÞ ;
ð14Þ i ¼ 1; 2; . . . ; n 1;
W 2;n ¼ 6ð1 þ aÞ; f 2;0 ¼ 2ðn 1Þ2þa ð3n þ aÞ þ n1þa 6n2 4nð3 þ aÞ þ ð2 þ aÞð3 þ aÞ ; W
ð15Þ ð16Þ ð17Þ
f 2;i ¼ 2ðn i 1Þaþ2 ð3i 3n aÞ 2ðn i þ 1Þaþ2 ð3i 3n þ aÞ ðn iÞaþ2 ð24 þ 8aÞ; W f 2;n ¼ 2a: W
i ¼ 1; 2; . . . ; n 1;
ð18Þ ð19Þ
So we can also get the other following numerical approach for fractional integral. Theorem 6. Suppose y(t) 2 C4[0, b], the fractional integral has the following numerical approach:
J a0;tn yðtÞ ¼
a aþ1 n n X X h h f y0 ðt Þ þ Oðh4 Þ; W W 2;i yðt i Þ þ Cða þ 4Þ i¼0 Cða þ 4Þ i¼0 2;i i
ð20Þ
f 2;i are defined as (14)–(19). where W2,i and W Proof. The numerical formula (20) is obvious. In effect, set
EðtÞ ¼ yðtÞ yn ðtÞ ¼
yð4Þ ðtÞ ðt ti Þ2 ðt t iþ1 Þ2 ; 4!
ð21Þ
where t 2 (ti,ti+1). Thus
kJ a0;tn yn ðtÞ J a0;tn yðtÞk1 ¼ kJ a0;tn ðyn ðtÞ yðtÞÞk1 6 kyn ðtÞ yðtÞk1
1 CðaÞ
Z
tn
ðt n sÞa1 ds ¼
0
t an ky ðtÞ yðtÞk1 Cða þ 1Þ n
a
tan b kyð4Þ k1 4 4 ¼ kyn ðtÞ yðtÞk1 ¼ Oðh Þ ¼ Oðh Þ: Cða þ 1Þ Cða þ 1Þ
In (20), the first-order derivative can be approximated by the back forward difference or the central difference. Since the spline function has good properties to fit the experimental data, in the following we use the spline function to approximate y. The cubic interpolation spline function reads as follows:
2 2 2 t ti t tiþ1 t tiþ1 t ti t t iþ1 yn ðtÞ ¼ 1 2 yðt i Þ þ 1 2 yðtiþ1 Þ þ ðt ti Þ mi ti tiþ1 ti tiþ1 t iþ1 ti t iþ1 ti ti tiþ1 2 t ti þ ðt tiþ1 Þ miþ1 ; t iþ1 t i
ð22Þ
where mi (i = 0, 1, 2, . . . , n 1) are the undetermined real numbers. Then one has
J a0;tn yðtÞ
Z
tn
0
ðt n sÞa1 yn ðsÞds: CðaÞ
ð23Þ
That is, a
J 0;tn yðtÞ
n1 Z X i¼0
þ þ
Z Z
t iþ1 ti tiþ1
ti
ti
t iþ1
! ( Z n1 t iþ1 X ðt n sÞa1 ðt n sÞa1 s ti s tiþ1 2 ds yðti Þ yn ðsÞds ¼ 12 CðaÞ CðaÞ ti t iþ1 t i t iþ1 ti i¼0 ! ! Z tiþ1 ðtn sÞa1 s tiþ1 s ti 2 ðt n sÞa1 s tiþ1 2 ds yðt iþ1 Þ þ ds mi 12 ðs ti Þ CðaÞ tiþ1 t i tiþ1 t i CðaÞ ti tiþ1 ti ) 2 ! a aþ1 a1 n n X X ðtn sÞ s ti h h f m; W ds miþ1 ¼ W 3;i yðt i Þ þ ð24Þ ðs t iþ1 Þ CðaÞ t iþ1 ti Cða þ 4Þ i¼0 Cða þ 4Þ i¼0 3;i i
f 3;i ¼ W f 2;i ; i ¼ 0; 1; 2; . . . ; n. where W 3;i ¼ W 2;i ; W For the undetermined real number mi, i = 1, 2, . . . , n 1, the expressions are as follows:
ki mi1 þ 2mi þ li miþ1 ¼ C i ;
i ¼ 1; 2; . . . ; n 1;
in which
t iþ1 t i ; t iþ1 ti1 t t li ¼ 1 ki ¼ i i1 ; t iþ1 ti1 C i ¼ 3½ki yðti1 ; ti Þ þ li yðti ; t iþ1 Þ: ki ¼
ð25Þ
3356
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368
The above formulas have only n 1 linear equations, but the unknown numbers are n + 1. In order to determine the unknown mi, here we present two ways: 1. Taking first-order derivative of spline interpolation function at boundary as m0 ¼ fn0 ðx0 Þ; mn ¼ fn0 ðxn Þ, then we get the following m1,i:
m1;i
8 0 i ¼ 0; > < f ðx0 Þ; ¼ f 0 ðxn Þ; i ¼ n; > : the solution mi of ð25Þ; i ¼ 1; 2; . . . ; n 1:
ð26Þ
In (25) mi is substituted by m1,i, C1 and Cn1 are updated by C1 k1m1,0 and Cn1 ln1m1,n respectively. Under these conditions, we get a system of Eq. (25) whose coefficient matrix is a triangular one. So the method of forward elimination and backward substitution can be used. 2. Taking fn00 ðx0 Þ ¼ f 00 ðx0 Þ; fn00 ðxn Þ ¼ f 00 ðxn Þ, then we get the following equations with respect to m1,i:
8 x1 x0 00 > < 2m1;0 þ m1;1 ¼ 3f ðx0 ; x1 Þ 2 f ðx0 Þ; ki m1;i1 þ 2m1;i þ li m1;iþ1 ¼ C i ; i ¼ 1; 2; . . . ; n 1; > : m1;n1 þ 2m1;n ¼ 3f ðxn1 ; xn Þ þ xn x2 n1 f 00 ðxn Þ:
ð27Þ
Thus we also get a system of equations whose coefficient matrix is also a triangular one. Thus, we have another numerical approach for fractional integral. Theorem 7. Suppose y(t) 2 C2[0, b]. Then the fractional integral has the following numerical approach:
J a0;tn yðtÞ ¼
a aþ1 n n X X h h f m þ Oðh4 Þ; W W 3;i yðt i Þ þ Cða þ 4Þ i¼0 Cða þ 4Þ i¼0 3;i 1;i
ð28Þ
f 3;i are defined as (14)–(19), and m1,i can be determined by the above suggestion (26) or (27). where W3,i and W
Proof. According to [25], we have
jyn ðtÞ yðtÞj 6
7 4 Lh ; 8
where L is a constant. Thus
kJ a0;tn yn ðtÞ J a0;tn yðtÞk1 ¼ kJ a0;tn ðyn ðtÞ yðtÞÞk1 6 kyn ðtÞ yðtÞk1
1 CðaÞ
Z
tn
ðt n sÞa1 ds ¼
0
t an ky ðtÞ yðtÞk1 Cða þ 1Þ n
a
tan 7b Lkyð4Þ k1 4 4 ¼ kyn ðtÞ yðtÞk1 ¼ Oðh Þ ¼ Oðh Þ: Cða þ 1Þ 8Cða þ 1Þ
2.2. Fractional derivative Due to a wide application of the Caputo derivative, we put our attention to the numerical approach of it. Oldham and Spanier [1] have done some researches on the numerical approach to Grünwald–Letnikov derivative and Riemann–Liouville derivative. Besides, Liu et al. used the finite difference to approximate Riemann–Liouville derivative based on the relationship between Grünwald–Letnikov derivative and Riemann–Liouville derivative. It shows that the approaches are very efficient. Here we approximate the Caputo derivative which is similar to the above discussions for the fractional integral. In the following, in order to compute the numerical solution y(m), we always assume that it makes sense in the small neighborhood of the two endpoints. Theorem 8. When y(t) 2 C2+m[0, b], m = [a] + 1, then the Caputo derivative has the following numerical approach: a
¼
ma n X h 2 W yðmÞ ðt i Þ þ Oðh Þ; Cðm a þ 2Þ i¼0 5;i
ð29Þ
a
¼
ma maþ1 n n X X h h f m þ Oðh4 Þ: W W 6;i yðmÞ ðt i Þ þ Cðm a þ 4Þ i¼0 Cðm a þ 4Þ i¼0 6;i 2;i
ð30Þ
C D0;t n yðtÞ
C D0;t n yðtÞ
And when y(t) 2 C4+m[0, b], m = [a] + 1, then the Caputo derivative has the following numerical approach: a
C D0;t n yðtÞ
¼
ma maþ1 n n X X h h f yðmþ1Þ ðt i Þ þ Oðh4 Þ; W W 7;i yðmÞ ðt i Þ þ Cðm a þ 4Þ i¼0 Cðm a þ 4Þ i¼0 7;i
ð31Þ
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368
3357
f 6;i are defined as follows: where W5,i, W6,i and W
8 > nma ðm a þ 1 nÞ þ ðn 1Þmaþ1 ; i ¼ 0; > > > < ðn i 1Þmaþ1 þ ðn i þ 1Þmaþ1 ; W 5;i ¼ > > 2ðn iÞmaþ1 ; i ¼ 1; . . . ; n 1; > > : 1; i ¼ n; 8 6ðn 1Þ2þma ð1 þ 2n þ m aÞ þ nma ð12n3 6n2 ð3 þ m aÞ > > > > > > < þð1 þ m aÞð2 þ m aÞð3 þ m aÞÞ; i ¼ 0; W 6;i ¼ 6ð4ðn iÞ3þma þ ðn i 1Þ2þma ð1 þ 2i 2n m þ aÞ > > > > þð1 i þ nÞ2þma ð1 þ 2i 2n þ m aÞÞ; i ¼ 1; . . . ; n 1; > > : 6ð1 þ m aÞ; i ¼ n; 8 2ðn 1Þ2þma ð3n þ m aÞ > > > > > þn1þma 6n2 4nð3 þ m aÞ þ ð2 þ m aÞð3 þ m aÞ ; i ¼ 0; > > > > maþ2 < ð3i 3n m þ aÞ f 6;i ¼ 2ðn i 1Þ W maþ2 > > ð3i 3n þ m aÞ 2ðn i þ 1Þ > > > > maþ2 > ð24 þ 8ðm aÞÞ; i ¼ 1; 2; . . . ; n 1; ðn iÞ > > : 2a; i ¼ n:
As to (29), y(m)(ti) have been given by [26,27]
8 yðt þhÞyðt hÞ 2 i i þ Oðh Þ; > 2h > > > 2 > yðt i þhÞ2yðt i Þþyðt i hÞ > þ Oðh Þ; > > h2 > > > 2n > P < 2 ð2nÞ! 1 ð1Þk k!ð2nkÞ! yðt i þ ðn kÞhÞ þ Oðh Þ; ym ðti Þ :¼ um ðhÞ ¼ h2n k¼0 > > > 2nþ1 > P > ð2nþ1Þ! 1 > > 2nþ1 ð1Þk k!ð2nþ1kÞ! > h > > k¼0 > > : 2 12 ðyðt i þ ðn kÞhÞ þ yðt i þ ðn þ 1 kÞhÞÞ þ Oðh Þ; And combining the above formula with Richardson extrapolation, y
m ¼ 1; m ¼ 2; m ¼ 2j;
m ¼ 2j þ 1; j ¼ 1; 2; . . .
(m)
(ti) to (30) and (31) are obtained as follows:
4 h 1 4 um ðhÞ þ Oðh Þ: yðmÞ ðt i Þ ¼ um 3 2 3 f 7;i ¼ W f 6;i , and m2,i are those values m1,i given in Theorem 7 in which y(t) is substituted by y(m)(t). And W7,i = W6,i, W Proof. The proof is simple. Once combining the definition of Caputo derivative and the already designed numerical approaches for the fractional integral, the remain thing is only the substitution. h
3. Numerical approach for the fractional differential equations In this section, we discuss a new numerical scheme for solving the following system: a
C D0;t yðtÞ
¼ f ðt; yðtÞÞ;
ðkÞ
yðkÞ ð0Þ ¼ y0 ;
k ¼ 0; . . . ; ½a:
ð32Þ
It is well known that the initial value problem (32) is equivalent to a Volterra integral equation
yðtÞ ¼
½a X
m¼0
ðmÞ
y0
tm 1 þ 2 tt ðt uÞa1 f ðu; yðuÞÞdu m! CðaÞ 0
ð33Þ
in the sense that a continuous function is a solution of (32) if and only if it is also a solution of (33). Here, we indicate the approach briefly. The main problem is to solve the integral on the right-hand side of (33) by numerical method. For the consequence of the non-local structure of the fractional-order differential operators, we use the compound quadrature formula with respect to the weight function (tk+1 )a1 to replace the integral naturally. The product trapezoidal quadrature formula is a good choice. The scheme constructed by using the trapezoidal quadrature formula to replace the integral is proved to work [28,6]. We may also construct the higher order of scheme by using the compound Simpson formula. We need not only nodes tj (j = 0,1,. . .,k + 1) but also t jþ1 ðj ¼ 0; 1; . . . ; kÞ. In other words, we apply the approximation 2
3358
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368
Z
t kþ1
ðtkþ1 zÞa1 gðzÞdz ¼
0
Z
tkþ1
ðt kþ1 zÞa1 g~kþ1 ðzÞdz;
ð34Þ
0
where g~kþ1 is the piecewise quadratic interpolation for g with nodes and knots chosen at the tj ; t jþ1 ðj ¼ 0; 1; 2; . . . ; k þ 1Þ. 2 Using the standard techniques of the Simpson theory, we write the integral on the right-hand side of (34) as
Z
t kþ1
ðtkþ1 zÞa1 g~kþ1 ðzÞdz ¼
0
kþ1 X
aj;kþ1 gðt j Þ þ
j¼0
k X
bj;kþ1 gðt jþ1 Þ;
ð35Þ
2
j¼0
where
aj;kþ1
h i 8 aþ2 aþ2 4ha > ðk þ 2Þ ðk þ 1Þ > a ð a þ1Þð a þ2Þ > > > > > h i > a a > > aðahþ1Þ ðk þ 1Þaþ1 3ðk þ 2Þaþ1 þ ha ðk þ 2Þa ; j ¼ 0; > > > > > < h i a aþ2 4ha ¼ ðk jÞaþ2 aðahþ1Þ aðaþ1Þðaþ2Þ ðk þ 2 jÞ > > > > > h i > > > ðk þ 2 jÞaþ1 þ 6ðk þ 1 jÞaþ1 þ ðk jÞaþ1 ; 1 6 j 6 k; > > > > > > > : 2ha a aðaþ1Þðaþ2Þ ð1 2Þ; j ¼ k þ 1; h i h i 8h 4h ðk þ 2 jÞaþ2 ðk þ 1 jÞaþ2 ðk þ 2 jÞaþ1 þ ðk þ 1 jÞaþ1 ; aða þ 1Þða þ 2Þ aða þ 1Þ a
bj;kþ1 ¼
ð36Þ
a
0 6 j 6 k: ð37Þ
In order to construct the explicit scheme avoiding iterations. We replace the integral on the right-hand side of Eq. (33) by the product rectangle rule to calculate the predictor of yk+1, then
ykþ1
" # ½a k k X X t jkþ1 ðjÞ 1 X P ¼ a f ðt ; y Þ þ akþ1;kþ1 f ðt kþ1 ; ykþ1 Þ þ bj;kþ1 f ðt jþ1 ; ykþ1 Þ ; y þ 2 2 j! 0 CðaÞ j¼0 j;kþ1 j j j¼0 j¼0
yPkþ1 ¼
ð38Þ
½a k X t jkþ1 ðjÞ 1 X c f ðt ; y Þ; y0 þ j! CðaÞ j¼0 j;kþ1 j j j¼0
ð39Þ
where aj,k+1, bj,k+1 are shown in (36) and (37) respectively, and
h
cj;kþ1 ¼
a
a
½ðk þ 1 jÞa ðk jÞa :
ð40Þ
The remaining problem is to determine the values ykþ1 . We use the product rectangle rule to calculate the ykþ1 . Thus 2
ykþ1 ¼ 2
j ½a t X kþ1 2
j¼0
j!
ðjÞ
y0 þ
2
k 1 X d f ðt ; y Þ; CðaÞ j¼0 j;kþ1 j j
ð41Þ
where
dj;kþ1 ¼
h
a
kþ
a
1 j 2
a
a 1 ; k j 2
0 6 j 6 k:
ð42Þ
Thus a new scheme is described by (38), (39) and (41) with the weights aj,k+1, bj,k+1, cj,k+1, dj,k+1 defined as above, respectively. 3.1. Improved algorithms In this section, we want to improve the above designed algorithm. Consider the predictor–corrector format to the intermediate nodes, we have
ykþ1 ¼ 2
yPkþ1 ¼ 2
j ½a t X kþ1 2
j¼0
j!
2
j!
þ
" k X
# P ej;kþ1 f ðt j ; yj Þ þ ekþ1;kþ1 f t kþ1 ; ykþ1 ;
j¼0
j ½a t 1 X kþ j¼0
ðjÞ y0
ðjÞ
y0 þ
k 1 X d f ðt ; y Þ; CðaÞ j¼0 j;kþ1 j j
2
2
ð43Þ
ð44Þ
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368
3359
where dj,k+1 are defined as (42) and ej,k+1 are given by
8 1a a h aþ1 i ðkþ2Þ h ha 1 aþ1 > > k þ 12 ; j ¼ 0; > Cðaþ1Þ þ Cðaþ2Þ k 2 > > h i > > a aþ1 aþ1 > 1 aþ1 a þ1 h 3 > 3ð12 Þ j ¼ k; þ Cðahþ1Þ 12 ; > > > Cðaþ2Þ 2 > > a > 1 aþ1 : 2h ; j ¼ k þ 1:
ej;kþ1
ð45Þ
Cðaþ2Þ 2
At the beginning, we apply the product rectangle rule to calculate it
a aha 1 ; Cða þ 2Þ 2
e0;1 ¼ e1;1
a a h 1 ¼ : Cða þ 2Þ 2
ð46Þ ð47Þ
We call the above algorithm Improved Algorithm I (descried by (38), (39), (43) and (44)), namely, using the intermediate nodes with the predictor–corrector method (similar to the classical Adams methods), which can be upgraded to ensure an overall accuracy. Improved accuracy of this algorithm found in the middle node is still inadequate, so we can use the formula of Simpson to calculate the intermediate nodes. Thus
ykþ1 ¼ 2
yPkþ1 ¼ 2
½a k k X X tjkþ1 ðjÞ X ej;kþ1 f ðt j1 ; yj1 Þ þ ekþ1;kþ1 f tkþ1 ; yPkþ1 þ fj;kþ1 f ðt j ; yj Þ þ fkþ1;kþ1 f t kþ1 ; yPkþ1 ; y0 þ 2 2 2 2 j! j¼1 j¼0 j¼0 j ½a t 1 X kþ 2
j¼0
j!
ðjÞ
y0 þ
k 1 X d f ðt ; y Þ; CðaÞ j¼0 j;kþ1 j j
ð48Þ
ð49Þ
where
dj;kþ1
8 h i < ha k þ 1 j a k 1 j a ; 0 6 j 6 k 1; 2 2 a ¼ : ha 1a ; j ¼ k; a 2
ð50Þ
ej;kþ1
8
a aþ1 aþ1 a aþ2 aþ2 > < 4 ðkþ32jÞ þðkjþ12Þ h ðkþ32jÞ ðkjþ12Þ 8h ; 1 6 j 6 k; Cðaþ2Þ Cðaþ3Þ ¼ > : 2ha 1 aþ1 ; j ¼ k þ 1; Cðaþ2Þ 2
ð51Þ
fj;kþ1
8 aþ1 aþ1 a aþ2 aþ2 a 3ðkþ12Þ þðk12Þ h 4 ðk12Þ ðkþ12Þ h > ha 1 a > k þ ; j ¼ 0; > > Cðaþ1Þ 2 Cðaþ2Þ Cðaþ3Þ > <
aþ2 aþ2 a aþ1 aþ1 aþ1 a 3 1 3 1 1 ¼ 4 ðkþ2jÞ ðkj2Þ h ðkþ2jÞ þ6ðkjþ2Þ þðkj2Þ h ; 1 6 j 6 k; > Cðaþ3Þ Cðaþ2Þ > > h h > 1aþ1 i aþ2 i > 3 aþ1 4ha 3 aþ2 : ha þ 5 12 þ ; j ¼ k þ 1: 2 2 2 Cðaþ2Þ Cðaþ3Þ
ð52Þ
The above algorithm is called Improved Algorithm II (described by (38), (39), (48) and (49)), namely, using the Simpson formats to calculate the intermediate nodes in order to enhance accuracy. In order to further improve the algorithm’s overall accuracy, there is still room for improvement in the predictor. Observing the above algorithm, one can find at each step whatever format we used in the predictor, its main format is not changed. If the value after correction continues to serve as a predictor, one can improve the accuracy of the whole algorithm, so as to enhance the overall accuracy of the format. However the computational complexity is almost not increased. Thus
yPkþ1 ¼
" # ½a k k X X t jkþ1 ðjÞ 1 X aj;kþ1 f ðtj ; yj Þ þ akþ1;kþ1 f ðt kþ1 ; e y Pkþ1 Þ þ bj;kþ1 f tjþ1 ; ykþ1 ; y0 þ 2 2 j! CðaÞ j¼0 j¼0 j¼0
ð53Þ
ykþ1 ¼
½a k k X X tjkþ1 ðjÞ X ~Pkþ1 Þ; ej;kþ1 f t j1 ; yk1 þ ekþ1;kþ1 f t kþ1 ; ypkþ1 þ fj;kþ1 f ðtj ; yj Þ þ fkþ1;kþ1 f ðt kþ1 ; y y0 þ 2 2 2 2 j! j¼1 j¼0 j¼0
ð54Þ
ypkþ1 ¼
½a k k X X tjkþ1 ðjÞ X ej;kþ1 f t j1 ; yk1 þ ekþ1;kþ1 f t kþ1 ; yPkþ1 þ fj;kþ1 f ðtj ; yj Þ þ fkþ1;kþ1 f ðt kþ1 ; e y Pkþ1 Þ; y0 þ 2 2 2 2 j! j¼1 j¼0 j¼0
ð55Þ
2
2
3360
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368
where e y Pkþ1 is updated by (39). The remaining format is the same as that of Improved Algorithm II. The above algorithm is called Improved Algorithm III (described by (38), (53), (54) and (55)), namely, based on the Improved Algorithm II, the value after correction continues to serve as the predictor value by using the Simpson formula to calculate the intermediate nodes. 3.2. Error analysis First of all, we focus on the errors of the compound Simpson formulas and a corresponding result for compound rectangle rule and compound trapezoidal formula. They are easy to be verified by using the Simpson theory, we omit the details. Lemma 1. Let z 2 C3[0,T], then
Z kzð3Þ k tkþ1 kþ1 k X X a1 1 a1 4 ðt kþ1 tÞ zðtÞdt aj;kþ1 zðt j Þ bj;kþ1 zðt jþ1 Þ 6 t kþ1 h : 2 0 4! j¼0 j¼0
ð56Þ
Lemma 2. Let z 2 C1[0,T], then
Z kz0 k tkþ1 k X a1 1 a ðt kþ1 tÞ zðtÞdt cj;kþ1 zðt j Þ 6 t h: 0 a kþ1 j¼0
ð57Þ
Lemma 3. Let z 2 C1[0,T], then
Z kz0 k tkþ1 k X a1 1 a ðt kþ1 tÞ zðtÞdt dj;kþ1 zðt j Þ 6 t h: 0 a kþ1 j¼0
ð58Þ
Then we give the error analysis for the (39). Lemma 4. Assume that the initial value problem solution y satisfies
Z tkþ1 k X d ðt kþ1 tÞa1 C Da yðtÞdt cj;kþ1C Da yðtj Þ 6 Ctmkþ1 h 0 j¼0
ð59Þ
with m P 0 and d > 0. Then we have
yðt kþ1 Þ yP 6 C 1 hd þ C 2 max yðt j Þ y ; j kþ1 06j6k
where C1, C2 > 0. Proof. By the construction of yPkþ1 , we find that
yðt kþ1 Þ yP j ¼ kþ1
Z k X 1 tkþ1 ðt kþ1 tÞa1 f ðt; yðtÞÞdt cj;kþ1 f ðt j ; yj Þj CðaÞ 0 j¼0 Z k k X 1 tkþ1 1 X a1 a a 6 ðt kþ1 tÞ CD yðtÞdt cj;kþ1CD yðtj Þ þ c jf ðtj ; yðt j ÞÞ f ðt j ; yj Þj CðaÞ j¼0 j;kþ1 CðaÞ 0 j¼0 6
k k X C 1 tmkþ1 d 1 X d d cj;kþ1 jyðt j Þ yj j 6 C 1 h þ C 2 max jyðt j Þ yj j cj;kþ1 6 C 1 h þ C 2 max jyðt j Þ yj j: h þ 06j6k 06j6k CðaÞ CðaÞ j¼0 j¼0
Based on the error estimate of the preceding section, we now present a general convergence result for this new scheme. Theorem 9. Assume that the initial value problem solution y satisfies
Z tkþ1 k X d a1 a a 1 ðtkþ1 tÞ C D yðtÞdt cj;kþ1C D yðt j Þ 6 C 1 tmkþ1 h 1; 0 j¼0 Z tkþ1 k X d a1 2 a a 2 ðtkþ1 tÞ C D yðtÞdt dj;kþ1C D yðtj Þ 6 C 2 tmkþ1 h 2; 0 j¼0 Z tkþ1 k k X X d 3 ðtkþ1 tÞa1 C Da yðtÞdt aj;kþ1C Da yðt j Þ akþ1;kþ1C Da yðt j Þ bj;kþ1C Da yðt jþ1 Þ 6 C 3 t mkþ1 h 3; 2 0 j¼0 j¼0
3361
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368
with m1, m2, m3 P 0 and d1, d2, d3 > 0. Then we have q
max jyðtj Þ yj j ¼ Oðh Þ;
06j6N
where q = min{d1 + a, d2, d3} and N = [T/h]. Proof. Firstly by Lemma 4, we have d
jyðtkþ1 Þ yPkþ1 j 6 C 1 h 1 þ C 2 max jyðtj Þ yj j: 06j6k
Secondly, from formula of ykþ1 and Lemma 4, we easily have 2
d
jyðtkþ1 Þ ykþ1 j 6 C 3 h 2 þ C 4 2
2
So
k X
d
dj;kþ1 jyðtj Þ yj j 6 C 3 h 2 þ C 4 max jyðt j Þ yj j: 06j6k
j¼0
Z k k X X 1 tkþ1 ðt kþ1 tÞa1 f ðt; yðtÞÞdt aj;kþ1 f ðt j ; yj Þ akþ1;kþ1 f ðtkþ1 ; yPkþ1 Þ bj;kþ1 f t jþ1 ; yjþ1 2 2 CðaÞ 0 j¼0 j¼0 Z kþ1 k k X X 1 tkþ1 1 X 6 ðt kþ1 tÞa1 C Da yðtÞdt aj;kþ1 C Da yðtj Þ bj;kþ1CD a y t jþ1 þ a ðf ðt j ; yðt j ÞÞ 2 CðaÞ 0 CðaÞ j¼0 j;kþ1 j¼0 j¼0
jyðtkþ1 Þ ykþ1 j ¼
k 1 1 X bj;kþ1 f tjþ1 ; yðt jþ1 Þ f t jþ1 ; yjþ1 akþ1;kþ1 ðf ðtkþ1 ; yðtkþ1 ÞÞ f ðtkþ1 ; yPkþ1 ÞÞ þ 2 2 2 2 CðaÞ CðaÞ j¼0 k X d d 6 C5h 3 þ C6 aj;kþ1 jyðtj Þ yj j þ C 7 akþ1;kþ1 C 1 h 1 þ C 2 max jyðtj Þ yj j f ðt j ; yj ÞÞ þ
06j6k
j¼0
þ C8
k X
d
bj;kþ1 C 3 h 2 þ C 4 max jyðt j Þ yj j
06j6k
j¼0 d1 þa
d
6 C5h 3 þ C7h
minfd3 ;d1 þa;d2 g
d
þ C 8 h 2 þ C 9 max jyðt j Þ yj j 6 Ch 06j6k
:
In the following, we mainly deal with the error analysis for the improved algorithms I–III respectively. Theorem 10. Assume that the initial value problem solution y satisfies
Z tkþ1 k X d 2 ðt kþ1 tÞa1 C Da yðtÞdt dj;kþ1C Da yðtj Þ 6 C 1 t mkþ1 1 h 1 ; 2 0 2 j¼0 Z tkþ1 k X d a1 a 2 ðt kþ1 tÞ C D yðtÞdt ej;kþ1C Da yðt j Þ ekþ1;kþ1C Da yðt jþ1 Þ 6 C 2 t mkþ2 1 h 2 ; 2 2 0 2 j¼0 with q
max jyðtj Þ yj j ¼ Oðh Þ;
06j6N
where m1, m2 P 0 and d1, d2 > 0 for some n, T > 0. Then we have p
max jyðt jþ1 Þ yjþ1 j ¼ Oðh Þ;
06j6N1
2
2
where p = min(q, d1 + a, d2) and N = [T/h]. Proof. Firstly, we find that
Z k X 1 tkþ12 a1 ðt kþ1 tÞ f ðt; yðtÞÞdt dj;kþ1 f ðtj ; yj Þ 2 CðaÞ 0 j¼0 Z k k X X 1 tkþ12 d 6 ðtkþ1 tÞa1 Da yðtÞdt dj;kþ1 Da yðtj Þ 6 C 1 h 1 þ C 2 max jyðtj Þ yj j dj;kþ1 2 06j6k CðaÞ 0 j¼0 j¼0
yðt kþ1 Þ yPkþ1 ¼ 2
2
d
6 C 1 h 1 þ C 2 max jyðt j Þ yj j: 06j6k
3362
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368
Finally, one gets that
Z k k a1 X 1 tkþ12 1 X tkþ1 t f ðt; yðtÞÞdt ej;kþ1 f ðt j ; yj Þ ekþ1;kþ1 f t kþ1 ; ykþ1 þ e jf ðt ; yðt j ÞÞ 2 2 2 CðaÞ 0 CðaÞ j¼0 j;kþ1 j j¼0 1 f ðt j ; yj Þj þ ekþ1;kþ1 f ðt kþ1 ; y t kþ1 Þ f t kþ1 ; yPkþ1 2 2 2 2 CðaÞ d a 6 C 3 h 2 þ C 4 max jyðtj Þ yj j þ C 5 h jy tkþ1 ykþ1 j
jyðtkþ1 Þ ykþ1 j ¼ 2
2
2
06j6k
d2
6 C 3 h þ C 4 max jyðtj Þ yj j þ C 5 h
aþd1
06j6k
p
6 Ch :
2
minfd1 þa;d2 g
þ C 1 max jyðt j Þ yj j 6 C 1 h 06j6k
þ C 2 max jyðt j Þ yj j 06j6k
Theorem 11. Assume that the initial value problem solution y satisfies
Z tkþ1 k X d a1 a a 2 ðtkþ1 tÞ C D yðtÞdt dj;kþ1C D yðtj Þ 6 C 1 t mkþ1 1 h 1 ; 2 0 2 j¼0 Z tkþ1 k X 2 ðtkþ1 tÞa1 C Da yðtÞdt ej;kþ1C Da yðt j1 Þ ekþ1;kþ1C Da y t kþ1 2 2 2 0 j¼1 k X d fj;kþ1C Da yðt j Þ fkþ1;kþ1 C Da yðt kþ1 Þ 6 C 2 tmkþ2 1 h 2 ; 2 j¼0 with q
max jyðtj Þ yj j ¼ Oðh 1 Þ;
06j6N
and q
max jyðtj Þ yPj j ¼ Oðh 2 Þ;
06j6N
where m1, m2 P 0 and d1, d2 > 0. Then for n, T > 0, we have
p max yðtjþ1 Þ yjþ1 ¼ Oðh Þ; 2 2
06j6N1
where p = min(q1, q1 + a, q2 + a, d1 + a, d2) and N = [T/h]. Proof. Based on the above discussion, we find that
d yðt kþ1 Þ yPkþ1 6 C 1 h 1 þ C 2 max jyðt j Þ yj j: 2
06j6k
2
In the following, one gets
yðt kþ1 Þ ykþ1 ¼ 2
2
Z k X 1 tkþ12 ðtkþ1 tÞa1 f ðt; yðtÞÞdt fj;kþ1 f ðt j ; yj Þ fkþ1;kþ1 f ðtkþ1 ; yPkþ1 Þ 2 CðaÞ 0 j¼0 Z k X 1 tkþ12 P ej;kþ1 f ðtj1 ; yj1 Þ ekþ1;kþ1 f t kþ1 ; ykþ1 6 ðtkþ1 tÞa1 f ðt; yðtÞÞdt 2 2 2 2 2 CðaÞ 0 j¼1
k X j¼0
fj;kþ1 f ðtj ; yðt j ÞÞ fkþ1;kþ1 f ðt kþ1 ; yðt kþ1 ÞÞ
k X j¼1
1 ej;kþ1 f ðt j1 ; yðt j1 ÞÞ ekþ1;kþ1 f ðtkþ1 ; yðt kþ1 ÞÞ þ 2 2 2 2 CðaÞ
k 1 X 1 ej;kþ1 f ðtj1 ; yj1 Þ f ðt j1 ; yðtj1 ÞÞ þ fj;kþ1 jf ðt j ; yj Þ f ðtj ; yðt j ÞÞj þ fkþ1;kþ1 f tkþ1 ; yPkþ1 2 2 2 2 C ð a Þ C ð a Þ j¼1 j¼0 1 d f ðt kþ1 ; yðt kþ1 ÞÞj þ ekþ1;kþ1 f t jþ1 ; yPjþ1 f ðt jþ1 ; yðt jþ1 ÞÞ 6 C 1 h 2 þ C 2 max yðt jþ1 Þ yjþ1 2 2 2 2 2 2 06j6k CðaÞ a a d þ C 3 max jyðt j Þ yj j þ C 4 yðtkþ1 Þ yPkþ1 h þ C 5 yðt kþ1 Þ yPkþ1 h 6 C 1 h 2 þ C 2 max yðt jþ1 Þ yjþ1 2 2 2 2 06j6k 06j6k q aþq aþd aþq p p þ C 3 h 1 þ C 4 h 2 þ C 5 h 1 þ C 6 h 1 6 Ch þ C 2 max yðt jþ1 Þ yjþ1 6 Ch : k X
06j6k
2
2
3363
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368
Theorem 12. Assume that the initial value problem solution y satisfies
Z tkþ1 k X d a1 a a 1 ðt kþ1 tÞ C D yðtÞdt cj;kþ1C D yðt j Þ 6 C 1 tmkþ1 h 1; 0 j¼0 Z tkþ1 k k X X d a1 a a a 2 ðt kþ1 tÞ C D yðtÞdt aj;kþ1C D yðtj Þ bj;kþ1 D yðt jþ1 Þ 6 C 3 tmkþ1 h 2; 2 0 j¼0 j¼0 q max yðt jþ1 Þ yjþ1 ¼ Oðh Þ; 06j6k
2
2
where 0 6 j 6 N, m1, m2, m3 P 0 and d1, d2, d3 > 0 for n, T > 0. Then we have
yðt k Þ yP ¼ Oðhp Þ; k
where p = min{d1 + a, q, d2}. Proof. From Lemma 4, we have
d yðt kþ1 Þ e y Pkþ1 6 C 1 h 1 þ C 2 max jyðtj Þ yj j: 06j6k
Next we focus on the error for yPkþ1 . It follows
Z k k X X 1 tkþ1 a1 P e ðtkþ1 tÞ f ðt; yðtÞÞdt aj;kþ1 f ðt j ; yj Þ akþ1;kþ1 f ðt kþ1 ; y kþ1 Þ bj;kþ1 f ðtjþ1 ; yjþ1 Þ 2 2 CðaÞ 0 j¼0 j¼0 Z kþ1 k k t kþ1 X X 1 1 X 6 ðt kþ1 tÞa1 C Da yðtÞdt aj;kþ1 C Da yðt j Þ bj;kþ1C Da yðtjþ1 Þ þ a jf ðt j ; yðt j ÞÞ 2 CðaÞ j¼0 j;kþ1 CðaÞ 0 j¼0 j¼0
yðt kþ1 Þ yP ¼ kþ1
k 1 1 X y Pkþ1 þ bj;kþ1 f ðt jþ1 ; yðt jþ1 ÞÞ f ðtjþ1 ; yjþ1 Þ akþ1;kþ1 f ðtkþ1 ; yðtkþ1 ÞÞ f t kþ1 ; e 2 2 2 2 CðaÞ CðaÞ j¼0 a d 6 C 1 h 2 þ C 2 max jyðt j Þ yj j þ C 3 yðt kþ1 Þ e y Pkþ1 h þ C 4 max yðt jþ1 Þ yjþ1 f ðt j ; yj Þj þ
06j6k
2
06j6k
aþd1
d2
6 C 1 h þ C 2 max jyðt j Þ yj j þ C 3 h 06j6k
q
þ C 4 max jyðt j Þ yj j þ C 5 h 6 Ch 06j6k
2
minfd1 þa;q;d2 g
:
3.3. Stability analysis An interpretation of stability of a numerical scheme is that for small errors in the initial conditions cause small errors in the solution [29,23]. In the following, we firstly deal with the stability analysis of the new scheme (38), (39), (41). Theorem 13. Suppose yk+1 and y0kþ1 are numerical solutions in (38), respectively. Then
ðjÞ
0ðjÞ
and the initial conditions are given by y0 and y0
jykþ1 y0kþ1 j 6 Kky0 y00 k1 for any k, i.e. the new scheme (38), (39) and (41) are numerically stable. Sketch of the proof. According to the expression of (38), we easily get
jykþ1
y0kþ1 j
" # ½a X k k X tjkþ1 ðjÞ 1 X P ¼ a f ðt ; y Þ þ akþ1;kþ1 f ðtkþ1 ; ykþ1 Þ þ bj;kþ1 f ðtjþ1 ; ykþ1 Þ y þ 2 2 j¼0 j! 0 CðaÞ j¼0 j;kþ1 j j j¼0 " #! j ½ a k k Xt X 1 X kþ1 0ðjÞ 0 0P 0 a f ðt ; y Þ þ akþ1;kþ1 f ðtkþ1 ; ykþ1 Þ þ bj;kþ1 f t jþ1 ; ykþ1 y þ 2 2 j! 0 CðaÞ j¼0 j;kþ1 j j j¼0 j¼0 6
½a k X X tjkþ1 1 aj;kþ1 jf ðtj ; yj Þ f ðt j ; y0j Þj þ akþ1;kþ1 f ðt kþ1 ; yPkþ1 Þ f ðt kþ1 ; y0P ky0 y00 k1 þ kþ1 Þ j! C ð a Þ j¼0 j¼0 ! k X bj;kþ1 f ðtjþ1 ; ykþ1 Þ f ðt jþ1 ; y0kþ1 Þ 6 Kky0 y00 k1 þ C 1 max jyj y0j j þ C 2 jyPkþ1 y0P þ kþ1 j þ C 3 max jykþ1
j¼0 y0kþ1 j 2
2
2
2
2
6 K 1 ky0 y00 k1 þ K 2 max jyj y0j j 6 Kky0 y00 k1 :
06j6k
06j6k
Next, we deal with the stability analysis for the Improved Algorithms I–III respectively.
06j6k
2
3364
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368
Theorem 14. Suppose ykþ1 and y0kþ1 are numerical solutions in (43), 2 2 respectively. Then
ðjÞ
0ðjÞ
ðjÞ
0ðjÞ
and the initial conditions are given by y0 and y0
ykþ1 y0kþ1 6 Kky0 y00 k1 2
2
for any k, i.e. the Improved Algorithm I is numerically stable. Proof. The proof is similar to that of Theorem 13, so is omitted here. h By almost the same reasoning, we can get Theorem 15. Theorem 15. Suppose ykþ1 and y0kþ1 are numerical solutions in (48), 2 2 respectively. Then
and the initial conditions are given by y0 and y0
ykþ1 y0kþ1 6 Kky0 y00 k1 2
2
for any k, i.e. the Improved Algorithm II is numerically stable. ðjÞ
0ðjÞ
Theorem 16. Suppose yPkþ1 and y0P kþ1 are numerical solutions in (53), and the initial conditions are given by y0 and y0 respectively. Then
P y y0P 6 Kky y0 k 0 kþ1 kþ1 0 1 for any k, i.e. the Improved Algorithm III is numerically stable. Sketch of the proof. According to the expression of (53), we get " # ½a j k k X P X tkþ1 ðjÞ 1 X P y y0P ¼ e a f ðt ; y Þ þ akþ1;kþ1 f ðtkþ1 ; y kþ1 Þ þ bj;kþ1 f ðtjþ1 ; ykþ1 Þ y þ kþ1 kþ1 2 2 j¼0 j! 0 CðaÞ j¼0 j;kþ1 j j j¼0
( ½a j Xt
kþ1
j¼0
6
½a j X t
kþ1
j¼0
þ
k X j¼0
j!
j!
0ðjÞ
y0 þ
" #) k k X 1 X 0 aj;kþ1 f ðtj ; y0j Þ þ akþ1;kþ1 f ðtkþ1 ; e y 0P Þ þ b f ðt ; y Þ 1 1 j;kþ1 jþ2 kþ kþ1 2 CðaÞ j¼0 j¼0
ky0 y00 k1 þ
k X 1 a jf ðt ; y Þ f ðtj ; y0j Þj þ akþ1;kþ1 f tkþ1 ; e y Pkþ1 f tkþ1 ; e y 0P kþ1 CðaÞ j¼0 j;kþ1 j j
! 0 bj;kþ1 f ðtjþ1 ; ykþ1 Þ f ðtjþ1 ; y0kþ1 Þ 6 Kky0 y00 k1 þ C 1 max jyj y0j j þ C 2 j e y Pkþ1 e y 0P kþ1 j þ C 3 max ykþ1 ykþ1 : 2
2
2
06j6k
2
06j6k
2
2
Because e y Pkþ1 and e y 0P kþ1 are updated by (39), so
P e y
kþ1
0 0 e y 0P kþ1 6 Kky0 y0 k1 þ C max jyj yj j: 06j6k
Besides, according to Theorem 15, we have
jykþ1 y0kþ1 j 6 Kky0 y00 k1 : 2
2
Thus
P y y0P 6 Kky y0 k þ C 1 max jy y0 j þ C 2 Kky y0 k þ C max jy y0 j þ C 3 ky y0 k 0 j 0 j 0 kþ1 kþ1 0 1 j 0 1 j 0 1 06j6k
6 Kky0 y00 k1 þ C max jyj y0j j 6 Kky0 y00 k1 :
06j6k
06j6k
The sum of coefficients in our designed algorithm are bounded. Besides, in the process of above proof, we do not distinguish the constants C1–C9 and K1–K2 for convenience. 4. Numerical examples In order to verify our constructed numerical schemes for the fractional calculus and fractional differential equations, we give the following numerical examples. And the numerical results show that the numerical computational algorithms are efficient.
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368
3365
Example 1. Consider the function f(x) = x4, x 2 [0, 1]. Then kf ð2Þ k1 ¼ 12; kf ð4Þ k1 ¼ 24. The numerical results are shown in Tables 1–4. Example 2. Consider the function f(x) = x4 and f(x) = x6, x 2 [0, 1]. The numerical results are shown in Tables 5–8. Example 3.
( a
C D0;t yðtÞ
¼
2
Cð3aÞ t
2a
2a 2 Cð3aÞ t
yðtÞ þ t 2 t;
1a 1 Cð2aÞ t
a > 1; yðtÞ þ t t; a 6 1:
ð60Þ
2
The numerical results are shown in Tables 9–14. Example 4. 0:5 C D0;t yðtÞ
¼
2
Cð2:5Þ
t 1:5 yðtÞ þ t2 :
ð61Þ
The initial value we chosen as y(0) = 0. Its true solution is y(t) = t2. The numerical results are shown in Table 15.
Table 1 The absolute error of J a0;1 x4 for verifying Theorem 5. h
a = 0.8
a=1
a = 1.6
0.1 0.05 0.025 0.0125 0.00625
4.161861E003 1.051866E003 2.644498E004 6.631336E005 1.660625E005
3.330000E003 8.331250E004 2.083203E004 5.208252E005 1.302078E005
1.505989E003 3.745245E004 9.348082E005 2.335869E005 5.838784E006
Table 2 The absolute error of J a0;1 x4 for verifying Theorem 6. h
a = 0.8
a=1
a = 1.6
0.1 0.05 0.025 0.0125 0.00625
3.542874E006 2.223854E007 1.393338E008 8.711287E010 5.008000E011
3.333333E006 2.083333E007 1.302083E008 8.138021E010 5.086262E011
2.334466E006 1.457874E007 9.110534E009 5.961440E010 5.476738E010
Table 3 The absolute error of J a0;1 x4 for verifying Theorem 7 (on the first suggestion). h
a = 0.8
a=1
a = 1.6
0.1 0.05 0.025 0.0125 0.00625
3.552355E006 2.225207E007 1.393541E008 8.711598E010 5.008050E011
3.333333E006 2.083333E007 1.302083E008 8.138021E010 5.086262E011
2.301434E006 1.452862E007 9.102809E009 5.960241E010 5.476719E010
Table 4 The absolute error of J a0;1 x4 for verifying Theorem 7 (on the second suggestion). h
a = 0.8
a=1
a = 1.6
0.1 0.05 0.025 0.0125 0.00625
9.953847E005 2.881589E005 8.580664E006 2.563213E006 7.636424E007
4.399626E005 1.031004E005 2.536443E006 6.316066E007 1.577471E007
1.048313E005 3.435521E006 9.687469E007 2.600681E007 6.763113E008
3366
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368 Table 5 The absolute error of C Da0;1 x4 for verifying (29). h
a = 0.8
a=1
a = 1.6
0.1 0.05 0.025 0.0125 0.00625
1.055717E002 2.892040E003 7.770508E004 2.059158E004 5.400317E005
0 0 0 0 0
1.994496E002 5.143291E003 1.315600E003 3.345430E004 8.470497E005
Table 6 The absolute error of C Da0;1 x4 for verifying (30) (on the first suggestion). h
a = 0.8
a=1
a = 1.6
0.1 0.05 0.025 0.0125 0.00625
1.441004E007 3.973507E009 1.180234E010 2.718403E011 9.706680E011
0 0 0 0 0
2.055812E006 1.159297E007 6.887261E009 4.459828E010 2.237222E009
Table 7 The absolute error of C Da0;1 x4 for verifying (30) (on the second suggestion). h
a = 0.8
a=1
a = 1.6
0.1 0.05 0.025 0.0125 0.00625
8.270386E004 3.710286E004 1.647685E004 7.272723E005 3.195694E005
0 0 0 0 0
1.189139E003 5.834701E004 2.408542E004 9.464361E005 3.653170E005
Table 8 The absolute error of C Da0;1 x6 for verifying (31). h
a = 0.8
a=1
a = 1.6
0.1 0.05 0.025 0.0125 0.00625
4.945389E005 3.431437E006 2.327757E007 1.554474E008 1.084208E009
0 0 0 0 0
9.759985E005 6.328401E006 4.063441E007 2.582014E008 2.267807E009
Table 9 Errors for Eq. (60) with t = 1, a < 1. h
a = 0.3
a = 0.5
a = 0.7
0.1 0.05 0.025 0.0125 0.00625 0.003125
0.0225 0.0124 0.0061 0.0029 0.0014 6.6099E004
0.0041 0.0031 0.0018 0.0010 5.2963E004 2.7256E004
0.0059 0.0018 4.9076E004 8.2064E005 2.2522E005 3.6356E005
Table 10 Errors for Eq. (60) with t = 1, a > 1. h
a = 1.25
a = 1.5
a = 1.85
0.1 0.05 0.025 0.0125 0.00625 0.003125
0.0255 0.0127 0.0063 0.0031 0.0016 7.8007E004
0.0200 0.0101 0.005 0.0025 0.0012 6.1623E004
0.0042 0.0031 0.0018 0.001 5.4928E004 2.8978E004
3367
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368 Table 11 (Improved Algorithm I) Errors for Eq. (60) with t = 1, a < 1. h
a = 0.3
a = 0.5
a = 0.7
0.1 0.05 0.025 0.0125 0.00625 0.003125
0.0021 4.5357E004 8.6128E005 1.1965E005 3.33131E007 1.2418E006
0.0018 7.1792E004 2.7651E004 1.0326E004 3.7853E005 1.3714E005
0.0066 0.0031 0.0014 5.7564E004 2.4134E004 1.002E004
Table 12 (Improved Algorithm I) Errors for Eq. (60) with t = 1, a > 1. h
a = 1.25
a = 1.5
a = 1.85
0.1 0.05 0.025 0.0125 0.00625 0.003125
0.0014 4.4287E004 1.3188E004 4.3188E005 1.3340E005 4.968E006
0.0018 7.1792E004 2.7651E004 1.0326E004 3.7853E005 1.3714E005
0.0073 0.0031 0.0013 5.9123E004 2.6345E004 1.18E004
Table 13 (Improved Algorithm II) Errors for Eq. (60) with t = 1, a < 1. h
a = 0.3
a = 0.5
a = 0.7
0.1 0.05 0.025 0.0125 0.00625 0.003125
0.0161 0.0059 0.0021 7.5212E004 2.7028E004 9.8093E005
0.0036 0.0011 3.402E004 1.0651E004 3.4101E005 1.1161E005
9.2062E005 1.1114E004 6.5809E005 3.2349E005 1.4721E005 6.4434E006
Table 14 (Improved Algorithm II) Errors for Eq. (60) with t = 1, a > 1. h
a = 1.25
a = 1.5
a = 1.85
0.1 0.05 0.025 0.0125 0.00625 0.003125
3.7015E004 9.3442E005 2.4233E005 6.4536E006 1.7602E006 4.8997E007
0.0015 5.2098E004 1.8259E004 6.4261E005 2.22667E005 8.0045E006
0.0130 0.0058 0.0026 0.0012 5.3348E004 2.4039E004
Table 15 Errors for Eq. (61) with t = 10. h
Improved Algorithm II
Improved Algorithm III
0.1 0.05 0.025 0.0125
0.0361 0.0117 0.0039 0.0013
0.0032 6.7576E004 1.5012E004 3.4362E005
References [1] K.B. Oldham, J. Spanier, The Fractional Calculus: Theory and Applications of Differentiation and Integration to Arbitrary Order, Academic Press, New York, 1974. [2] I. Podlubny, Fractional Differential Equations, Academic Press, 1999. [3] S.G. Samko, A.A. Kilbas, O.I. Marichev, Fractional Integrals and Derivatives Theory and Applications, Gordon and Breach, New York, 1993. [4] A.A. Kilbas, H.M. Srivastava, J.J. Trujillo, Theory and Application of Fractional Differential Equations, Elsevier, 2006. [5] Z. Odibat, Approximations of fractional integral and Caputo fractional derivatives, Appl. Math. Comput. 178 (2006) 527–533. [6] K. Diethelm, N.J. Ford, A.D. Freed, Y. Luchko, Algorithms for the fractional calculus: a selection of numerical methods, Comput. Methods Appl. Mech. Eng. 194 (6–8) (2005) 743–773. [7] C. Lubich, Discretized fractional calculus, SIAM J. Math. Anal. 17 (3) (1986) 704–719. [8] W. Deng, Short memory principle and a predictor–corrector approach for fractional differential equations, J. Comput. Appl. Math. 206 (2007) 174–188. [9] L. Yuan, O.P. Agrawal, A numerical scheme for dynamic systems containing fractional derivatives, ASME J. Vibr. Acoust. 124 (2002) 321–324. [10] K. Diethelm, An improvement of a nonclassical numerical method for the computation of fractional derivatives, Numer. Algorithms 47 (2008) 190– 361.
3368
C. Li et al. / Journal of Computational Physics 230 (2011) 3352–3368
[11] D.A. Murio, On the stable numerical evaluation of Caputo fractional derivatives, Comput. Math. Appl. 51 (2006) 1539–1550. [12] T. Miyakoda, Discretized fractional calculus with a series of chebyshev polynomial, Electron. Notes Theor. Comput. Sci. 225 (2009) 239–244. [13] K. Dithelm, N.J. Ford, A.D. Freed, A predictor–corrector approach for the numerical solution of fractional differential equations, Nonlinear Dyn. 29 (2002) 3–22. [14] K. Dithelm, N.J. Ford, A.D. Freed, Detailed error analysis for a fractional adams method, Numer. Algorithms 36 (2004) 31–52. [15] C.P. Li, C.X. Tao, On the fractional adams method, Comput. Math. Appl. 58 (8) (2009) 1573–1588. [16] C.P. Li, Y.H. Wang, Numerical algorithm based on adomian decomposition for fractional differential equations, Comput. Math. Appl. 57 (2009) 1672– 1681. [17] O.P. Agrawal, A general finite element formulation for fractional variational problems, J. Math. Anal. Appl. 337 (2008) 1–12. [18] T. Tang, A finite difference scheme for partial integro-differential equations with weakly singular kernel, Appl. Numer. Math. 11 (1993) 309–319. [19] Q. Sheng, T. Tang, Optimal convergence of an euler and finite difference method for nonlinear partial integro-differential equations, Math. Comput. Model. 21 (1995) 1–11. [20] I. Podlubny, Matrix approach to discrete fractional calculus, Int. J. Theor. Appl. 3 (2000) 359–386. [21] I. Podlubny, A. Chechkin, T. Skovranek, Y. Chen, B.M.V. Jara, Matrix approach to discrete fractional calculus. II: Partial fractional differential equations, J. Comput. Phys. 228 (2009) 3137–3153. [22] Q. Yang, F. Liu, I. Turner, Numerical methods for fractional partial differential equations with riesz space fractional derivatives, Appl. Math. Model. 34 (2010) 200–218. [23] R. Lin, F. Liu, Fractional high order methods for the nonlinear fractional ordinary differential equation, Nonlinear Anal. 66 (2007) 856–869. [24] X. Li, C. Xu, Existence and uniqueness of the weak solution of the space–time fractional diffusion equation and a spectral method approximation, Commun. Comput. Phys. 8 (2010) 1016–1051. [25] M. Xi, Numerical Analysis Method, University of Science and Technology of China Press, 2004. [26] Z. Odibat, Computational algorithms for computing the fractional derivatives of functions, Math. Comput. Simul. 79 (7) (2009) 2013–2020. [27] C. Gerald, P. Wheatley, Applied Numerical Analysis, Addison-Wesley, USA, 2004. [28] K. Diethelm, An algorithm for the numerical solution of differential equations of fractional order, Electron. Trans. Numer. Anal. 5 (1997) 1–6. [29] J.W. Thomas, Numerical Partial Differential Equations: Finite Difference Methods, Springer, 1995.
