images
Computer Modeling in Engineering & Sciences
images

DOI: 10.32604/cmes.2023.026009

ARTICLE

A Study on the Nonlinear Caputo-Type Snakebite Envenoming Model with Memory

Pushpendra Kumar1,*, Vedat Suat Erturk2, V. Govindaraj1 and Dumitru Baleanu3,4,5

1Department of Mathematics, National Institute of Technology Puducherry, Karaikal, 609609, India
2Department of Mathematics, Faculty of Arts and Sciences, Ondokuz Mayis University, Atakum, Samsun, 55200, Turkey
3Department of Mathematics, Cankaya University, Ankara, 06530, Turkey
4Institute of Space Sciences, Magurele, Bucharest, R76900, Romania
5Department of Computer Science and Mathematics, Lebanese American University, Beirut, 11022801, Lebanon
*Corresponding Author: Pushpendra Kumar. Email: kumarsaraswatpk@gmail.com
Received: 09 August 2022; Accepted: 10 November 2022

Abstract: In this article, we introduce a nonlinear Caputo-type snakebite envenoming model with memory. The well-known Caputo fractional derivative is used to generalize the previously presented integer-order model into a fractional-order sense. The numerical solution of the model is derived from a novel implementation of a finite-difference predictor-corrector (L1-PC) scheme with error estimation and stability analysis. The proof of the existence and positivity of the solution is given by using the fixed point theory. From the necessary simulations, we justify that the first-time implementation of the proposed method on an epidemic model shows that the scheme is fully suitable and time-efficient for solving epidemic models. This work aims to show the novel application of the given scheme as well as to check how the proposed snakebite envenoming model behaves in the presence of the Caputo fractional derivative, including memory effects.

Keywords: Mathematical model; Caputo fractional derivative; L1-predictor-corrector method; error estimation; stability; graphical simulations

1  Introduction

Nowadays, fractional calculus [1–3] is being applied to solve various real-world problems in terms of mathematical modeling. Different fractional derivatives [4] have been successfully used to model various problems. More specifically, several deadly epidemics have been modeled by using mathematical models in a fractional-order sense. It is a well-known fact that the fractional-order operators are non-local in nature and may be more effective for modeling history (memory)-dependent systems. Moreover, a fractional order can be fixed as any positive real number that better fits a real-data. So, by using such an operator, an accurate adjustment can be done in a model to fit with real data for better predicting the outbreaks of an epidemic. Recently, several applications of fractional derivatives have been recorded in epidemiology. In [5–7], the authors have studied the dynamics of the COVID-19 epidemic by using fractional-order mathematical models. In [8], the authors used non-singular Caputo-Fabrizio and Atangana-Baleanu fractional derivatives to study the dynamical nature of a malaria epidemic model. Kumar et al. in [9] have solved canine distemper virus and rabies epidemic models in the sense of generalized Caputo derivative. Kumar et al. [10] defined two different types of fractional-order models to study the dynamics of mosaic disease. A novel application of the generalized Caputo derivative in environmental infection dynamics can be seen in [11]. Sinan et al. [12] proposed the mathematical modeling of typhoid fever in terms of fractional derivatives. In [13], some novel analyses on the numerical modeling of biological systems using fractional derivatives have been given. Rihan et al. [14] studied fractional-order predator-prey models, including delay with Holling type-II functional response. In [15], some novel applications of delay differential equations were proposed. Work on the application of fractional derivatives in ecology and psychology can be understood from [16,17], respectively.

To analyze the various types of fractional-order systems, several computational schemes have been proposed by researchers. Odibat et al. [18] derived a new generalized form of the predictor-corrector (PC) scheme to investigate fractional initial value problems. Kumar et al. [19] introduced a new method to simulate fractional-order systems with various examples. In [20], the PC method was derived to simulate delay fractional differential equations. A modified form of the PC scheme in terms of the generalized Caputo derivative to solve delay-type systems has been introduced in reference [21]. Odibat et al. [22] have derived the generalized differential transform method for solving fractional impulsive differential equations. In [23], some computational schemes to solve delayed parabolic and time-fractional partial differential equations have been proposed. Shah et al. [24] simulated the dynamics of some important fractional order differential equations. In [25], some novel analyses of the Cauchy-type fractional-order dynamical system in terms of piecewise equations have been given. Some more applications of fractional derivatives can be seen in [26–28].

Jhinga et al. [29] introduced a novel finite-difference predictor-corrector (L1-PC) scheme to solve fractional-order systems in the sense of the Caputo derivative. That L1-PC scheme is not yet applied to any epidemic model to check how it will work and show the disease dynamics. In this paper, we fill this research gap by implementing the given L1-PC method on a Caputo-type snakebite envenoming (SBE) mathematical model, which was previously proposed in the integer-order sense in reference [30]. Currently, SBE is a deathly neglected disease, mainly in developing nations. The World Health Organization (WHO) recognized SBE as a fatal disease that generally results from the injection of a combination of various toxins (venom) following a venomous snakebite. Mainly poor, rural societies in tropical and subtropical nations all over the world are typically affected by SBE, which is a big threat to the health and well-being of about 5.8 billion population [31]. More information on SBE can be collected from the references [32–34]. To date, only a few research studies have gone into the mathematical modeling of SBE. In [35], the authors introduced a model using the law of mass action to analyze snakebite incidence. In [36], a mathematical model considering the socio-demographic components that influence the death risk from SBE in India was proposed.

