Loading [MathJax]/jax/output/SVG/jax.js
Research article Special Issues

Within-host delay differential model for SARS-CoV-2 kinetics with saturated antiviral responses

  • The present study discussed a model to describe the SARS-CoV-2 viral kinetics in the presence of saturated antiviral responses. A discrete-time delay was introduced due to the time required for uninfected epithelial cells to activate a suitable antiviral response by generating immune cytokines and chemokines. We examined the system's stability at each equilibrium point. A threshold value was obtained for which the system switched from stability to instability via a Hopf bifurcation. The length of the time delay has been computed, for which the system has preserved its stability. Numerical results show that the system was stable for the faster antiviral responses of epithelial cells to the virus concentration, i.e., quick antiviral responses stabilized patients' bodies by neutralizing the virus. However, if the antiviral response of epithelial cells to the virus increased, the system became unstable, and the virus occupied the whole body, which caused patients' deaths.

    Citation: Kaushik Dehingia, Anusmita Das, Evren Hincal, Kamyar Hosseini, Sayed M. El Din. Within-host delay differential model for SARS-CoV-2 kinetics with saturated antiviral responses[J]. Mathematical Biosciences and Engineering, 2023, 20(11): 20025-20049. doi: 10.3934/mbe.2023887

    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
  • The present study discussed a model to describe the SARS-CoV-2 viral kinetics in the presence of saturated antiviral responses. A discrete-time delay was introduced due to the time required for uninfected epithelial cells to activate a suitable antiviral response by generating immune cytokines and chemokines. We examined the system's stability at each equilibrium point. A threshold value was obtained for which the system switched from stability to instability via a Hopf bifurcation. The length of the time delay has been computed, for which the system has preserved its stability. Numerical results show that the system was stable for the faster antiviral responses of epithelial cells to the virus concentration, i.e., quick antiviral responses stabilized patients' bodies by neutralizing the virus. However, if the antiviral response of epithelial cells to the virus increased, the system became unstable, and the virus occupied the whole body, which caused patients' deaths.



    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

    (e1e2)=ΔtθD1(L1e1r1(Δt)L2e2r2(Δt))ΔtθD1L1e1r1(Δt)L2e2r2(Δt)ΔtθD1(Lmax{|e1|,|e2|}+max{|r1(Δt)|,|r2(Δt)|}), (4.11)

    where D=(ˉD11ˉD12D11D12). Since ˉD11D12ˉD12D11 is bounded, then

    max{|e1|,|e2|}C1Δtθmax{|r1(Δt)|,|r2(Δt)|}, (4.12)

    therefore,

    {|e1|CΔt3θ,|e2|CΔt3θ, (4.13)

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

    When i3, let's set ˉei=eiρei1,i=1,...,j, and ˉe0=e0 for (4.13), we have

    {|ˉe1|=|e1ρe0|CΔt3θ,|ˉe2|=|e2ρe1|53CΔt3θ. (4.14)

    Combining with (2.20) and (4.6), ˉej, j3, satisfy

    kj=1ˉAkkjΔˉej=f(tk,u(tk))f(tk,uk)rk(Δt). (4.15)

    Because of ej=jn=1ρjnˉen, then

    kj=3ˉAkkjΔˉej=Lkkj=3ρkjˉejˉrk(Δt), (4.16)

    where

    ˉrk(Δt)=rk(Δt)Lk2j=1ρkjˉej+2j=1ˉAkkjΔˉej, (4.17)

    by multiplying 2ˉek on both sides of (4.16), from Lemma A.1 in Liao [3], we can obtain that

    2ˉekkj=3ˉAkkjΔˉejkj=3ˉAkkjΔ(|ˉej|2),2ˉek[Lkkj=3ρkjˉejˉrk(Δt)]Lkj=3ρkj(|ˉej|2+|ˉek|2)+2|ˉek||ˉrk(Δt)|4Lkj=3ρkj|ˉej|2+2|ˉek||ˉrk(Δt)|.

    To sum up, we can get the following conclusions

    kj=3ˉAkkjΔ(|ˉej|2)4Lkj=3ρkj|ˉej|2+2|ˉek||ˉrk(Δt)|. (4.18)

    We change k to i, and then multiply both sides of (4.18) by a ki=3Pkki, then

    ki=3Pkkiij=3ˉAiijΔ(|ˉej|2)ki=3Pkkiij=34Lρij|ˉej|2+2ki=3Pkki|ˉei||ˉri(Δt)|. (4.19)

    From (4.2) and Δuj=ujuj1, we can get the following results

    ki=3Pkkiij=3ˉAiijΔ(|ˉej|2)=kj=3Δ(|ˉej|2)ki=jPkkiˉAiij=kj=3Δ(|ˉej|2)=|ˉek|2|ˉe2|2, (4.20)

    substituting (4.20) into (4.19) yields

    |ˉek|24Lki=3Pkkiij=3ρij|ˉej|2+|e2|2+2ki=3Pkki|ˉei||ˉri(Δt)|. (4.21)

    Next, we will prove it by mathematical induction

    |ˉek|2Eθ(2πAΛtθk)(|ˉe2|+2max3jkki=3pjji|ˉri(Δt)|)FkGk,fork2, (4.22)

    where Λ=12L, Eθ stands for the Mittag-Leffler function, we usually define it as: Eθ(z)=j=0ZjΓ(1+jθ), and Eθ(0)=1,ZR, and Z>0,Eθ(Z)>0, so FkFk12 for all k2. where Gk=|ˉe2|+2max3jkki=3Pjji|ˉri(Δt)|.

    When k=3, (i) if |ˉe3||e2|,G3=|e2|+2P30|ˉr3(Δt)||ˉe2|, we have |ˉe3||ˉe2|G3F3G3,

    (ii) if |ˉe3|>|e2|, from (4.21), it can be concluded that

    |ˉe3|24LP30|ˉe3|2+|ˉe2|2+2P30|ˉe3||ˉr3(Δt)|.

    According to (4.4), it can be concluded that

    4LP30πAΓ(2θ)ΔtθΛ12Δt1θ2πAΓ(2θ)Λ,

    so |ˉe3|2(|ˉe2|+2P30|ˉr3(Δt)|)2G3F3G3.

    When 4kK, we assume that: |ˉej|FjGj, for 3jk1, let |ˉej(k)|=max2ik1|ˉei|.

    (i) When |ˉek||ˉej(k)|, because Fj and Gj monotonically increase with j, so

    |ˉek||ˉej(k)|Fj(k)Gj(k)FkGk.

    (ii) When |ˉek||ˉej(k)|, from (4.21) we can obtain that

    |ˉek|2|ˉek|k1i=3Pkkiij=34Lρij|ˉej|+|ˉek|2Pk0kj=34Lρij+|ˉek||ˉe2|+2|ˉek|ki=3Pkki|ˉri(Δt)|.

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

    Pk0kj=34LρijπAΓ(2θ)ΔtθΛ<12Δt1θ2πAΓ(2θ)Λ.

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

    k1i=3PkkiEθ(μtθi)πAEθ(μtθk)1μ,for3kK.

    We have

    |ˉek|2k1i=3Pkkiij=34Lρij|ˉej|+2|ˉe2|+4ki=3Pkki|ˉri(Δt)|2Λk1i=3PkkiFiGk+2Gk2k1i=3Pkkiij=34LθijFjGj+2Gk[4Λk1i=3PkkiEθ(2πAΛtθi)+2]Gk{4ΛπA2πAΛ[Eθ(2πAtθk)1]+2}Gk=2Eθ(2πAΛtθk)Gk=FkGk.

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

    Next, we carefully estimate |ˉri(Δt)|. According to (4.17), we have

    |ˉrk(Δt)||rk(Δt)|+L2j=1ρkj|ˉej|+(ˉAkk1+ˉAkk2)|ˉe1|+ˉAkk2|ˉe2|. (4.23)

    Next, we estimate the upper bounds of ˉA(k)k1 and ˉA(k)k2. According to (2.37), Lemma 2.1 and (4) in Lemma A.1, we can obtain

    ˉAkk2=k2j=2Akkjρj2+ρk3ˉAk1Ak2k2j=2ρj2+ρk3ˉAk13Ak2+1ΔtθΓ(3θ)[θ41+(3θ2)2θ]=1ΔtθΓ(3θ)[54(4θ)+92(8θ)3θ+112(6+θ)2θ]1ΔtθΓ(3θ)[5+928+112(6+1)12]=1ΔtθΓ(3θ)1094.

    Therefore, ˉAkk21094ΔtθΓ(3θ). According to Lemma 2.2, we can get ˉAkk1<ˉAkk21094ΔtθΓ(3θ). Therefore, we can obtain from (4.23),

    |ˉri(Δt)|CΔt3θ+5LCΔt3+1092Γ(3θ)Δt3θ+54512Γ(3θ)CΔt3θ[C+5LCΔtθ+100Γ(3θ)C]Δt3θ.

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

    ki=3Pjji|ˉri(Δt)|=ki=3Pjjiω1θ  (ti)max3ik|ˉri(Δt)|ω1θ  (ti), (4.24)

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

    ki=3Pjjiω1θ  (ti)ki=1Pjjiω1θ  (ti)πA,1kK.

    Because of 1ω1θ  (ti)=Γ(1θ)tθi, therefore (4.24) becomes

    ki=3Pjji|ˉri(Δt)|πAmax3ik|ˉri(Δt)ω1θ  (ti)πA[C+5LCΔtθ+100Γ(3θ)C]Δt3θmax3ikΓ(1θ)tθi. (4.25)

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

    |ˉek|2Eθ(2πAΛtθk){53CΔtθ+2πA[C+5LCΔtθ+100Γ(3θ)C]Γ(1θ)tθk}Δt3θ.

    According to ek=ˉekρek1=kn=0ρknˉen, we have |ek|3max0nk|ˉen|. 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

    f(t,u(t))=Γ(4+θ)6t3+t3+θu(t), (5.1)
    f(t,u(t))=Γ(4+θ)6t3+t6+2θu2(t). (5.2)

    The exact solution of Eq (2.1) is u(t)=t3+θ. We take T=1 and choose the step size to be Δt=2l,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.

    Table 1.  Maximum errors eΔt and decay rate as functions of Δt with the right hand side (5.1).
    Δt θ=0.3 Rate θ=0.6 Rate θ=0.9 Rate
    1128 5.1545E-6 - 5.3113E-5 - 3.8546E-4 -
    1256 8.1757E-7 2.6564 1.0192E-5 2.3815 9.0841E-5 2.0851
    1512 1.2874E-7 2.6668 1.9457E-6 2.3891 2.1297E-5 2.0926
    11024 2.0164E-8 2.6746 3.7032E-7 2.3934 4.9804E-6 2.0963
    12048 3.1714E-9 2.6686 7.0366E-8 2.3958 1.1632E-6 2.0981

     | Show Table
    DownLoad: CSV
    Table 2.  Maximum errors eΔt and decay rate as Δt=0.01,0.99 with the right hand side (5.1).
    Δt θ=0.01 Rate θ=0.99 Rate
    1128 3.7469E-8 - 6.6342E-4 -
    1256 5.1473E-9 2.8638 1.6645E-4 1.9947
    1512 7.0275E-10 2.8727 4.1540E-5 2.0025
    11024 8.9656E-11 2.9705 1.0339E-5 2.0063
    12048 1.1231E-11 2.9968 2.5703E-6 2.0081

     | Show Table
    DownLoad: CSV
    Table 3.  Maximum errors eΔt and decay rate as functions of Δt with the right hand side (5.2).
    Δt θ=0.3 Rate θ=0.6 Rate θ=0.9 Rate
    1128 1.0154E-6 - 1.0167E-5 - 9.1822E-5 -
    1256 1.4941E-7 2.7647 1.9494E-6 2.3828 2.1549E-5 2.0911
    1512 2.3628E-8 2.6607 3.7176E-7 2.3906 5.0416E-6 2.0957
    11024 3.7049E-9 2.6729 7.0707E-8 2.3944 1.1777E-6 2.0978
    12048 5.6279E-10 2.7187 1.3497E-8 2.3891 2.7492E-7 2.0989

     | Show Table
    DownLoad: CSV
    Table 4.  Maximum errors eΔt and decay rate as Δt=0.01,0.99 with the right hand side (5.2).
    Δt θ=0.01 Rate θ=0.99 Rate
    1128 3.2116E-6 - 1.7072E-4 -
    1256 3.9867E-7 3.0099 4.2645E-5 2.0012
    1512 4.9490E-8 3.0099 1.0618E-5 2.0058
    11024 6.1435E-9 3.0099 2.6400E-6 2.0079
    12048 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

    f(t,u(t))=tθSin1,1α(t)+sin(t)u(t), (5.3)
    f(t,u(t))=tθSin1,1α(t)+sin2(t)u2(t). (5.4)

    The corresponding exact solution is u(t)=sin(t), and Sinα,β(t)=k=1(1)k+1t2k1Γ(α(2k1)+β). 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.

    Table 5.  Maximum errors eΔt and decay rate as functions of Δt with the right hand side (5.3).
    Δt θ=0.3 Rate θ=0.6 Rate θ=0.9 Rate
    1128 4.4720E-7 - 3.3763E-6 - 2.1961E-5 -
    1256 7.1175E-8 2.6515 6.4755E-7 2.3823 5.1721E-6 2.0861
    1512 1.1232E-8 2.6636 1.2353E-7 2.3901 1.2121E-6 2.0931
    11024 1.7643E-9 2.6704 2.3497E-8 2.3942 2.8338E-7 2.0967
    12048 2.6806E-10 2.7185 4.4390E-9 2.4041 6.6105E-8 2.0999

     | Show Table
    DownLoad: CSV
    Table 6.  Maximum errors eΔt and decay rate as functions of Δt with the right hand side (5.4).
    Δt θ=0.3 Rate θ=0.6 Rate θ=0.9 Rate
    1128 5.2408E-7 - 3.5338E-6 - 2.0939E-5 -
    1256 8.5515E-8 2.6155 6.8794E-7 2.3608 4.9822E-6 2.0713
    1512 1.3708E-8 2.6411 1.3230E-7 2.3784 1.1738E-6 2.0855
    11024 2.1730E-9 2.6573 2.5279E-8 2.3878 2.7517E-7 2.0928
    12048 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θ 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].

    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)F1=32(2θ)(k2)1θ+12(2θ)(k3)1θ2(2θ)k1θ+(k3)2θ3(k2)2θ+2k2θ>0,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θ>0,(2)32(4θ)(4+θ)31θ+(6θ)2θ>0,(3)34+5θ42(4θ)31θ+(4+θ)(6θ)(4θ)231θ2θ(6θ)2(4θ)2(2θ)2>0,(4)0<ρ=12(3+θ64θ21θ)<23,(5)˜f(θ)=3(4θ)2+6(4θ)(6θ)21θ4(4θ)(8θ)31θ+(6θ)241θ>0,
    (1)F1=32(2θ)(k2)1θ+12(2θ)(k3)1θ2(2θ)k1θ+(k3)2θ3(k2)2θ+2k2θ>0,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θ>0,(2)32(4θ)(4+θ)31θ+(6θ)2θ>0,(3)34+5θ42(4θ)31θ+(4+θ)(6θ)(4θ)231θ2θ(6θ)2(4θ)2(2θ)2>0,(4)0<ρ=12(3+θ64θ21θ)<23,(5)˜f(θ)=3(4θ)2+6(4θ)(6θ)21θ4(4θ)(8θ)31θ+(6θ)241θ>0,

    Proof. (1) Firstly, we will prove F1>0. let k2=x, Taylor formula is used to obtain the following results

    F1=(2θ)(1θ)xθ1+n=0ni=0(θi)(1x)nan,

    where an=12(n+1)(1)n8(n+1)2n(n+3)!=12(n+1)[(1)n162n](n+3)!<0, so +n=0ni=0(θi)(1x)nan is an alternating series with positive first term, so satisfies +n=0ni=0(θi)(1x)nan>0, so F1>0.

    Next, we will prove F2>0. Let k2=ˉx, Taylor formula is used to obtain the following results

    F2=(2θ)(1θ)ˉxθ1+n=0ni=0(θi)(1ˉx)n1(n+3)!bn,

    where bn=32(n+1)(1)n+1+2(n1)2n[1+(1)n],b0=112,b1=3, when n2, n is odd, bn>0; n is even, bn=(n1)[42n32]3>[832]3>0, so +n=1ni=0(θi)(1ˉx)n1(n+3)!bn is an alternating series with positive first term. So +n=0ni=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(θ)=3231θ[1+(4+θ)(ln3)]+2θ[1+(6θ)(ln2)],

    f(θ)=3θg(θ), where g(θ)=3(ln3)[2(4+θ)ln3]+(32)θ(ln2)[2+(6θ)ln2], because

    g(θ)=3(ln3)2+(32)θ(ln2){[2+(6θ)(ln2)]ln(32)ln2},g(θ)=(32)θ(ln2)ln(32){[2+(6θ)(ln2)]ln(32)2ln2}{32}θ(ln2)ln(32)H(θ),

    We can get H(θ)=ln2ln(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+5ln352ln2>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θ)231θ2θ(6θ)2(4θ)2(2θ)2-34+a131θ+a231θ2θ+a3(2θ)2, where

    a1=5θ42(4θ),a2=(4+θ)(6θ)(4θ)2,a3=(6θ)2(4θ)2.

    Next, we just need to prove that f1(θ)>0, using a Taylor expansion yields

    f1(θ)=34+[a121θ+a2212θ](1+12)1θ+a3(2θ)2=34+a121θ+a2212θ+a322θ+[a121θ+a2212θ][1θ1!12+(1θ)(θ)2!(12)2+(1θ)(θ)(θ1)3!(12)3+],

    where

    a121θ+a2212θ=2θ1(4θ)2[(5θ4)(4θ)+2(4+θ)(6θ)2θ]2θ1(4θ)2[(5θ4)(4θ)+(4+θ)(6θ)]0.

    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

    f1(θ)>34+a322θ+[a121θ+a2212θ][1+1θ1!12+(1θ)(θ)2!(12)2]=34+2θ8(4θ)2[5θ4+49θ3196θ2+368θ192+(θ4+7θ32θ248θ+144)21θ]34+2θ8(4θ)2[f2(θ)+f3(θ)].

    We wish prove f1(θ)>0, only need 34+2θ8(4θ)2[f2(θ)+f3(θ)]>0f2(θ)+f3(θ)>348(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

    ˜f(θ)>3(4θ)2+6(4θ)(6θ)21θ+[4(4θ)(8θ)+4(6θ)23+θ]31θ. (A.1)

    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

    ˜f(θ)>a1+a221θ+a331θ=a1+21θ[a2+a3(1+12)1θ]=a1+21θ{a2+a3[1+1θ1!12+(1θ)(θ)2(12)2]+a3[(1θ)(θ)(θ1)3!(12)3+(1θ)(θ)(θ1)(θ2)4!(12)4+]}=a1+21θ{a2+a3a4+a3(1θ)θ+k=0Πki=0(θ1i)1(k+3)!(12)k+3}a1+21θ{a2+a3a4+a3(1θ)θ+k=0bk},

    where bk=Πki=0(θ1i)1(k+3)!(12)k+3, and b0=(θ1)3!(12)3>0, |bk+1bk|=θ+2+k2(k+4)<1, so +k=0bk is a convergent alternating series, that is 0<+k=0bk<b0. Because of a3<0, so we get

    ˜f(θ)>a1+21θ{a2+a3a4+a3(1θ)θb0}=a1+21θ{a2+a3[a4+(1θ)θ(θ1)3!(12)3]}a1+21θ{a2+a3a5}˜f1(θ),

    where a5=a4+(1θ)θ(θ1)3!(12)3=48+24(1θ)6(θθ2)+(θθ3)48, by careful calculation, we can get

    ˜f1(θ)=3612(3+θ)[θ35θ28θ+48]+21θ112(3+θ)[θ616θ5+97θ4278θ3+88θ2+732θ+864]=112(3+θ){36(θ35θ28θ+48)+21θ(θ616θ5+97θ4278θ3+88θ2+732θ+864)}112(3+θ)˜f2(θ),

    because θ(0,1), so θ3<θ2,θ6>0.

    ˜f2(θ)>36(θ25θ28θ+48)+21θ(016θ5+97θ4278θ3+88θ2+732θ+864)=144(θ2+2θ12)+21θ(16θ5+97θ4278θ3+88θ2+732θ+864)=144(θ2+2θ12)+21θ˜f4(θ)˜f3(θ),

    where ˜f4(θ)=16θ5+97θ4278θ3+88θ2+732θ+864, by careful calculation, we can get ˜f3(θ)<0, that is ˜f3(θ) monotone decreasing, so ˜f3(1)<˜f3(θ)<˜f3(0), by direct calculation, we can get ˜f3(θ)=144(2θ+2)+21θ[˜f4(θ)(ln2)+˜f4(θ)],˜f3(0)>0,˜f3(0)<0, so ˜f3(0) changes from positive to negative, ˜f3(0) increases first and then decreases, ˜f3(0)=0,˜f3(1)=223>0, that is ˜f3(θ)>0 established, so ˜f(θ)>0.



    [1] K. Liang, Mathematical model of infection kinetics and its analysis for COVID-19, SARS and MERS, Infect. Genet. Evol., 82 (2020), 104306. https://doi.org/10.1016/j.meegid.2020.104306 doi: 10.1016/j.meegid.2020.104306
    [2] C. Yang, J. Wang, A mathematical model for the novel coronavirus epidemic in Wuhan, China, Math. Biosci. Eng., 17 (2020), 2708–2724. https://doi.org/10.3934/mbe.2020148 doi: 10.3934/mbe.2020148
    [3] S. K. Biswas, J. K. Ghosh, S. Sarkar, U. Ghosh, COVID-19 pandemic in India: a mathematical model study, Nonlinear Dyn., 102 (2020), 537–553. https://doi.org/10.1007/s11071-020-05958-z doi: 10.1007/s11071-020-05958-z
    [4] A. A. Arjani, Md. T. Nasseef, S. M. Kamal, B. V. S. Rao, M. Mahmud, Md. S. Uddin, Application of mathematical modeling in prediction of COVID-19 transmission dynamics, Arab. J. Sci. Eng., 47 (2022), 10163–10186. https://doi.org/10.1007/s13369-021-06419-4 doi: 10.1007/s13369-021-06419-4
    [5] A. J. Kucharski, T. W. Russell, C. Diamond, Early dynamics of transmission and control of COVID-19: a mathematical modelling study, Lancet Infect. Dis., 20 (2020), 553–558. https://doi.org/10.1016/S1473-3099(20)30144-4 doi: 10.1016/S1473-3099(20)30144-4
    [6] M. Zamir, F. Nadeem, M. A. Alqudah, T. Abdeljawad, Future implications of COVID-19 through Mathematical modeling, Results Phys., 33 (2022), 105097. https://doi.org/10.1016/j.rinp.2021.105097 doi: 10.1016/j.rinp.2021.105097
    [7] F. A. Rihan, H. J. Alsakaji, C. Rajivganthi, Stochastic SIRC epidemic model with time-delay for COVID-19, Adv. Differ. Equ., 502 (2020), 502(2020). https://doi.org/10.1186/s13662-020-02964-8 doi: 10.1186/s13662-020-02964-8
    [8] S. R. Bandekar, M. Ghosh, C. Rajivganthi, Impact of vaccination on the dynamics of COVID-19: A mathematical study using fractional derivatives, Int. J. Biomath., 17 (2024), 2350018. https://doi.org/10.1142/S1793524523500183 doi: 10.1142/S1793524523500183
    [9] H. J. Alsakaji, F. A. Rihan, A. Hashish, Dynamics of a stochastic epidemic model with vaccination and multiple time-delays for COVID-19 in the UAE, Complexity, 2022 (2022), 4247800. https://doi.org/10.1155/2022/4247800 doi: 10.1155/2022/4247800
    [10] C. Yang, Y. Yang, Z. Li, L. Zhang, Modeling and analysis of COVID-19 based on a time delay dynamic model, Math. Biosci. Eng., 18 (2020), 154–165. https://doi.org/10.3934/mbe.2021008 doi: 10.3934/mbe.2021008
    [11] M. Radha, S. Balamuralitharan, A study on COVID-19 transmission dynamics: stability analysis of SEIR model with Hopf bifurcation for effect of time delay, Adv. Differ. Equ., 2020 (2020), 523. https://doi.org/10.1186/s13662-020-02958-6 doi: 10.1186/s13662-020-02958-6
    [12] A. Raza, A. Ahmadian, M. Rafiq, M. C. Ang, S. Salahshour, M. Pakdaman, The impact of delay strategies on the dynamics of coronavirus pandemic model with nonlinear incidence rate, Fractals, 30 (2022), 2240121. https://doi.org/10.1142/S0218348X22401211 doi: 10.1142/S0218348X22401211
    [13] O. Babasola, O. Kayode, O. J. Peter, F. C. Onwuegbuche, F. A. Oguntolu, Time-delayed modelling of the COVID-19 dynamics with a convex incidence rate, Inform. Med. Unlocked., 35 (2022), 101124. https://doi.org/10.1016/j.imu.2022.101124 doi: 10.1016/j.imu.2022.101124
    [14] S. M. Al‑Tuwairqi, S. K. Al‑Harbi, A time‑delayed model for the spread of COVID‑19 with vaccination, Sci. Rep., 12 (2022), 19435. https://doi.org/10.1038/s41598-022-23822-5 doi: 10.1038/s41598-022-23822-5
    [15] S. Q. Du, W. Yuan, Mathematical modeling of interaction between innate and adaptive immune responses in COVID-19 and implications for viral pathogenesis, J. Med. Virol., 92 (2020), 1615–1628. https://doi.org/10.1002/jmv.25866 doi: 10.1002/jmv.25866
    [16] E. A. Hernandez-Vargas, J. X. Velasco-Hernandez, In-host Mathematical Modelling of COVID-19 in Humans, Annu. Rev. Control, 50 (2020), 448–456. https://doi.org/10.1016/j.arcontrol.2020.09.006 doi: 10.1016/j.arcontrol.2020.09.006
    [17] D. B. Prakash, D. K. K. Vamsi, D. B. Rajesh, C. B. Sanjeevi, Control Intervention Strategies for within-host, between-host and their efficacy in the treatment, spread of COVID-19 : a multi scale modeling approach, Comput. Math. Biophys., 8 (2020), 198–210. https://doi.org/10.1515/cmb-2020-0111 doi: 10.1515/cmb-2020-0111
    [18] B. Chhetri, V. M. Bhagat, D. K. K. Vamsi, V. S. Ananth, B. Prakash, R. Mandale, et al., Within-host mathematical modeling on crucial inflammatory mediators and drug interventions in COVID-19 identifies combination therapy to be most effective and optimal, Alex. Eng. J., 60 (2021), 2491–2512. https://doi.org/10.1016/j.aej.2020.12.011 doi: 10.1016/j.aej.2020.12.011
    [19] A. E. S. Almocera, G. Quiroz, E. A. Hernandez-Vargas, Stability analysis in COVID-19 within-host model with immune response. Commun. Nonlinear Sci. Numer. Simul., 95 (2021), 105584. https://doi.org/10.1016/j.cnsns.2020.105584 doi: 10.1016/j.cnsns.2020.105584
    [20] R. Ghostine, M. Gharamti, S. Hassrouny, I. Hoteit, Mathematical modeling of immune responses against SARS-CoV-2 using an ensemble Kalman Filter, Mathematics, 9 (2021), 2427. https://doi.org/10.3390/math9192427 doi: 10.3390/math9192427
    [21] S. M. E. K. Chowdhury, J. T. Chowdhury, S. F. Ahmed, P. Agarwal, I. A. Badruddin, S. Kamangar, Mathematical modelling of COVID-19 disease dynamics: interaction between immune system and SARS-CoV-2 within host, AIMS Math., 7 (2021), 2618–2633. https://doi.org/10.3934/math.2022147 doi: 10.3934/math.2022147
    [22] C. Li, J. Xu, J. Liu, Y. Zhou, The within-host viral kinetics of SARS-CoV-2, Math. Biosci. Eng., 17 (2020), 2853–2861. https://doi.org/10.3934/mbe.2020159 doi: 10.3934/mbe.2020159
    [23] B. J. Nath, K. Dehingia, V. N. Mishra, Y. M. Chu, H. K. Sarmah, Mathematical analysis of a within-host model of SARS-CoV-2, Adv. Differ. Equ., 2021 (2021), 113. https://doi.org/10.1186/s13662-021-03276-1 doi: 10.1186/s13662-021-03276-1
    [24] I. Ghosh, Within host dynamics of SARS‑CoV‑2 in humans: modeling immune responses and antiviral treatments, SN comput. sci., 2 (2021), 482. https://doi.org/10.1007/s42979-021-00919-8 doi: 10.1007/s42979-021-00919-8
    [25] I. M. Elbaz, M. A. Sohaly, H. El‑Metwally, Modeling the stochastic within‑host dynamics SARS‑CoV‑2 infection with discrete delay, Theor. Biosci., 141 (2020), 365–374. https://doi.org/10.1007/s12064-022-00379-5 doi: 10.1007/s12064-022-00379-5
    [26] V. Staroverov, S. Nersisyan, A. Galatenko, D. Alekseev, S. Lukashevich, F. Polyakov, et al., Development of a novel mathematical model that explains SARS-CoV-2 infection dynamics in Caco-2 cells, PeerJ, 11 (2023), e14828. https://doi.org/10.7717/peerj.14828 doi: 10.7717/peerj.14828
    [27] L. G. de Pillis, R. Caffrey, G. Chen, et al., A mathematical model of the within-host kinetics of SARS-CoV-2 neutralizing antibodies following COVID-19 vaccination, J. Theor. Biol., 556 (2023), 111280. doi:10.1016/j.jtbi.2022.111280 doi: 10.1016/j.jtbi.2022.111280
    [28] T. A. Miura, Respiratory epithelial cells as master communicators during viral infections, Curr. Clin. Microbiol., 6 (2019), 10–17. https://doi.org/10.1007/s40588-019-0111-8 doi: 10.1007/s40588-019-0111-8
    [29] A. N. Chatterjee, F. A. Basir, M. A. Almuqrin, J. Mondal, I. Khan, SARS-CoV-2 infection with lytic and non-lytic immune responses: a fractional order optimal control theoretical study, Results in Physics, 26 (2021), 104260. https://doi.org/10.1016/j.rinp.2021.104260 doi: 10.1016/j.rinp.2021.104260
    [30] H. J. Alsakaji, S. Kundu, F. A. Rihan, Delay differential model of one-predator two-prey system with Monod-Haldane and holling type Ⅱ functional responses, Appl. Math. Comput., 397 (2021). https://doi.org/10.1016/j.amc.2020.125919 doi: 10.1016/j.amc.2020.125919
    [31] X. Yang, L. Chen, J. Chen, Permanence and positive periodic solution for the single-species non-autonomous delay diffusive models, Comput. Math. Applic., 32 (1996), 109–116. https://doi.org/10.1016/0898-1221(96)00129-0 doi: 10.1016/0898-1221(96)00129-0
    [32] P. V. Driessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci., 180 (2002), 29–48. https://doi.org/10.1016/S0025-5564(02)00108-6 doi: 10.1016/S0025-5564(02)00108-6
    [33] J. P. LaSalle, The stability of dynamical systems, SIAM, Philadelphia, Pa, USA, 1976. https://doi.org/10.21236/ADA031020
    [34] Y. Kuang, Delay differential equations with applications in population dynamics, Academic Press, Inc.: Boston, MA, USA, 1993.
    [35] H. I. Freedman, V. S. H. Rao, Stability criteria for a system involving two-time delays, SIAM J. Appl. Math., 46 (1986), 552–560. https://doi.org/10.1137/0146037 doi: 10.1137/0146037
    [36] L. H. Erbe, H. I. Freedman, V. S. H. Rao, Three species food chain models with mutual interference and time delays, Math. Biosci., 80 (1986), 57–80. https://doi.org/10.1016/0025-5564(86)90067-2 doi: 10.1016/0025-5564(86)90067-2
    [37] H. I. Freedman, V. S. H. Rao, The trade-off between mutual interference and time lags in predator-prey systems, Bull. Math. Biol., 45 (2019), 1983,991–1004. https://doi.org/10.1016/S0092-8240(83)80073-1 doi: 10.1016/S0092-8240(83)80073-1
    [38] Y. Fadaei, F. A. Rihan, C. Rajivganthi, Immunokinetic model for COVID-19 patients, Complexity, 2022 (2022), 8321848. https://doi.org/10.1155/2022/8321848 doi: 10.1155/2022/8321848
    [39] G. Li, Y. Fan, Y. Lai, T. Han, Z. Li, P. Zhou, et al., Coronavirus infections and immune responses, J. Med. Virol., 92 (2020), 424–432. https://doi.org/10.1002/jmv.25685 doi: 10.1002/jmv.25685
  • 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
  • © 2023 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(1807) PDF downloads(74) Cited by(5)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog