- Email: [email protected]
A robust numerical method for a fractional differential equation
Applied Mathematics and Computation 315 (2017) 445–452
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
A robust numerical method for a fractional differential equation Zhongdi Cen a,b, Anbo Le b,∗, Aimin Xu b a b
College of Finance and Trade, Ningbo Dahongying University, Ningbo, China Institute of Mathematics, Zhejiang Wanli University, Ningbo, China
a r t i c l e
i n f o
a b s t r a c t
MSC: 65L05 65R20 65L20
This paper is devoted to giving a rigorous numerical analysis for a fractional differential equation with order α ∈ (0, 1). First the fractional differential equation is transformed into an equivalent Volterra integral equation of the second kind with a weakly singular kernel. Based on the a priori information about the exact solution, an integral discretization scheme on an a priori chosen adapted mesh is proposed. By applying the truncation error estimate techniques and a discrete analogue of Gronwall’s inequality, it is proved that the numerical method is first-order convergent in the discrete maximum norm. Numerical results indicate that this method is more accurate and robust than finite difference methods when α is close to 0.
Keywords: Fractional differential equation Caputo fractional derivative Volterra integral equation Adapted mesh Convergence analysis
© 2017 Elsevier Inc. All rights reserved.
1. Introduction Fractional calculus has played a significant role in a variety of scientific and engineering fields, such as finance [12,18], control [11], viscoelasticity [6], and hydrology [1,13]. Those fractional models, described in the form of fractional differential equations, have been proved to be more appropriate for modeling many physical phenomena than the traditional integer order models, because fractional calculus enables the description of the memory properties of various materials and processes [17, Chapter 10]. The analytical solutions of most fractional differential equations can not be obtained, so approximate and numerical techniques must be used. In this paper, we consider the following fractional differential equation
Dα∗ u(x ) + f (x, u ) = 0,
x ∈ := (0, 1],
(1.1)
u (0 ) = γ ,
(1.2)
where Dα ∗ denotes the Caputo fractional derivative defined by
Dα∗ u(x ) =
1 (1 − α )
x 0
(x − t )−α u (t )dt ,
0 < α < 1,
(1.3)
the constant γ and function f are given. We assume that function f satisfies
¯ × R ∩ C 1 ( × R ). f (x, u ) ∈ C ∗
Corresponding author. E-mail address: [email protected] (A. Le).
http://dx.doi.org/10.1016/j.amc.2017.08.011 0 096-30 03/© 2017 Elsevier Inc. All rights reserved.
(1.4)
446
Z. Cen et al. / Applied Mathematics and Computation 315 (2017) 445–452
¯ (see e.g., [3, Under this hypothesis, the fractional differential equation (1.1) and a unique solution u(x ) ∈ C () (1.2) exists ¯ × R ∩ C q ( × R ) for some integer q ≥ 1 there exists a Theorem 6.8]). In [3, Theorem 6.28] it is shown that if f (x, u ) ∈ C positive constant C0 such that
( j) u (x ) ≤ C0 ,
if j = 0,
C0 xα − j ,
if j = 1, 2, . . . , q + 1,
(1.5)
for all x ∈ . These bounds indicate that the derivatives of the solution u may blow up at the end point x = 0. This singular behavior complicates the construction of the discretization scheme and the convergence analysis of the numerical method. Since analytic solutions for fractional differential equations are generally impossible to attain, various numerical methods for solving the fractional differential equations have been proposed. There are some papers which take account of the possibly singular behavior of solutions of the fractional differential equations. Numerical schemes based on the collocation methods are developed in [4,7–10,14–16], finite difference schemes are proposed in [2,5,19,20]. Because of the presence of the singular behavior, standard numerical methods fail to give accurate approximation, even for high-order methods. In this paper we first transform the fractional order initial value problem (1.1) and (1.2) into an equivalent Volterra integral equation of the second kind with a weakly singular kernel. Then we discretize the Volterra integral equation on an apriori chosen adapted mesh. A rigorous analysis about the convergence of the discretization scheme is given by taking account of the possibly singular behavior of the solution. Applying the truncation error estimate techniques and a discrete analogue of Gronwall’s inequality [21], we will prove that the scheme is first-order convergent in the discrete maximum norm. Numerical results are given to display that this method is more accurate and robust than finite difference methods when α is close to 0. The rest of the paper is organized as follows. The discretization scheme is described in Section 2. Convergence analysis of the scheme is given in Section 3. Finally, numerical experiments are presented in Section 4. Notation. Throughout the paper, C will denote a generic positive constant that is independent of the mesh. Note that C ¯ . is not necessarily the same at each occurrence. To simplify the notation we set gi = g(xi ) for any function g ∈ C () 2. Discretization scheme Based on the properties of the exact solution u(x) we construct an apriori mesh N ≡ {0 = x0 < x1 < · · · < xN = 1} with the mesh points
⎧ r i ⎪ ⎪ , ⎨ N xi = 1 ⎪ ⎪ ⎩ 1 + 2 1 − 2r ( i − N ) , r 2
N
0≤i≤
N , 2
(2.1)
N < i ≤ N, 2
2
where the discretization parameter N is a positive even integer, r = α1 > 1. The mesh points (2.1) are more densely at the region near x = 0 since the exact solution of the fractional differential equations (1.1) and (1.2) may be singular near x = 0; r away from x = 0 a uniform mesh is used. Pedas et al. [7,14–16] used the mesh xi = Ni for 0 ≤ i ≤ N, which may be coarse when i close to N. Our mesh is a slight modification of that in [7,14–16]. It is well known that the initial value problem (1.1) and (1.2) can be written as the following equivalent Volterra integral equation of the second kind with a weakly singular kernel [3, Lemma 6.2]
u (x ) = u (0 ) −
1
(α )
x 0
(x − t )α−1 f (t , u(t ) )dt,
u (0 ) = γ .
(2.2) (2.3)
In the following we will discrete this integral equation instead of the differential equations (1.1) and (1.2). An approximation to the integral can be obtained by a quadrature formula
xi 0
≈
(xi − t )α−1 f (t , u(t ) )dt =
i
f ( xk , uk )
k=1
=
i 1
α
i
k=1
xk xk−1
xk
xk−1
(xi − t )α−1 f (t , u(t ) )dt
(xi − t )α−1 dt
f (xk , uk ) (xi − xk−1 )α − (xi − xk )α .
k=1
Then, we have the following discretization scheme for the problem (2.2) and (2.3):
uNi = uN0 − uN0 = γ .
i
1 f xk , uNk (xi − xk−1 )α − (xi − xk )α , 1 ≤ i ≤ N, (α + 1 )
(2.4)
k=1
(2.5)
Z. Cen et al. / Applied Mathematics and Computation 315 (2017) 445–452
447
3. Convergence analysis In this section we give an apriori error estimate for the discretization scheme (2.4) and (2.5) on the apriori mesh (2.1). Let ziN = uN − ui , where uN is the solution of problem (2.4), (2.5) and ui is the solution of problem (2.2) and (2.3) at mesh i i point xi . Then, for the error ziN we have i
1 ak zkN (xi − xk−1 )α − (xi − xk )α = Ri , (α + 1 )
ziN +
1 ≤ i ≤ N,
(3.1)
k=1
z0N = 0,
(3.2)
where i
1
Ri =
(α ) −
xk xk−1
k=1
f (t , u(t ) )(xi − t )α −1 dt
i
1 f (xk , uk ) (xi − xk−1 )α − (xi − xk )α , (α + 1 )
(3.3)
k=1
∂f xk , uk + λk uNk − uk , ∂u
ak =
0 < λk < 1.
(3.4)
First we state a useful technical result, which is similar to the result in [19, Lemma 4.5], but under the nonuniform mesh
N . Lemma 3.1. For 2 ≤ i ≤ N there exists a constant C independent of N such that i
k=2
−1 ≤ CN −1 , (xi − xk−1 )α − (xi − xk )α hk xαk−1
(3.5)
where hk = xk − xk−1 is the mesh size. Proof. Let t denote the smallest integer that is greater than or equal to t for any t ∈ R. For the analysis of the estimates we distinguish two cases. Case I: 2 ≤ i ≤ N/2. By applying the mean value theorem we have
i/2 k=2
≤
−1 (xi − xk−1 )α − (xi − xk )α hk xαk−1
i/2 k=2
i/2 α−1
−1 −1 α (xi − xk )α−1 h2k xαk−1 ≤ α xi − x i/2 h2k xαk−1
r
i N
=α
−
1−r
i ≤C N
≤ CN −2
k=2
k i
N
i/2 1 N2
i/2 k=2
i/2
k=2
r α−1 i/2 r k=2
2r−2 k N
k−1 N
k N
−
k−1 N
r 2
k−1 N
r (α−1)
1−r
r−1 ≤ CN −1 ,
(3.6)
where we have used rα = 1 and k ≤ 2(k − 1 ) for k ≥ 2. Furthermore, we have
i
k= i/2 +1
≤
−1 (xi − xk−1 )α − (xi − xk )α hk xαk−1
max
i/2 +1≤k≤i
−1 hk xαk−1
i
(xi − xk−1 )α − (xi − xk )α
k= i/2 +1 α α −1 ≤ x i/2 xi − x i/2 max hk
i/2 +1≤k≤i
448
Z. Cen et al. / Applied Mathematics and Computation 315 (2017) 445–452
≤
≤
i/2
r (α−1) rα i N
N i 2N
1−r r−1 i N
r max
i N
i/2 +1≤k≤i
k N
k−1 N
−
r
1 N
≤ CN −1 ,
(3.7)
where we also have used rα = 1. Therefore, from (3.6) and (3.7) we conclude that the lemma holds true for Case I. Case II: N/2 < i ≤ N. By using the same method as that in Case I we also can obtain
i/2 k=2
−1 ≤ CN −1 . (xi − xk−1 )α − (xi − xk )α hk xαk−1
(3.8)
Furthermore, we have
i
k= i/2 +1
≤
−1 (xi − xk−1 )α − (xi − xk )α hk xαk−1
max
i/2 +1≤k≤i
−1 hk xαk−1
= xi − x i/2 ≤ xαi ·
α
·
max
i/2 +1≤k≤i
i
max
i/2 +1≤k≤i
−1 hk xαk−1
−1 hk xαk−1 .
(3.9)
rα r
≤r
i N
i N
k N
r−1 k N
k= i/2 +1
For i/2 + 1 ≤ k ≤ N/2 we have −1 xαi hk xαk−1 ≤
(xi − xk−1 )α − (xi − xk )α
1 N
− k 2N
k−1 N
r
1−r
k−1 N
r (α−1)
≤ CN −1 ,
(3.10)
where we also have used rα = 1 and k ≤ 2(k − 1 ) for k ≥ 2. For N/2 < k ≤ i we have
N/2 −1 xαi hk xαk−1 ≤ xαi 2(1 − 2−r )N −1 N
r (α−1)
≤ CN −1 ,
(3.11)
where we also have used rα = 1. Therefore, from (3.8)–(3.11) we conclude that the lemma also holds true for Case II.
The next lemma gives us a useful formula for the truncation error. Lemma 3.2. Let f satisfy the hypothesis (1.4) and
∂ f (x, u ) ≤ Cxα−1 for x ∈ . ∂x
(3.12)
Then the truncation errors of the discrete scheme (2.4) and (2.5) satisfy
|Ri | ≤ CN−1 , 1 ≤ i ≤ N,
(3.13)
where C is a positive constant independently of N. Proof. For the analysis of the truncation errors we distinguish two cases. Case 1: i = 1. From (3.3) we have
|R1 | ≤
1
(α )
x1 0
(x1 − t )α−1 | f (t , u(t ) )|dt +
≤ C xα1 = C N −1 ,
1
(α + 1 )
| f (x1 , u1 )|xα1 (3.14)
where we have used the hypothesis (1.4), the bound on u given by (1.5) and rα = 1. From this we know that the lemma holds true for Case 1.
Z. Cen et al. / Applied Mathematics and Computation 315 (2017) 445–452
449
Case 2: 1 < i ≤ N. From (3.3) we have i
1
Ri =
(α ) −
xk xk−1
k=1
f (t , u(t ) )(xi − t )α −1 dt
i
1 f (xk , uk ) (xi − xk−1 )α − (xi − xk )α (α + 1 ) k=1
= R1i + R2i , where
R1i =
R2i =
1
(α ) −
(3.15) x1
0
f (t , u(t ) )(xi − t )α −1 dt
1
i
1
(α )
k=2
f (x1 , u1 ) xαi − (xi − x1 )α ,
(α + 1 )
xk xk−1
(3.16)
f (t , u(t ) )(xi − t )α −1 dt
i
1 − f (xk , uk ) (xi − xk−1 )α − (xi − xk )α . (α + 1 )
(3.17)
k=2
By using the similar method as that in Case 1, from (3.16) we have
1 Ri ≤
1
(α )
x1 0
| f (t , u(t ) )|(xi − t )α−1 dt
1 | f (x1 , u1 )| xαi − (xi − x1 )α (α + 1 ) ≤ C xαi − (xi − x1 )α ≤ Cx1 (xi − x1 )α −1 +
≤ C xα1 = C N −1 ,
(3.18)
where we also have used the hypothesis (1.4), the bound on u given by (1.5) and rα = 1. From (3.17) we have
2 Ri = ≤
i xk α −1 dt [ f (t , u(t ) ) − f (xk , uk )](xi − t ) (α ) xk−1 k=2 1
i
1
(α )
| f (ξk , u(ξk ) ) − f (xk , uk )|
k=2
xk
xk−1
(xi − t )α−1 dt
by the mean value theorem for integrals with ξk ∈ (xk−1 , xk ). Hence, we can obtain
2 Ri ≤
i
1 | f (ξk , u(ξk ) ) − f (xk , uk )| · (xi − xk−1 )α − (xi − xk )α (α + 1 ) k=2
1
i
{[| f (ξk , u(ξk ) ) − f (ξk , uk )| + | f (ξk , uk ) − f (xk , uk )|] (α + 1 ) k=2 · (xi − xk−1 )α − (xi − xk )α i ∂ f
1 (ξk , μk ) · |u(ξk ) − uk | ≤ (α + 1 ) ∂u k=2 ∂ f +hk (η , u ) (xi − xk−1 )α − (xi − xk )α ∂x k k i ∂ f
α α ≤C (xi − xk−1 ) − (xi − xk ) hk u (ζk ) + (ηk , uk ) ∂x ≤
k=2
≤C
i
k=2
−1 , (xi − xk−1 )α − (xi − xk )α hk xαk−1
450
Z. Cen et al. / Applied Mathematics and Computation 315 (2017) 445–452
where we again have used the mean value theorem with ηk , ζk ∈ (xk−1 , xk ) and μk ∈ uk−1 , uk , the bound on u given by (1.5) and the hypotheses (1.4) and (3.12). Now, an appeal to Lemma 3.1 yields
2 Ri ≤ CN−1 .
(3.19)
Therefore, combining (3.18) and (3.19) with (3.15) we conclude that the lemma also holds true for Case 2.
Now we can get the following a priori error estimate for the discrete scheme (2.4) and (2.5). Theorem 3.3. Let u be the solution of (2.2), (2.3) and uN be the solution of the discrete scheme (2.4) and (2.5) on the mesh N . Then, under the hypotheses (1.4) and (3.12) we have
N ui − u(xi ) ≤ CN−1 , 1 ≤ i ≤ N
for sufficiently large N, where C is a positive constant independent of N. Proof. From (3.1) we have i−1 N
zi ≤ pi + wi vk zkN , 1 ≤ i ≤ N,
(3.20)
k=1
where
(α + 1 )|Rk | 1 , wk = , (α + 1 ) + ak hα (α + 1 ) + ak hα k k vk = |ak | (xi − xk−1 )α − (xi − xk )α , 1 ≤ k ≤ i. pk =
By applying the discrete analogue of Gronwall’s inequality [21, Theorem 3], we have
i−1 i−1 k N −1 zi ≤ pi + wi 1 + v jw j · vk pk 1 + v jw j , 1 ≤ i ≤ N. k=1
j=1
(3.21)
j=1
It is easy to conclude that for sufficiently large N
0 < pk ≤ CN −1
0 < wk ≤ C for 1 ≤ k ≤ i,
and
(3.22)
where we have used Lemma 3.2. Hence we can obtain i−1
k=1
vk pk
k
1 + v jw j
−1
≤ CN −1
i−1
vk
k=1
j=1
≤ CN −1
i−1
(xi − xk−1 )α − (xi − xk )α
k=1
= CN −1 xαi − (xi − xi−1 )α
≤ CN −1 , and i−1
1 + v j w j ≤ exp
j=1
i−1
(3.23)
v jw j
j=1
≤ exp C
i−1
xi − x j−1
α
j=1
= exp C xαi − (xi − xi−1 )α ≤ C,
− xi − x j
α
(3.24)
where we have used the inequalities (3.22) and (1 + y ) ≤ ey for y ≥ −1. Therefore, substituting (3.22)–(3.24) into (3.21) we have
N zi ≤ CN−1 , 1 ≤ i ≤ N,
which completes the proof.
Z. Cen et al. / Applied Mathematics and Computation 315 (2017) 445–452
451
4. Numerical experiments In this section, we verify experimentally the theoretical results obtained in the preceding section. Error estimates and convergence rates for the discrete scheme are presented for the following example which has been presented in [2]. We will give the comparisons of numerical solution on various schemes. Example. We consider the following fractional differential equation
Dα u(x ) + f (x, u(x )) = 0,
x ∈ (0, 1],
u ( 0 ) = 1, where 0 < α < 1, f (x, u(x )) = 2u(x ) + sin(u(x )) + 0.1u2 (x ) + s(x ), and s(x) is chosen such that the exact solution is
u ( x ) = xα + 2x + 1.
Table 1 Error estimates and convergence rates on an a priori adapted mesh for example.
α = 0.1 α = 0.2 α = 0.3 α = 0.4 α = 0.5 α = 0.6 α = 0.7 α = 0.8 α = 0.9
N = 64
N = 128
N = 256
N = 512
N = 1024
N = 2048
4.6534e−3 0.857 7.2937e−3 0.896 8.0308e−3 0.927 7.5963e−3 0.950 6.6375e−3 0.968 5.5351e−3 0.981 4.4872e−3 0.985 3.8578e−3 0.930 3.0207e−3 0.860
2.5688e−3 0.878 3.9183e−3 0.916 4.2241e−3 0.943 3.9311e−3 0.964 3.3922e−3 0.978 2.8051e−3 0.988 2.2665e−3 0.991 2.0243e−3 0.949 1.6646e−3 0.884
1.3976e−3 0.895 2.0770e−3 0.930 2.1966e−3 0.955 2.0155e−3 0.973 1.7216e−3 0.985 1.4144e−3 0.992 1.1400e−3 0.995 1.0488e−3 0.963 9.0217e−4 0.903
7.5173e−4 0.908 1.0899e−3 0.942 1.1328e−3 0.965 1.0267e−3 0.980 8.6974e−4 0.990 7.1100e−4 0.995 5.7204e−4 0.997 5.3808e−4 0.973 4.8234e−4 0.919
4.0062e−4 0.919 5.6737e−4 0.951 5.8048e−4 0.972 5.2059e−4 0.985 4.3801e−4 0.993 3.5672e−4 0.997 2.8663e−4 0.998 2.7404e−4 0.981 2.5504e−4 0.932
2.1189e−4 − 2.9350e−4 − 2.9602e−4 − 2.6305e−4 − 2.2011e−4 − 1.7875e−4 − 1.4350e−4 − 1.3883e−4 − 1.3364e−4 −
Table 2 Comparisons of our numerical method with three finite difference schemes on error estimates and convergence rates for example.
α = 0.05
Our method FDM1 FDM2 FDM3
α = 0.1
Our method FDM1 FDM2 FDM3
α = 0.15
Our method FDM1 FDM2 FDM3
N = 64
N = 128
N = 256
N = 512
N = 1024
N = 2048
2.5766e−3 0.834 1.3992e−2 0.033 6.2556e−2 −0.015 1.1971e−2 −0.051 4.6534e−3 0.857 2.3243e−2 0.060 5.8221e−3 −0.538 2.0039e−2 0.105 6.2386e−3 0.878 2.8934e−2 0.088 2.1720e−3 1.080 2.3506e−2 0.163
1.4450e−3 0.856 1.3676e−2 0.030 6.3214e−2 0.007 1.2401e−2 0.053 2.5688e−3 0.878 2.2302e−2 0.057 8.4544e−3 −0.080 1.8637e−2 0.107 3.3946e−3 0.898 2.7224e−2 0.087 1.0272e−3 −0.717 2.0995e−2 0.172
7.9855e−4 0.873 1.3393e−2 0.029 6.2905e−2 −0.035 1.1955e−2 0.053 1.3976e−3 0.895 2.1437e−2 0.056 8.9343e−3 −0.043 1.7300e−2 0.111 1.8213e−3 0.914 2.5635e−2 0.087 1.6887e−3 0.658 1.8635e−2 0.182
4.3613e−4 0.886 1.3129e−2 0.028 6.4469e−2 −0.134 1.1524e−2 0.053 7.5173e−4 0.908 2.0619e−2 0.056 9.2016e−3 0.010 1.6023e−2 0.114 9.6656e−4 0.926 2.4132e−2 0.089 1.0706e−3 −0.191 1.6431e−2 0.191
2.3593e−4 0.898 1.2876e−2 0.028 7.0739e−2 −1.480 1.1107e−2 0.054 4.0062e−4 0.919 1.9837e−2 0.056 9.1365e−3 0.002 1.4801e−2 0.118 5.0855e−4 0.937 2.2695e−2 0.090 1.2218e−3 −0.036 1.4389e−2 0.201
1.2661e−4 − 1.2632e−2 − 1.9734e−1 − 1.0703e−2 − 2.1189e−4 − 1.9081e−2 − 9.1248e−3 − 1.3635e−2 − 2.6570e−4 − 2.1316e−2 − 1.2529e−3 − 1.2518e−2 −
452
Z. Cen et al. / Applied Mathematics and Computation 315 (2017) 445–452
0 The nonlinear problem (2.4) and (2.5) is solved by using Newton’s method with uN, = 1.5(i = 1, 2, . . . , N ) as an initial i guess. We iteratively compute uN, j until
j j−1 max uN, − uN, < 10−6 . i i
1≤i≤N
We measure the accuracy in the discrete maximum norm
eN = max uNi − ui , 1≤i≤N
and the convergence rate
N
r = log2
eN e 2N
.
The numerical results for example are tabulated in Table 1. Table 1 shows that the computed solution converges to the exact solution with first-order accuracy for sufficiently large N and the numerical results do not depend strongly on the value of α , which supports the convergence estimate of Theorem 3.3. For comparison we also use the numerical methods based on the finite difference method on the uniform mesh (FDM1), the finite difference method on the a priori mesh N (FDM2) and the finite difference method on an a posteriori adapted mesh [2] (FDM3) to compute this example in Table 2. From Table 2 we confirm that our method proposed in this paper is more accurate and robust than finite difference methods when α is close to 0. Acknowledgment We would like to thank the anonymous reviewers for some suggestions for the improvement of this paper. The work was supported by Humanities and Social Sciences Planning Fund of Ministry of Education of China (Grant No. 14YJC790 0 06), Project of Philosophy and Social Science Research in Zhejiang Province (Grant No. 15NDJC243YB), Zhejiang Province Natural Science Foundation (Grant No. Y17D010024), Ningbo Municipal Natural Science Foundation (Grant Nos. 2017A610131, 2017A610140) and Ningbo Municipal Soft Science Foundation (Grant No. 2015A10045). References [1] P. Becker-Kern, M.M. Meerschaert, H.P. Scheffler, Limit theorem for continuous-time random walks with two scales, J. Appl. Prob. 41 (2) (2004) 455–466. [2] Z. Cen, A. Le, A. Xu, A posteriori error analysis for a fractional differential equation, Int. J. Comput. Math. 94 (6) (2017) 1185–1195. [3] K. Diethlm, The analysis of fractional differential equations, in: Lecture Notes in Mathematics, 2004, Springer, Berlin, 2010. [4] N.J. Ford, M.L. Morgado, M. Rebelo, A nonpolynomial collocation method for fractional terminal value problems, J. Comput. Appl. Math. 275 (2015) 392–402. [5] J.L. Gracia, M. Stynes, Central difference approximation of convection in caputo fractional derivative two-point boundary value problems, J. Comput. Appl. Math. 273 (2015) 103–115. [6] R.C. Koeller, Application of fractional calculus to the theory of viscoelasticity, J. Appl. Mech. 51 (2) (1984) 229–307. [7] M. Kolk, A. Pedas, E. Tamme, Modified spline collocation for linear fractional differential equations, J. Comput. Appl. Math. 283 (2015) 28–40. [8] N. Kopteva, M. Stynes, An efficient collocation method for a caputo two-point boundary value problem, BIT Numer. Math. 55 (2015) 1105–1123. [9] H. Liang, M. Stynes, Collocation methods for general caputo two-point boundary value problems, https://www.researchgate.net/publication/313878934. [10] X. Ma, C. Huang, Numerical solution of fractional integro-differential equations by a hybrid collocation method, Appl. Math. Comput. 219 (12) (2013) 6750–6760. [11] J.T. Machado, Discrete time fractional-order controllers, Fract. Calc. Appl. Anal. 4 (1) (2001) 47–66. [12] M.M. Meerschaert, E. Scalas, Coupled continuous time random walks in finance, Phys. A 370 (1) (2006) 114–118. [13] M.M. Meerschaert, Y. Zhang, B. Baeumer, Particle tracking for fractional diffusion with two time scales, Comput. Math. Appl. 59 (3) (2010) 1078–1086. [14] A. Pedas, E. Tamme, Piecewise polynomial collocation for linear boundary value problems of fractional differential equations, J. Comput. Appl. Math. 236 (13) (2012) 3349–3359. [15] A. Pedas, E. Tamme, Numerical solution of nonlinear fractional differential equations by spline collocation methods, J. Comput. Appl. Math. 255 (2014) 216–230. [16] A. Pedas, E. Tamme, Spline collocation for nonlinear fractional boundary value problems, Appl. Math. Comput. 244 (2014) 502–513. [17] I. Podlubny, Fractional Differential Equations, Academic Press, New York, 1999. [18] M. Raberto, E. Scalas, F. Mainardi, Waiting-times and returns in high-frequency financial data: an empirical study, Phys. A 314 (1–4) (2002) 749–755. [19] M. Stynes, J.L. Gracia, A finite difference method for a two-point boundary value problem with a Caputo fractional derivative, IMA J. Numer. Anal. 35 (2) (2015) 689–721. [20] M. Stynes, J.L. Gracia, Boundary layers in a two-point boundary value problem with a Caputo fractional derivative, Comput. Method Appl. Math. 15 (1) (2015) 79–95. [21] D. Willett, J.S.W. Wong, On the discrete analogues of some generalizations of Gronwall’s inequality, Monatshefte für Mathematik 69 (4) (1965) 362–367.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
A robust numerical method for a fractional differential equation Zhongdi Cen a,b, Anbo Le b,∗, Aimin Xu b a b
College of Finance and Trade, Ningbo Dahongying University, Ningbo, China Institute of Mathematics, Zhejiang Wanli University, Ningbo, China
a r t i c l e
i n f o
a b s t r a c t
MSC: 65L05 65R20 65L20
This paper is devoted to giving a rigorous numerical analysis for a fractional differential equation with order α ∈ (0, 1). First the fractional differential equation is transformed into an equivalent Volterra integral equation of the second kind with a weakly singular kernel. Based on the a priori information about the exact solution, an integral discretization scheme on an a priori chosen adapted mesh is proposed. By applying the truncation error estimate techniques and a discrete analogue of Gronwall’s inequality, it is proved that the numerical method is first-order convergent in the discrete maximum norm. Numerical results indicate that this method is more accurate and robust than finite difference methods when α is close to 0.
Keywords: Fractional differential equation Caputo fractional derivative Volterra integral equation Adapted mesh Convergence analysis
© 2017 Elsevier Inc. All rights reserved.
1. Introduction Fractional calculus has played a significant role in a variety of scientific and engineering fields, such as finance [12,18], control [11], viscoelasticity [6], and hydrology [1,13]. Those fractional models, described in the form of fractional differential equations, have been proved to be more appropriate for modeling many physical phenomena than the traditional integer order models, because fractional calculus enables the description of the memory properties of various materials and processes [17, Chapter 10]. The analytical solutions of most fractional differential equations can not be obtained, so approximate and numerical techniques must be used. In this paper, we consider the following fractional differential equation
Dα∗ u(x ) + f (x, u ) = 0,
x ∈ := (0, 1],
(1.1)
u (0 ) = γ ,
(1.2)
where Dα ∗ denotes the Caputo fractional derivative defined by
Dα∗ u(x ) =
1 (1 − α )
x 0
(x − t )−α u (t )dt ,
0 < α < 1,
(1.3)
the constant γ and function f are given. We assume that function f satisfies
¯ × R ∩ C 1 ( × R ). f (x, u ) ∈ C ∗
Corresponding author. E-mail address: [email protected] (A. Le).
http://dx.doi.org/10.1016/j.amc.2017.08.011 0 096-30 03/© 2017 Elsevier Inc. All rights reserved.
(1.4)
446
Z. Cen et al. / Applied Mathematics and Computation 315 (2017) 445–452
¯ (see e.g., [3, Under this hypothesis, the fractional differential equation (1.1) and a unique solution u(x ) ∈ C () (1.2) exists ¯ × R ∩ C q ( × R ) for some integer q ≥ 1 there exists a Theorem 6.8]). In [3, Theorem 6.28] it is shown that if f (x, u ) ∈ C positive constant C0 such that
( j) u (x ) ≤ C0 ,
if j = 0,
C0 xα − j ,
if j = 1, 2, . . . , q + 1,
(1.5)
for all x ∈ . These bounds indicate that the derivatives of the solution u may blow up at the end point x = 0. This singular behavior complicates the construction of the discretization scheme and the convergence analysis of the numerical method. Since analytic solutions for fractional differential equations are generally impossible to attain, various numerical methods for solving the fractional differential equations have been proposed. There are some papers which take account of the possibly singular behavior of solutions of the fractional differential equations. Numerical schemes based on the collocation methods are developed in [4,7–10,14–16], finite difference schemes are proposed in [2,5,19,20]. Because of the presence of the singular behavior, standard numerical methods fail to give accurate approximation, even for high-order methods. In this paper we first transform the fractional order initial value problem (1.1) and (1.2) into an equivalent Volterra integral equation of the second kind with a weakly singular kernel. Then we discretize the Volterra integral equation on an apriori chosen adapted mesh. A rigorous analysis about the convergence of the discretization scheme is given by taking account of the possibly singular behavior of the solution. Applying the truncation error estimate techniques and a discrete analogue of Gronwall’s inequality [21], we will prove that the scheme is first-order convergent in the discrete maximum norm. Numerical results are given to display that this method is more accurate and robust than finite difference methods when α is close to 0. The rest of the paper is organized as follows. The discretization scheme is described in Section 2. Convergence analysis of the scheme is given in Section 3. Finally, numerical experiments are presented in Section 4. Notation. Throughout the paper, C will denote a generic positive constant that is independent of the mesh. Note that C ¯ . is not necessarily the same at each occurrence. To simplify the notation we set gi = g(xi ) for any function g ∈ C () 2. Discretization scheme Based on the properties of the exact solution u(x) we construct an apriori mesh N ≡ {0 = x0 < x1 < · · · < xN = 1} with the mesh points
⎧ r i ⎪ ⎪ , ⎨ N xi = 1 ⎪ ⎪ ⎩ 1 + 2 1 − 2r ( i − N ) , r 2
N
0≤i≤
N , 2
(2.1)
N < i ≤ N, 2
2
where the discretization parameter N is a positive even integer, r = α1 > 1. The mesh points (2.1) are more densely at the region near x = 0 since the exact solution of the fractional differential equations (1.1) and (1.2) may be singular near x = 0; r away from x = 0 a uniform mesh is used. Pedas et al. [7,14–16] used the mesh xi = Ni for 0 ≤ i ≤ N, which may be coarse when i close to N. Our mesh is a slight modification of that in [7,14–16]. It is well known that the initial value problem (1.1) and (1.2) can be written as the following equivalent Volterra integral equation of the second kind with a weakly singular kernel [3, Lemma 6.2]
u (x ) = u (0 ) −
1
(α )
x 0
(x − t )α−1 f (t , u(t ) )dt,
u (0 ) = γ .
(2.2) (2.3)
In the following we will discrete this integral equation instead of the differential equations (1.1) and (1.2). An approximation to the integral can be obtained by a quadrature formula
xi 0
≈
(xi − t )α−1 f (t , u(t ) )dt =
i
f ( xk , uk )
k=1
=
i 1
α
i
k=1
xk xk−1
xk
xk−1
(xi − t )α−1 f (t , u(t ) )dt
(xi − t )α−1 dt
f (xk , uk ) (xi − xk−1 )α − (xi − xk )α .
k=1
Then, we have the following discretization scheme for the problem (2.2) and (2.3):
uNi = uN0 − uN0 = γ .
i
1 f xk , uNk (xi − xk−1 )α − (xi − xk )α , 1 ≤ i ≤ N, (α + 1 )
(2.4)
k=1
(2.5)
Z. Cen et al. / Applied Mathematics and Computation 315 (2017) 445–452
447
3. Convergence analysis In this section we give an apriori error estimate for the discretization scheme (2.4) and (2.5) on the apriori mesh (2.1). Let ziN = uN − ui , where uN is the solution of problem (2.4), (2.5) and ui is the solution of problem (2.2) and (2.3) at mesh i i point xi . Then, for the error ziN we have i
1 ak zkN (xi − xk−1 )α − (xi − xk )α = Ri , (α + 1 )
ziN +
1 ≤ i ≤ N,
(3.1)
k=1
z0N = 0,
(3.2)
where i
1
Ri =
(α ) −
xk xk−1
k=1
f (t , u(t ) )(xi − t )α −1 dt
i
1 f (xk , uk ) (xi − xk−1 )α − (xi − xk )α , (α + 1 )
(3.3)
k=1
∂f xk , uk + λk uNk − uk , ∂u
ak =
0 < λk < 1.
(3.4)
First we state a useful technical result, which is similar to the result in [19, Lemma 4.5], but under the nonuniform mesh
N . Lemma 3.1. For 2 ≤ i ≤ N there exists a constant C independent of N such that i
k=2
−1 ≤ CN −1 , (xi − xk−1 )α − (xi − xk )α hk xαk−1
(3.5)
where hk = xk − xk−1 is the mesh size. Proof. Let t denote the smallest integer that is greater than or equal to t for any t ∈ R. For the analysis of the estimates we distinguish two cases. Case I: 2 ≤ i ≤ N/2. By applying the mean value theorem we have
i/2 k=2
≤
−1 (xi − xk−1 )α − (xi − xk )α hk xαk−1
i/2 k=2
i/2 α−1
−1 −1 α (xi − xk )α−1 h2k xαk−1 ≤ α xi − x i/2 h2k xαk−1
r
i N
=α
−
1−r
i ≤C N
≤ CN −2
k=2
k i
N
i/2 1 N2
i/2 k=2
i/2
k=2
r α−1 i/2 r k=2
2r−2 k N
k−1 N
k N
−
k−1 N
r 2
k−1 N
r (α−1)
1−r
r−1 ≤ CN −1 ,
(3.6)
where we have used rα = 1 and k ≤ 2(k − 1 ) for k ≥ 2. Furthermore, we have
i
k= i/2 +1
≤
−1 (xi − xk−1 )α − (xi − xk )α hk xαk−1
max
i/2 +1≤k≤i
−1 hk xαk−1
i
(xi − xk−1 )α − (xi − xk )α
k= i/2 +1 α α −1 ≤ x i/2 xi − x i/2 max hk
i/2 +1≤k≤i
448
Z. Cen et al. / Applied Mathematics and Computation 315 (2017) 445–452
≤
≤
i/2
r (α−1) rα i N
N i 2N
1−r r−1 i N
r max
i N
i/2 +1≤k≤i
k N
k−1 N
−
r
1 N
≤ CN −1 ,
(3.7)
where we also have used rα = 1. Therefore, from (3.6) and (3.7) we conclude that the lemma holds true for Case I. Case II: N/2 < i ≤ N. By using the same method as that in Case I we also can obtain
i/2 k=2
−1 ≤ CN −1 . (xi − xk−1 )α − (xi − xk )α hk xαk−1
(3.8)
Furthermore, we have
i
k= i/2 +1
≤
−1 (xi − xk−1 )α − (xi − xk )α hk xαk−1
max
i/2 +1≤k≤i
−1 hk xαk−1
= xi − x i/2 ≤ xαi ·
α
·
max
i/2 +1≤k≤i
i
max
i/2 +1≤k≤i
−1 hk xαk−1
−1 hk xαk−1 .
(3.9)
rα r
≤r
i N
i N
k N
r−1 k N
k= i/2 +1
For i/2 + 1 ≤ k ≤ N/2 we have −1 xαi hk xαk−1 ≤
(xi − xk−1 )α − (xi − xk )α
1 N
− k 2N
k−1 N
r
1−r
k−1 N
r (α−1)
≤ CN −1 ,
(3.10)
where we also have used rα = 1 and k ≤ 2(k − 1 ) for k ≥ 2. For N/2 < k ≤ i we have
N/2 −1 xαi hk xαk−1 ≤ xαi 2(1 − 2−r )N −1 N
r (α−1)
≤ CN −1 ,
(3.11)
where we also have used rα = 1. Therefore, from (3.8)–(3.11) we conclude that the lemma also holds true for Case II.
The next lemma gives us a useful formula for the truncation error. Lemma 3.2. Let f satisfy the hypothesis (1.4) and
∂ f (x, u ) ≤ Cxα−1 for x ∈ . ∂x
(3.12)
Then the truncation errors of the discrete scheme (2.4) and (2.5) satisfy
|Ri | ≤ CN−1 , 1 ≤ i ≤ N,
(3.13)
where C is a positive constant independently of N. Proof. For the analysis of the truncation errors we distinguish two cases. Case 1: i = 1. From (3.3) we have
|R1 | ≤
1
(α )
x1 0
(x1 − t )α−1 | f (t , u(t ) )|dt +
≤ C xα1 = C N −1 ,
1
(α + 1 )
| f (x1 , u1 )|xα1 (3.14)
where we have used the hypothesis (1.4), the bound on u given by (1.5) and rα = 1. From this we know that the lemma holds true for Case 1.
Z. Cen et al. / Applied Mathematics and Computation 315 (2017) 445–452
449
Case 2: 1 < i ≤ N. From (3.3) we have i
1
Ri =
(α ) −
xk xk−1
k=1
f (t , u(t ) )(xi − t )α −1 dt
i
1 f (xk , uk ) (xi − xk−1 )α − (xi − xk )α (α + 1 ) k=1
= R1i + R2i , where
R1i =
R2i =
1
(α ) −
(3.15) x1
0
f (t , u(t ) )(xi − t )α −1 dt
1
i
1
(α )
k=2
f (x1 , u1 ) xαi − (xi − x1 )α ,
(α + 1 )
xk xk−1
(3.16)
f (t , u(t ) )(xi − t )α −1 dt
i
1 − f (xk , uk ) (xi − xk−1 )α − (xi − xk )α . (α + 1 )
(3.17)
k=2
By using the similar method as that in Case 1, from (3.16) we have
1 Ri ≤
1
(α )
x1 0
| f (t , u(t ) )|(xi − t )α−1 dt
1 | f (x1 , u1 )| xαi − (xi − x1 )α (α + 1 ) ≤ C xαi − (xi − x1 )α ≤ Cx1 (xi − x1 )α −1 +
≤ C xα1 = C N −1 ,
(3.18)
where we also have used the hypothesis (1.4), the bound on u given by (1.5) and rα = 1. From (3.17) we have
2 Ri = ≤
i xk α −1 dt [ f (t , u(t ) ) − f (xk , uk )](xi − t ) (α ) xk−1 k=2 1
i
1
(α )
| f (ξk , u(ξk ) ) − f (xk , uk )|
k=2
xk
xk−1
(xi − t )α−1 dt
by the mean value theorem for integrals with ξk ∈ (xk−1 , xk ). Hence, we can obtain
2 Ri ≤
i
1 | f (ξk , u(ξk ) ) − f (xk , uk )| · (xi − xk−1 )α − (xi − xk )α (α + 1 ) k=2
1
i
{[| f (ξk , u(ξk ) ) − f (ξk , uk )| + | f (ξk , uk ) − f (xk , uk )|] (α + 1 ) k=2 · (xi − xk−1 )α − (xi − xk )α i ∂ f
1 (ξk , μk ) · |u(ξk ) − uk | ≤ (α + 1 ) ∂u k=2 ∂ f +hk (η , u ) (xi − xk−1 )α − (xi − xk )α ∂x k k i ∂ f
α α ≤C (xi − xk−1 ) − (xi − xk ) hk u (ζk ) + (ηk , uk ) ∂x ≤
k=2
≤C
i
k=2
−1 , (xi − xk−1 )α − (xi − xk )α hk xαk−1
450
Z. Cen et al. / Applied Mathematics and Computation 315 (2017) 445–452
where we again have used the mean value theorem with ηk , ζk ∈ (xk−1 , xk ) and μk ∈ uk−1 , uk , the bound on u given by (1.5) and the hypotheses (1.4) and (3.12). Now, an appeal to Lemma 3.1 yields
2 Ri ≤ CN−1 .
(3.19)
Therefore, combining (3.18) and (3.19) with (3.15) we conclude that the lemma also holds true for Case 2.
Now we can get the following a priori error estimate for the discrete scheme (2.4) and (2.5). Theorem 3.3. Let u be the solution of (2.2), (2.3) and uN be the solution of the discrete scheme (2.4) and (2.5) on the mesh N . Then, under the hypotheses (1.4) and (3.12) we have
N ui − u(xi ) ≤ CN−1 , 1 ≤ i ≤ N
for sufficiently large N, where C is a positive constant independent of N. Proof. From (3.1) we have i−1 N
zi ≤ pi + wi vk zkN , 1 ≤ i ≤ N,
(3.20)
k=1
where
(α + 1 )|Rk | 1 , wk = , (α + 1 ) + ak hα (α + 1 ) + ak hα k k vk = |ak | (xi − xk−1 )α − (xi − xk )α , 1 ≤ k ≤ i. pk =
By applying the discrete analogue of Gronwall’s inequality [21, Theorem 3], we have
i−1 i−1 k N −1 zi ≤ pi + wi 1 + v jw j · vk pk 1 + v jw j , 1 ≤ i ≤ N. k=1
j=1
(3.21)
j=1
It is easy to conclude that for sufficiently large N
0 < pk ≤ CN −1
0 < wk ≤ C for 1 ≤ k ≤ i,
and
(3.22)
where we have used Lemma 3.2. Hence we can obtain i−1
k=1
vk pk
k
1 + v jw j
−1
≤ CN −1
i−1
vk
k=1
j=1
≤ CN −1
i−1
(xi − xk−1 )α − (xi − xk )α
k=1
= CN −1 xαi − (xi − xi−1 )α
≤ CN −1 , and i−1
1 + v j w j ≤ exp
j=1
i−1
(3.23)
v jw j
j=1
≤ exp C
i−1
xi − x j−1
α
j=1
= exp C xαi − (xi − xi−1 )α ≤ C,
− xi − x j
α
(3.24)
where we have used the inequalities (3.22) and (1 + y ) ≤ ey for y ≥ −1. Therefore, substituting (3.22)–(3.24) into (3.21) we have
N zi ≤ CN−1 , 1 ≤ i ≤ N,
which completes the proof.
Z. Cen et al. / Applied Mathematics and Computation 315 (2017) 445–452
451
4. Numerical experiments In this section, we verify experimentally the theoretical results obtained in the preceding section. Error estimates and convergence rates for the discrete scheme are presented for the following example which has been presented in [2]. We will give the comparisons of numerical solution on various schemes. Example. We consider the following fractional differential equation
Dα u(x ) + f (x, u(x )) = 0,
x ∈ (0, 1],
u ( 0 ) = 1, where 0 < α < 1, f (x, u(x )) = 2u(x ) + sin(u(x )) + 0.1u2 (x ) + s(x ), and s(x) is chosen such that the exact solution is
u ( x ) = xα + 2x + 1.
Table 1 Error estimates and convergence rates on an a priori adapted mesh for example.
α = 0.1 α = 0.2 α = 0.3 α = 0.4 α = 0.5 α = 0.6 α = 0.7 α = 0.8 α = 0.9
N = 64
N = 128
N = 256
N = 512
N = 1024
N = 2048
4.6534e−3 0.857 7.2937e−3 0.896 8.0308e−3 0.927 7.5963e−3 0.950 6.6375e−3 0.968 5.5351e−3 0.981 4.4872e−3 0.985 3.8578e−3 0.930 3.0207e−3 0.860
2.5688e−3 0.878 3.9183e−3 0.916 4.2241e−3 0.943 3.9311e−3 0.964 3.3922e−3 0.978 2.8051e−3 0.988 2.2665e−3 0.991 2.0243e−3 0.949 1.6646e−3 0.884
1.3976e−3 0.895 2.0770e−3 0.930 2.1966e−3 0.955 2.0155e−3 0.973 1.7216e−3 0.985 1.4144e−3 0.992 1.1400e−3 0.995 1.0488e−3 0.963 9.0217e−4 0.903
7.5173e−4 0.908 1.0899e−3 0.942 1.1328e−3 0.965 1.0267e−3 0.980 8.6974e−4 0.990 7.1100e−4 0.995 5.7204e−4 0.997 5.3808e−4 0.973 4.8234e−4 0.919
4.0062e−4 0.919 5.6737e−4 0.951 5.8048e−4 0.972 5.2059e−4 0.985 4.3801e−4 0.993 3.5672e−4 0.997 2.8663e−4 0.998 2.7404e−4 0.981 2.5504e−4 0.932
2.1189e−4 − 2.9350e−4 − 2.9602e−4 − 2.6305e−4 − 2.2011e−4 − 1.7875e−4 − 1.4350e−4 − 1.3883e−4 − 1.3364e−4 −
Table 2 Comparisons of our numerical method with three finite difference schemes on error estimates and convergence rates for example.
α = 0.05
Our method FDM1 FDM2 FDM3
α = 0.1
Our method FDM1 FDM2 FDM3
α = 0.15
Our method FDM1 FDM2 FDM3
N = 64
N = 128
N = 256
N = 512
N = 1024
N = 2048
2.5766e−3 0.834 1.3992e−2 0.033 6.2556e−2 −0.015 1.1971e−2 −0.051 4.6534e−3 0.857 2.3243e−2 0.060 5.8221e−3 −0.538 2.0039e−2 0.105 6.2386e−3 0.878 2.8934e−2 0.088 2.1720e−3 1.080 2.3506e−2 0.163
1.4450e−3 0.856 1.3676e−2 0.030 6.3214e−2 0.007 1.2401e−2 0.053 2.5688e−3 0.878 2.2302e−2 0.057 8.4544e−3 −0.080 1.8637e−2 0.107 3.3946e−3 0.898 2.7224e−2 0.087 1.0272e−3 −0.717 2.0995e−2 0.172
7.9855e−4 0.873 1.3393e−2 0.029 6.2905e−2 −0.035 1.1955e−2 0.053 1.3976e−3 0.895 2.1437e−2 0.056 8.9343e−3 −0.043 1.7300e−2 0.111 1.8213e−3 0.914 2.5635e−2 0.087 1.6887e−3 0.658 1.8635e−2 0.182
4.3613e−4 0.886 1.3129e−2 0.028 6.4469e−2 −0.134 1.1524e−2 0.053 7.5173e−4 0.908 2.0619e−2 0.056 9.2016e−3 0.010 1.6023e−2 0.114 9.6656e−4 0.926 2.4132e−2 0.089 1.0706e−3 −0.191 1.6431e−2 0.191
2.3593e−4 0.898 1.2876e−2 0.028 7.0739e−2 −1.480 1.1107e−2 0.054 4.0062e−4 0.919 1.9837e−2 0.056 9.1365e−3 0.002 1.4801e−2 0.118 5.0855e−4 0.937 2.2695e−2 0.090 1.2218e−3 −0.036 1.4389e−2 0.201
1.2661e−4 − 1.2632e−2 − 1.9734e−1 − 1.0703e−2 − 2.1189e−4 − 1.9081e−2 − 9.1248e−3 − 1.3635e−2 − 2.6570e−4 − 2.1316e−2 − 1.2529e−3 − 1.2518e−2 −
452
Z. Cen et al. / Applied Mathematics and Computation 315 (2017) 445–452
0 The nonlinear problem (2.4) and (2.5) is solved by using Newton’s method with uN, = 1.5(i = 1, 2, . . . , N ) as an initial i guess. We iteratively compute uN, j until
j j−1 max uN, − uN, < 10−6 . i i
1≤i≤N
We measure the accuracy in the discrete maximum norm
eN = max uNi − ui , 1≤i≤N
and the convergence rate
N
r = log2
eN e 2N
.
The numerical results for example are tabulated in Table 1. Table 1 shows that the computed solution converges to the exact solution with first-order accuracy for sufficiently large N and the numerical results do not depend strongly on the value of α , which supports the convergence estimate of Theorem 3.3. For comparison we also use the numerical methods based on the finite difference method on the uniform mesh (FDM1), the finite difference method on the a priori mesh N (FDM2) and the finite difference method on an a posteriori adapted mesh [2] (FDM3) to compute this example in Table 2. From Table 2 we confirm that our method proposed in this paper is more accurate and robust than finite difference methods when α is close to 0. Acknowledgment We would like to thank the anonymous reviewers for some suggestions for the improvement of this paper. The work was supported by Humanities and Social Sciences Planning Fund of Ministry of Education of China (Grant No. 14YJC790 0 06), Project of Philosophy and Social Science Research in Zhejiang Province (Grant No. 15NDJC243YB), Zhejiang Province Natural Science Foundation (Grant No. Y17D010024), Ningbo Municipal Natural Science Foundation (Grant Nos. 2017A610131, 2017A610140) and Ningbo Municipal Soft Science Foundation (Grant No. 2015A10045). References [1] P. Becker-Kern, M.M. Meerschaert, H.P. Scheffler, Limit theorem for continuous-time random walks with two scales, J. Appl. Prob. 41 (2) (2004) 455–466. [2] Z. Cen, A. Le, A. Xu, A posteriori error analysis for a fractional differential equation, Int. J. Comput. Math. 94 (6) (2017) 1185–1195. [3] K. Diethlm, The analysis of fractional differential equations, in: Lecture Notes in Mathematics, 2004, Springer, Berlin, 2010. [4] N.J. Ford, M.L. Morgado, M. Rebelo, A nonpolynomial collocation method for fractional terminal value problems, J. Comput. Appl. Math. 275 (2015) 392–402. [5] J.L. Gracia, M. Stynes, Central difference approximation of convection in caputo fractional derivative two-point boundary value problems, J. Comput. Appl. Math. 273 (2015) 103–115. [6] R.C. Koeller, Application of fractional calculus to the theory of viscoelasticity, J. Appl. Mech. 51 (2) (1984) 229–307. [7] M. Kolk, A. Pedas, E. Tamme, Modified spline collocation for linear fractional differential equations, J. Comput. Appl. Math. 283 (2015) 28–40. [8] N. Kopteva, M. Stynes, An efficient collocation method for a caputo two-point boundary value problem, BIT Numer. Math. 55 (2015) 1105–1123. [9] H. Liang, M. Stynes, Collocation methods for general caputo two-point boundary value problems, https://www.researchgate.net/publication/313878934. [10] X. Ma, C. Huang, Numerical solution of fractional integro-differential equations by a hybrid collocation method, Appl. Math. Comput. 219 (12) (2013) 6750–6760. [11] J.T. Machado, Discrete time fractional-order controllers, Fract. Calc. Appl. Anal. 4 (1) (2001) 47–66. [12] M.M. Meerschaert, E. Scalas, Coupled continuous time random walks in finance, Phys. A 370 (1) (2006) 114–118. [13] M.M. Meerschaert, Y. Zhang, B. Baeumer, Particle tracking for fractional diffusion with two time scales, Comput. Math. Appl. 59 (3) (2010) 1078–1086. [14] A. Pedas, E. Tamme, Piecewise polynomial collocation for linear boundary value problems of fractional differential equations, J. Comput. Appl. Math. 236 (13) (2012) 3349–3359. [15] A. Pedas, E. Tamme, Numerical solution of nonlinear fractional differential equations by spline collocation methods, J. Comput. Appl. Math. 255 (2014) 216–230. [16] A. Pedas, E. Tamme, Spline collocation for nonlinear fractional boundary value problems, Appl. Math. Comput. 244 (2014) 502–513. [17] I. Podlubny, Fractional Differential Equations, Academic Press, New York, 1999. [18] M. Raberto, E. Scalas, F. Mainardi, Waiting-times and returns in high-frequency financial data: an empirical study, Phys. A 314 (1–4) (2002) 749–755. [19] M. Stynes, J.L. Gracia, A finite difference method for a two-point boundary value problem with a Caputo fractional derivative, IMA J. Numer. Anal. 35 (2) (2015) 689–721. [20] M. Stynes, J.L. Gracia, Boundary layers in a two-point boundary value problem with a Caputo fractional derivative, Comput. Method Appl. Math. 15 (1) (2015) 79–95. [21] D. Willett, J.S.W. Wong, On the discrete analogues of some generalizations of Gronwall’s inequality, Monatshefte für Mathematik 69 (4) (1965) 362–367.