The motivation of this study is to justify the application of the aforementioned L1-PC scheme in epidemiology by solving a mathematical model of SBE. Also, the aim is to check how the proposed SBE model will behave in the presence of the Caputo fractional derivative, which is a non-local differential operator that allows memory effects in the system. The given study is designed using the following systems: Section 2 recalls some preliminaries of fractional calculus. The description of the considered model of SBE in terms of the Caputo derivative is given in Section 3. The necessary numerical analysis, like the solution algorithm, error estimation, and method stability, is performed in Section 4. A brief discussion of the proposed methodology, its novelty, and outputs is given in Section 5 with supporting conclusions.

2  Preliminaries

Here we recall some preliminaries of fractional calculus.

Definition 1. A real function f(s), s>0 belongs to the space

a) Cη,η∈R if there exists a real number q>η, such that f(s) = sqf1(s), f1∈C[0,∞). Clearly,   Cη⊂Cα if α≤η.

b) Cηm,m∈N∪{0} if fm∈Cη.

Definition 2. [2] The Riemann-Liouville fractional integral of f(t)∈Cη(η≥−1) is defined as follows:

Jγf(t)=1Γ(γ)∫0t(t−s)γ−1f(s)ds,J0f(t)=f(t).

Definition 3. [2] The Caputo fractional derivative of f∈C−1m is written by

