Loading [MathJax]/jax/element/mml/optable/GeneralPunctuation.js
Research article

A high order numerical method for solving Caputo nonlinear fractional ordinary differential equations

  • Received: 01 August 2021 Accepted: 09 September 2021 Published: 15 September 2021
  • MSC : 26A33, 65L12, 65M06, 65M12

  • In this paper, we construct a high order numerical scheme for Caputo nonlinear fractional ordinary differential equations. Firstly, we use the piecewise Quadratic Lagrange interpolation method to construct a high order numerical scheme for Caputo nonlinear fractional ordinary differential equations, and then analyze the local truncation error of the high order numerical scheme. Secondly, based on the local truncation error, the convergence order of 3θ order is obtained. And the convergence are strictly analyzed. Finally, the numerical simulation of the high order numerical scheme is carried out. Through the calculation of typical problems, the effectiveness of the numerical algorithm and the correctness of theoretical analysis are verified.

    Citation: Xumei Zhang, Junying Cao. A high order numerical method for solving Caputo nonlinear fractional ordinary differential equations[J]. AIMS Mathematics, 2021, 6(12): 13187-13209. doi: 10.3934/math.2021762

    Related Papers:

    [1] Alessandra Jannelli, Maria Paola Speciale . On the numerical solutions of coupled nonlinear time-fractional reaction-diffusion equations. AIMS Mathematics, 2021, 6(8): 9109-9125. doi: 10.3934/math.2021529
    [2] Sara S. Alzaid, Pawan Kumar Shaw, Sunil Kumar . A study of Ralston's cubic convergence with the application of population growth model. AIMS Mathematics, 2022, 7(6): 11320-11344. doi: 10.3934/math.2022632
    [3] Xiaopeng Yi, Chongyang Liu, Huey Tyng Cheong, Kok Lay Teo, Song Wang . A third-order numerical method for solving fractional ordinary differential equations. AIMS Mathematics, 2024, 9(8): 21125-21143. doi: 10.3934/math.20241026
    [4] Najat Almutairi, Sayed Saber . Chaos control and numerical solution of time-varying fractional Newton-Leipnik system using fractional Atangana-Baleanu derivatives. AIMS Mathematics, 2023, 8(11): 25863-25887. doi: 10.3934/math.20231319
    [5] Sara S. Alzaid, Pawan Kumar Shaw, Sunil Kumar . A numerical study of fractional population growth and nuclear decay model. AIMS Mathematics, 2022, 7(6): 11417-11442. doi: 10.3934/math.2022637
    [6] Ali Turab, Hozan Hilmi, Juan L.G. Guirao, Shabaz Jalil, Nejmeddine Chorfi, Pshtiwan Othman Mohammed . The Rishi Transform method for solving multi-high order fractional differential equations with constant coefficients. AIMS Mathematics, 2024, 9(2): 3798-3809. doi: 10.3934/math.2024187
    [7] Din Prathumwan, Thipsuda Khonwai, Narisara Phoochalong, Inthira Chaiya, Kamonchat Trachoo . An improved approximate method for solving two-dimensional time-fractional-order Black-Scholes model: a finite difference approach. AIMS Mathematics, 2024, 9(7): 17205-17233. doi: 10.3934/math.2024836
    [8] Naimi Abdellouahab, Keltum Bouhali, Loay Alkhalifa, Khaled Zennir . Existence and stability analysis of a problem of the Caputo fractional derivative with mixed conditions. AIMS Mathematics, 2025, 10(3): 6805-6826. doi: 10.3934/math.2025312
    [9] Abdon Atangana, Seda İğret Araz . A successive midpoint method for nonlinear differential equations with classical and Caputo-Fabrizio derivatives. AIMS Mathematics, 2023, 8(11): 27309-27327. doi: 10.3934/math.20231397
    [10] Khalid K. Ali, Mohamed A. Abd El Salam, Mohamed S. Mohamed . Chebyshev fifth-kind series approximation for generalized space fractional partial differential equations. AIMS Mathematics, 2022, 7(5): 7759-7780. doi: 10.3934/math.2022436
  • In this paper, we construct a high order numerical scheme for Caputo nonlinear fractional ordinary differential equations. Firstly, we use the piecewise Quadratic Lagrange interpolation method to construct a high order numerical scheme for Caputo nonlinear fractional ordinary differential equations, and then analyze the local truncation error of the high order numerical scheme. Secondly, based on the local truncation error, the convergence order of 3θ order is obtained. And the convergence are strictly analyzed. Finally, the numerical simulation of the high order numerical scheme is carried out. Through the calculation of typical problems, the effectiveness of the numerical algorithm and the correctness of theoretical analysis are verified.



    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 k2 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 k2 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.

    We consider the following initial value problem of nonlinear fractional ordinary differential equations

    {0Dθtu(t)=f(t,u(t)),0tT,0<θ<1,(2.1)u(0)=u0,(2.2)

    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

    0Dθtu(t)=t0ω1θ  (tτ)u(τ)dτ, (2.3)

    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

    J[t0,t2]u(t)=ϖ0,0(t)u0+ϖ1,0(t)u1+ϖ2,0(t)u2, (2.4)

    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

    ϖ0,0(t)=(tt1)(tt2)2Δt2, ϖ1,0(t)=(tt0)(tt2)Δt2, ϖ2,0(t)=(tt0)(tt1)2Δt2.

    Let Δuk=ukuk1,Δu0=u0,k1. 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

    0Dθtu(t1)=t10ω1θ  (t1τ)u(τ)dτt10ω1θ  (t1τ)[J[t0,t2]u(τ)]dτ=1Δt2t10ω1θ  (t1τ)[(τt12)Δu2(τt32)Δu1]dτ=B11Δu1+B21Δu2 0DθΔtu1, (2.5)
    0Dθtu(t2)=t20ω1θ  (t2τ)u(τ)dτt20ω1θ  (t2τ)[J[t0,t2]u(τ)]dτ=1Δt2t20ω1θ  (t2τ)[(τt12)Δu2(τt32)Δu1]dτ=B12Δu1+B22Δu2 0DθΔtu2, (2.6)

    where

    B1j=1Δt2tj0ω1θ  (tjτ)(τt32)dτ,  B2j=1Δt2tj0ω1θ  (tjτ)(τt12)dτ,j=1,2. (2.7)

    When k3, we have

    0Dθtu(tk)=t10ω1θ  (tkτ)u(τ)dτ+k1j=1tj+1tjω1θ  (tkτ)u(τ)dτt10ω1θ  (tkτ)[J[t0,t2]u(τ)]dτ+k1j=1tj+1tjω1θ  (tkτ)[J[tj1,tj+1]u(τ)]dτ=t10ω1θ  (tkτ)[2τt1t22Δt2u0+2τt0t2Δt2u1+2τt0t12Δt2u2]dτ   +k1j=1tj+1tjω1θ  (tkτ)[2τtjtj+12Δt2uj1+2τtj1tj+1Δt2uj+2τtj1tj2Δt2uj+1]dτ=1Δt2t10ω1θ  (tkτ)[(τt12)Δu2(τt32)Δu1]dτ   +1Δt2k1j=1tj+1tjω1θ  (tkτ)[(τtj12)Δuj+1(τtj+12)Δuj]dτ=1Δt2[t20ω1θ  (tkτ)(τt32)dτ]Δu1+1Δt2[t20ω1θ  (tkτ)(τt12)dτ]Δu2   +1Δt2k1j=2tj+1tjω1θ  (tkτ)[(τtj12)Δuj+1(τtj+12)Δuj]dτ 0DθΔtuk,

    we can conclude that

    0DθΔtuk=Akk1Δu1+Akk2Δu2+k1j=3AkkjΔuj+Ak0Δuk,3kK, (2.8)

    where

    Ak0=1Δt2tktk1ω1θ  (tkτ)(τtk32)dτ,  Akk1=1Δt2t20ω1θ  (tkτ)(τt32)dτ,Akk2=1Δt2[t20ω1θ  (tkτ)(τt12)dτt3t2ω1θ  (tkτ)(τt52)dτ],Akkj=1Δt2[tjtj1ω1θ  (tkτ)(τtj32)dτtj+1tjω1θ  (tkτ)(τtj+12)dτ],j=3,,k1. (2.9)

    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

    0DθΔtuk=fk,k=1,2,,K. (2.10)

    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

    1Γ(1θ)t0u(s)(ts)θds=u(0)t1θ  Γ(1θ)(1θ)+t0f(s)ds. (2.11)

    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 Akkj and judge its size relation.

    Lemma 2.1. Ak0>0, the sign of Ak1 is uncertain, Ak2>Ak3>>Akk1>0.

    Proof. Through calculation, we can get the following results

    Ak0=1Δt2Γ(2θ)tktk1(τtk32)d(tkτ)1θ  =1ΔtθΓ(3θ)(2θ2)>0, (2.12)
    Ak1=1Δt2[tk1tk2ω1θ  (tkτ)(τtk52)dτtktk1ω1θ  (tkτ)(τtk12)dτ]=1ΔtθΓ(3θ)[(6θ)2θ4+θ]1ΔtθΓ(3θ)M(θ). (2.13)

    The sign of Ak1 is determined by the sign of M(θ), which satisfies

    M(θ)=2θ(log2)[2+(6θ)log2]>0,  M(1)=12(15log2)<0.

    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,,k2, we have

    Akkj=1ΔtθΓ(3θ){12(2θ)[(kj1)1θ  2(kj)1θ  +(kj+1)1θ  ]+(kj1)2θ2(kj)2θ+(kj+1)2θ}, (2.14)

    let ˉj=kj, ˉj2, then using Taylor expansion, we can get the following results

    Akkj=1ΔtθΓ(2θ)ˉjθ+i=0aiˉji, (2.15)

    where ai=Πin=0(1θn)1(i+2)![(1)i(i2)+i+42]. Next, we will prove ˉjθ+i=0aiˉji is decrease function on ˉj. Because +i=0aiˉji is a convergence series, and the radius of convergence exists, it can exchange the derivative and the sum.

    (ˉjθ+i=0aiˉji)=ˉjθ[θˉj1+i=0aiˉji++i=1iaiˉji1]=ˉjθ[θˉj1a0+(1+θ)ˉj2a1++i=2(i+θ)aiˉji1].

    We can find +i=2(i+θ)aiˉji1 is an alternating series with positive first term. So +i=2(i+θ)aiˉji1>0, we have

    θˉj1a0+(1+θ)ˉj2a1++i=2(i+θ)aiˉji1>θˉj1a0+(1+θ)ˉj2a1=θˉj1(1θ)+(1+θ)ˉj2(12)(1θ)θ=(1θ)θ[ˉj112(1+θ)ˉj2]>0.

    Therefore, we have already proved ˉjθ+i=0aiˉji is decrease function on ˉj. According to (2.15), we obtain

    Ak2>Ak3>>Akk3. (2.16)

    Next, let's prove Akk3>Akk2>Akk1>0. Through careful calculation, we can get

    Akk2Akk1=1Γ(3θ)ΔtθF1,Akk3Akk2=1Γ(3θ)ΔtθF2,

    where

    F1=32(2θ)(k2)1θ  +12(2θ)(k3)1θ  2(2θ)k1θ  +(k3)2θ3(k2)2θ+2k2θ,F2=12(2θ)k1θ  +32(2θ)(k2)1θ  32(2θ)(k3)1θ  +12(2θ)(k4)1θ  k2θ+3(k2)2θ3(k3)2θ+(k4)2θ,

    According to (1) in Lemma A.1, we get F1>0,F2>0. Therefore, we have

    Akk2Akk1>0,   Akk3Akk2>0. (2.17)

    The following step, we will prove Akk1>0. By carefully calculation, we have

    Akk1=1ΔtθΓ(3θ){2θ2[(k2)1θ  +3k1θ  ]+(k2)2θk2θ}, (2.18)

    According to (2.33), we have

    Akk1>0. (2.19)

    Combining (2.12), (2.13), (2.16), (2.17) with (2.19), we already proved Lemma 2.1.

    We let Δˉuj=Δ(ujρuj1), Δˉu0=Δu0. Let ρ=12(1Ak1Ak0)=12(3+θ64θ21θ  ). Therefore, from (2.8), we can get that

    0DθΔtuk=Ak0Δuk+Ak1Δuk1+Ak2Δuk2++Akk2Δu2+Akk1Δu1=Ak0Δ(ukρuk1)+(Ak1+ρAk0)Δ(uk1ρuk2)++(Akkj+ρAkkj1++ρkjAk0)Δ(ujρuj1)++(Akk1+ρAkk2++ρk1Ak0)Δ(u1ρu0)+(ρAkk1+ρ2Akk2++ρkAk0)Δu0=kj=1ˉAkkjΔˉuj+Δu0ni=1ρiAkki, (2.20)

    where

    ˉAkkj=ki=jρijAkki. (2.21)

    Next, we discuss the relationship between the size of ˉAkkj.

    Lemma 2.2. ˉAk0>ˉAk1>ˉAk2>>ˉAkk1>0.

    Proof. When k=3, we only need to prove ˉA30>ˉA31>ˉA32>0. By careful calculation, we can get

    ˉA31ˉA30=12ΔtθΓ(3θ)[32(4θ)(4+θ)31θ  +(6θ)2θ].

    According to (2) in Lemma A.1, we get

    ˉA30>ˉA31. (2.22)

    By direct calculation, we have

    ˉA32ˉA31=A32A31+ρ(ˉA31ˉA30)=4θ2ΔtθΓ(3θ)[34+5θ42(4θ)31θ  +(4+θ)(6θ)(4θ)231θ  2θ(6θ)2(4θ)2(2θ)2].

    According to (3) in Lemma A.1, we have

    ˉA31>ˉA32. (2.23)

    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,

    ˉA32>0. (2.24)

    Combining (2.22), (2.23) and (2.24), when k=3, we know that Lemma 2.2 holds.

    When k4, through careful calculation, we can draw the following conclusions

    ˉAk1ˉAk0=2θ2ΔtθΓ(3θ)[32+(θ6)2θ4θ]=2θ2ΔtθΓ(3θ)ρ. (2.25)

    According to (4) in Lemma A.1, we can get ρ>0, we have ˉAk0>ˉAk1. Because

    ˉAk2ˉAk1=Ak2Ak1+ρ(ˉAk1ˉAk0)=18(4θ)ΔtθΓ(3θ)˜f(θ),

    where ˜f(θ)=3(4θ)2+6(4θ)(6θ)21θ  4(4θ)(8θ)31θ  +(6θ)241θ  . According to (5) in Lemma A.1, we can get ˜f(θ)>0. It can be concluded that

    ˉAk1>ˉAk2. (2.26)

    By careful calculation, we can get ˉAk3ˉAk2=Ak3Ak2+ρ(ˉAk2ˉAk1). According to ρ>0, (2.26) and Lemma 2.1, we can get

    ˉAk2>ˉAk3. (2.27)

    By mathematical induction, and Lemma 2.1 we can conclude that

    ˉAk1jˉAk2j=Akk1jAkk2j+ρ(ˉAkk2jˉAkk3j)>0,j=0,,k4.

    So we have

    ˉAk0>ˉAk1>ˉAk2>>ˉAkk1. (2.28)

    Next we prove that ˉAkk1>0. According to (2.21), we have

    ˉAkk1=Akk1+k2i=2ρijAkki+ρk2ˉAk1.

    According to (2.37), we get ˉAk1>0. Using Lemma 2.1, ρ>0, we have

    ˉAkk1>0. (2.29)

    Combining (2.28) with (2.29), we already proved Lemma 2.2.

    Lemma 2.3. There is a constant πA6 such that the discrete kernel satisfies the lower bound.

    ˉAkj1πAΔttkjtkj1ω1θ  (tkτ)dτ=1πAΔtθΓ(2θ)[(j+1)1θ  j1θ  ]. (2.30)

    Proof. When k4, According to (2.21), Lemma 2.1 and Lemma 2.2, we can get

    ˉAkkj=Akkj+k2i=j+1Akkiρij+ρk1jˉAk1Akkj.

    For j=1, we have

    ˉAkk1Akk1=1ΔtθΓ(3θ){2θ2[(k2)1θ  +3k1θ  ]+(k2)2θk2θ}, (2.31)

    let ˜j=k1,˜j3, using Taylor formula to expand the calculation, we can get the following results

    Akk1=˜jθΔtθΓ(2θ)+i=0ai, (2.32)

    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=2aia2, so

    Akk1=˜jθΔtθΓ(2θ)(a0+a1++i=2ai)˜jθΔtθΓ(1θ)(1836)=˜jθΔtθΓ(1θ)97. (2.33)

    We have

    πA97. (2.34)

    For j=2, we have

    ˉAkk2Akk2=1ΔtθΓ(3θ){2θ2[(k3)1θ  2(k2)1θ  k1θ  ]+(k3)2θ2(k2)2θ+k2θ},

    let ˆj=k2,ˆj2, and then using Taylor expansion, similar to j=1, it can be obtained by calculation

    Akk2=ˆjθΔtθΓ(3θ)+i=0ai,

    where ai=Πin=0(1θn)1(i+2)!(1ˆj)i[i2(1)i+2i22i+1], because +i=3ai is a convergent staggered series, and a3>0, 0<+i=3ai<a3, we get

    ˉAkk2ˆjθΔtθΓ(2θ)[a0+a1+a2]ˆjθΔtθΓ(2θ)[151221242]=ˆjθΔtθΓ(2θ)3748,

    therefore, we have

    πA4837. (2.35)

    For j=3,4,,k2, we have

    ˉAkkjAkkj=1ΔtθΓ(3θ){12(2θ)[(kj1)1θ  2(kj)1θ  +(kj+1)1θ  ]+(kj1)2θ2(kj)2θ+(kj+1)2θ},

    let ˉj=kj, ˉj2, we have

    Akkj=ˉj1θ  ΔtθΓ(3θ){12(2θ)[(11ˉj)1θ  2+(1+1ˉj)1θ  ]+ˉj[(11ˉj)2θ2+(1+1ˉj)2θ]},

    then using Taylor expansion, we can get the following results

    Akkj=ˉjθΔtθΓ(2θ)+i=0ai,

    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=2aia2, so

    ˉAkkjˉjθΔtθΓ(1θ)11θ  [1θ+33!(1θ)(θ)1ˉj]ˉjθΔtθΓ(1θ)(114),

    so we can get

    πA43. (2.36)

    When j=k1,

    ˉAk1=ρˉAk0+Ak1=1ΔtθΓ(3θ)[θ41+(3θ2)2θ]12ΔtθΓ(1θ)[θ41+(3θ2)2θ]14ΔtθΓ(1θ). (2.37)

    Therefore, we have

    πA4. (2.38)

    When j=k, according to (2.30), we can get

    ˉAk0=Ak0=1ΔtθΓ(3θ)(2θ2)12ΔtθΓ(2θ)(2θ2)1ΔtθΓ(2θ). (2.39)

    We have

    πA1. (2.40)

    When k=3, ˉA30 is in (2.39), so 1πA1. According to (2.30), we only need

    ˉA31=ρˉA30+A31=1ΔtθΓ(3θ)[θ41+(θ23)2θ+(θ2+2)31θ  ]=21θ  1ΔtθΓ(2θ)θ41+(θ23)2θ+(θ2+2)31θ  (21θ  1)(2θ)21θ  12ΔtθΓ(2θ)1(21θ  1)(2θ)21θ  14ΔtθΓ(2θ). (2.41)

    Therefore, we have

    πA4. (2.42)

    Next, we calculate ˉA32,

    ˉA32=A32+ρˉA31=1ΔtθΓ(3θ)18(4θ)[(θ4)2+4(4θ)(θ6)2θ+4(6θ)2(2θ)2   +[6(4θ)24(6θ)(4+θ)2θ]31θ  ]31θ  21θ  ΔtθΓ(2θ)18(4θ)(31θ  21θ  )(2θ)(2+θ)   [θ36θ2+32+4(θ3+17θ276θ+96)2θ+4(θ34θ272)22θ]31θ  21θ  ΔtθΓ(2θ)116(4θ)(2+θ)[θ36θ2+32+4(θ3+17θ276θ+96)2θ   +4(θ34θ272)22θ]31θ  21θ  ΔtθΓ(2θ)2416(4θ)(2+θ)31θ  21θ  6ΔtθΓ(2θ).

    Therefore, we have

    πA6. (2.43)

    Combining (2.33)–(2.38) with (2.12)–(2.43), we have already proved Lemma 2.3.

    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

    rk(Δt)=0Dθtu(tk)0DθΔtuk,k1. (3.1)

    Then there exists a constant Cu depending only on the function u, such that for all Δt>0,

    |rk(Δt)|CuΔt3θ. (3.2)

    Proof. Our error estimation will be established on the following Taylor theorem

    u(t)J[t0,t2]u(t)=u(3)(ξ(τ))6(τt0)(τt1)(τt2),τ[t0,t2], (3.3)

    where ξ(τ) is a function defined on [t0,t2] with range (t0,t2). We first estimate r1(Δt)

    |r1(Δt)|=|1Γ(1θ){t10u(τ)(t1τ)θdτt10[J[t0,t2]u(τ)](t1τ)θdτ}|=θΓ(1θ)|t10[u(τ)J[t0,t2]u(τ)](t1τ)θ1dτ|=θΓ(1θ)|t10u(3)(ξ(τ))6(τt0)(τt2)(t1τ)θdτ|θΓ(1θ)maxt[0,T]|u(3)(t)|6t10τ(t2τ)(t1τ)θdτ=θ3(3θ)Γ(2θ)maxt[0,T]|u(3)(t)|Δt3θCuΔt3θ.

    This proves (3.2) for k=1. The case k=2 can be similarly proven, and here we omit the details.

    When k3, we have

    |rk(Δt)|=|0Dθtu(tk)0DθΔtu(tk)|=1Γ(1θ)|t10(tkτ)θ[u(τ)J[t0,t2]u(τ)]dτ+k1j=1tj+1tj(tkτ)θ[u(τ)J[tj1,tj+1]u(τ)]dτ||M+N|.

    Similar to the proof of |r1(Δt)|, we have

    |M|CuΔt3θ. (3.4)

    For N, we have

    |N|=1Γ(1θ)|k1j=1tj+1tj(tkτ)θd[u(τ)J[tj1,tj+1]u(τ)]|=1Γ(1θ)|k2j=1tj+1tj(tkτ)θd[u(τ)J[tj1,tj+1]u(τ)]+tktk1(tkτ)θd[u(τ)J[tj1,tj+1]u(τ)]|=1Γ(1θ){|k2j=1tj+1tju(3)(ξ)6(τtj1)(τtj)(τtj+1)d(tkτ)θtktk1u(3)(ξ)6(τtk2)(τtk1)(τtk)d(tkτ)θ|}3θ27Γ(1θ)maxt[0,T]|u(3)(t)|Δt3tj+1tj(tkτ)θ1dτ+θ6Γ(1θ)maxt[0,T]|u(3)(t)||Δt2tktk1(tkτ)θdτ3θ27Γ(1θ)maxt[0,T]|u(3)(t)|Δt3θ+θ6Γ(1θ)maxt[0,T]|u(3)(t)|Δt3θ.

    In the above inequality, we use |(τtj1)(τtj)(τtj+1)|239Δt3.

    We can get the following conclusions

    |rk(Δt)|=|M+N|CuΔt3θ. (3.5)

    Complete the proof of Theorem 3.1.

    In order to obtain the convergence analysis, we now introduce an important tool: complementary discrete convolution kernel. Because of ωθωβ=ωθ+β, therefore

    tsωθ(tμ)ω1θ  (μs)dμ=ω1(ts)=1,for0<s<t<. (4.1)

    A class of complementary discrete convolution kernels Pkki with the same properties are given

    ki=mPkkiˉAiim1,for3mkK. (4.2)

    According to [3], we have

    Pk0=1ˉAk0,Pki=1ˉAki0i1j=0(ˉAkjij1ˉAkjij)Pkj,for1ik3, (4.3)

    where Pki>0,i=0,1,...,k3 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

    0PkkiπAΓ(2θ)Δtθ,1ikN, (4.4)
    kj=1Pkkiω1θ  (ti)πA,1kK. (4.5)

    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

    f(tk,u(tk))f(tk,uk)=Lk(u(tk)uk)=Lkek,k=1,2,...,K, (4.6)

    and |Lk|L,k1.

    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

    Δt1θ2πAΓ(2θ)Λ, (4.7)

    where πA6, and Λ=12L, then the following error estimates hold,

    |u(tk)uk|CΔt3θ, (4.8)

    here C is only dependent on the function u.

    Proof. Because Δu1=u1u0,Δu2=u2u1, through calculation, we can get

    {0Dθtu(t1)=B11u0+(B11B21)u1+B21u2,0Dθtu(t2)=B12u0+(B12B22)u1+B22u2, (4.9)

    let B11=ΔtθˉD10, B11B21=ΔtθˉD11, B21=ΔtθˉD12; B12=ΔtθD10, B12B22=ΔtθD11, B22=ΔtθD12. It is easy to obtain by calculation

    {ΔtθˉD11e1+ΔtθˉD12e2=f(t1,u(t1))f(t1,u1)r1(Δt),ΔtθD11e1+ΔtθD12e2=f(t2,u(t2))f(t2,u2)r2(Δt). (4.10)

    Among (4.6), then we can get

    {ΔtθˉD11e1+ΔtθˉD12e2=L1e1r1(Δt),ΔtθD11e1+ΔtθD12e2=L2e2r2(Δt).

    That is

    (4.11)

    where D = \begin{pmatrix} \bar D^{1}_{1} & \bar D^{1}_{2}\\ D^{1}_{1} & D^{1}_{2}\\ \end{pmatrix}. Since \bar D^{1}_{1}D^{1}_{2}-\bar D^{1}_{2}D^{1}_{1} is bounded, then

    \begin{eqnarray} \max \{ \vert e_{1} \vert , \vert e_{2} \vert \} \leq C_{1}\Delta t^{\theta} \max \{ \vert r_{1} (\Delta t) \vert , \vert r_{2}(\Delta t) \vert \}, \end{eqnarray} (4.12)

    therefore,

    \begin{eqnarray} \begin{array}{lll}{ } \begin{cases} \vert e_{1} \vert \leq C\Delta t^{3-\theta}, \\ \vert e_{2} \vert \leq C\Delta t^{3-\theta}, \\ \end{cases} \end{array} \end{eqnarray} (4.13)

    when j = 1, 2, (4.8) is proved.

    When i\geq 3 , let's set \bar e_{i} = e_{i}-\rho e_{i-1}, i = 1, ..., j , and \bar e_{0} = e_{0} for (4.13), we have

    \begin{eqnarray} \begin{array}{llll}{ } \begin{cases} \vert \bar e_{1}\vert = \vert e_{1}-\rho e_{0}\vert \leq C\Delta t^{3-\theta}, \\ \vert \bar e_{2}\vert = \vert e_{2}-\rho e_{1}\vert \leq \frac{5}{3}C\Delta t^{3-\theta}.\\ \end{cases} \end{array} \end{eqnarray} (4.14)

    Combining with (2.20) and (4.6), \bar e_{j} , j\geq 3 , satisfy

    \begin{eqnarray} \begin{array}{llll}{ } \sum\limits_{j = 1}^{k}\bar A_{k-j}^{k}\Delta \bar e_{j} = f(t_{k}, u(t_{k}))-f(t_{k}, u_{k})-r_{k}(\Delta t). \end{array} \end{eqnarray} (4.15)

    Because of e_{j} = \sum_{n = 1}^{j}\rho^{j-n}\bar e_{n} , then

    \begin{eqnarray} \begin{array}{lll}{ } \sum\limits_{j = 3}^{k} \bar A_{k-j}^{k}\Delta \bar e_{j} = L_{k}\sum\limits_{j = 3}^{k}\rho^{k-j}\bar e_{j}-\bar r_{k}(\Delta t), \end{array} \end{eqnarray} (4.16)

    where

    \begin{eqnarray} \begin{array}{lll}{ } \bar r_{k}(\Delta t) = r_{k}(\Delta t)-L_{k}\sum\limits_{j = 1}^{2}\rho^{k-j}\bar e_{j}+\sum\limits_{j = 1}^{2}\bar A_{k-j}^{k}\Delta \bar e_{j}, \end{array} \end{eqnarray} (4.17)

    by multiplying 2\bar e_{k} on both sides of (4.16), from Lemma A.1 in Liao [3], we can obtain that

    \begin{eqnarray*} &\!\!\!&\!\!\!2\bar e_{k}\sum\limits_{j = 3}^{k}\bar A_{k-j}^{k}\Delta \bar e_{j}\geq \sum\limits_{j = 3}^{k}\bar A_{k-j}^{k}\Delta(\vert \bar e_{j}\vert^{2}), \\ &\!\!\!&\!\!\! 2\bar e_{k}[L_{k}\sum\limits_{j = 3}^{k}\rho^{k-j}\bar e_{j}-\bar r_{k}(\Delta t)] \leq L \sum\limits_{j = 3}^{k}\rho^{k-j}(\vert \bar e_{j}\vert^{2}+\vert \bar e_{k}\vert^{2})+2\vert \bar e_{k}\vert \vert \bar r_{k}(\Delta t)\vert\\ &\!\!\!&\!\!\! \leq 4L \sum\limits_{j = 3}^{k}\rho^{k-j}\vert \bar e_{j}\vert^{2}+2\vert \bar e_{k}\vert \vert \bar r_{k}(\Delta t)\vert. \end{eqnarray*}

    To sum up, we can get the following conclusions

    \begin{eqnarray} \begin{array}{lll}{ } \sum\limits_{j = 3}^{k}\bar A_{k-j}^{k}\Delta(\vert \bar e_{j}\vert^{2})\leq 4L\sum\limits_{j = 3}^{k}\rho^{k-j}\vert \bar e_{j}\vert^{2}+2\vert \bar e_{k}\vert \cdot\vert \bar r_{k}(\Delta t)\vert. \end{array} \end{eqnarray} (4.18)

    We change k to i , and then multiply both sides of (4.18) by a \sum_{i = 3}^{k}P_{k-i}^{k} , then

    \begin{eqnarray} \begin{array}{llll}{ } \sum\limits_{i = 3}^{k}P_{k-i}^{k}\sum\limits_{j = 3}^{i}\bar A_{i-j}^{i}\Delta (\vert \bar e_{j}\vert^{2})\leq \sum\limits_{i = 3}^{k}P_{k-i}^{k}\sum\limits_{j = 3}^{i}4L\rho^{i-j}\vert \bar e_{j}\vert^{2}+2\sum\limits_{i = 3}^{k}P_{k-i}^{k}\vert \bar e_{i}\vert \vert \bar r_{i}(\Delta t)\vert . \end{array} \end{eqnarray} (4.19)

    From (4.2) and \Delta u_{j} = u_{j}-u_{j-1} , we can get the following results

    \begin{eqnarray} \begin{array}{lll}{ } \sum\limits_{i = 3}^{k}P_{k-i}^{k}\sum\limits_{j = 3}^{i}\bar A_{i-j}^{i}\Delta (\vert \bar e_{j}\vert^{2}) = \sum\limits_{j = 3}^{k}\Delta(\vert\bar e_{j}\vert^{2})\sum\limits_{i = j}^{k}P_{k-i}^{k}\bar{A}_{i-j}^{i} = \sum\limits_{j = 3}^{k}\Delta(\vert\bar e_{j}\vert^{2}) = \vert \bar e_{k}\vert^{2}-\vert \bar e_{2}\vert^{2}, \end{array} \end{eqnarray} (4.20)

    substituting (4.20) into (4.19) yields

    \begin{eqnarray} \begin{array}{lll}{ } \vert \bar e_{k}\vert^{2}\leq 4L\sum\limits_{i = 3}^{k}P_{k-i}^{k}\sum\limits_{j = 3}^{i}\rho^{i-j}\vert \bar e_{j}\vert^{2}+\vert e_{2}\vert^{2}+2\sum\limits_{i = 3}^{k}P_{k-i}^{k}\vert \bar e_{i}\vert \vert \bar r_{i}(\Delta t)\vert. \end{array} \end{eqnarray} (4.21)

    Next, we will prove it by mathematical induction

    \begin{eqnarray} \begin{array}{llll}{ } \vert \bar e_{k}\vert \leq 2E_{\theta}(2\pi_{A}\Lambda t_{k}^{\theta})(\vert \bar e_{2}\vert+2 \max \limits_{3 \leq j \leq k} \sum\limits_{i = 3}^{k} p_{j-i}^{j}\left|\bar{r}_{i}(\Delta t)\right|)\doteq F_{k}G_{k}, \quad for \quad k\geq 2, \end{array} \end{eqnarray} (4.22)

    where \Lambda = 12L , E_{\theta} stands for the Mittag-Leffler function, we usually define it as: E_{\theta}(z) = \sum_{j = 0}^{\infty}\frac{Z^{j}}{\Gamma(1+j\theta)} , and E_{\theta}(0) = 1, \quad \forall Z\in R , and Z > 0, E_{\theta}^{'}(Z) > 0 , so F_{k}\geq F_{k-1}\geq2 for all k\geq2 . where G_{k} = |\bar{e}_{2}|+2 {\max_{3 \leq j\leq k}} \sum_{i = 3}^{k} P_{j-i}^{j}\left|\bar{r}_{i}(\Delta t)\right| .

    When k = 3 , (i) if |\bar e_{3}| \leq | e_{2}|, G_{3} = |e_{2}|+2P_{0}^{3}|\bar r_{3}(\Delta t)| \geq |\bar e_{2}| , we have |\bar e_{3}|\leq |\bar e_{2}| \leq G_{3}\leq F_{3}G_{3},

    (ii) if |\bar e_{3}| > |e_{2}| , from (4.21), it can be concluded that

    \begin{eqnarray*} |\bar e_{3}|^{2}\leq4 L P_{0}^{3} |\bar e_{3}|^{2}+ |\bar e_{2}|^{2}+2P_{0}^{3}|\bar e_{3}| |\bar r_{3}(\Delta t)|. \end{eqnarray*}

    According to (4.4), it can be concluded that

    \begin{eqnarray*} \begin{array}{lll}{ } 4LP_{0}^{3}\leq \pi_{A}\Gamma(2-\theta)\Delta t^{\theta}\Lambda \leq \frac{1}{2} \quad \Longrightarrow \Delta t\leq\frac{1}{\sqrt[\theta]{2\pi_{A}\Gamma(2-\theta)\Lambda}}, \end{array} \end{eqnarray*}

    so |\bar e_{3}|\leq 2(|\bar{e}_{2}|+2P_{0}^{3}|\bar{r}_{3}(\Delta t)|)\leq 2G_{3}\leq F_{3}G_{3}.

    When 4\leq k\leq K , we assume that: |\bar e_{j}|\leq F_{j}G_{j} , for 3 \leq j\leq k-1 , let | \bar e_{j(k)}| = \max_{2\leq i\leq k-1}|\bar e_{i}|.

    (i) When |\bar e_{k}|\leq |\bar e_{j(k)}| , because F_{j} and G_{j} monotonically increase with j , so

    \begin{eqnarray*} |\bar e_{k}|\leq |\bar e_{j(k)}|\leq F_{j(k)}G_{j(k)}\leq F_{k}G_{k}. \end{eqnarray*}

    (ii) When |\bar e_{k}|\geq |\bar e_{j(k)}| , from (4.21) we can obtain that

    \begin{eqnarray*} |\bar e_{k}|^{2}\leq |\bar e_{k}|\sum\limits_{i = 3}^{k-1}P_{k-i}^{k}\sum\limits_{j = 3}^{i}4L\rho^{i-j}|\bar e_{j}|+|\bar e_{k}|^{2}P_{0}^{k}\sum\limits_{j = 3}^{k}4L\rho^{i-j} +{ } |\bar e_{k}||\bar e_{2}|+2|\bar e_{k}|\sum\limits_{i = 3}^{k}P_{k-i}^{k}|\bar r_{i}(\Delta t)|. \end{eqnarray*}

    From (4.4), we give the limit of the maximum step size, which means that

    \begin{eqnarray*} \begin{array}{lll}{ } P_{0}^{k}\sum\limits_{j = 3}^{k}4L\rho^{i-j}\leq \pi_{A}\Gamma(2-\theta)\Delta t^{\theta}\Lambda < \frac{1}{2} \Longrightarrow \Delta t\leq \frac{1}{\sqrt[\theta]{2\pi_{A}\Gamma(2-\theta)\Lambda}}. \end{array} \end{eqnarray*}

    Using of Lemma 2.3 in [3] for any real \mu > 0 ,

    \begin{eqnarray*} \begin{array}{lll}{ } \sum\limits_{i = 3}^{k-1}P_{k-i}^{k}E_{\theta}(\mu t_{i}^{\theta})\leq \pi_{A}\frac{E_{\theta}(\mu t_{k}^{\theta})-1}{\mu}, \quad for \quad 3\leq k\leq K. \end{array} \end{eqnarray*}

    We have

    \begin{eqnarray*} |\bar e_{k}| &\leq& 2\sum\limits_{i = 3}^{k-1}P_{k-i}^{k}\sum\limits_{j = 3}^{i}4L\rho^{i-j}|\bar e_{j}|+2|\bar e_{2}|+4\sum\limits_{i = 3}^{k}P_{k-i}^{k}|\bar r_{i}(\Delta t)| \leq 2\Lambda \sum\limits_{i = 3}^{k-1}P_{k-i}^{k}F_{i}G_{k}+2G_{k}\\ & \leq & 2\sum\limits_{i = 3}^{k-1}P_{k-i}^{k}\sum\limits_{j = 3}^{i}4L \theta^{i-j}F_{j}G_{j}+2G_{k} \leq[4\Lambda\sum\limits_{i = 3}^{k-1}P_{k-i}^{k}E_{\theta}(2\pi_{A}\Lambda t_{i}^{\theta})+2]G_{k}\\ &\leq & \{4\Lambda \cdot \frac{\pi_{A}}{2\pi_{A}\Lambda}[E_{\theta}(2\pi_{A}t_{k}^{\theta})-1]+2\}G_{k} = 2E_{\theta}(2\pi_{A}\Lambda t_{k}^{\theta})G_{k} = F_{k}G_{k}. \end{eqnarray*}

    To sum up, the proof of (4.22) is complete.

    Next, we carefully estimate |\bar r_{i}(\Delta t)| . According to (4.17), we have

    \begin{eqnarray} |\bar r_{k}(\Delta t)|\leq { } |r_{k}(\Delta t)|+L\sum\limits_{j = 1}^{2}\rho^{k-j}|\bar e_{j}|+(\bar A_{k-1}^{k}+\bar A_{k-2}^{k})|\bar e_{1}|+\bar A_{k-2}^{k}|\bar e_{2}|. \end{eqnarray} (4.23)

    Next, we estimate the upper bounds of \bar A_{k-1}^{(k)} and \bar A_{k-2}^{(k)} . According to (2.37), Lemma 2.1 and (4) in Lemma A.1, we can obtain

    \bar A_{k-2}^{k} \\ = \sum\limits_{j = 2}^{k-2}A_{k-j}^{k}\rho^{j-2}+\rho^{k-3}\bar A_{1}^{k} \leq A_{2}^{k}\sum\limits_{j = 2}^{k-2}\rho^{j-2}+\rho^{k-3}\bar A_{1}^{k} \leq 3A_{2}^{k}+\frac{1}{\Delta t^{\theta}\Gamma(3-\theta)}[\frac{\theta}{4}-1+(3-\frac{\theta}{2})2^{-\theta}]\\ = \frac{1}{\Delta t^{\theta}\Gamma(3-\theta)}\Big[\frac{5}{4}(4-\theta)+\frac{9}{2}(8-\theta)3^{-\theta} +\frac{11}{2}(-6+\theta)\cdot2^{-\theta}\Big]\\ \leq\frac{1}{\Delta t^{\theta}\Gamma(3-\theta)}\Big[5+\frac{9}{2}\cdot 8+\frac{11}{2}(-6+1)\cdot\frac{1}{2}\Big] = \frac{1}{\Delta t^{\theta}\Gamma(3-\theta)}\frac{109}{4}.

    Therefore, \bar A_{k-2}^{k}\leq \frac{109}{4\Delta t^{\theta}\Gamma(3-\theta)} . According to Lemma 2.2, we can get \bar A_{k-1}^{k} < \bar{A}_{k-2}^{k}\leq \frac{109}{4\Delta t^{\theta}\Gamma(3-\theta)} . Therefore, we can obtain from (4.23),

    \begin{eqnarray*} |\bar r_{i}(\Delta t)| &\leq& C\Delta t^{3-\theta}+5LC\Delta t^{3}+\frac{109}{2\Gamma(3-\theta)}\Delta t^{3-\theta}+\frac{545}{12\Gamma(3-\theta)}C\Delta t^{3-\theta}\\ &\leq& \Big[C+5LC\Delta t^{\theta}+\frac{100}{\Gamma(3-\theta)}C\Big]\Delta t^{3-\theta}. \end{eqnarray*}

    In the expression on the right side of (4.22),

    \begin{eqnarray} \begin{array}{lll}{ } \sum\limits_{i = 3}^{k}P_{j-i}^{j}|\bar r_{i}(\Delta t)| = \sum\limits_{i = 3}^{k}P_{j-i}^{j}\omega_{1-\theta}~~(t_{i})\cdot \max\limits_{3\leq i\leq k}\frac{|\bar r_{i}(\Delta t)|}{\omega_{1-\theta}~~(t_{i})}, \end{array} \end{eqnarray} (4.24)

    according to (4.5) in Lemma 4.1, we can obtain that

    \begin{eqnarray*} \begin{array}{lll}{ } \sum\limits_{i = 3}^{k}P_{j-i}^{j}\omega_{1-\theta}~~(t_{i})\leq \sum\limits_{i = 1}^{k}P_{j-i}^{j}\omega_{1-\theta}~~(t_{i})\leq \pi_{A}, \quad 1\leq k \leq K. \end{array} \end{eqnarray*}

    Because of \frac{1}{\omega_{1-\theta}~~(t_{i})} = \Gamma(1-\theta)t_{i}^{\theta} , therefore (4.24) becomes

    \sum\limits_{i = 3}^{k}P_{j-i}^{j}|\bar r_{i}(\Delta t)|\leq \pi_{A}\max\limits_{3\leq i\leq k}\frac{|\bar r_{i}(\Delta t)}{\omega_{1-\theta}~~(t_{i})} \\\leq { } \pi_{A}\Big[C+5LC\Delta t^{\theta}+\frac{100}{\Gamma(3-\theta)}C\Big]\Delta t^{3-\theta}\max\limits_{3\leq i\leq k}\Gamma(1-\theta)t_{i}^{\theta}. (4.25)

    Based on (4.13), (4.25) and (4.22), we obtain

    \begin{eqnarray*} \begin{array}{lll}{ } |\bar e_{k}|\leq 2E_{\theta}(2\pi_{A}\Lambda t_{k}^{\theta})\Big\{\frac{5}{3}C\Delta t^{\theta}+2\pi_{A}\Big[C+5LC\Delta t^{\theta}+\frac{100}{\Gamma(3-\theta)}C\Big]\Gamma(1-\theta)t_{k}^{\theta}\Big\}\Delta t^{3-\theta}. \end{array} \end{eqnarray*}

    According to e_{k} = \bar e_{k}-\rho e_{k-1} = \sum_{n = 0}^{k}\rho^{k-n}\bar e_{n} , we have |e_{k}|\leq 3 \max_{0\leq n \leq k } |\bar e_{n}| . Then the proof is completed.

    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

    \begin{eqnarray} &&f(t, u(t)) = \frac{\Gamma(4+\theta)}{6}t^{3}+t^{3+\theta}-u(t), \end{eqnarray} (5.1)
    \begin{eqnarray} &&f(t, u(t)) = \frac{\Gamma(4+\theta)}{6}t^{3}+t^{6+2\theta}-u^{2}(t). \end{eqnarray} (5.2)

    The exact solution of Eq (2.1) is u(t) = t^{3+\theta} . We take T = 1 and choose the step size to be \Delta t = 2^{-l}, l = 7, 8, \ldots, 11 . The error we will display is defined by e_{\Delta t} = { }\max_{k = 1, 2, \cdots, K}|u(t_{k})-u_{k}|, K = T/\Delta t.

    In (5.1), f is a linear case of u . From Table 1, the convergence order is computed by \log_2(e_{2\Delta t}/e_{\Delta t}) . it is observed that for \theta = 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 \theta = 0.01 and 0.99 and obtain that the convergence rates are close to 2.99 and 2.01 . That tells us that as \theta\rightarrow0 or 1 , the convergence rate still 3-\theta . 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-\theta for 0 < \theta < 1 .

    Table 1.  Maximum errors e_{\Delta t} and decay rate as functions of \Delta t with the right hand side (5.1).
    \Delta t \theta=0.3 Rate \theta=0.6 Rate \theta=0.9 Rate
    \frac{1}{128} 5.1545E-6 - 5.3113E-5 - 3.8546E-4 -
    \frac{1}{256} 8.1757E-7 2.6564 1.0192E-5 2.3815 9.0841E-5 2.0851
    \frac{1}{512} 1.2874E-7 2.6668 1.9457E-6 2.3891 2.1297E-5 2.0926
    \frac{1}{1024} 2.0164E-8 2.6746 3.7032E-7 2.3934 4.9804E-6 2.0963
    \frac{1}{2048} 3.1714E-9 2.6686 7.0366E-8 2.3958 1.1632E-6 2.0981

     | Show Table
    DownLoad: CSV
    Table 2.  Maximum errors e_{\Delta t} and decay rate as \Delta t = 0.01, 0.99 with the right hand side (5.1).
    \Delta t \theta=0.01 Rate \theta=0.99 Rate
    \frac{1}{128} 3.7469E-8 - 6.6342E-4 -
    \frac{1}{256} 5.1473E-9 2.8638 1.6645E-4 1.9947
    \frac{1}{512} 7.0275E-10 2.8727 4.1540E-5 2.0025
    \frac{1}{1024} 8.9656E-11 2.9705 1.0339E-5 2.0063
    \frac{1}{2048} 1.1231E-11 2.9968 2.5703E-6 2.0081

     | Show Table
    DownLoad: CSV
    Table 3.  Maximum errors e_{\Delta t} and decay rate as functions of \Delta t with the right hand side (5.2).
    \Delta t \theta=0.3 Rate \theta=0.6 Rate \theta=0.9 Rate
    \frac{1}{128} 1.0154E-6 - 1.0167E-5 - 9.1822E-5 -
    \frac{1}{256} 1.4941E-7 2.7647 1.9494E-6 2.3828 2.1549E-5 2.0911
    \frac{1}{512} 2.3628E-8 2.6607 3.7176E-7 2.3906 5.0416E-6 2.0957
    \frac{1}{1024} 3.7049E-9 2.6729 7.0707E-8 2.3944 1.1777E-6 2.0978
    \frac{1}{2048} 5.6279E-10 2.7187 1.3497E-8 2.3891 2.7492E-7 2.0989

     | Show Table
    DownLoad: CSV
    Table 4.  Maximum errors e_{\Delta t} and decay rate as \Delta t = 0.01, 0.99 with the right hand side (5.2).
    \Delta t \theta=0.01 Rate \theta=0.99 Rate
    \frac{1}{128} 3.2116E-6 - 1.7072E-4 -
    \frac{1}{256} 3.9867E-7 3.0099 4.2645E-5 2.0012
    \frac{1}{512} 4.9490E-8 3.0099 1.0618E-5 2.0058
    \frac{1}{1024} 6.1435E-9 3.0099 2.6400E-6 2.0079
    \frac{1}{2048} 7.6264E-10 3.0099 6.5589E-7 2.0090

     | Show Table
    DownLoad: CSV

    Example 5.2. we consider the problem (2.1) and (2.2) with the following right hand side function

    \begin{eqnarray} &&f(t, u(t)) = t^{-\theta} Sin_{1, 1-\alpha}(t)+\sin(t)-u(t), \end{eqnarray} (5.3)
    \begin{eqnarray} &&f(t, u(t)) = t^{-\theta} Sin_{1, 1-\alpha}(t)+\sin^2(t)-u^2(t). \end{eqnarray} (5.4)

    The corresponding exact solution is u(t) = \sin(t) , and Sin_{\alpha, \beta}(t) = \sum_{k = 1}^\infty (-1)^{k+1}\frac{t^{2k-1}}{\Gamma(\alpha(2k-1)+\beta)} . 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-\theta for 0 < \theta < 1 .

    Table 5.  Maximum errors e_{\Delta t} and decay rate as functions of \Delta t with the right hand side (5.3).
    \Delta t \theta=0.3 Rate \theta=0.6 Rate \theta=0.9 Rate
    \frac{1}{128} 4.4720E-7 - 3.3763E-6 - 2.1961E-5 -
    \frac{1}{256} 7.1175E-8 2.6515 6.4755E-7 2.3823 5.1721E-6 2.0861
    \frac{1}{512} 1.1232E-8 2.6636 1.2353E-7 2.3901 1.2121E-6 2.0931
    \frac{1}{1024} 1.7643E-9 2.6704 2.3497E-8 2.3942 2.8338E-7 2.0967
    \frac{1}{2048} 2.6806E-10 2.7185 4.4390E-9 2.4041 6.6105E-8 2.0999

     | Show Table
    DownLoad: CSV
    Table 6.  Maximum errors e_{\Delta t} and decay rate as functions of \Delta t with the right hand side (5.4).
    \Delta t \theta=0.3 Rate \theta=0.6 Rate \theta=0.9 Rate
    \frac{1}{128} 5.2408E-7 - 3.5338E-6 - 2.0939E-5 -
    \frac{1}{256} 8.5515E-8 2.6155 6.8794E-7 2.3608 4.9822E-6 2.0713
    \frac{1}{512} 1.3708E-8 2.6411 1.3230E-7 2.3784 1.1738E-6 2.0855
    \frac{1}{1024} 2.1730E-9 2.6573 2.5279E-8 2.3878 2.7517E-7 2.0928
    \frac{1}{2048} 3.4105E-10 2.6716 4.8110E-9 2.3935 6.4316E-8 2.0970

     | Show Table
    DownLoad: CSV

    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-\theta for 0 < \theta < 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].

    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.

    The authors declare that there is no conflict of interests regarding the publication of this article.

    Lemma A.1.

    (1)\; F_{1} = { } -\frac{3}{2}(2-\theta)(k-2)^{1-\theta}+\frac{1}{2}(2-\theta)(k-3)^{1-\theta}\\-2(2-\theta)k^{1-\theta} + { } (k-3)^{2-\theta}-3(k-2)^{2-\theta}+2k^{2-\theta} > 0,\\ \; F_{2} = { } \frac{1}{2}(2-\theta)k^{1-\theta}+\frac{3}{2}(2-\theta)(k-2)^{1-\theta}-\frac{3}{2}(2-\theta)(k-3)^{1-\theta} +\frac{1}{2}(2-\theta)(k-4)^{1-\theta}\\ \; - { } k^{2-\theta}+3(k-2)^{2-\theta}-3(k-3)^{2-\theta}+(k-4)^{2-\theta} > 0,\\ (2)\; \frac{3}{2}(4-\theta)-(4+\theta)3^{1-\theta}+(6-\theta)2^{-\theta} > 0,\\ (3)\; -\frac{3}{4}+\frac{5\theta-4}{2(4-\theta)}\cdot3^{1-\theta} +\frac{(4+\theta)(6-\theta)}{(4-\theta)^2}\cdot3^{1-\theta}\cdot2^{-\theta} -\frac{(6-\theta)^2}{(4-\theta)^2}(2^{-\theta})^2 > 0,\\ (4)\; 0 < \rho = \frac{1}{2}(3+\frac{\theta-6}{4-\theta}\cdot2^{1-\theta}) < \frac{2}{3},\\ (5)\; \tilde f(\theta) = -3(4-\theta)^2+6(4-\theta)(6-\theta)2^{1-\theta}-4(4-\theta)(8-\theta)3^{1-\theta} +(6-\theta)^2\cdot4^{1-\theta} > 0,
    (1)\; F_{1} = { } -\frac{3}{2}(2-\theta)(k-2)^{1-\theta}+\frac{1}{2}(2-\theta)(k-3)^{1-\theta}\\-2(2-\theta)k^{1-\theta} + { } (k-3)^{2-\theta}-3(k-2)^{2-\theta}+2k^{2-\theta} > 0,\\ \; F_{2} = { } \frac{1}{2}(2-\theta)k^{1-\theta}+\frac{3}{2}(2-\theta)(k-2)^{1-\theta}-\frac{3}{2}(2-\theta)(k-3)^{1-\theta} +\frac{1}{2}(2-\theta)(k-4)^{1-\theta}\\ \; - { } k^{2-\theta}+3(k-2)^{2-\theta}-3(k-3)^{2-\theta}+(k-4)^{2-\theta} > 0,\\ (2)\; \frac{3}{2}(4-\theta)-(4+\theta)3^{1-\theta}+(6-\theta)2^{-\theta} > 0,\\ (3)\; -\frac{3}{4}+\frac{5\theta-4}{2(4-\theta)}\cdot3^{1-\theta} +\frac{(4+\theta)(6-\theta)}{(4-\theta)^2}\cdot3^{1-\theta}\cdot2^{-\theta} -\frac{(6-\theta)^2}{(4-\theta)^2}(2^{-\theta})^2 > 0,\\ (4)\; 0 < \rho = \frac{1}{2}(3+\frac{\theta-6}{4-\theta}\cdot2^{1-\theta}) < \frac{2}{3},\\ (5)\; \tilde f(\theta) = -3(4-\theta)^2+6(4-\theta)(6-\theta)2^{1-\theta}-4(4-\theta)(8-\theta)3^{1-\theta} +(6-\theta)^2\cdot4^{1-\theta} > 0,

    Proof. (1) Firstly, we will prove F_{1} > 0 . let k-2 = x , Taylor formula is used to obtain the following results

    \begin{eqnarray*} \begin{array}{llll}{ } F_{1} = { } (2-\theta)(1-\theta)x^{-\theta-1}\sum\limits_{n = 0}^{+\infty}\prod\limits_{i = 0}^{n}(-\theta-i)(\frac{1}{x})^{n}\cdot a_{n}, \end{array} \end{eqnarray*}

    where a_{n} = \frac{\frac{1}{2}(n+1)\cdot(-1)^{n}-8(n+1)\cdot2^{n}}{(n+3)!} = \frac{\frac{1}{2}(n+1)[(-1)^{n}-16\cdot2^{n}]}{(n+3)!} < 0, so \sum_{n = 0}^{+\infty}\prod_{i = 0}^{n}(-\theta-i)(\frac{1}{x})^{n}\cdot a_{n} is an alternating series with positive first term, so satisfies \sum_{n = 0}^{+\infty}\prod_{i = 0}^{n}(-\theta-i)(\frac{1}{x})^{n}\cdot a_{n} > 0, so F_{1} > 0.

    Next, we will prove F_2 > 0 . Let k-2 = \bar{x} , Taylor formula is used to obtain the following results

    \begin{eqnarray*} F_{2} = (2-\theta)(1-\theta)\bar{x}^{-\theta-1}\sum\limits_{n = 0}^{+\infty}\prod\limits_{i = 0}^{n}(-\theta-i) (\frac{1}{\bar{x}})^{n}\frac{1}{(n+3)!}b_{n}, \end{eqnarray*}

    where b_{n} = \frac{3}{2}(n+1)(-1)^{n+1}+2(n-1)\cdot2^{n}[1+(-1)^{n}], b_{0} = -\frac{11}{2}, b_{1} = 3, when n\geq 2, n is odd, b_{n} > 0 ; n is even, b_{n} = (n-1)[4\cdot2^{n}-\frac{3}{2}]-3 > [8-\frac{3}{2}]-3 > 0 , so \sum_{n = 1}^{+\infty}\prod_{i = 0}^{n}(-\theta-i)(\frac{1}{\bar{x}})^{n}\frac{1}{(n+3)!}b_{n} is an alternating series with positive first term. So \sum_{n = 0}^{+\infty}\prod_{i = 0}^{n}(-\theta-i)(\frac{1}{\bar{x}})^{n}\frac{1}{(n+3)!}b_{n} > (-\theta)\frac{1}{3!}b_{0}+0 > 0. We get F_2 > 0 .

    (2) Let f(\theta) = \frac{3}{2}(4-\theta)-(4+\theta)3^{1-\theta}+(6-\theta)2^{-\theta}, through careful calculation, we get

    f^{\prime}(\theta) = -\frac{3}{2}-3^{1-\theta}[1+(4+\theta)(-\ln3)]+2^{-\theta}[-1+(6-\theta)(-\ln2)],

    f^{\prime\prime}(\theta) = 3^{-\theta}\cdot g(\theta), where g(\theta) = 3(\ln3)[2-(4+\theta)\ln3]+(\frac{3}{2})^{\theta}(\ln2)[2+(6-\theta)\ln2], because

    \begin{eqnarray*} g^{\prime}(\theta)&\!\!\! = &\!\!\!-3(\ln3)^2+(\frac{3}{2})^{\theta}(\ln2)\{[2+(6-\theta)(\ln2)]\ln({\frac{3}{2}})-\ln2\},\\ g^{\prime\prime}(\theta)&\!\!\! = &\!\!\!(\frac{3}{2})^{\theta}(\ln2)\ln(\frac{3}{2}) \{[2+(6-\theta)(\ln2)]\ln(\frac{3}{2})-2\ln2\}\doteq \{\frac{3}{2}\}^{\theta}(\ln2)\ln(\frac{3}{2})H(\theta), \end{eqnarray*}

    We can get H^{\prime}(\theta) = -\ln2\cdot \ln(\frac{3}{2}) < 0, so we get H(\theta) > H(1) > 0 , therefore g^{\prime}(\theta) < g^{\prime}(1) < 0, g(\theta) monotone decreasing, g(1) < g(\theta) < g(0) = -3.6227 < 0 , that is f^{\prime}(\theta) monotone decreasing, 0 < f^{\prime}(1) < f^{\prime}(\theta) < f^{\prime}(0) , where f^{\prime}(1) = -3+5\ln3-\frac{5}{2}\ln2 > 0 , that is f(\theta) monotone increasing, f(0) < f(\theta) < f(1) , where f(0) = 0 , to sum up, we can get f(\theta) > 0 .

    (3) Let f_{1}(\theta) = -\frac{3}{4}+\frac{5\theta-4}{2(4-\theta)}\cdot3^{1-\theta} +\frac{(4+\theta)(6-\theta)}{(4-\theta)^2}\cdot3^{1-\theta}\cdot2^{-\theta} -\frac{(6-\theta)^2}{(4-\theta)^2}(2^{-\theta})^2 \doteq -\frac{3}{4}+a_{1}\cdot3^{1-\theta}+a_{2}\cdot3^{1-\theta}2^{-\theta}+a_{3}\cdot(2^{-\theta})^{2}, where

    a_{1} = \frac{5\theta-4}{2(4-\theta)},a_{2}\\ = \frac{(4+\theta)(6-\theta)}{(4-\theta)^2},a_{3} = -\frac{(6-\theta)^2}{(4-\theta)^2}.

    Next, we just need to prove that f_{1}(\theta) > 0, using a Taylor expansion yields

    \begin{array}{lll}{ } f_{1}(\theta) = { } -\frac{3}{4}+[a_{1}\cdot2^{1-\theta}+a_{2}\cdot2^{1-2\theta}](1+\frac{1}{2})^{1-\theta}+a_{3}\cdot(2^{-\theta})^{2} =\\ { } -\frac{3}{4}+a_{1}\cdot2^{1-\theta}+a_{2}\cdot2^{1-2\theta}+a_{3}\cdot2^{-2\theta}\\ { }+[a_{1}\cdot2^{1-\theta}+a_{2}\cdot2^{1-2\theta}][\frac{1-\theta}{1!}\frac{1}{2}+\frac{(1-\theta)(-\theta)}{2!}(\frac{1}{2})^{2}+\frac{(1-\theta)(-\theta)(-\theta-1)}{3!}(\frac{1}{2})^{3}+\cdots], \end{array}

    where

    \begin{eqnarray*} a_{1}2^{1-\theta}+a_{2}2^{1-2\theta}& = & 2^{-\theta}\cdot\frac{1}{(4-\theta)^{2}}[(5\theta-4)(4-\theta)+2(4+\theta)(6-\theta)2^{-\theta}]\\ &\geq&2^{-\theta}\cdot\frac{1}{(4-\theta)^{2}}[(5\theta-4)(4-\theta)+(4+\theta)(6-\theta)] \geq0. \end{eqnarray*}

    Because of \frac{1-\theta}{1!}\cdot\frac{1}{2}+\frac{(1-\theta)(-\theta)}{2!}\cdot(\frac{1}{2})^{2} +\frac{(1-\theta)(-\theta)(-\theta-1)}{3!}(\frac{1}{2})^{3}+\ldots\doteq\sum_{k = 0}^{+\infty}a_{k} is an alternating series with positive first term, and \sum_{k = 0}^{+\infty}a_{k} = a_{0}+a_{1}+\sum_{k = 2}^{+\infty}a_{k} , where \sum_{k = 2}^{+\infty}a_{k} is an alternating series with positive first term, so 0 < \sum_{k = 2}^{+\infty}a_{k} < a_{2} , we have

    \begin{eqnarray*} \begin{array}{lll}{ } &f_{1}(\theta) > { } -\frac{3}{4}+a_{3}\cdot2^{-2\theta}+[a_{1}\cdot2^{1-\theta}+a_{2}\cdot2^{1-2\theta}][1+\frac{1-\theta}{1!}\cdot\frac{1}{2}+\frac{(1-\theta)(-\theta)}{2!}(\frac{1}{2})^{2}]\\ & = { } -\frac{3}{4}+\frac{2^{-\theta}}{8(4-\theta)^{2}}[-5\theta^4+49\theta^3-196\theta^2+368\theta-192+(-\theta^4+7\theta^3-2\theta^2-48\theta+144)2^{1-\theta}]\\ &\doteq { } -\frac{3}{4}+\frac{2^{-\theta}}{8(4-\theta)^{2}}[f_{2}(\theta)+f_{3}(\theta)]. \end{array} \end{eqnarray*}

    We wish prove f_{1}(\theta) > 0 , only need -\frac{3}{4}+\frac{2^{-\theta}}{8(4-\theta)^{2}}[f_{2}(\theta)+f_{3}(\theta)] > 0\Longleftrightarrow f_{2}(\theta)+f_{3}(\theta) > \frac{3}{4}\cdot8(4-\theta)^{2}2^{\theta} = 6(4-\theta)^{2}2^{\theta} , that is f_{2}(\theta)+f_{3}(\theta) > 6(4-\theta)^{2}2^{\theta}\doteq6f_{4}(\theta) , so to prove f_{1}(\theta) > 0 , just prove f_{2}(\theta)+f_{3}(\theta)-6f_{4}(\theta) > 0 , let's remember \bar{f}(\theta) = f_{2}(\theta)+f_{3}(\theta)-6f_{4}(\theta) , since \bar{f}(\theta) first increases and then decreases, and the two endpoints are \bar{f}(0) = 0, \bar{f}(1) = 16 > 0 , \bar{f}(\theta) > 0 is always true, so f_{1}(\theta) > 0 .

    (4) Firstly, we will prove \rho > 0 . Because \rho = \frac{3(4-\theta)+(\theta-6)2^{1-\theta}}{2(4-\theta)}, let \bar{g}(\theta) = 3(4-\theta)+(\theta-6)2^{1-\theta} , by calculation, we can get \bar{g}(\theta) is monotonically increasing function on \theta , that is \bar{g}(\theta) > \bar{g}(0) = 0 . therefore, we have \rho > 0.

    In addition, \rho-\frac{2}{3} = \frac{5(4-\theta)+3(\theta-6)2^{1-\theta}}{6(4-\theta)}\doteq\frac{\bar{g}_{1}(\theta)}{6(4-\theta)}. we can directly find \bar{g}_{1}(\theta) < \bar{g}_{1}(1) = 0 , so 0 < \rho < \frac{2}{3} .

    (5) We use Lagrange's mean value theorem 4^{1-\theta} > \frac{4}{3+\theta}\cdot 3^{1-\theta} , and we have

    \begin{eqnarray} \tilde f(\theta) > { } -3(4-\theta)^{2}+6(4-\theta)(6-\theta)2^{1-\theta} +[-4(4-\theta)(8-\theta)+\frac{4(6-\theta)^{2}}{3+\theta}] 3^{1-\theta}. \end{eqnarray} (A.1)

    Let a_{1} = -3(4-\theta)^{2}, a_{2} = 6(4-\theta)(6-\theta), a_{3} = \frac{4}{3+\theta}[-(4-\theta)(8-\theta)(3+\theta)+(6-\theta)^{2}], a_{4} = 1+\frac{1-\theta}{1!}\frac{1}{2}+\frac{(1-\theta)(-\theta)}{2}(\frac{1}{2})^{2} . (A.1) is becomes

    \begin{eqnarray*} \tilde f(\theta)& > & a_{1}+a_{2}\cdot 2^{1-\theta}+a_{3}\cdot 3^{1-\theta} = a_{1}+2^{1-\theta}[a_{2}+a_{3}(1+\frac{1}{2})^{1-\theta}]\\ & = & a_{1}+2^{1-\theta}\{a_{2}+a_{3}[1+\frac{1-\theta}{1!}\frac{1}{2}+\frac{(1-\theta)(-\theta)}{2}(\frac{1}{2})^{2}]\\ &&+ a_{3}\cdot[\frac{(1-\theta)(-\theta)(-\theta-1)}{3!}(\frac{1}{2})^{3} +\frac{(1-\theta)(-\theta)(-\theta-1)(-\theta-2)}{4!}(\frac{1}{2})^{4}+\ldots]\}\\ & = & a_{1}+2^{1-\theta}\{a_{2}+a_{3}\cdot a_{4}+a_{3}(1-\theta)\theta\sum\limits_{k = 0}^{+\infty}\Pi_{i = 0}^{k}(-\theta-1-i)\cdot \frac{-1}{(k+3)!}(\frac{1}{2})^{k+3}\}\\ &\doteq & a_{1}+2^{1-\theta}\{a_{2}+a_{3}\cdot a_{4}+a_{3}(1-\theta)\theta\sum\limits_{k = 0}^{+\infty}b_{k}\}, \end{eqnarray*}

    where b_{k} = \Pi_{i = 0}^{k}(-\theta-1-i)\cdot \frac{-1}{(k+3)!}(\frac{1}{2})^{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

    \begin{eqnarray*} \tilde f(\theta)& > & a_{1}+2^{1-\theta}\{a_{2}+a_{3}\cdot a_{4}+a_{3}(1-\theta)\theta b_{0}\} = a_{1}+2^{1-\theta}\{a_{2}+a_{3}[a_{4}+(1-\theta)\theta\cdot \frac{-(-\theta-1)}{3!}(\frac{1}{2})^{3}]\}\\ &\doteq& a_{1}+2^{1-\theta}\{a_{2}+a_{3}\cdot a_{5}\} \doteq \tilde f_{1}(\theta), \end{eqnarray*}

    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

    \begin{eqnarray*} &\!\!\!\!\!\!&\!\!\!\!\!\!\tilde f_{1}(\theta) = \frac{-36}{12(3+\theta)}[\theta^{3}-5\theta^{2}-8\theta+48]\\ &\!\!\!\!\!\!&\!\!\!\!\!\!+ 2^{1-\theta}\frac{1}{12(3+\theta)}[\theta^{6}-16\theta^{5}+97\theta^{4}-278\theta^{3}+88\theta^{2}+732\theta+864]\\ &\!\!\!\!\!\!&\!\!\!\!\!\! = \frac{1}{12(3+\theta)}\{-36(\theta^{3}-5\theta^{2}-8\theta+48) +2^{1-\theta}(\theta^{6}-16\theta^{5}+97\theta^{4}-278\theta^{3}+88\theta^{2}+732\theta+864)\}\\ &\!\!\!\!\!\!&\!\!\!\!\!\!\doteq \frac{1}{12(3+\theta)}\tilde {f}_{2}(\theta), \end{eqnarray*}

    because \theta\in(0, 1), so \theta^{3} < \theta^{2}, \theta^{6} > 0.

    \begin{eqnarray*} &\!\!\!&\!\!\!\tilde {f}_{2}(\theta) > -36(\theta^{2}-5\theta^{2}-8\theta+48)+2^{1-\theta}(0-16\theta^{5}+97\theta^{4} -278\theta^{3}+88\theta^{2}+732\theta+864)\\ &\!\!\!&\!\!\! = 144(\theta^{2}+2\theta-12)+2^{1-\theta}(-16\theta^{5} +97\theta^{4}-278\theta^{3}+88\theta^{2}+732\theta+864)\\ &\!\!\!&\!\!\! = 144(\theta^{2}+2\theta-12)+2^{1-\theta}\tilde f_{4}(\theta) \doteq \tilde f_{3}(\theta), \end{eqnarray*}

    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 .



    [1] C. W. Lv, C. J. Xu, Error analysis of a high order method for time-fractional diffusion equations, SIAM J. Sci. Comput., 38 (2016), A2699–A2724. doi: 10.1137/15M102664X
    [2] J. Y. Cao, C. J. Xu, Z. Q. Wang, A high order finite difference/spectral approximations to the time fractional diffusion equations, Adv. Mater. Res., 875 (2014), 781–785.
    [3] J. Y. Cao, C. J. Xu, A high order schema for the nmerical solution of the fractional ordinary differential equations, J. Comput. Phys., 238 (2013), 154–168. doi: 10.1016/j.jcp.2012.12.013
    [4] H. L. Liao, W. McLean, J. W. Zhang, A discrete gronwall inequality with applications to numerical schemes for subdiffusion problems, SIAM J. Numer. Anal., 57 (2019), 218–237. doi: 10.1137/16M1175742
    [5] W. P. Bu, L. Ji, Y. F. Tang, J. Zhou, Space-time finite element method for the distributed-order time fractional reaction diffusion equations, Appl. Numer. Math., 152 (2020), 446–465. doi: 10.1016/j.apnum.2019.11.010
    [6] G. H. Gao, Z. Z. Sun, H. W. Zhang, A new fractional numerical differentiation formula to approximate the Caputo fractional derivative and its applications, J. Comput. Phys., 259 (2014), 33–50. doi: 10.1016/j.jcp.2013.11.017
    [7] C. P. Xie, S. M. Fang, Finite difference scheme for time fractional diffusion equation with fractional boundary conditions, Math. Method. Appl. Sci., 43 (2020), 3473–3487. doi: 10.1002/mma.6132
    [8] C. P. Li, M. Cai, High-order approximation to Caputo derivatives and Caputo-type advection-diffusion equations: revisited, Numer. Funct. Anal. Optim., 38 (2017), 861–890. doi: 10.1080/01630563.2017.1291521
    [9] J. Alessandra, S. M. Paola, On the numerical solutions of coupled nonlinear time-fractional reaction-diffusion equations, AIMS Mathematics, 6 (2021), 9109–9125. doi: 10.3934/math.2021529
    [10] Y. Zhang, X. B. Bao, L. B. Liu, Analysis of a finite difference scheme for a nonlinear Caputo fractional differential equation on an adaptive grid, AIMS Mathematics, 6 (2021), 8611–8624. doi: 10.3934/math.2021500
    [11] Z. J. Meng, M. X. Yi, J. Huang, L. Song, Numerical solutions of nonlinear fractional differential equations by alternative Legendre polynomials, Appl. Math. Comput., 336 (2018), 454–464.
    [12] S. Abbas, D. Mehdi, A new operational matrix for solving fractional-order differential equations, Comput. Math. Appl., 59 (2010), 1326–1336. doi: 10.1016/j.camwa.2009.07.006
    [13] M. R. Eslahchi, M. Dehghan, M. Parvizi, Application of the collocation method for solving nonlinear fractional integro-differential equations, J. Comput. Appl. Math., 257 (2014), 105–128. doi: 10.1016/j.cam.2013.07.044
    [14] Y. N. Zhang, Z. Z. Sun, X. Zhao, Compact alternating direction implicit scheme for the two-dimensional fractional diffusion-wave equation, SIAM J. Numer. Anal., 50 (2012), 1535–1555. doi: 10.1137/110840959
    [15] W. H. Luo, C. P. Li, T. Z. Huang, X. M. Gu, G. C. Wu, A high-order accurate numerical scheme for the Caputo derivative with applications to fractional diffusion problems, Numer. Funct. Anal. Optim., 39 (2018), 600–622. doi: 10.1080/01630563.2017.1402346
    [16] S. Lee, H. Kim, B. Jang, A fast and high-order numerical method for nonlinear fractional-order differential equations with non-singular kernel, Appl. Numer. Math., 163 (2021), 57–76. doi: 10.1016/j.apnum.2021.01.013
    [17] N. Khadijeh, D. Raziyeh, Galerkin finite element method for nonlinear fractional differential equations, Numer. Algorithms, 88 (2021), 113–141. doi: 10.1007/s11075-020-01032-2
    [18] Y. L. Guo, Z. Q. Wang, An hp-version Chebyshev collocation method for nonlinear fractional differential equations, Appl. Numer. Math., 158 (2020), 194–211. doi: 10.1016/j.apnum.2020.08.003
    [19] Z. D. Gu, Spectral collocation method for nonlinear Riemann-Liouville fractional differential equations, Appl. Numer. Math., 157 (2020), 654–669. doi: 10.1016/j.apnum.2020.07.003
    [20] I. Podlubny, Fractional differential equations, New York: Academic Press, 1999.
    [21] H. Xi, Y. Gu, Generalized finite difference method for electroelastic analysis of three-dimensional piezoelectric structures, Appl. Math. Lett., 117 (2021), 107084. doi: 10.1016/j.aml.2021.107084
    [22] Y. Gu, C. M. Fan, Z. J. Fu, Localized method of fundamental solutions for three-dimensional elasticity problems: theory, Adv. Appl. Math. Mech., 13 (2021), 1520–1534. doi: 10.4208/aamm.OA-2020-0134
  • This article has been cited by:

    1. Xiaopeng Yi, Chongyang Liu, Huey Tyng Cheong, Kok Lay Teo, Song Wang, A third-order numerical method for solving fractional ordinary differential equations, 2024, 9, 2473-6988, 21125, 10.3934/math.20241026
  • Reader Comments
  • © 2021 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(2353) PDF downloads(132) Cited by(1)

Figures and Tables

Tables(6)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog