1.
Introduction
Fractional calculus is a branch of study that of any order integral or derivative. It is obtained by replacing the integer order of differential equations with the fractional order. Because fractional derivative has heredity or memory, the fractional differential equation is widely used in various fields. In the past few decades, fractional derivative has been studied by a large number of scholars. They firstly proved a (3−θ)-order convergence and stability scheme for time-fractional derivative θ under the time step k≥2 in [1]. The similar (3−θ)-order scheme to efficiently solve the time-fractional diffusion equation in time was first constructed in [2]. Based on the idea of block-by-block approach method, the reference [3] presented a general technique to construct high order schemes for the numerical solution of the fractional ordinary differential equations. The reference [4] established a fractional Gronwall inequality to prove some stability and convergence estimates of schemes for fractional reaction-subdiffusion problems. In [5], a space-time finite element method was used to solve the distributed-order time fractional reaction diffusion equations. The scheme with (3−θ) convergence order for time step k≥2 and (2−θ) convergence order for time step k=1 was constructed for time-fractional derivative in [6]. They proved the (2−θ)-order convergence and stability of the time scheme, where θ is the order of the time-fractional derivative in [7]. A series of high order numerical approximations to θ-Caputo derivatives (0<θ<2) and Riemann-Liouville derivatives were given in [8]. They investigated numerical solutions of coupled nonlinear time-fractional reaction-diffusion equations obtained by the Lie symmetry analysis in [9]. They built a robust finite difference scheme for a nonlinear Caputo fractional differential equation on an adaptive grid with first-order convergent result in [10]. Spectral methods was popularly used to solve fractional derivative in [11,12,13]. An implicit numerical method was constructed to solve the two-dimensional θ(1<θ<2) time fractional diffusion-wave equation with the convergence order 3−θ in [14]. They used the piecewise linear and quadratic Lagrange interpolation functions to construct a (3−θ)-order numerical approximate method for the Caputo fractional derivative in [15]. In [16], a fast high order numerical method for nonlinear fractional-order differential equations with non-singular kernel was constructed. The numerical scheme for nonlinear fractional differential equations was constructed by the Galerkin finite element method in [17]. The reference [18] established a hp-version Chebyshev collocation method for nonlinear fractional differential equations. The Spectral collocation method for nonlinear Riemann-Liouville fractional differential equations was given in [19]. The above high order schemes are either implicit or low-order schemes are used in the first step or obtained by discretizing the equivalent form of the original equation. Therefore, we will give a high order numerical scheme with uniform convergence order in this paper using the direct numerical discretization of original equations based on the idea of [1,2,3].
This paper is arranged as follows: In Section 2 the higher order numerical scheme is proposed. And the local truncation error of the constructed higher order numerical scheme is given in Section 3. The convergence analysis is studied in Section 4. In Section 5 some numerical examples are given. Finally, some conclusion are given in Section 6.
2.
A high order finite difference scheme for Caputo derivative
We consider the following initial value problem of nonlinear fractional ordinary differential equations
where θ is the order of the fractional derivative, 0Dθtu(t) in (2.1) is defined as the Caputo fractional derivatives of order θ given by
with ω1−θ is defined by ω1−θ (t)=1tθΓ(1−θ), and Γ(⋅) denotes Gamma function [20].
First, the finite difference scheme is introduced to discretize the fractional derivative. Dividing the interval [0,T] into K equal subintervals, set tk=kΔt,k=0,1,…,K, with Δt=TK, the numerical solution of (2.1) at tk is uk, and 0Dθtu(tk)=0Dθtuk, fk=f(tk,uk).
Next, we start discretizing the fractional derivative (2.1). Firstly, the values of u(t) on t1 and t2 are determined. Using Quadratic Lagrange interpolation, the approximation formula of u(t) on the interval [t0, t2] is
where ϖi,0(t),i=0,1,2, are the Quadratic Lagrangian interpolation basis functions on point t0, t1 and t2, are defined as following
Let Δuk=uk−uk−1,Δu0=u0,∀k≥1. when k=1,2, we use 0Dθt(J[t0,t2]u)(t1), 0Dθt(J[t0,t2]u)(t2) to approximate 0Dθtu(t1), 0Dθtu(t2), respectively, then we get
where
When k≥3, we have
we can conclude that
where
According to approximations (2.5), (2.6) and (2.8), the higher order numerical scheme of Eq (2.1) corresponding to the initial value (2.2) is as follows
Remark 2.1. The numerical algorithm in this paper is different from [15]. In [15], in order to obtain (3−θ)-order numerical scheme for the linear form of (2.1), they changed (2.1) into the following equation
And using the linear Lagrange interpolation in [t0,t1] to u(s), they obtain a (3−θ)-order numerical scheme. In this paper, we use the Quadratic Lagrange interpolation in [t0,t2] to u(s) in (2.1) and directly obtain the (3−θ)-order numerical scheme (2.5) and (2.6).
Next, we will introduce a lemma to prove the sign of Akk−j and judge its size relation.
Lemma 2.1. Ak0>0, the sign of Ak1 is uncertain, Ak2>Ak3>…>Akk−1>0.
Proof. Through calculation, we can get the following results
The sign of Ak1 is determined by the sign of M(θ), which satisfies
Therefore M(θ) is a decrease function on θ∈(0,1). By M(0)=2>0 and M(1)=−12<0, we know that M(θ) has only one zero θ0 in (0,1), and M(θ)>0 if θ∈(0,θ0) and M(θ)<0 if θ∈(θ0,1), which agrees with the conclusion of Lemma 2.1.
For j=3,…,k−2, we have
let ˉj=k−j, ˉj≥2, then using Taylor expansion, we can get the following results
where ai=Πin=0(1−θ−n)1(i+2)![(−1)i(−i2)+i+42]. Next, we will prove ˉj−θ∑+∞i=0aiˉj−i is decrease function on ˉj. Because ∑+∞i=0aiˉj−i is a convergence series, and the radius of convergence exists, it can exchange the derivative and the sum.
We can find ∑+∞i=2(i+θ)aiˉj−i−1 is an alternating series with positive first term. So ∑+∞i=2(i+θ)aiˉj−i−1>0, we have
Therefore, we have already proved ˉj−θ∑+∞i=0aiˉj−i is decrease function on ˉj. According to (2.15), we obtain
Next, let's prove Akk−3>Akk−2>Akk−1>0. Through careful calculation, we can get
where
According to (1) in Lemma A.1, we get F1>0,F2>0. Therefore, we have
The following step, we will prove Akk−1>0. By carefully calculation, we have
According to (2.33), we have
Combining (2.12), (2.13), (2.16), (2.17) with (2.19), we already proved Lemma 2.1.
We let Δˉuj=Δ(uj−ρuj−1), Δˉu0=Δu0. Let ρ=12(1−Ak1Ak0)=12(3+θ−64−θ⋅21−θ ). Therefore, from (2.8), we can get that
where
Next, we discuss the relationship between the size of ˉAkk−j.
Lemma 2.2. ˉAk0>ˉAk1>ˉAk2>…>ˉAkk−1>0.
Proof. When k=3, we only need to prove ˉA30>ˉA31>ˉA32>0. By careful calculation, we can get
According to (2) in Lemma A.1, we get
By direct calculation, we have
According to (3) in Lemma A.1, we have
Next, we will prove ˉA32=A32+ρˉA31>0. According to (2.41), we get ˉA31>0. Hence, we just need to prove A32>0. By calculation, we have A32=12ΔtθΓ(3−θ)[4−θ−θ⋅32−θ]>0, therefore,
Combining (2.22), (2.23) and (2.24), when k=3, we know that Lemma 2.2 holds.
When k≥4, through careful calculation, we can draw the following conclusions
According to (4) in Lemma A.1, we can get ρ>0, we have ˉAk0>ˉAk1. Because
where ˜f(θ)=−3(4−θ)2+6(4−θ)(6−θ)21−θ −4(4−θ)(8−θ)31−θ +(6−θ)2⋅41−θ . According to (5) in Lemma A.1, we can get ˜f(θ)>0. It can be concluded that
By careful calculation, we can get ˉAk3−ˉAk2=Ak3−Ak2+ρ(ˉAk2−ˉAk1). According to ρ>0, (2.26) and Lemma 2.1, we can get
By mathematical induction, and Lemma 2.1 we can conclude that
So we have
Next we prove that ˉAkk−1>0. According to (2.21), we have
According to (2.37), we get ˉAk1>0. Using Lemma 2.1, ρ>0, we have
Combining (2.28) with (2.29), we already proved Lemma 2.2.
Lemma 2.3. There is a constant πA≥6 such that the discrete kernel satisfies the lower bound.
Proof. When k≥4, According to (2.21), Lemma 2.1 and Lemma 2.2, we can get
For j=1, we have
let ˜j=k−1,˜j≥3, using Taylor formula to expand the calculation, we can get the following results
where ai=Πii=0(1−θ−i)1(i+2)!(1˜j)i[(−1)i+1i2+3i+42], a0=(1−θ), a1=8(1−θ)(−θ)12˜j, ∑+∞i=2ai is a convergent alternating series, and a2>0, therefore, we have 0≤∑+∞i=2ai≤a2, so
We have
For j=2, we have
let ˆj=k−2,ˆj≥2, and then using Taylor expansion, similar to j=1, it can be obtained by calculation
where ai=Πin=0(1−θ−n)1(i+2)!(1ˆj)i[i−2⋅(−1)i+2−i22i+1], because ∑+∞i=3ai is a convergent staggered series, and a3>0, 0<∑+∞i=3ai<a3, we get
therefore, we have
For j=3,4,…,k−2, we have
let ˉj=k−j, ˉj≥2, we have
then using Taylor expansion, we can get the following results
where ai=Πin=0(1−θ−n)1(i+2)!(1ˉj)i[(−1)i(−i2)+i+42]. Because ∑+∞i=2ai is a convergent alternating series and a2>0,0≤∑+∞i=2ai≤a2, so
so we can get
When j=k−1,
Therefore, we have
When j=k, according to (2.30), we can get
We have
When k=3, ˉA30 is in (2.39), so 1πA≤1. According to (2.30), we only need
Therefore, we have
Next, we calculate ˉA32,
Therefore, we have
Combining (2.33)–(2.38) with (2.12)–(2.43), we have already proved Lemma 2.3.
3.
Estimation of the truncation errors
Now we turn to derive an estimate for the truncation errors of the scheme (2.10). We start with deriving an error estimate for the finite difference operator 0DθΔtuk.
Theorem 3.1. Assume u(t)∈C3[0,T]. Let
Then there exists a constant Cu depending only on the function u, such that for all Δt>0,
Proof. Our error estimation will be established on the following Taylor theorem
where ξ(τ) is a function defined on [t0,t2] with range (t0,t2). We first estimate r1(Δt)
This proves (3.2) for k=1. The case k=2 can be similarly proven, and here we omit the details.
When k≥3, we have
Similar to the proof of |r1(Δt)|, we have
For N, we have
In the above inequality, we use |(τ−tj−1)(τ−tj)(τ−tj+1)|≤2√39Δt3.
We can get the following conclusions
Complete the proof of Theorem 3.1.
4.
Convergence analysis
In order to obtain the convergence analysis, we now introduce an important tool: complementary discrete convolution kernel. Because of ωθ∗ωβ=ωθ+β, therefore
A class of complementary discrete convolution kernels Pkk−i with the same properties are given
According to [3], we have
where Pki>0,i=0,1,...,k−3 was obtained from Lemma 2.2.
Next, we introduce Lemma 2.1 of [3], as follows.
Lemma 4.1. In (4.3), the discrete kernel Pki has the property of (4.2) and satisfies the following conditions
Now we can analyze the convergence of 0DθΔtuj. First, we consider f(t,u) and assume that it satisfies the following differential mean value theorem and that there exists a constant L, such that
and |Lk|≤L,∀k≥1.
Theorem 4.1. Assumes that u is the exact solution of (2.1) and (2.2), and {uk}Kk=1 is the numerical solution of (2.10). If θ∈(0,1) and step Δt satisfy
where πA≥6, and Λ=12L, then the following error estimates hold,
here C is only dependent on the function u.
Proof. Because Δu1=u1−u0,Δu2=u2−u1, through calculation, we can get
let −B11=Δt−θˉD10, B11−B21=Δt−θˉD11, B21=Δt−θˉD12; −B12=Δt−θD10, B12−B22=Δt−θD11, B22=Δt−θD12. It is easy to obtain by calculation
Among (4.6), then we can get
That is
where D=(ˉD11ˉD12D11D12). Since ˉD11D12−ˉD12D11 is bounded, then
therefore,
when j=1,2, (4.8) is proved.
When i≥3, let's set ˉei=ei−ρei−1,i=1,...,j, and ˉe0=e0 for (4.13), we have
Combining with (2.20) and (4.6), ˉej, j≥3, satisfy
Because of ej=∑jn=1ρj−nˉen, then
where
by multiplying 2ˉek on both sides of (4.16), from Lemma A.1 in Liao [3], we can obtain that
To sum up, we can get the following conclusions
We change k to i, and then multiply both sides of (4.18) by a ∑ki=3Pkk−i, then
From (4.2) and Δuj=uj−uj−1, we can get the following results
substituting (4.20) into (4.19) yields
Next, we will prove it by mathematical induction
where Λ=12L, Eθ stands for the Mittag-Leffler function, we usually define it as: Eθ(z)=∑∞j=0ZjΓ(1+jθ), and Eθ(0)=1,∀Z∈R, and Z>0,E′θ(Z)>0, so Fk≥Fk−1≥2 for all k≥2. where Gk=|ˉe2|+2max3≤j≤k∑ki=3Pjj−i|ˉri(Δt)|.
When k=3, (i) if |ˉe3|≤|e2|,G3=|e2|+2P30|ˉr3(Δt)|≥|ˉe2|, we have |ˉe3|≤|ˉe2|≤G3≤F3G3,
(ii) if |ˉe3|>|e2|, from (4.21), it can be concluded that
According to (4.4), it can be concluded that
so |ˉe3|≤2(|ˉe2|+2P30|ˉr3(Δt)|)≤2G3≤F3G3.
When 4≤k≤K, we assume that: |ˉej|≤FjGj, for 3≤j≤k−1, let |ˉej(k)|=max2≤i≤k−1|ˉei|.
(i) When |ˉek|≤|ˉej(k)|, because Fj and Gj monotonically increase with j, so
(ii) When |ˉek|≥|ˉej(k)|, from (4.21) we can obtain that
From (4.4), we give the limit of the maximum step size, which means that
Using of Lemma 2.3 in [3] for any real μ>0,
We have
To sum up, the proof of (4.22) is complete.
Next, we carefully estimate |ˉri(Δt)|. According to (4.17), we have
Next, we estimate the upper bounds of ˉA(k)k−1 and ˉA(k)k−2. According to (2.37), Lemma 2.1 and (4) in Lemma A.1, we can obtain
Therefore, ˉAkk−2≤1094ΔtθΓ(3−θ). According to Lemma 2.2, we can get ˉAkk−1<ˉAkk−2≤1094ΔtθΓ(3−θ). Therefore, we can obtain from (4.23),
In the expression on the right side of (4.22),
according to (4.5) in Lemma 4.1, we can obtain that
Because of 1ω1−θ (ti)=Γ(1−θ)tθi, therefore (4.24) becomes
Based on (4.13), (4.25) and (4.22), we obtain
According to ek=ˉek−ρek−1=∑kn=0ρk−nˉen, we have |ek|≤3max0≤n≤k|ˉen|. Then the proof is completed.
5.
Numerical results
In this section, two examples are presented to verify the effectiveness of our numerical method (2.10).
Example 5.1. we consider the problem (2.1) with u(0)=0, and
The exact solution of Eq (2.1) is u(t)=t3+θ. We take T=1 and choose the step size to be Δt=2−l,l=7,8,…,11. The error we will display is defined by eΔt=maxk=1,2,⋯,K|u(tk)−uk|,K=T/Δt.
In (5.1), f is a linear case of u. From Table 1, the convergence order is computed by log2(e2Δt/eΔt). it is observed that for θ=0.3,0.6 and 0.9, the convergence rates are close to 2.7,2.4 and 2.1, respectively. This is in a good agreement with the theoretical prediction. In Table 2, we can take θ=0.01 and 0.99 and obtain that the convergence rates are close to 2.99 and 2.01. That tells us that as θ→0 or 1, the convergence rate still 3−θ. In (5.2), f is the nonlinear case of u. From Tables 3 and 4, once again these results confirm that the convergence of the numerical solution is close to of order 3−θ for 0<θ<1.
Example 5.2. we consider the problem (2.1) and (2.2) with the following right hand side function
The corresponding exact solution is u(t)=sin(t), and Sinα,β(t)=∑∞k=1(−1)k+1t2k−1Γ(α(2k−1)+β). f is a linear and nonlinear case of u in (5.3) and (5.4), respectively.
We repeat the same calculation as in Example 5.1 using the proposed numerical scheme. Tables 5 and 6 show the maximum errors and decay rates of the step size for several a ranging from 0.3 to 0.9. This is consistent with the theoretical prediction. The convergence rate is close to 3−θ for 0<θ<1.
6.
Conclusions
In this paper, we presented a high order numerical method for solving Caputo nonlinear fractional ordinary differential equations. The numerical method was constructed by using the Quadratic Lagrange interpolation. By careful error estimation, the proposed scheme is of order 3−θ for 0<θ<1. Finally, numerical experiments are carried out to verify the theoretical prediction. In the future, we plan to apply this kind of methods to 3D fractional partial differential equations with time derivatives based the idea of [21] and [22].
Acknowledgments
This research was supported by National Natural Science Foundation of China (Grant No. 11901135, 11961009), Foundation of Guizhou Science and Technology Department (Grant No. [2020]1Y015, [2017]1086), Foundation for graduate students of Guizhou Provincial Department of Education(Grant No.YJSCXJH[2020]136). The authors thank the anonymous referees for their valuable suggestions to improve the quality of this work significantly.
Conflict of interest
The authors declare that there is no conflict of interests regarding the publication of this article.
A.
Proof of some inequalities
Lemma A.1.
Proof. (1) Firstly, we will prove F1>0. let k−2=x, Taylor formula is used to obtain the following results
where an=12(n+1)⋅(−1)n−8(n+1)⋅2n(n+3)!=12(n+1)[(−1)n−16⋅2n](n+3)!<0, so ∑+∞n=0∏ni=0(−θ−i)(1x)n⋅an is an alternating series with positive first term, so satisfies ∑+∞n=0∏ni=0(−θ−i)(1x)n⋅an>0, so F1>0.
Next, we will prove F2>0. Let k−2=ˉx, Taylor formula is used to obtain the following results
where bn=32(n+1)(−1)n+1+2(n−1)⋅2n[1+(−1)n],b0=−112,b1=3, when n≥2, n is odd, bn>0; n is even, bn=(n−1)[4⋅2n−32]−3>[8−32]−3>0, so ∑+∞n=1∏ni=0(−θ−i)(1ˉx)n1(n+3)!bn is an alternating series with positive first term. So ∑+∞n=0∏ni=0(−θ−i)(1ˉx)n1(n+3)!bn>(−θ)13!b0+0>0. We get F2>0.
(2) Let f(θ)=32(4−θ)−(4+θ)31−θ+(6−θ)2−θ, through careful calculation, we get
f′′(θ)=3−θ⋅g(θ), where g(θ)=3(ln3)[2−(4+θ)ln3]+(32)θ(ln2)[2+(6−θ)ln2], because
We can get H′(θ)=−ln2⋅ln(32)<0, so we get H(θ)>H(1)>0, therefore g′(θ)<g′(1)<0, g(θ) monotone decreasing, g(1)<g(θ)<g(0)=−3.6227<0, that is f′(θ) monotone decreasing, 0<f′(1)<f′(θ)<f′(0), where f′(1)=−3+5ln3−52ln2>0, that is f(θ) monotone increasing, f(0)<f(θ)<f(1), where f(0)=0, to sum up, we can get f(θ)>0.
(3) Let f1(θ)=−34+5θ−42(4−θ)⋅31−θ+(4+θ)(6−θ)(4−θ)2⋅31−θ⋅2−θ−(6−θ)2(4−θ)2(2−θ)2≐-34+a1⋅31−θ+a2⋅31−θ2−θ+a3⋅(2−θ)2, where
Next, we just need to prove that f1(θ)>0, using a Taylor expansion yields
where
Because of 1−θ1!⋅12+(1−θ)(−θ)2!⋅(12)2+(1−θ)(−θ)(−θ−1)3!(12)3+…≐∑+∞k=0ak is an alternating series with positive first term, and ∑+∞k=0ak=a0+a1+∑+∞k=2ak, where ∑+∞k=2ak is an alternating series with positive first term, so 0<∑+∞k=2ak<a2, we have
We wish prove f1(θ)>0, only need −34+2−θ8(4−θ)2[f2(θ)+f3(θ)]>0⟺f2(θ)+f3(θ)>34⋅8(4−θ)22θ=6(4−θ)22θ, that is f2(θ)+f3(θ)>6(4−θ)22θ≐6f4(θ), so to prove f1(θ)>0, just prove f2(θ)+f3(θ)−6f4(θ)>0, let's remember ˉf(θ)=f2(θ)+f3(θ)−6f4(θ), since ˉf(θ) first increases and then decreases, and the two endpoints are ˉf(0)=0,ˉf(1)=16>0, ˉf(θ)>0 is always true, so f1(θ)>0.
(4) Firstly, we will prove ρ>0. Because ρ=3(4−θ)+(θ−6)21−θ2(4−θ), let ˉg(θ)=3(4−θ)+(θ−6)21−θ, by calculation, we can get ˉg(θ) is monotonically increasing function on θ, that is ˉg(θ)>ˉg(0)=0. therefore, we have ρ>0.
In addition, ρ−23=5(4−θ)+3(θ−6)21−θ6(4−θ)≐ˉg1(θ)6(4−θ). we can directly find ˉg1(θ)<ˉg1(1)=0, so 0<ρ<23.
(5) We use Lagrange's mean value theorem 41−θ>43+θ⋅31−θ, and we have
Let a1=−3(4−θ)2,a2=6(4−θ)(6−θ),a3=43+θ[−(4−θ)(8−θ)(3+θ)+(6−θ)2],a4=1+1−θ1!12+(1−θ)(−θ)2(12)2. (A.1) is becomes
where bk=Πki=0(−θ−1−i)⋅−1(k+3)!(12)k+3, and b_{0} = \frac{-(-\theta-1)}{3!}(\frac{1}{2})^{3} > 0, |\frac{b_{k+1}}{b_{k}}| = \frac{\theta+2+k}{2(k+4)} < 1, so \sum_{k = 0}^{+\infty}b_{k} is a convergent alternating series, that is 0 < \sum_{k = 0}^{+\infty}b_{k} < b_{0}. Because of a_3 < 0 , so we get
where a_{5} = a_{4}+(1-\theta)\theta\cdot \frac{-(-\theta-1)}{3!}(\frac{1}{2})^{3} = \frac{48+24(1-\theta)-6(\theta-\theta^{2})+(\theta-\theta^{3})}{48}, by careful calculation, we can get
because \theta\in(0, 1), so \theta^{3} < \theta^{2}, \theta^{6} > 0.
where \tilde f_{4}(\theta) = -16\theta^{5}+97\theta^{4}-278\theta^{3}+88\theta^{2}+732\theta+864, by careful calculation, we can get \tilde f^{\prime\prime}_{3}(\theta) < 0 , that is \tilde f^{\prime}_{3}(\theta) monotone decreasing, so \tilde f^{\prime}_{3}(1) < \tilde f^{\prime}_{3}(\theta) < \tilde f^{\prime}_{3}(0) , by direct calculation, we can get \tilde f^{\prime}_{3}(\theta) = 144(2\theta+2)+2^{1-\theta}[\tilde f_{4}(\theta)(-\ln2)+\tilde f^{\prime}_{4}(\theta)], \tilde f^{\prime}_{3}(0) > 0, \tilde f^{\prime}_{3}(0) < 0, so \tilde f^{\prime}_{3}(0) changes from positive to negative, \tilde f_{3}(0) increases first and then decreases, \tilde f_{3}(0) = 0, \tilde f_{3}(1) = 223 > 0, that is \tilde f_{3}(\theta) > 0 established, so \tilde f(\theta) > 0 .