Dtγf(t)={dmf(t)dtm,ifγ=m∈N1Γ(m−γ)∫0t(t−ϑ)m−γ−1f(m)(ϑ)dϑ,ifm−1<γ<m,m∈N.(1)

Theorem 1. [37] Consider the function f such that f and CDc+γf are continuous for γ∈(0,1]. Then, ∀t∈(c,d], there exists some k∈(c,t) following the constraint

f(t)=f(c)+1Γ(γ+1)CDc+γf(k)(t−c)γ.

Therefore, from Theorem 1, we say that if CDc+γf(t)>0, for allt∈[c,d], then f is strictly increasing, and ifCDc+γf(t)<0, for all t∈[c,d], then the function f is strictly decreasing.

Lemma 1. [38] Let γ∈(0,1),n∈N and define the vectors X:=(x1,x2,…,xn) and Y:=(y1,y2,…,yn). For each i=1,2,…,n, let us consider Gi:[c,d]×Rn→R be a continuous function fulfils the Lipschitz condition with respect to the second variable such that

|Gi(t,X)−Gi(t,Y)|≤Li‖X−Y‖,

whereLi is a constant. Let us write G:=(G1,G2,…,Gn) and define the two fractional differential equations

CDc+γX(t)=G(t,X)+1jandCDc+γX(t)=G(t,X),(2)

with the same initial constraints, where j is a positive integer. If jX∗:=(jx1∗,…,jxn∗) andX∗:=(x1∗,x2∗,…,xn∗) are the solutions of (2), simultaneously, thenjX∗(t)→X∗(t) asj→∞ for all t∈[c,d].

3  Model Description

Here we introduce the generalized form of the previously published integer-order model [30] into Caputo-type fractional-order sense (CDγ). In our model, the total human population NH(t) at time t is split out into the following mutual exclusive classes: unaware susceptible humans, (SU(t)), aware susceptible humans, (SE(t)), SBE population, (I(t)), humans getting early remedy with antivenom, (TE(t)), peoples getting late remedy with antivenom, (TL(t)), peoples suffering from early adverse reaction (EAR) at the time of early remedy, (VE(t)), peoples facing EAR at the time of late remedy, (VL(t)), recovered peoples with disabilities, (RD(t)), and recovered peoples without disabilities, (RW(t)). Therefore, the total population size is defined as

N(t)=SU(t)+SE(t)+I(t)+TE(t)+TL(t)+VE(t)+VL(t)+RD(t)+RW(t).

The total snake population is defined by (NS(t)) and D(t) shows the total deaths caused by snakebite. For taking the equal time-dimension day−γ on the both sides of the fractional-order model, the power γ is applied to the parameters in time unit day−1 in the classical case. Therefore, the fractional-order nonlinear model for SBE is given as follows:

CDγSU=ΛHγ−(λ+K1)SU,CDγSE=ϵγSU+ϕ1γRD+ϕ2γRW−(Π1λ+K2)SE,CDγI=(Π1SE+SU)λ−K3I,CDγTE=τγkI−K4TE,CDγTL=τγΠ2I−K5TL,CDγVE=α1γTE−K6VE,CDγVL=α2γTL−K7VL,CDγRD=σ1γρ1TL+σ2γρ2VL−K8RD,CDγRW=r1γTE+r2γVE+σ1γΠ3TL+σ2γΠ4VL−K9RW,CDγNS=ΛSγNS(1−NSKS)−μSγNS,CDγD=δ1γI+(TL+VL)δ2γ,(3)

λ(t)=βγNSNH+NS,

where,

K1=ϵγ+μHγ,K2=μHγ,K3=τγ+δ1γ+μHγ,K4=α1γ+r1γ+μHγ,K5=α2γ+σ1γ+δ2γ+μHγ,K6=r2γ+μHγ,K7=σ2γ+δ2γ+μHγ,K8=ϕ1γ+μHγ,K9=ϕ2γ+μHγ,Π1=1−θ,Π2=1−k,Π3=1−ρ1,Π4=1−ρ2, with the initial conditions

SU(0)>0,SE(0)≥0,I(0)≥0,TE(0)≥0,TL(0)≥0,VE(0)≥0,VL(0)≥0,RD(0)≥0,RW(0)≥0,NS(0)≥0,D(0)≥0.(4)

The definitions of the model parameters alongwith their numerical values are given in Table 1.

Theorem 2. There exists a unique solution for the model (3) and (4) which belongs to(R0+)11:={(SU∗,SE∗,I∗,TE∗,TL∗,VE∗,VL∗,RD∗,RW∗,NS∗,D∗)∈R+11}.

Proof. By using the Theorem 3.1 and Remark 3.2 of [39], the global existence of the unique solution can easily be proved. Now to prove the non-negativity of the solution, let us write the following auxiliary system of fractional differential equations:

CDγSU=ΛHγ−(λ+K1)SU+1k,CDγSE=ϵγSU+ϕ1γRD+ϕ2γRW−(Π1λ+K2)SE+1k,CDγI=(Π1SE+SU)λ−K3I+1k,CDγTE=τγkI−K4TE+1k,CDγTL=τγΠ2I−K5TL+1k,CDγVE=α1γTE−K6VE+1k,CDγVL=α2γTL−K7VL+1k,CDγRD=σ1γρ1TL+σ2γρ2VL−K8RD+1k,CDγRW=r1γTE+r2γVE+σ1γΠ3TL+σ2γΠ4VL−K9RW+1k,CDγNS=ΛSγNS(1−NSKS)−μSγNS+1k,CDγD=δ1γI+(TL+VL)δ2γ+1k,(5)

with k∈N. We will show that solution of (5) (SUk∗(t),SEk∗(t),Ik∗(t),TEk∗(t),TLk∗(t),VEk∗(t), VLk∗(t),RDk∗(t),RWk∗(t),NSk∗(t),Dk∗(t)) belongs to (R0+)11, ∀t≥0. For obtaining a contradiction, we opine that there exists a point of time where the condition fails. Let t0:=inf{t>0|(SUk∗(t),SEk∗(t),Ik∗(t),TEk∗(t),TLk∗(t),VEk∗(t), VLk∗(t),RDk∗(t),RWk∗(t),NSk∗(t), Dk∗(t))∈(R0+)11}. Thus, (SUk∗(t0),SEk∗(t0),Ik∗(t0),TEk∗(t0),TLk∗(t0),VEk∗(t0),VLk∗(t0),RDk∗(t0), RWk∗(t0),NSk∗(t0), Dk∗(t0))∈(R0+)11 and one of the quantities (SUk∗(0),SEk∗(0),Ik∗(0),TEk∗(0),TLk∗(0),VEk∗(0), VLk∗(0),RDk∗(0), RWk∗(0),NSk∗(0),Dk∗(0)) is zero. Suppose that Ik∗(t0)=0. Since

CDγIk∗(t0)=(Π1SUk∗(t0)+SEk∗(t0))λ+1k>0

by continuity of CD0+γIk∗, we conclude that CD0+γIk∗([t0,t0+ξ)⊂R+, for some ξ>0. By Theorem 1, Ik∗([t0,t0+ξ)⊂R0+ and so Ik∗ is non negative. In an analogous way we can justify that the remaining functions SUk∗,SEk∗,TEk∗,TLk∗,VEk∗, VLk∗,RDk∗,RWk∗,NSk∗, and Dk∗ are non-negative, establishing a contradiction. Using Lemma 1 as k→∞, we get that the solution (SUk∗(t),SEk∗(t),Ik∗(t),TEk∗(t),TLk∗(t),VEk∗(t), VLk∗(t),RDk∗(t),RWk∗(t),NSk∗(t),Dk∗(t)) of (5) belongs to (R0+)11, ∀t≥0, giving the required result.

Now before performing the further numerical simulations on the proposed fractional-order model (3), we rewrite the model into a compact form by representing it in terms of an initial value problem, defined as follows:

Let us consider

{f1(t,SU,…,D)=ΛHγ−(λ+K1)SU,f2(t,SU,…,D)=ϵγSU+ϕ1γRD+ϕ2γRW−(Π1λ+K2)SE,f3(t,SU,…,D)=(Π1SE+SU)λ−K3I,f4(t,SU,…,D)=τγkI−K4TE,f5(t,SU,…,D)=τγΠ2I−K5TL,f6(t,SU,…,D)=α1γTE−K6VE,f7(t,SU,…,D)=α2γTL−K7VL,f8(t,SU,…,D)=σ1γρ1TL+σ2γρ2VL−K8RD,f9(t,SU,…,D)=r1γTE+r2γVE+σ1γΠ3TL+σ2γΠ4VL−K9RW,f10(t,SU,…,D)=ΛSγNS(1−NSKS)−μSγNS,f11(t,SU,…,D)=δ1γI+(TL+VL)δ2γ.(6)

By using (6), we have

CDγ𝒜(t)=Φ(t,𝒜(t)),t∈[0,T],0<γ≤1,𝒜(0)=𝒜0,(7)

where

𝒜(t)={SU(t)SE(t)I(t)TE(t)TL(t)VE(t)VL(t)RD(t)RW(t)NS(t)D(t),𝒜0(t)={SU0(t)SE0(t)I0(t)TE0(t)TL0(t)VE0(t)VL0(t)RD0(t)RW0(t)NS0(t)D0(t),Φ(t,𝒜(t))={f1(t,SU,…,D)f2(t,SU,…,D)f3(t,SU,…,D)f4(t,SU,…,D)f5(t,SU,…,D)f6(t,SU,…,D)f7(t,SU,…,D)f8(t,SU,…,D)f9(t,SU,…,D)f10(t,SU,…,D)f11(t,SU,…,D).(8)

4  Numerical Analysis on the Model

In this section, we perform the necessary numerical simulations (solution derivation, error estimation and stability) to derive the solution of the proposed fractional-order model (3) by using the L1-predictor-corrector scheme [29].

Consider the above given initial value problem (IVP) for 0<γ<1,

CDγ𝒜(t)=Φ(t,𝒜(t)),t∈[0,T],𝒜(0)=𝒜0.(9)

where CDγ represents the Caputo derivatives and Φ:[0,T]×D→R,D⊂R. Split the time span [0, T] into N subintervals. Take an uniform grid with step size of h=TN with tk=kh,k=0,1,…,N.

4.1 Derivation of the Solution

According to the L1-PC method, the Caputo fractional derivative is numerically defined by

[CDγ𝒜(t)]t=tn=1Γ(1−γ)∫0tk(tn−s)−γ𝒜′(s)ds=1Γ(1−γ)∑k=0n−1∫tktk+1(tn−s)−γ𝒜′(s)ds≈1Γ(1−γ)∑k=0n−1∫tktk+1(tn−s)−γ𝒜(tk+1)−𝒜(tk)hds=∑k=0n−1bn−k−1(𝒜(tk+1)−𝒜(tk)),(10)

where,

bk=h−γΓ(2−γ)[(k+1)1−γ−k1−γ].

We approximate CDγ𝒜(t) by the Eq. (10) and put it into (9) to get

[CDγ𝒜(t)]t=tn=∑k=0n−1bn−k−1(𝒜(tk+1)−𝒜(tk))=Φ(tn,𝒜n),(11)

where 𝒜k defines the approximate value of the solution of (9) at t=tk and

bn−k−1=h−γΓ(2−γ)[(n−k)1−γ−(n−k−)1−γ].

(11) can be rewritten as

bn−1(𝒜1−𝒜0)+bn−2(𝒜2−𝒜1)+⋯+b0(𝒜n−𝒜n−1)=Φ(tn,𝒜n).(12)

After rewriting the terms (12), we get the following from:

b0𝒜n=b0𝒜n−1−∑k=0n−2bk+1𝒜n−1−k+∑k=1n−1bk𝒜n−1−k+Φ(tn,𝒜n).(13)

Substituting

b0=h−γΓ(2−γ) and bk=h−γΓ(2−γ)[(n−k)1−γ−(n−k−)1−γ].

in (12), we get

𝒜n=𝒜n−1−(21−γ−11−γ)𝒜n−1−∑k=1n−2((2+k)1−γ−(1+k)(1−γ))𝒜n−1−k+∑k=1n−2((1+k)1−γ−(k)(1−γ))𝒜n−1−k+(n1−γ−(n−1)1−γ)𝒜0+Γ(2−γ)hγΦ(tn,𝒜n)=(n1−γ−(n−1)1−γ)𝒜0+∑k=1n−1[2(n−k)1−γ−(n+1−k)(1−γ)−(n−1−k)(1−γ)]𝒜k+Γ(2−γ)hγΦ(tn,𝒜n).(14)

Define

ak:=k+1(1−γ)−k(1−γ).(15)

Remark that ak′s have the following characteristics:

•   ak>0,k=0,1,…,n−1.

•   a0=1>a1>⋯>ak, and ak→0 as k→∞.

•   ∑k=0n−1(ak−ak+1)+an=(1−a1)+∑k=1n−2(ak−ak+1)+an−1=1.

Given Eqs. (14) and (15), the following form can be obtained:

𝒜n=an−1𝒜0+∑k=1n−1(an−1−k−an−k)𝒜k+Γ(2−γ)hγΦ(tn,𝒜n).(16)

We can see that Eq. (16) is of the form 𝒜n=g+N(𝒜n), if we identify

g=an−1𝒜0+∑k=1n−1(an−1−k−an−k)𝒜k

and

N(𝒜n)=Γ(2−γ)hγΦ(tn,𝒜n).

Hence using the scheme of the DGJ method gives approximate value of 𝒜n, given by

𝒜n,0=g=an−1𝒜0+∑k=1n−1(an−1−k−an−k)𝒜k,𝒜n,0=N(𝒜n,0)=Γ(2−γ)hγΦ(tn,𝒜n),𝒜n,2=N(𝒜n,0+𝒜n,1−N(𝒜n,0).

The three term approximation of 𝒜n≈𝒜n,0+𝒜n,0+𝒜n,2. Therefore, this approximated solution of DGJ scheme gives the following predictor-corrector algorithm called as L1-PCM.

𝒜np=an−1𝒜0+∑k=1n−1(an−1−k−an−k)𝒜k,znp=N(𝒜n)=Γ(2−γ)hγΦ(tn,𝒜np),𝒜nc=𝒜np+Γ(2−γ)hγΦ(tn,𝒜np+znp),(17)

where 𝒜np and znp are predictors and 𝒜nc is the corrector.

Using the above given methodology, the approximation equations of the proposed model (3) in terms of L1-PC method are derived as follows:

{SUnc=SUnp+Γ(2−γ)hγf1(tn,SUnp+z1np,…,Dnp+z11np),SEnc=SEnp+Γ(2−γ)hγf2(tn,SUnp+z1np,…,Dnp+z11np),Inc=Inp+Γ(2−γ)hγf3(tn,SUnp+z1np,…,Dnp+z11np),TEnc=TEnp+Γ(2−γ)hγf4(tn,SUnp+z1np,…,Dnp+z11np),TLnc=TLnp+Γ(2−γ)hγf5(tn,SUnp+z1np,…,Dnp+z11np),VEnc=VEnp+Γ(2−γ)hγf6(tn,SUnp+z1np,…,Dnp+z11np),VLnc=VLnp+Γ(2−γ)hγf7(tn,SUnp+z1np,…,Dnp+z11np),RDnc=RDnp+Γ(2−γ)hγf8(tn,SUnp+z1np,…,Dnp+z11np),RWnc=RWnp+Γ(2−γ)hγf9(tn,SUnp+z1np,…,Dnp+z11np),NSnc=NSnp+Γ(2−γ)hγf10(tn,SUnp+z1np,…,Dnp+z11np),Dnc=Dnp+Γ(2−γ)hγf11(tn,SUnp+z1np,…,Dnp+z11np),(18)

where,

{SUnp=an−1SU0+∑k=1n−1(an−1−k−an−k)SUk,SEnp=an−1SE0+∑k=1n−1(an−1−k−an−k)SEk,Inp=an−1I0+∑k=1n−1(an−1−k−an−k)Ik,TEnp=an−1TE0+∑k=1n−1(an−1−k−an−k)TEk,TLnp=an−1TL0+∑k=1n−1(an−1−k−an−k)TLk,VEnp=an−1VE0+∑k=1n−1(an−1−k−an−k)VEk,VLnp=an−1VL0+∑k=1n−1(an−1−k−an−k)VLk,RDnp=an−1RD0+∑k=1n−1(an−1−k−an−k)RDk,RWnp=an−1RW0+∑k=1n−1(an−1−k−an−k)RWk,NSnp=an−1NS0+∑k=1n−1(an−1−k−an−k)NSk,Dnp=an−1D0+∑k=1n−1(an−1−k−an−k)Dk,(19)

and

{z1np=N(SUn)=Γ(2−γ)hγf1(tn,SUnp,…,Dnp),z2np=N(SEn)=Γ(2−γ)hγf2(tn,SUnp,…,Dnp),z3np=N(In)=Γ(2−γ)hγf3(tn,SUnp,…,Dnp),z4np=N(TEn)=Γ(2−γ)hγf4(tn,SUnp,…,Dnp),z5np=N(TLn)=Γ(2−γ)hγf5(tn,SUnp,…,Dnp),z6np=N(VEn)=Γ(2−γ)hγf6(tn,SUnp,…,Dnp),z7np=N(VLn)=Γ(2−γ)hγf7(tn,SUnp,…,Dnp),z8np=N(RDn)=Γ(2−γ)hγf8(tn,SUnp,…,Dnp),z9np=N(RWn)=Γ(2−γ)hγf9(tn,SUnp,…,Dnp),z10np=N(NSn)=Γ(2−γ)hγf10(tn,SUnp,…,Dnp),z11np=N(Dn)=Γ(2−γ)hγf11(tn,SUnp,…,Dnp).(20)

The above given algorithm is coded in Mathematica and the bunch of Figs. 1 and 2 plotted to explore the outputs of the scheme at various fractional-order values. From the Figs. 1c–1f, we identify when the fractional order decreases, the population of the SBE humans; I(t), peoples getting early and late treatment with antivenom; TE(t) and TL(t), respectively, and the humans facing EAR at the time of early treatment; VE(t) increases, respectively. The time-respective variations in the rest of the model classes can be seen from the Figs. 2a–2e. Overall, we can say that the proposed L1-PC scheme performed well to explore the proposed model dynamics in the sense of the Caputo fractional derivative. The accuracy and stability of the given scheme are now discussed in our further analysis.

Figure 1: Dynamics of the model classes SU,SE,I,TE,TL,VE at fractional orders γ=1 (solid), γ=0.95 (dashed), γ=0.90 (dot-dashed), and γ=0.85 (dotted)

Figure 2: Dynamics of the model classes VL,RD,RW,NS,D at fractional orders γ=1 (solid), γ=0.95 (dashed), γ=0.90 (dot-dashed), and γ=0.85 (dotted)

4.2 Error Analysis

The brief analysis on the error estimation of L1-PC scheme has been given in the studies [29,40,41] and now investigated below. The error estimate is given by

|[CDγ𝒜(t)]t=tn−∑k=0n−1bn−k−1(𝒜k+1−𝒜k)|≤Ch2−γ,(21)

here C is a positive constant depends on γ and 𝒜.

Derive rn by

rn:=Γ(2−γ)hγ[[CDγ𝒜(t)]t=tn−∑k=0n−1bn−k−1(𝒜k+1−𝒜k)‖].(22)

In view of (21)

|rn|=Γ(2−γ)hγ|[CDγ𝒜(t)]t=tn−∑k=0n−1bn−k−1(𝒜k+1−𝒜k)‖|≤Γ(2−γ)Ch2.(23)

To derive the error estimation, we will use the lemmas given below:

Lemma 2. [42] For 0<γ<1 and ak′s (as given in Eq. (15)), we have

k−γak−1−1≤11−γ,k=1,2,…,N.

Lemma 3. [29] Consider 𝒜(tk) as exact solution of the proposed IVP and 𝒜kp be the approximate solution calculated from the algorithm (17). Then for 0<γ<1, we have

|𝒜(tk)−𝒜kp|≤Cak−1−1,k=1,2,…N,

where ak′s are given in Eq. (15).

Lemma 4. [29] Consider 𝒜(tk) as exact solution of the proposed IVP and 𝒜kp be the approximate value evaluated from Eq. (17). Then for 0<γ<1, we have

|𝒜(tk)−𝒜kp|≤CγTγh2−γ,k=1,2,⋯N,

where Cγ=C/(1−γ).

Theorem 3. Consider𝒜(t) as exact solution of the proposed IVP (9), Φ(t,𝒜(t)) satisfies the Lipschitz property respect to the variable𝒜 with a constant L, andΦ(t,𝒜(t)),𝒜(t)∈C1[0,T]. Also,𝒜kc defines the approximate solutions att=tk calculated by using L1-PC method. Then for0<γ<1, we have

|𝒜(tk)−𝒜kc|≤C1Tγh2−γ,k=1,2,…N,(24)

where C1=d/(1−γ) and d is a constant.

Proof. Let ek=𝒜(tk)−𝒜kc and ekp=𝒜(tk)−𝒜kp. Using Eqs. (9), (17), and (22), we get

en=enp+Γ(2−γ)hγ(Φ(tn,𝒜(tn))+N(𝒜(tn)))−Φ(tn,𝒜np+N((𝒜np))).

Further observe that

|en|≤|enp|+Γ(2−γ)hγ|Φ(tn,𝒜(tn))+N(𝒜(tn)))−Φ(tn,𝒜np+N((𝒜np)))|≤|enp|+LΓ(2−γ)hγ|𝒜(tn)−𝒜np+N(𝒜(tn)))−N((𝒜np))|≤|enp|+LΓ(2−γ)hγ|enp|+L2(Γ(2−γ))2h2γ|𝒜(tn)−𝒜np|≤|enp|+LΓ(2−γ)hγ|enp|+L2(Γ(2−γ))2h2γ|enp|≤[1+LΓ(2−γ)hγ+L2(Γ(2−γ))2h2γ]|enp|.(25)

Using Lemma 4 in Eq. (25), we get

|en|≤[1+LΓ(2−γ)hγ+L2(Γ(2−γ))2h2γ]CγTγh2−γ≤[1+LΓ(2−γ)+L2(Γ(2−γ))2]CγTγh2−γ.

Therefore

|en|≤CγTγh2−γ,

where C1 is a constant.

The behaviour of the absolute remainder error is plotted in the bundle of Figs. 3 and 4 at the previously used fractional-order values γ=1,0.95,0.90,0.85. For all classes, when the fractional order γ increases the error decreases as we move away from the left endpoint, namely 0. For SU,SE,VL,RD, and D classes, very small errors are obtained at the other points of the considered interval while the errors deteriorate at points very close to the origin for TE,RW and NS classes. Moreover, it is obvious that quite small errors are obtained for the integer-order case.

Figure 3: Absolute remainder error in the solution of the model classes SU,SE,I,TE,TL,VE at fractional orders γ=1 (solid), γ=0.95 (dashed), γ=0.90 (dot-dashed), and γ=0.85 (dotted)

Figure 4: Absolute remainder error in the solution of the model classes VL,RD,RW,NS,D at fractional orders γ=1 (solid), γ=0.95 (dashed), γ=0.90 (dot-dashed), and γ=0.85 (dotted)

4.3 Stability Analysis

The term stability means that small deviations in the initial values do not make the big changes in the numerical solutions. Consider that 𝒜nc and vnc(n=1,2…,N) are two solutions calculated by the numerical scheme (17). For δ0=|𝒜0−v0|, there exists two positive quantities k and h', such that

|𝒜nc−vnc|≤kδ0forh∈(0,h′),1≤n≤N,

here h is the step size given in Eq. (9).

Theorem 4. SupposeΦ(t,𝒜) follows the Lipschitz property respect to the variable𝒜 with a constant L and𝒜nc(n=1,2…,N) are the solutions established from the scheme (17), then the scheme (17) is stable.

Proof. We have to prove that

|𝒜nc−vnc|≤C|𝒜0−v0|.

Denote by η0:=(1+(LΓ(2−γ))+L2(Γ(2−γ))2hγ). Note that

|𝒜nc−vnc|≤|𝒜np−vnp|+LΓ(2−γ)hγ(|𝒜np−vnp|+|N(𝒜np)−N(vnp)|).(26)

Further observe that

|𝒜nc−vnc|=|an−1(𝒜0−v0)+∑k=1n−1(an−1−k−an−k)(𝒜k−vk)|≤an−1|𝒜0−v0|+∑k=1n−1(an−1−k−an−k)|𝒜k−vk|≤|𝒜0−v0|+∑k=1n−1(an−1−k−an−k)|𝒜k−vk|.

Using discrete form of Gronwalls inequality and Eq. (15), we obtain

|𝒜np−vnp|≤c|𝒜0−v0|,(27)

where c is a constant and

|N(𝒜np)−N(vnp)|=|Γ(2−γ)hγ(Φ(tn,𝒜np))−Φ(tn,vnp))|≤LΓ(2−γ)hγ|𝒜np−vnp|.(28)

Using (27) and (28) in (26), we get

|𝒜nc−vnc|≤|𝒜np−vnp|+LΓ(2−γ)hγ|𝒜np−vnp|+L2(Γ(2−γ))2h2γ|𝒜np−vnp|≤|𝒜np−vnp|+LΓ(2−γ)hγ|𝒜np−vnp|+L2(Γ(2−γ))2hγ|𝒜np−vnp|≤(1+LΓ(2−γ)+L2(Γ(2−γ))2)hγ|𝒜np−vnp|≤η0c|𝒜0−v0|≤C|𝒜0−v0|,

where C is a constant.

5  Conclusion

In this research work, we have performed a novel implementation of a finite-difference predictor-corrector scheme on a fractional-order nonlinear snakebite envenoming model in terms of the Caputo fractional derivative. The numerical solution of the model has been plotted by using several graphs to justify the behavior of the model at various fractional-order values. The analysis related to the error estimation and stability of the scheme has also been derived from exploring method’s accuracy. In the error estimation, the fractional order γ increases and the error decreases which justifies that our fractional-order analysis should not be performed at small values of order γ. From the given observations, it is clear that the proposed method can easily be implemented on various epidemic models, and the algorithm is time-efficient for more accurate solutions. In the future, the given scheme can be widely used to solve various nonlinear mathematical models related to real-life phenomena.

Funding Statement: The authors received no specific funding for this study.

Author Contributions: P. Kumar: Conceptualization, Visualization, Resources, Formal analysis, Investigation, Writing-original draft. V.S. Erturk: Conceptualization, Investigation, Visualization, Software, Writing-review & editing. V. Govindaraj: Investigation, Resources, Formal analysis, Writing-review & editing. D. Baleanu: Investigation, Resources, Formal analysis, Writing-review & editing.

Availability of Data and Materials: The data used in this research is available/mentioned in the manuscript.

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

References

1. Kilbas, A., Srivastava, H. M., Trujillo, J. J. (2006). Theory and applications of fractional differential equations. Netherlands: Elsevier Science.

2. Podlubny, I. (1998). Fractional differential equations: An introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications. USA: Elsevier.

3. Oldham, K., Spanier, J. (1974). The fractional calculus theory and applications of differentiation and integration to arbitrary order. USA: Elsevier.

4. Caputo, M., Fabrizio, M. (2015). A new definition of fractional derivative without singular kernel. Progress in Fractional Differentiation & Applications, 1(2), 73–85.

5. Erturk, V. S., Kumar, P. (2020). Solution of a COVID-19 model via new generalized caputo-type fractional derivatives. Chaos, Solitons & Fractals, 139, 110280. https://doi.org/10.1016/j.chaos.2020.110280

6. Ahmad, S., Ullah, A., Al-Mdallal, Q. M., Khan, H., Shah, K. et al. (2020). Fractional order mathematical modeling of COVID-19 transmission. Chaos, Solitons & Fractals, 139, 110256. https://doi.org/10.1016/j.chaos.2020.110256

7. Higazy, M. (2020). Novel fractional order SIDARTHE mathematical model of COVID-19 pandemic. Chaos, Solitons & Fractals, 138, 110007. https://doi.org/10.1016/j.chaos.2020.110007

8. Abboubakar, H., Kumar, P., Rangaig, N. A., Kumar, S. (2020). A malaria model with Caputo-Fabrizio and Atangana-Baleanu derivatives. International Journal of Modeling, Simulation, and Scientific Computing, 12(2), 2150013. https://doi.org/10.1142/S1793962321500136

9. Kumar, P., Erturk, V. S., Yusuf, A., Nisar, K. S., Abdelwahab, S. F. (2021). A study on canine distemper virus (CDV) and rabies epidemics in the red fox population via fractional derivatives. Results in Physics, 25, 104281. https://doi.org/10.1016/j.rinp.2021.104281

10. Kumar, P., Erturk, V. S., Almusawa, H. (2021). Mathematical structure of mosaic disease using microbial biostimulants via Caputo and Atangana-Baleanu derivatives. Results in Physics, 24, 104186. https://doi.org/10.1016/j.rinp.2021.104186

11. Kumar, P., Erturk, V. S. (2021). Environmental persistence influences infection dynamics for a butterfly pathogen via new generalised Caputo type fractional derivative. Chaos, Solitons & Fractals, 144, 110672. https://doi.org/10.1016/j.chaos.2021.110672

12. Sinan, M., Shah, K., Kumam, P., Mahariq, I., Ansari, K. J. et al. (2022). Fractional order mathematical modeling of typhoid fever disease. Results in Physics, 32, 105044. https://doi.org/10.1016/j.rinp.2021.105044

13. Rihan, F. A. (2013). Numerical modeling of fractional-order biological systems. Abstract and Applied Analysis, 2013. https://doi.org/10.1155/2013/816803

14. Rihan, F. A., Lakshmanan, S., Hashish, A. H., Rakkiyappan, R., Ahmed, E. (2015). Fractional-order delayed predator–prey systems with holling type-II functional response. Nonlinear Dynamics, 80(1), 777–789. https://doi.org/10.1007/s11071-015-1905-8

15. Rihan, F. A. (2021). Delay differential equations and applications to biology. Singapore: Springer.

16. Kumar, P., Erturk, V. S., Banerjee, R., Yavuz, M., Govindaraj, V. (2021). Fractional modeling of plankton-oxygen dynamics under climate change by the application of a recent numerical algorithm. Physica Scripta, 96(12), 124044. https://doi.org/10.1088/1402-4896/ac2da7

17. Kumar, P., Erturk, V. S., Murillo-Arcila, M. (2021). A complex fractional mathematical modeling for the love story of Layla and Majnun. Chaos, Solitons & Fractals, 150, 111091. https://doi.org/10.1016/j.chaos.2021.111091

18. Odibat, Z. Baleanu, D. (2020). Numerical simulation of initial value problems with generalized caputo-type fractional derivatives. Applied Numerical Mathematics, 156, 94–105. https://doi.org/10.1016/j.apnum.2020.04.015

19. Kumar, P., Erturk, V. S., Kumar, A. (2021). A new technique to solve generalized Caputo type fractional differential equations with the example of computer virus model. Journal of Mathematical Extension, 15, 1–23.

20. Bhalekar, S. Daftardar-Gejji, V. (2011). A predictor-corrector scheme for solving nonlinear delay differential equations of fractional order. Journal of Fractional Calculus and Applications, 1(5), 1–9.

21. Odibat, Z., Erturk, V. S., Kumar, P., Govindaraj, V. (2021). Dynamics of generalized Caputo type delay fractional differential equations using a modified predictor-corrector scheme. Physica Scripta, 96(12), 125213. https://doi.org/10.1088/1402-4896/ac2085

22. Odibat, Z., Erturk, V. S., Kumar, P., Ben Makhlouf, A., Govindaraj, V. (2022). An implementation of the generalized differential transform scheme for simulating impulsive fractional differential equations. Mathematical Problems in Engineering, 2022. https://doi.org/10.1155/2022/8280203

23. Rihan, F. A. (2010). Computational methods for delay parabolic and time-fractional partial differential equations. Numerical Methods for Partial Differential Equations, 26(6), 1556–1571. https://doi.org/10.1002/num.20504

24. Shah, K., Arfan, M., Ullah, A., Al-Mdallal, Q., Ansari, K. J. et al. (2022). Computational study on the dynamics of fractional order differential equations with applications. Chaos, Solitons & Fractals, 157, 111955. https://doi.org/10.1016/j.chaos.2022.111955

25. Shah, K., Abdeljawad, T., Ali, A. (2022). Mathematical analysis of the Cauchy type dynamical system under piecewise equations with caputo fractional derivative. Chaos, Solitons & Fractals, 161, 112356. https://doi.org/10.1016/j.chaos.2022.112356

26. Amin, M., Abbas, M., Baleanu, D., Iqbal, M. K., Riaz, M. B. (2021). Redefined extended cubic B-spline functions for numerical solution of time-fractional telegraph equation. Computer Modeling in Engineering & Sciences, 127(1), 361–384. https://doi.org/10.32604/cmes.2021.012720

27. Wang, Y., Veeresha, P., Prakasha, D. G., Baskonus, H. M., Gao, W. (2022). Regarding deeper properties of the fractional order Kundu-Eckhaus equation and massive thirring model. Computer Modeling in Engineering & Sciences, 133(3), 697–717. https://doi.org/10.32604/cmes.2022.021865

28. Dai, D., Li, X., Li, Z., Zhang, W., Wang, Y. (2023). Numerical simulation of the fractional-order lorenz chaotic systems with caputo fractional derivative. Computer Modeling in Engineering & Sciences, 135(2), 1371–1392. https://doi.org/10.32604/cmes.2022.022323

29. Jhinga, A., Daftardar-Gejji, V. (2018). A new finite-difference predictor-corrector method for fractional differential equations. Applied Mathematics and Computation, 336, 418–432. https://doi.org/10.1016/j.amc.2018.05.003

30. Abdullahi, S. A., Habib, A. G., Hussaini, N. (2021). Control of snakebite envenoming: A mathematical modeling study. PLoS Neglected Tropical Diseases, 15(8), e0009711. https://doi.org/10.1371/journal.pntd.0009711

31. World Health Organization. Snakebite envenoming: A strategy for prevention and control. Geneva. License: CC BY-NC-SA 3.0 IGO.2019.

32. Harrison, R. A., Hargreaves, A., Wagstaff, S. C., Faragher, B., Lalloo, D. G. (2009). Snake envenoming: A disease of poverty. PLoS Neglected Tropical Diseases, 3(12), e569. https://doi.org/10.1371/journal.pntd.0000569

33. Lancet, T. (2017). Snake-bite envenoming: A priority neglected tropical disease. Lancet, 390(10089), 2.

34. Chippaux, J. P. (2017). Snakebite envenomation turns again into a neglected tropical disease! Journal of Venomous Animals and Toxins Including Tropical Diseases, 23(1), 1–2. https://doi.org/10.1186/s40409-017-0127-6

35. Bravo-Vega, C. A., Cordovez, J. M., Renjifo-Ibáñez, C., Santos-Vega, M., Sasa, M. (2019). Estimating snakebite incidence from mathematical models: A test in Costa Rica. PLoS Neglected Tropical Diseases, 13(12), e0007914. https://doi.org/10.1371/journal.pntd.0007914

36. Kim, S. (2020). Introduction of a mathematical model to characterize relative risk of snakebite envenoming, a neglected tropical disease. Oregon State University, USA.

37. Diethelm, K. (2012). The mean value theorems and a nagumo-type uniqueness theorem for Caputo’s fractional calculus. Fractional Calculus and Applied Analysis, 15(2), 304–313. https://doi.org/10.2478/s13540-012-0022-3

38. Almeida, R., Brito da Cruz, A., Martins, N., Monteiro, M. T. T. (2019). An epidemiological MSEIR model described by the Caputo fractional derivative. International Journal of Dynamics and Control, 7(2), 776–784. https://doi.org/10.1007/s40435-018-0492-1

39. Lin, W. (2007). Global existence theory and chaos control of fractional differential equations. Journal of Mathematical Analysis and Applications, 332(1), 709–726. https://doi.org/10.1016/j.jmaa.2006.10.040

40. Langlands, T. A. M., Henry, B. I. (2005). The accuracy and stability of an implicit solution method for the fractional diffusion equation. Journal of Computational Physics, 205(2), 719–736. https://doi.org/10.1016/j.jcp.2004.11.025

41. Sun, Z. Z., Wu, X. (2006). A fully discrete difference scheme for a diffusion-wave system. Applied Numerical Mathematics, 56(2), 193–209. https://doi.org/10.1016/j.apnum.2005.03.003

42. Lin, Y., Xu, C. (2007). Finite difference/spectral approximations for the time-fractional diffusion equation. Journal of Computational Physics, 225(2), 1533–1552. https://doi.org/10.1016/j.jcp.2007.02.001

images
This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.