Computer Systems Science & Engineering DOI:10.32604/csse.2022.019655
Article

Deterministic and Stochastic Fractional Order Model for Lesser Date Moth

Moustafa El-shahed1,* and Asmaa M. Al-Dububan2

1Department of Mathematics, Unaizah College of Sciences and Arts, Qassim University, P.O. Box 3771, Unaizah, 51911, Saudi Arabia
2Department of Mathematics, Buraydah College of Sciences and Arts, Qassim University, P.O. Box 1162, Buraydah, 51431, Saudi Arabia
*Corresponding Author: Moustafa El-shahed. Email: elshahedm@yahoo.com
Received: 21 April 2021; Accepted: 23 May 2021

Abstract: In this paper, a deterministic and stochastic fractional order model for lesser date moth (LDM) using mating disruption and natural enemies is proposed and analysed. The interaction between LDM larvae, fertilized LDM female, unfertilized LDM female, LDM male and the natural enemy is investigated. In order to clarify the characteristics of the proposed deterministic fractional order model, the analysis of existence, uniqueness, non-negativity and boundedness of the solutions of the proposed fractional-order model are examined. In addition, some sufficient conditions are obtained to ensure the local and global stability of equilibrium points. The occurrence of local bifurcation near the equilibrium points is investigated with the help of Sotomayor’s theorem. Numerical simulations are conducted to illustrate the properties of the proposed fractional order model with respect to the intrinsic growth rate of the LDM larvae, natural enemy’s mortality rate, predation rate, sex pheromone trap parameter, fractional order and environmental noise. The impact of mating disruption on lesser date moth is demonstrated. Also, a numerical approximation method is developed for the proposed stochastic fractional-order model.

Keywords: Lesser date moth; stochastic; stability; natural enemies; Sotomayor’s theorem; mating disruption

1  Introduction

The palm tree is considered one of the oldest major and basic crops in Southwest Asia, North Africa and many other places of the world. Palm trees are affected by many agricultural pests that cause significant losses to the palm trees and their fruits, as well as affecting the age and growth of the palm tree if it is left without control [1]. The LDM is one of the most dangerous pests of young palm trees and immature palm fruits. The damage is mainly due to the way in which the LDM larvae feed, as soon as they leave the egg until the date of their entry the virgin dwelling develops feeding and causing tunnels in the affected part of the palm in all directions and depths without any early signs showing the infection [2–5]. One of the most promising strategies for controlling LDM is the use of a mating disruption using the sex pheromone traps [5]. Natural enemies’ use to stop pest infestation has long been recommended as a clean and environmentally friendly way to protect crops. The main natural enemies that are used in agricultural pest control are larval predators. Goniozus swirskiana can be considered as one of the most important natural enemies that can attack the LDM [6,7]. In the real world, plant diseases models are always affected by the environmental noise. Thus, the stochastic models may be a more appropriate way of modelling agricultural pests in many circumstances [8]. Recently, fractional calculus has been applied to describe different mathematical models, and it has been shown to be more accurate in some cases compared to the classical models [9]. The main objective of this paper is to propose and analyse a deterministic and stochastic fractional order LDM model with mating disruption and sex pheromone traps taking into consideration the effect of the natural enemy on LDM.

The paper is arranged as follows: In Section 2, the deterministic mathematical model is described as well as the existence, uniqueness, non-negativity, and boundedness of the solutions of LDM system are verified. The local and global stability of equilibrium points of the LDM system is analyzed in Section 3. Using Sotomayor’s theorem, the local bifurcation conditions are derived in Section 4. In Section 5, we extend the deterministic fractional-order LDM model to the stochastic case and a numerical approximation method developed for the proposed stochastic fractional-order case. Some numerical examples are given in Section 6 to illustrate the theoretical findings. Finally, the discussion and conclusion are given in Section 7.

2  Dynamic of the Deterministic Fractional Order Model

Following [10–14], the model of lesser date moth with mating disruption and sex pheromone

trap is describing by the following system

cαDαx=rx(1−xk)−α1x−βxza+x−μ1x,cαDαy=ϵα1x−α2y+δf−μ2y,cαDαf=α2y−δf−μ3f,cαDαm=(1−ϵ)α1x−γpmy+p−μ4m,cαDαz=eβxza+x−μ5z, (1)

where CDα is the Caputo fractional derivative of order α and 0<α<1 . The detailed explanation of the variables and parameters for system (1) are listed in Tab. 1.

2.1 Existence and Uniqueness

In this section, we investigate the existence and uniqueness of the solutions of the fractional order system (1) in the region Ω×(0,T] where

Theorem 1. For each X0=(x0,y0,f0,m0,z0)∈Ω , there exists a unique solution X(t)∈Ω of the fractional order system (1), which is defined for all t≥0 .

Proof. Define a mapping F(X)=(F1(X),F2(X),F3(X),F4(X),F5(X)) , in which

F1(X)=rx(1−xk)−α1x−βxza+x−μ1x,F2(X)=ϵα1x−α2y+δf−μ2y,F3(X)=α2y−δf−μ3f,F4(X)=(1−ϵ)α1x−γpmy+p−μ4m,F5(X)=eβxza+x−μ5z. (2)

For any X,X¯∈Ω, it follows from (2) that

where

Hence, F(X) satisfies the Lipschitz condition with respect to X . According to Cresson et al. [15], as F(X) locally Lipschitz. Then there exists unique local solution to the fractional order system (1).

2.2 Non-Negativity and Boundedness

The following results show the non-negativity of the solutions of the fractional order system (1). According to Cresson et al. [15], a model of the form dXdt=F(X) satisfies the positivity property if and only if for all i=1,….,5 , Fi(X)≥0 for all X∈R+5 such that Xi=0 . Thus, the solution of the integer-order model (1), with nonnegative initial conditions remains nonnegative. Also, the solution satisfies the Lipschitz condition, as stated in Theorem 1. By Theorem 5 and Theorem 6 in Cresson et al. [15], the solution of the fractional-order model (1) also satisfies the non-negativity. The boundedness of the solutions of model (1) are given in the following theorem.

Theorem 2. All the solutions of the fractional-order LDM Model (1) starting in R+5 are uniformly bounded.

Proof. The approach of [16,17] is utilized. Let (x(t),y(t),f(t),m(t),z(t)) to be any solution of the system (1) with non-negative initial conditions. Let M(t)=x(t)+y(t)+f(t)+m(t)+z(t), then

where, ν<min{μ1,μ2,μ3,μ4,μ5} , thus, cαDαM(t)+νM≤rk4 . In accordance with Lemma 9 in Choi et al. [18], it follows that, 0≤M(t)≤M(0)Eα(−νtα)+rk4νtαEα,α+1(−νtα),

where Eα is the Mittag-Leffler function. According to Lemma 5 and Corollary 6 in Choi et al. [18], it follows

Hence all the solutions of fractional-order LDM model (1) that start in R+5 are uniformly bounded in the region, H={(x,y,f,m,z)∈R+5:M(t)≤rk4ν+ξ,for\ anyξ>0} .

3  Equilibria and Stability

The LDM model (1) has the following three equilibrium points:

1) E0=(0,0,0,0,0) , which always exists.

2) The free natural enemy equilibrium point E1=(x1,y1,f1,m1,0) , where

The free natural enemy equilibrium point exists positively if R0>1, where, R0=rα1+μ1 is the basic offspring number obtained by using the next generation method [19].

3) The coexistence equilibrium point E2=(x2,y2,f2,m2,z2) , where

4) The coexistence equilibrium point E2 exists if eβ>μ5 and R0>1+raμ5k(α1+μ1)(eβ−μ5).

The locally and globally asymptotically stable of equilibrium points of LDM system (1) are now investigated. The stability analysis of the equilibrium point E0=(0,0,0,0,0) is not considered because in this case all the population will go to extinction.

The stability of free natural enemy equilibrium point E1=(x1,y1,f1,m1,0) is investigated as follows.

Theorem 3. If R0<1+raμ5k(α1+μ1)(eβ−μ5), then the equilibrium point E1 is locally asymptotically stable.

Proof. The first three eigenvalues of J(E1) are λ1=r−2rx1k−(α1+μ1)=(α1+μ1)(1−R0) , λ2=−(γpy1+p+μ4) and λ3=eβx1a+x1−μ5 . The other roots are determined by

λ2+(α2+δ+μ2+μ3)λ+δμ2+μ3(α2+μ2)=0. (3)

The roots of Eq. (3) have negative real parts. It can be observed that λ3<0, when eβx1a+x1<μ5 which equivalent R0<1+raμ5k(α1+μ1)(eβ−μ5). So, E1 is locally asymptotically stable if 1<R0<1+raμ5k(α1+μ1)(eβ−μ5).

Theorem 4. If R0<1+raμ5k(α1+μ1)(eβ−μ5) , then the natural enemy extinction equilibrium point E1 is globally asymptotically stable.

Proof. The following positive definite Lyapunov function is considered.

By calculating the time derivative of L1 along the solution of system (1), one obtains,

Then cαDαL1≤0 , when eβx1a+x1<μ5 which equivalent R0<1+raμ5k(α1+μ1)(eβ−μ5). According to generalized Lyapunov–Lasalle’s invariance principle [20], the natural enemy extinction equilibrium point E1 is globally asymptotically stable when R0<1+raμ5k(α1+μ1)(eβ−μ5).

The stability of the coexistence equilibrium point E2=(x2,y2,f2,m2,z2) is investigated as follows.

The first eigenvalue of J(E2) is λ1=−γpp+y2−μ4 . The other eigenvalues are determined by

λ4+θ1λ3+θ2λ2+θ3λ+θ4=0, (4)

where

where B11=rx2k−βz2x2(a+x2)2 and B44=γpp+y2+μ4 . Then, the proposition proposed in Matouk [21] can be used to determine the stability conditions of the equilibrium point E2 . When B11>0 , then θi>0,i=1,2,3,4. Also, θ1θ2θ3−θ32−θ12θ4=(B11(a+x2)3Ω4+Ω2Ω3)(B11δμ2(a+x2)3+Ω1(α2+μ2)+eaβ2w2x2(δ+μ3))2(a+x2)9, where

when Ω2Ω3>0 , then θ1θ2θ3−θ32−θ12θ4>0 , therefore all the eigenvalues of the Jacobian matrix J(E2) near the equilibrium point E2 have negative real parts. Thus, due to the Routh-Hurwitz criterion the equilibrium point E2 is locally asymptotically stable. The local stability of the coexistence equilibrium point E2 is given in the following theorem.

Theorem 5. If Ω2Ω3>0 and B11>0 , then the coexistence equilibrium point E2 is locally asymptotically stable.

The global stability of the coexistence equilibrium point is investigated in the following theorem.

Theorem 6. If βz2a(a+x2)<rk , then the coexistence equilibrium point E2 is globally asymptotically stable.

Proof. The following positive definite Lyapunov function is considered.

By calculating the time derivative of L2 along the solution of system (1), one obtains,

Choosing H=eaa+x2, then cαDαL2<0 . According to generalized Lyapunov–Lasalle’s invariance principle [20], the coexistence equilibrium point E2 is globally asymptotically stable when βz2a(a+x2)<rk .

4  Bifurcation Analysis

In this section the local bifurcations near the equilibrium points of LDM model (1) are inves-tigated with the help of Sotomayor’s theorem [22]. The Hopf bifurcation theorem given in Liu [23] is also presented to discuss the bifurcation analysis of the underlying system. One can compute

D2F(x,φ)(U,U)=((2aβz(a+x)3−2rk)ζ12−2aβζ1ζ5(a+x)200−2γpmζ22(p+y)3+2γpζ2ζ4(p+y)2eaβζ1ζ5(a+x)2−2eaβzζ12(a+x)3), (5)

where φ is any bifurcation parameter and U=(ζ1,ζ2,ζ3,ζ4,ζ5)T is any eigenvector.

Theorem 7. The LDM system (1) undergoes a transcritical bifurcation with respect to the bifurcation parameter r around E0 if r=r∗=α1+μ1 .

Proof. Let V1=(ν1,ν2,ν3,ν4,ν5)T be the eigenvector corresponding to the zero eigenvalue of the matrix (E0) r∗=α1+μ1 , hence J(E0)V1=0 , gives

where ν1 is any non zero real number. Similarly, suppose V2=(τ1,τ2,τ3,τ4,τ5)T be the eigenvector corresponding to the zero eigenvalue of the matrices J(E0) , thus J(E0)TV2=0 gives V2=(1,0,0,0,0)T . Consider, ∂F∂r=Fr(X,r)=(x(1−xk),0,0,0,0)T , thus, V2TFr(E0,r∗)=0. Therefore, according to Sotomayor’s theorem for local bifurcation, the LDM model (1) has no saddle-node bifurcation near E0 at r∗=α1+μ1 . One can note that r=α1+μ1 is equivalent to R0=1 . Now,

then V2TDFr(E0,r∗)V1=ν1≠0. Using (5), one obtains

Thus, according to Sotomayor’s theorem, the LDM system (1) has a transcritical bifurcation at r∗=α1+μ1 as the parameter r passes through the value r∗ , thus the proof is complete.

Theorem 8. The LDM system (1) undergoes a transcritical bifurcation with respect to the bifurcation parameter μ5 around E1=(x1,y1,f1,m1,0) if μ5=μ5∗=eβx1a+x1 .

Proof. The Jacobian matrix of the LDM system (1) at the free enemy equilibrium point E1 with μ5=μ5∗ has zero eigenvalue takes the form

where A11=r−2rx1k−(α1+μ1),A22=(α2+μ2),A33=(δ+μ3) and A44=γpy1+p+μ4. J(E1)V3=0 , gives the eigenvector corresponding to the zero eigenvalue of the matrix J(E1) , hence

where ν5 is any non zero real number. Similarly, J(E1)TV4=0 gives the eigenvector corresponding to the zero eigenvalue of the matrix J(E1)TV4=0 , hence V4=(0,0,0,0,1)T. Consider, ∂F∂μ5=Fμ5(X,μ5)=(0,0,0,0,−z)T , thus, V4TFμ5(E1,μ5∗)=0. Therefore, according to Sotomayor’s theorem for local bifurcation, the LDM model (1) has no saddle-node bifurcation near E1 at μ5=μ5∗ . Now,

then, V4TDFμ5(E1,μ5∗)V3=−ν5≠0. Using (5), one obtains

Thus, according to Sotomayor’s theorem, the LDM system (1) has a transcritical bifurcation at μ5∗=eβx1a+x1 as the parameter μ5 passes through the value μ5∗ , thus the proof is complete.

In this part, we shall show that as the coexistence equilibrium loses stability, periodic solutions can bifurcate from the positive equilibrium. We first give the following lemma.

Lemma 9 . The characteristic Eq.(4) has a pair of purely imaginary roots and the remaining roots have negative real parts if and only if z2=r(a+x2)2kβ and θ1θ2θ3−θ32−θ12θ4=0.

Suppose (4) has two eigenvalues which have negative real parts and two complex conjugates eigenvalues (call them λ=m(φ)±in(φ) ) such that m(φ∗)=0,n(φ∗)>0, dmdφ|φ=φ∗≠0. Substituting λ=m(φ)±in(φ) into (4), and separating the real and imaginary, we get

m4+θ1m3+θ2m3+θ3m+θ4−(6m2+3θ1m+θ2)n2+n4=0, (6)

4m3+3θ1m2+2θ2m+θ3−(4m+θ1)n2=0, (7)

Substituting (6) into (7), differentiating with respect to φ and utilizing m(φ∗)=0 and n(φ∗)≠0, we have

Theorem 10. For the coexistence equilibrium point E2 of the integer order LDM system (1), the system around E2 enters into the Hopf bifurcation when φ passes φ∗ if the coefficients θj(φ)(j=1,2,3,4) at φ=φ∗ satisfying the following condition:

1.     Φ(φ∗)=[θ1(φ)θ2(φ)θ3(φ)−θ32(φ)−θ12(φ)θ4(φ)]|φ=φ∗=0,

2.     (4θ4−θ1θ3−θ22)|φ=φ∗≠0 ,

3.     dΦ(φ)dφ|φ=φ∗≠0 .

According to Theorem 10, there exists a Hopf bifurcation in the LDM model (1) with α=1 , where the Hopf bifurcation is controlled by φ .

5  Stochastic Fractional Order Model

This section extends the deterministic fractional-order LDM model (1) to the following stochastic fractional-order model.

CαDαx=rx(1−xk)−α1x−βxza+x−μ1x+σ1xdW1dt,CαDαy=ϵα1x−α2y+δf−μ2y+σ2ydW2dt,CαDαf=α2y−δf−μ3f+σ3fdW3dt,CαDαm=(1−ϵ)α1x−γpmy+p−μ4m+σ4mdW4dt,CαDαz=eβxza+x−μ5z+σ5zdW5dt, (8)

where Wi(i=1,2,3,4,5) are independent standard Brownian motions with Wi(0)=0 and σi>0 denote the intensities of the white noise. The stochastic fractional-order LDM model (8) can be written in the general form:

CαDαX(t)=F(X)+g(X)dWdt, (9)

where F(x) is given in (2), g(x)=(σ1x,σ2y,σ3f,σ4m,σ5z) and dWdt=(dW1dt,dW2dt,dW3dt,dW4dt,dW5dt)T. Applying Riemann–Liouville integral to both sides of (9), one can obtain the following stochastic Volterra integral equation.

X(t)=X0+∫0tF(X)(t−s)α−1Γ(α)ds+∫0tg(X)(t−s)α−1dW(s)Γ(α)ds. (10)

According to Wang et al. [24,25], under some conditions on the coefficient functions, the global existence and uniqueness of solutions for the stochastic fractional-order system (8) can be investigated. Because Grunwald-Letnikov’s definition is the most straightforward from the point of view of numerical implementation, so we will use it to solve the LDM system of fractional order stochastic differential equations. Grunwald-Letnikov ( GLDα ) fractional derivative of order α defined by Aminikhah et al. [26,27]

GLαDαf(t)=Limh→0⁡h−α∑j=0[t−ah](−1)j(αj)f(t−jh), (11)

where [t−ah] means the integer part of t−ah . This formula (11) can be reduced to

GLαDαf(tn)≈h−α∑j=0nwjαf(tn−j), (12)

where h is the time step, tn=nh and wjα are the Grunwald-Letnikov coefficients satisfy the following recurrence relationship

If f(t) is continuous function and f′(t) is integrable function in the interval [0,T] , then the relation between Caputo and Grunwald-Letnikov fractional derivative takes the form [28,29]

CαDαf(t)=GLαDαf(t)−f(0)t−αΓ(1−α)=GLαDαf(t)−Limn→∞(−1)nhα(α−1n)f(0)≈1hα∑j=0nwjα(f(tn−j)−f(0)). (13)

Now, the fractional order stochastic LDM model (8) can be written as

xn=x0+hα(rxn−1(1−xn−1k)−(α1+μ1)xn−1−βxn−1zn−1a+xn−1+σ1xn−1hζ1n)−∑j=1nwjα(xn−j−x0)yn=y0+hα(ϵα1xn−(α2+μ2)yn−1+δfn−1+σ2yn−1hζ2n)−∑j=1nwjα(yn−j−y0)fn=f0+hα(α2yn−(δ+μ3)fn−1+σ3fn−1hζ3n)−∑j=1nwjα(fn−j−f0)mn=m0+hα((1−ϵ)α1xn−γpmn−1yn+p−μ4mn−1+σ4mn−1hζ4n)−∑j=1nwjα(mn−j−m0)zn=z0+hα(eβxnzn−1a+xn−μ5zn−1+σ5zn−1hζ5n)−∑j=1nwjα(zn−j−z0), (14)

where, σi and ζin represent real constants and a 5D Gaussian white noise processes, respectively, i=1,2,3,4,5. ζi satisfy the follows:

δij is Kronecker delta and δ(t1−tj) is the Dirac delta function.