# Stability and Bifurcation Analysis of a Discrete Singular Bioeconomic System

Authors: Qamar Din; A. M. Yousef; A. A. Elsadany
Journal: Discrete Dynamics in Nature and Society
 (2021)

---

## Abstract

The main concern of this paper is to discuss stability and bifurcation analysis for a class of discrete predator-prey interaction with Holling type II functional response and harvesting effort. Firstly, we establish a discrete singular bioeconomic system, which is based on the discretization of a system of differential algebraic equations. It is shown that the discretized system exhibits much richer dynamical behaviors than its corresponding continuous counterpart. Our investigation reveals that, in the discretized system, two types of bifurcations (i.e., period-doubling and Neimark–Sacker bifurcations) can be studied; however, the dynamics of the continuous model includes only Hopf bifurcation. Moreover, the state delayed feedback control method is implemented for controlling the chaotic behavior of the bioeconomic model. Numerical simulations are presented to illustrate the theoretical analysis. The maximal Lyapunov exponents (MLE) are computed numerically to ensure further dynamical behaviors and complexity of the model.


---

## Body

## 1. Introduction

Bioeconomics is linked closely to the early development of ideas in fisheries economics due to the pioneering work of Canadian economists Gordon [1] and Anthony Scott (in 1955). Their basic theories used recent developments in modeling of biological fisheries, initially the contributions made by Schaefer in 1954 and 1957 on introducing a systematic connection between fishing mechanism and growth of biological type through the implementation of mathematical modeling verified by experimental studies, and also associated itself to resource protection, ecology, and the environment [2]. These concepts were developed from the multifishing science environment in Canada at that time. Modeling and fishing science developed rapidly during an innovative and productive period, especially among Canadian fishing researchers of various disciplines. Fishing mortality and population modeling were launched for economists, and novel interdisciplinary methods of modeling became obtainable for the economists, which made it feasible to measure the economical and biological impacts of various fisheries management decisions and fishing activities. Modern bioeconomics related to fisheries science can furnish perception into developing methods to deal with the overexploitation and complexities of overcapacity in marine fisheries where most are affected by lack of solid governance, changing coastal ecosystem dynamics, and natural fluctuations [3].Moreover, Gordon in [1] suggested economic theory keeping in view the common property of resource, which was based on the effect of the harvest effort on an ecosystem by taking into account an economic perspective, assuming that xt and et denote the density of harvested population and the harvest effort in an ecosystem, respectively; then the total cost is equal to cet, and the total revenue is equal to petxt, where c denotes the cost of harvest effort, and p is used for the unit price of harvested population. Then, the economic interest μ for the harvest effort by the harvested population is given by(1)μ=etpxt−c.Taking into account Gordon [1] theory, Zhang et al. [4] studied a class of bioeconomic system with implementation of theory for singular systems. Their study was consisted of bifurcation analysis and chaos control for the proposed bioeconomic model. Later on, Liu et al. [5] reported stability, bifurcation analysis, and state feedback control for a class of predator-prey interaction with harvest effort on predator and stage structure for prey. Chakraborty et al. [6] studied a bioeconomic system with implementation of theory of differential algebraic equations. They investigated stability, Hopf bifurcation, and state feedback control for a class of predator-prey interaction with time-delayed effect. Zhang et al. [7] investigated a singular bioeconomic model for prey-predator interaction with diffusion and time delay. Zhang et al. [8] carried out comprehensive study related to theory, applications, complexity, and control of singular bioeconomic systems. Zhang et al. [9] explored the Hopf bifurcation for a predator-prey type bioeconomic system with two delays and predator harvesting. Meng and Zhang [10] discussed the qualitative behavior of a delayed singular bioeconomic predator-prey model without and with stochastic fluctuation. Liu et al. [11] analyzed the local dynamics and Hopf bifurcation for a biological economic model with Holling type II functional response and harvesting effort on prey. Liu et al. [12] proposed a singular fishery model for a class of prey-predator interaction with gestation delay for predator and maturation delay for prey. Liu et al. [13] formulated and discussed a singular predator-prey model by implementing commercial harvesting on predator with gestation delay for predator and maturation delay for prey. In [14], Li et al. studied a singular bioeconomic predator-prey system with Holling type II functional response and nonlinear harvesting on prey. Meng and Wu [15] discussed a singular prey-predator system with two delays, nonlinear predator harvesting and Beddington–DeAngelis functional response. Babaei and Shafiee [16] reported stability analysis, bifurcation, chaotic behavior, and control for a singular bioeconomic model of prey-predator interaction governed by an algebraic equation and 3-dimensional differential equations.In case of mathematical modeling of predator-prey interaction, the research concerning interspecific interactions has been numerously based on continuous predator-prey systems of two variables. On the other hand, particular species, covering several classes of insects and seasonal plants, have nonoverlapping generations successively, and consequently, their population undergoes in discrete time steps. Populations with nonoverlapping generations can be modeled suitably with difference equations, otherwise known as iterative maps or discrete dynamical systems. Several authors have shown that nonoverlapping generations governed by iterative maps reveal complex and chaotic behavior, and the dynamics in such cases may yield a much richer set of patterns than those examined in continuous-time systems (cf. [17–24]). Furthermore, discrete-time models also have been used for rich dynamics of bioeconomic systems for some classes of predator-prey interaction. For example, in [25], the authors analyzed complex dynamics of a discrete-time bioeconomic system for predator-prey interaction with the implementation of Euler approximations. Wu and Chen [26] implemented the Poincaré scheme for discretization of a singular bioeconomic model and they analyzed period-doubling bifurcation, Neimark–Sacker bifurcation, and stability behavior. Liu et al. [27] studied the chaotic behavior of a discrete singular system related to the bioeconomic model of the prey-predator type.Taking into account predator interaction with logistic growth and Holling type II functional response for prey population, we have the following system [28]:(2)dxdt=xd−kx−ya+x,dydt=ybxa+x−r,where x=xt and y=yt denote state variables for the densities of prey and predator at time t, respectively. Moreover, d is the intrinsic growth rate of prey, r represents the natural death rate of the predator, a is used for half capturing saturation constant, and b represents the maximal growth rate of the predator. Furthermore, d/k denotes the environmental carrying capacity for the prey population.Keeping in view (1) and (2), we obtain the following predator-prey biological economic model with Holling type II functional response with harvest effort:(3)dxdt=xd−kx−ya+x−e,dydt=ybxa+x−r,μ=epx−c.Applying the forward Euler scheme to system (3), we obtain the discrete-time predator-prey biological economic model with Holling type II functional response as follows:(4)xn+1=xn+hxnd−kxn−yna+xn−en,yn+1=yn+hynbxna+xn−r,μ=enpxn−c,where h is the integral step size for the Euler approximation. In this paper, we discuss some dynamical aspects of the discrete singular model (4). For this, the first existence of interior (positive) fixed point and local dynamics of system (4) about biologically feasible equilibrium are carried out. Secondly, it is proved that system (4) undergoes period-doubling bifurcation and Neimark–Sacker bifurcation by varying the economic profit μ as the bifurcation parameter. Thirdly, a state delayed feedback control strategy is applied to avoid bifurcating and chaotic behavior of bioeconomic model (4). At the end, numerical examples are presented for verification and illustration of our theoretical discussion.


## 2. Fixed Points and Stability Analysis

In order to study the qualitative behavior of the solutions of the nonlinear model (4), we study the existence of fixed points and their stability properties. From system (4), we can see that there exists a fixed point X0≔x0,y0,e0 in R3+ if and only if X0 is a solution of the following equations:(5)x=x+hxd−kx−ya+x−e,y=y+hy−r+bxa+x,μ=epx−c.Through a simple calculation, we obtain(6)X0=x0,y0,e0=arb−r,a+x0d−kx0−e0,μpx0−c.For biological considerations, we focus on the dynamics of the positive fixed point of the system (4). Thus, throughout the paper, we assume the conditions for the existence of a unique positive fixed point of system (4) as follows:(7)b>r,bc<rap+c,bd−rak+db−r2−μrap+c−bc>0.Inμk-plane, the existence of a unique positive fixed point of system (4) is depicted in Figure 1.Figure 1
Coexistence region for system (4) with a=6.76, b=9.4, c=1.58, r=0.74, p=9.2, and d=8.83.The generalized Jacobian matrixJx0,y0,e0 of system (4) about interior (positive) fixed point x0,y0,e0 is computed as follows:(8)Jx0,y0,e0=1−hkx0+hx0y0a+x02−hx0a+x0−hx0abhy0a+x0210e0p0px0−c.Then, it is easy to see that the generalized characteristic equation of the Jacobian matrixJx0,y0,e0 can be written as(9)det1−hkx0+hx0y0a+x02−λ−hx0a+x0−hx0abhy0a+x021−λ0e0p0px0−c=0,which on simplification yields (10)λ2+Pλ+Q=0,where(11)P=hkx0−hx0y0a+x02−he0px0px0−c−2=1arp/b−r−c2R1μ−S1+2arpb−r−c2=H−Q−1,Q=1−hkx0+hx0y0a+x02+he0px0px0−c+abh2x0y0a+x03=1arp/b−r−c2R2μ−S2−1arpb−r−c2,where(12)H=abx0y0a+x03>0,R1=−hrbap+c,S1=−hrbd−akb+rb−r,R2=hrbap+c1−hr+hcb,S2=hrbd−akb+rb−r+hdb−r−hark.LetFλ=λ2+Pλ+Q, then(13)F1=1+P+Q=H>0,F−1=1−P+Q=2arp/b−r−c2R3μ−S3−2arpb−r−c2,where(14)R3=hrbap+c1−hr2+hcb2,S3=hrbd−akb+rb−r+hdb−r2−hark2.In order to discuss the stability of the fixed point ofx0,y0,e0, we need the following lemma.Lemma 1 (see [23]).
ConsiderSζ=ζ2+Pζ+Q. Moreover, assuming that S1>0 and ζ1 and ζ2 are two roots of Sζ=0, then the following hold true:(i)
ζ1<1andζ2<1if and only ifS−1>0, andS0=Q<1(ii)
ζ1<1andζ2>1(orζ1>1andζ2<1) if and only ifS−1<0(iii)
ζ1>1andζ2>1if and only ifS−1>0andS0>1(iv)
ζ1=−1andζ2≠1if and only ifS−1=0andP≠0,2(v)
ζ1andζ2are complex andζ1=1andζ2=1if and only ifP2−4Q<0andQ=1Assume thatζ1 and ζ2 are roots for the characteristic equation of the variational matrix Jx0,y0,e0 about interior (positive) fixed point x0,y0,e0 which are known as eigenvalues for the equilibrium point x0,y0,e0. Taking into account the topological types related to the fixed point x0,y0,e0 of system (4), we say that the fixed point x0,y0,e0 is a sink (asymptotically stable) if ζ1<1 and ζ2<1; x0,y0,e0 is called a source (repeller) if ζ1>1 and ζ2>1; x0,y0,e0 is called a saddle if ζ1>1 and ζ2<1 (or ζ1<1 and ζ2>1; and x0,y0,e0 is called nonhyperbolic if either ζ1=1 or ζ2=1, on the other hand, if ζ1=−1 and ζ2≠1 are necessary conditions for the emergence of period-doubling bifurcation and −2<P<2 with Q=1 are necessary conditions for the occurrence of Neimark–Sacker bifurcation. Moreover, if S1=1+P+Q>0, then all cases of Lemma 1 are depicted in Figure 2 in PQ-plane.Figure 2
Dynamical classification of planar system withS1>0.Keeping in view Lemma1, the following theorem is presented for local dynamics of system (4) about its positive fixed point.Theorem 1.
Assume thatb>r, bc<rap+c, and bd−rak+d/b−r2−μ/rap+c−bc>0; then, there exists unique interior (positive) fixed point x0,y0,e0 for system (4) satisfying the following conditions:(i)
x0,y0,e0is a sink if and only ifR2μ<S2arp/b−r−c2andR3μ>S3−2arp/b−r−c2(ii)
x0,y0,e0is a source if and onlyR2μ>S2arp/b−r−c2andR3μ>S3−2arp/b−r−c2(iii)
x0,y0,e0is a saddle if and only ifR3μ<S3−2arp/b−r−c2(iv)
x0,y0,e0is nonhyperbolic if one of the following conditions hold true:(1)
R3μ=S3−2arp/b−r−c2,R1μ≠S1+2arp/b−r−c2,S1+4arp/b−r−c2(2)
R2μ=S2arp/b−r−c2andHμ<4Next, fora=3.9, b=14, c=3.7, r=3.4, p=5.8, d=21.1, and h=0.25, the dynamical classification of interior equilibrium of system (4) is depicted in μk-plane (see Figure 3).Figure 3
Dynamical classification of interior equilibrium of system (4) in μk-plane.Taking into account part (iv.1) of Theorem1, it is easily to observe that the eigenvalues about the equilibrium point x0,y0,e0 are given by λ1=−1 and λ2=1−P with λ2≠−1. On the other hand, if part (iv.2) of Theorem 1 holds true, then one can obtain that the eigenvalues for the equilibrium point x0,y0,e0 are conjugate complex numbers with modulus one.Next, we consider the following set:(15)FB=0<μ=1R3S3−2arpb−r−c2<d−karb−rparb−r−ca,b,c,d,k,p,r,μ,h:R1μ≠S1+2arpb−r−c2,S1+4arpb−r−c2.It is easy to see that the steady statex0,y0,e0 undergoes period-doubling (flip) bifurcation whenever parameters vary in a small neighborhood of FB.Next, we consider the following curve:(16)HB=a,b,c,d,k,p,r,μ,h:0<μ=S2R2arpb−r−c2<d−karb−rparb−r−c,Hμ<4.On the other hand,x0,y0,e0 undergoes the Neimark–Sacker (Hopf) bifurcation whenever parameters vary in a small neighborhood of HB. In Section 3, we discuss the occurrence of period-doubling bifurcation around the interior fixed point x0,y0,e0 with the variation of parameters in a small neighborhood of FB, and the emergence of the Neimark–Sacker bifurcation about x0,y0,e0 with varying the parameters in a small neighborhood of HB.


## 3. Bifurcation Analysis

Keeping in view the analysis of Section3, we discuss the period-doubling bifurcation and Neimark–Sacker bifurcation about the positive fixed point x0,y0,e0 in this section. For this, choose parameter μ as bifurcation parameter for investigating the period-doubling bifurcation and Neimark–Sacker bifurcation about x0,y0,e0 by implementing the novel normal form theory of discrete singular systems, the center manifold theorem, and the bifurcation theory of discrete systems [29–35].


### 3.1. Period-Doubling Bifurcation

First, we start our investigation related to period-doubling bifurcation for system (4) about its equilibrium x0,y0,e0 with a variation of parameters in a small neighborhood of FB. For this, we choose parameters a,b,c,d,k,p,r,μ1,h arbitrarily from FB, taking into account system (4) with a,b,c,d,k,p,r,μ1,h∈FB. In this case, system (4) is described by the following 2-dimensional map:(17)x⟶x+hxd−kx−ya+x−e,y⟶y+hy−r+bxa+x,μ1=epx−c.It is easy to see that system (17) has a unique positive fixed point x0,y0,e0 such that the eigenvalues are given by λ1=−1 and λ2=1−P with λ2≠−1. Assume that a,b,c,d,k,p,r,μ1,h∈FB with μ1=1/R1S1−2arp/b−r−c2. Next, we take μ∗ as a bifurcation parameter, considering a perturbation for system (17) as follows:(18)x⟶x+hxd−kx−ya+x−e,y⟶y+hy−r+bxa+x,μ1+μ∗=epx−c,where μ∗≪1 is taken as a small perturbation parameter.Then, system (17) can be described in the following way:(19)x,yT⟶fx,y,e,gx,y,e=0,where(20)fx,y,e≔f1x,y,e,f2x,y,eT,f1x,y,e=x+hxd−kx−ya+x−e,f2x,y,e=y+hy−r+bxa+x,gx,y,e=epx−c−μ1.Then, it is easy to see thatDgx0,y0,e=e0p,0,px0−c such that the rank of Dgx0,y0,e0=1. On the other hand, a local parameterization ψ for the 2-dimensional smooth manifold defined by Mg=x,y,e∈R3;gx,y,e=0,and rankDgx,y,e=1 about x,y,e∈Bx0,y0,e0⊂Mg is given as follows: for any x,y,e∈Bx0,y0,e0, there is Z∈A⊂R2 such that(21)X=ψZ=X0+U0Z+V0HZ,gψZ=0,where(22)U0=100100,V0=001,X=x,y,e, X0=x0,y0,e0, Z=z1,z2, and H:R2⟶R is a smooth mapping. For further details, the interested reader is referred to [36]. Taking into account the definition of ψ, we obtain the following:(23)Dψz=DgXU0T−1001001,for arbitrarily chosen x,y,e∈Bx0,yo,e0. Then, system (19) is written as follows:(24)Z⟶fZ,Z∈A⊂R2,where A=ψ−1Bx0,yo,e0.Taking into account (24), we have(25)z1⟶a1z1+a2z2+a13z1μ∗+a11z12+a12z1z2+a113z12μ∗+a111z13+a112z12z2+Oz1+z2+μ∗4,z2⟶b1z1+b2z2+b13z1μ∗+b11z12+b12z1z2+b113z12μ∗+b111z13+b112z12z2+Oz1+z2+μ∗4.In order to compute the coefficients related the normal form, we need some further computation. From the aforementioned computation, one can easily get(26)Df1X=1+hd−2kx−aya+x2−e,−hxa+x,−hxDf2X=abhya+x2,1+h−r+bxa+x,0DgX=ep,0,px−cDψz=DgXU0T−1001001=ep0px−c100010−1001001=1001−eppx−c0=Dz1ψz,Dz2ψz,and from this, it follows that(27)Dz1ψz=1,0,−eppx−cT,Dz2ψz=0,1,0T.Assume thatfiziψZ represents the derivative of fiψZ with respect to zi and taking fiziX=fiziψZ. In a similar way, one can adopt notations for fizizjψZ and fizizjzkψZ. Next, it is easy to see that(28)f1z1X=Df1XDz1ψz=1+hd−2kx−aya+x2−e+epxpx−c,f1z2X=Df1XDz2ψz=−hxa+x,f2z1X=Df2XDz1ψz=abhya+x2,f2z2X=Df2XDz2ψz=1+h−r+bxa+x.PuttingX0 into (28), one has the following:(29)f1z1X0=1+h−kx0+e0px0px0−c+x0d−kx0−e0a+x0=1+hrbd+akb+rb−r+hrbb−r2ap+capr−cb+cr2μ1+μ∗,f1z2X0=−x0a+x0=−hrb,f2z1X0=abhy0a+x02=abhd−kx0−e0a+x0=hb−rd−karb−r−μ1+μ∗b−rapr−cb+cr,f2z2X0=1+h−r+bx0a+x0=1.Then, from (28), it is easy to see that(30)Df1z1X=h−2k+2aya+x3−epcpx−c2,−aha+x2,h−1+pxpx−c,Df1z2X=−aha+x2,0,0,Df2z1X=−2abhya+x3,abha+x2,0,Df2z2X=abha+x2,0,0.Consequently, one has the following:(31)f1z1z1X=Df1z1XDz1ψz=h−2k+2aya+x3−2epcpx−c2,f1z1z2X=Df1z1XDz2ψz=−aha+x2,f2z1z1X=Df2z1XDz1ψz=−2abhya+x3,f2z1z2X=Df2z1XDz2ψz=abha+x2.Putting the value ofX0 in (31), one has the following:(32)f1z1z1X0=h−2k+2ad−kx0−e0a+x02−2e0pcpx0−c2=h−2k+2db−r2ab2−2krb−rb2−2b−r3apr−cb+cr1ab2+pcapr−cb+cr2μ1+μ∗,f1z1z2X0=−aha+x02=−hb−r2ab2,f2z1z1X0=−2abhd−kx0−e0a+x02=2hb−r2abakrb−r−d+2hb−r3abapr−cb+crμ1+μ∗,f2z1z2X0=abha+x02=hb−r2ab.Then, from (31), it follows that(33)Df1z1z1X=−6ahya+x4+4hep2cpx−c3,2aha+x3,−2hpcpx−c2,Df1z1z2X=2aha+x3,0,0,Df2z1z1X=6abhya+x4,−2abha+x3,0,Df2z1z2X=−2abha+x3,0,0.Then, we get(34)f1z1z1z1X=Df1z1z1XDz1ψz=−6ahya+x4+6hep2cpx−c3,f1z1z1z2X=Df1z1z1XDz2ψz=2aha+x3,f2z1z1z1X=Df2z1z1XDz1ψz=6abhya+x4,f2z1z1z2X=Df2z1z1XDz2ψz=−2abha+x3.Putting the value ofX0 into (34), we obtain that(35)f1z1z1z1X0=−6ahd−kx0−e0a+x03+6he0p2cpx0−c3=6hkrb−r2ab3−6hdb−r3a2b3+6hb−r4apr−cb+cr1a2b3+p2capr−cb+cr3μ1+μ∗,f1z1z1z2X0=2aha+x03=2hb−r3a2b3,f2z1z1z1X0=6abhd−kx0−e0a+x03=6hb−r3a2b2d−akrb−r−6hb−r4a2b2apr−cb+crμ1+μ∗,f2z1z1z2X0=−2abha+x03=−2hb−r3a2b2.Consequently, it follows from (29), (32), and (35) that(36)a1=1+hrbd+akb+rb−r+hrbb−r2ap+capr−cb+cr2μ1,a2=−hrb,a13=hrbb−r2ap+capr−cb+cr2,a11=h−2k+2db−r2ab2−2krb−rb2−2b−r3apr−cb+cr1ab2+pcapr−cb+cr2μ1,a12=−hb−r2ab2,a113=−2hb−r3apr−cb+cr1ab2+pcapr−cb+cr2a111=6hkrb−r2ab3−6hdb−r3a2b3+6hb−r4apr−cb+cr1a2b3+p2capr−cb+cr3μ1,a112=2hb−r3a2b3,b1=hb−rd−b−rapr−cb+crμ1−akr,b2=1,b13=−hb−r2apr−cb+cr,b11=2hb−r2abakrb−r−d+2hb−r3abapr−cb+crμ1,b12=hb−r2ab,b113=2hb−r3abapr−cb+cr,b111=6hb−r3a2b2d−akrb−r−6hb−r4a2b2apr−cb+crμ1,b112=−2hb−r3a2b2.Next, we consider the following nonsingular matrix:(37)T=a2a2−1−a1λ2−a1,implementing the following translation:(38)z1z2=Tuv.Then, it is easy to see that system (24) can be written as follows:(39)uv⟶−100λ2uv+f1u,v,μ∗f2u,v,μ∗,where(40)f1u,v,μ∗=a11λ2−a1−a2b11a21+λ2z12+a12λ2−a1−a2b12a21+λ2z1z2+a13λ2−a1−a2b13a21+λ2z1μ∗+a113λ2−a1−a2b113a21+λ2z12μ∗+a111λ2−a1−a2b111a21+λ2z13+a112λ2−a1−a2b112a21+λ2z12z2+Oz1+z2+μ∗4f2u,v,μ∗=a111+a1+a2b11a21+λ2z12+a121+a1+a2b12a21+λ2z1z2+a131+a1+a2b13a21+λ2z1μ∗+a1131+a1+a2b113a21+λ2z12μ∗+a1111+a1+a2b111a21+λ2z13+a1121+a1+a2b112a21+λ2z12z2+Oz1+z2+μ∗4z1=a2u+v,z2=−1+a1u+λ2−a1v,z1z2=a2−1+a1u2+λ2−2a1−1uv+a2λ2−a1v2,z12=a22u2+uv+v2,z13=a23u3+3u2v+3uv2+v3,z12z2=a22−1+a1u3+λ2−3a1−2u2v+2λ2−3a1−1uv2+λ2−1v3.In order to determine the center manifoldWc0,0,0 of (39) about the equilibrium point 0,0 in a small neighborhood of μ∗, we implement an application of the center manifold theorem [29], and it is easy to see that there exists a center manifold of the following form:(41)Wc0,0,0=u,v,μ∗∈R3,v=hu,μ∗,h0,0=0,Dh0,0=0,foru and μ∗ sufficiently small. Furthermore, we assume that(42)hu,μ∗=c0u2+c1uμ∗+c2μ∗2+Ou+μ∗3.Then, it is easy to see that the center manifold satisfies the following:(43)h−u+f1u,hu,μ∗,μ∗,μ∗=λ2hu,μ∗+f2u,hu,μ∗,μ∗.Taking into account (42) and (43) and then comparing coefficients for (43), one can easily obtain that(44)c0=11−λ221+a1a2a11−a121+a1−a2b12+a22b11,c1=a131+a1+a2b131+λ22,c2=0.Furthermore, taking into account (39), it is restricted to the center manifold Wc0,0,0 as follows:(45)F:u⟶−u+d11u2+d12uμ∗+d112u2μ∗+d122uμ∗2+d111u3+Ou+μ∗4,where(46)d11=1λ2+1λ2−a1a11a2−a121+a1+a113a2−b11a22+1+a1b12a2−b113a22,d12=1λ2+1a13λ2−1−b13a2,d112=1λ2+1λ2−a1a11a2c1+a12c1λ2−2a1−1+a13c0−b11a22c1−b12a2c1λ2−2a1−1−b13a2c0,d122=1λ2+1a13c1λ2−1−b13a2c1,d111=1λ2+1λ2−a1a11a2c0+a111a22+a12c0λ2−2a1−1−a112a21+a1−b11a22c0−b12a2c0λ2−2a1−1−b111a23+b112a221+a1.In order for system (45) to undergo a flip bifurcation, we require that two discriminatory quantities α1and α2 are not zero, where(47)α1=2∂2F∂μ∗∂u+∂F∂μ∗∂F∂u0,0=2d12,α2=12∂2F∂u22+13∂3F∂u30,0=2d111+d112.Keeping in view the aforementioned computation and bifurcation theory given in [36, 37], we have the following theorem.Theorem 2.
Assume thata,b,c,d,k,p,r,μ1∈FB. If α2≠0, then system (17) undergoes a flip bifurcation at the fixed point x0,y0,e0 when the parameter μ varies in a small neighborhood of μ1. Moreover, if α2>0 (resp., α2<0), then the period-2 orbits that bifurcate from x0,y0,e0 are stable (resp., unstable).In the last Section5 related to numerical simulation, we choose some parametric values for system (4) such that it undergoes period-doubling bifurcation about positive equilibrium as μ varies in the suitable interval (see Figure 4).Figure 4
Bifurcation diagrams and LLE for system (3) with a=5.8, b=8.8, c=2.1, d=3.6, k=2.5, r=1.6, p=2.2, h=0.98, μ∈0.19,0.22, and x0,y0,e0=1.28,0.6428,0.287. (a) Bifurcation diagram for xn. (b) Bifurcation diagram for yn. (c) Bifurcation diagram for en. (d) Largest Lyapunov exponents.
(a)(b)(c)(d)


### 3.2. Neimark–Sacker Bifurcation

Finally, we discuss the Neimark–Sacker bifurcation about equilibriumx0,y0,e0 with variation of parameters a,b,c,d,k,p,r,μ2,h in a small neighborhood of HB. For this, we choose parameters a,b,c,d,k,p,r,μ2,h arbitrarily from HB, taking into account system (4) with a,b,c,d,k,p,r,μ2,h∈HB. In this case, system (4) is described by the following 2-dimensional map:(48)x⟶x+xd−kx−ya+x−e,y⟶y+y−r+bxa+x,μ2=epx−c.Then, it is easy to observe that map (48) has a unique positive fixed point x0,y0,e0. Next, we take μ∗ as a bifurcation parameter and consider a perturbation corresponding to map (48) given as follows:(49)x⟶x+xd−kx−ya+x−e,y⟶y+y−r+bxa+x,μ2+μ∗=epx−c,where μ∗≪1 is taken as a small perturbation parameter.Furthermore, from (24), we have(50)z1⟶a1z1+a2z2+a11z12+a12z1z2+a111z13+a112z12z2+Oz1+z24,z2⟶b1z1+b2z2+b11z12+b12z1z2+b111z13+b112z12z2+Oz1+z24,where coefficients a1,a2,a11,a12,a111,a112 and b1,b2,b11,b12,b111,b112 given in (24) are similar to the coefficients in (36) by replacing μ1 with μ2+μ∗. Moreover, it must be noted that the characteristic equation corresponding to the linearization of (50) about z1,z2=0,0 is given as follows:(51)λ2+pμ∗λ+qμ∗=0,where(52)pμ∗=hμ∗−qμ∗−1,qμ∗=1arp/b−r−c2R2μ−S2−1arpb−r−c2,hμ∗=abx0y0a+x03>0.Suppose thata,b,c,d,k,p,r,μ2,h∈HB; then eigenvalues about 0,0 are conjugate complex numbers denoted by λ and λ¯ with modulus 1 such that(53)λ,λ¯=−pμ∗2±i24qμ∗−p2μ∗,and, consequently, one has(54)λμ∗=0=q0=1,l=dλdμ∗|μ∗=0=R22r1−c2≠0.On the other hand, it is necessary that, atμ∗=0, one must have λm,λ¯m≠1,m=1,2,3,4, or equivalently one has p0≠−2,0,1,2. Next, it must be noted that p0=h0−q0−1=Hμ2−2 and Hμ2<4 because a,b,c,d,k,p,r,μ2,h∈HB. As a result, one has p0≠−2,2. Consequently, it is required that p0≠0,1, which is equivalently written as follows:(55)Hμ2≠2,3.Consequently, the eigenvalues of system (50) which are given by λ,λ¯ about fixed point 0,0 do not lie in the intersection of the unit circle with the coordinate axes at μ∗=0 whenever conditions (55) hold true. Moreover, in order to discuss the normal form for system (50) about μ∗=0, we choose μ∗=0, σ=1−Hμ2/2, ω=1/2Hμ24−Hμ2,(56)T=a20σ−a1−ω,taking into account the following translation:(57)z1z2=Tuv.Then, it is easy to see that system (50) takes the following form:(58)uv⟶σ−ωωσuv+f1¯u,vf2¯u,v,where(59)f1¯u,v=a11a2z12+a12a2z1z2+a111a2z13+a112a2z12z2+Oz1+z24f2¯u,v=σ−a1a11−a2b11a2ωz12+σ−a1a12−a2b12a2ωz1z2+σ−a1a111−a2b111a2ωz13+σ−a1a112−a2b112a2ωz12z2+Oz1+z24z1=a2u,z2=σ−a1u−ωv,z1z2=a2σ−a1u2+a2ωuv,z12=a22u2,z13=a23u3,z12z2=a22σ−a1u3−ωu2v.Therefore, at point0,0, we have(60)f¯1uu=2a11+2a12σ−a1,f¯1uv=a12ω,f¯1uuu=6a111a22,f¯1uuv=−a112a2ω,f¯1vv=f¯1uvv=f¯1vvv=0,f¯2uu=2ωσ−a1a11a2+a12σ−a1−a2b12−a22b11,f¯2uv=σ−a1a12−a2b12,f¯2uuu=6a2ωσ−a1a2a111+σ−a1a112−a2b112−a2b111,f¯2uuv=2a2a2b112−σ−a1a112,f¯2vv=f¯2uvv=f¯2vvv=0.Then, system (58) undergoes a Neimark–Sacker bifurcation, if we have the following nonzero discriminatory quantity [30, 32]:(61)θ=−Re1−2λλ¯21−λL11L12−12L112−L212+Reλ¯L22μ∗=0,where(62)L11=14f¯1uu+f¯1vv+if¯2uu+f¯2vv,L12=18f¯1uu−f¯1vv+2f¯2uv+if¯2uu−f¯2vv−2f¯1uv,L21=18f¯1uu−f¯1vv−2f¯2uv+if¯2uu−f¯2vv+2f¯1uv,L22=116f¯1uuu+f¯1uvv+f¯2uuv+f¯2vvv+if¯2uuu+f¯2uvv−f¯2uuv−f¯2vvv,Keeping in mind the aforementioned computation, one has the following result.Theorem 3.
Assume that condition (55) is satisfied and θ≠0, then the Neimark–Sacker bifurcation exists at the fixed point x0,y0,e0 in a small neighborhood of μ2. Further, if θ<0 (resp., θ>0 ), then an attracting (resp., repelling) closed invariant curve bifurcates from the fixed point for μ>μ2 (resp., μ<μ2 ).


## 3.1. Period-Doubling Bifurcation

First, we start our investigation related to period-doubling bifurcation for system (4) about its equilibrium x0,y0,e0 with a variation of parameters in a small neighborhood of FB. For this, we choose parameters a,b,c,d,k,p,r,μ1,h arbitrarily from FB, taking into account system (4) with a,b,c,d,k,p,r,μ1,h∈FB. In this case, system (4) is described by the following 2-dimensional map:(17)x⟶x+hxd−kx−ya+x−e,y⟶y+hy−r+bxa+x,μ1=epx−c.It is easy to see that system (17) has a unique positive fixed point x0,y0,e0 such that the eigenvalues are given by λ1=−1 and λ2=1−P with λ2≠−1. Assume that a,b,c,d,k,p,r,μ1,h∈FB with μ1=1/R1S1−2arp/b−r−c2. Next, we take μ∗ as a bifurcation parameter, considering a perturbation for system (17) as follows:(18)x⟶x+hxd−kx−ya+x−e,y⟶y+hy−r+bxa+x,μ1+μ∗=epx−c,where μ∗≪1 is taken as a small perturbation parameter.Then, system (17) can be described in the following way:(19)x,yT⟶fx,y,e,gx,y,e=0,where(20)fx,y,e≔f1x,y,e,f2x,y,eT,f1x,y,e=x+hxd−kx−ya+x−e,f2x,y,e=y+hy−r+bxa+x,gx,y,e=epx−c−μ1.Then, it is easy to see thatDgx0,y0,e=e0p,0,px0−c such that the rank of Dgx0,y0,e0=1. On the other hand, a local parameterization ψ for the 2-dimensional smooth manifold defined by Mg=x,y,e∈R3;gx,y,e=0,and rankDgx,y,e=1 about x,y,e∈Bx0,y0,e0⊂Mg is given as follows: for any x,y,e∈Bx0,y0,e0, there is Z∈A⊂R2 such that(21)X=ψZ=X0+U0Z+V0HZ,gψZ=0,where(22)U0=100100,V0=001,X=x,y,e, X0=x0,y0,e0, Z=z1,z2, and H:R2⟶R is a smooth mapping. For further details, the interested reader is referred to [36]. Taking into account the definition of ψ, we obtain the following:(23)Dψz=DgXU0T−1001001,for arbitrarily chosen x,y,e∈Bx0,yo,e0. Then, system (19) is written as follows:(24)Z⟶fZ,Z∈A⊂R2,where A=ψ−1Bx0,yo,e0.Taking into account (24), we have(25)z1⟶a1z1+a2z2+a13z1μ∗+a11z12+a12z1z2+a113z12μ∗+a111z13+a112z12z2+Oz1+z2+μ∗4,z2⟶b1z1+b2z2+b13z1μ∗+b11z12+b12z1z2+b113z12μ∗+b111z13+b112z12z2+Oz1+z2+μ∗4.In order to compute the coefficients related the normal form, we need some further computation. From the aforementioned computation, one can easily get(26)Df1X=1+hd−2kx−aya+x2−e,−hxa+x,−hxDf2X=abhya+x2,1+h−r+bxa+x,0DgX=ep,0,px−cDψz=DgXU0T−1001001=ep0px−c100010−1001001=1001−eppx−c0=Dz1ψz,Dz2ψz,and from this, it follows that(27)Dz1ψz=1,0,−eppx−cT,Dz2ψz=0,1,0T.Assume thatfiziψZ represents the derivative of fiψZ with respect to zi and taking fiziX=fiziψZ. In a similar way, one can adopt notations for fizizjψZ and fizizjzkψZ. Next, it is easy to see that(28)f1z1X=Df1XDz1ψz=1+hd−2kx−aya+x2−e+epxpx−c,f1z2X=Df1XDz2ψz=−hxa+x,f2z1X=Df2XDz1ψz=abhya+x2,f2z2X=Df2XDz2ψz=1+h−r+bxa+x.PuttingX0 into (28), one has the following:(29)f1z1X0=1+h−kx0+e0px0px0−c+x0d−kx0−e0a+x0=1+hrbd+akb+rb−r+hrbb−r2ap+capr−cb+cr2μ1+μ∗,f1z2X0=−x0a+x0=−hrb,f2z1X0=abhy0a+x02=abhd−kx0−e0a+x0=hb−rd−karb−r−μ1+μ∗b−rapr−cb+cr,f2z2X0=1+h−r+bx0a+x0=1.Then, from (28), it is easy to see that(30)Df1z1X=h−2k+2aya+x3−epcpx−c2,−aha+x2,h−1+pxpx−c,Df1z2X=−aha+x2,0,0,Df2z1X=−2abhya+x3,abha+x2,0,Df2z2X=abha+x2,0,0.Consequently, one has the following:(31)f1z1z1X=Df1z1XDz1ψz=h−2k+2aya+x3−2epcpx−c2,f1z1z2X=Df1z1XDz2ψz=−aha+x2,f2z1z1X=Df2z1XDz1ψz=−2abhya+x3,f2z1z2X=Df2z1XDz2ψz=abha+x2.Putting the value ofX0 in (31), one has the following:(32)f1z1z1X0=h−2k+2ad−kx0−e0a+x02−2e0pcpx0−c2=h−2k+2db−r2ab2−2krb−rb2−2b−r3apr−cb+cr1ab2+pcapr−cb+cr2μ1+μ∗,f1z1z2X0=−aha+x02=−hb−r2ab2,f2z1z1X0=−2abhd−kx0−e0a+x02=2hb−r2abakrb−r−d+2hb−r3abapr−cb+crμ1+μ∗,f2z1z2X0=abha+x02=hb−r2ab.Then, from (31), it follows that(33)Df1z1z1X=−6ahya+x4+4hep2cpx−c3,2aha+x3,−2hpcpx−c2,Df1z1z2X=2aha+x3,0,0,Df2z1z1X=6abhya+x4,−2abha+x3,0,Df2z1z2X=−2abha+x3,0,0.Then, we get(34)f1z1z1z1X=Df1z1z1XDz1ψz=−6ahya+x4+6hep2cpx−c3,f1z1z1z2X=Df1z1z1XDz2ψz=2aha+x3,f2z1z1z1X=Df2z1z1XDz1ψz=6abhya+x4,f2z1z1z2X=Df2z1z1XDz2ψz=−2abha+x3.Putting the value ofX0 into (34), we obtain that(35)f1z1z1z1X0=−6ahd−kx0−e0a+x03+6he0p2cpx0−c3=6hkrb−r2ab3−6hdb−r3a2b3+6hb−r4apr−cb+cr1a2b3+p2capr−cb+cr3μ1+μ∗,f1z1z1z2X0=2aha+x03=2hb−r3a2b3,f2z1z1z1X0=6abhd−kx0−e0a+x03=6hb−r3a2b2d−akrb−r−6hb−r4a2b2apr−cb+crμ1+μ∗,f2z1z1z2X0=−2abha+x03=−2hb−r3a2b2.Consequently, it follows from (29), (32), and (35) that(36)a1=1+hrbd+akb+rb−r+hrbb−r2ap+capr−cb+cr2μ1,a2=−hrb,a13=hrbb−r2ap+capr−cb+cr2,a11=h−2k+2db−r2ab2−2krb−rb2−2b−r3apr−cb+cr1ab2+pcapr−cb+cr2μ1,a12=−hb−r2ab2,a113=−2hb−r3apr−cb+cr1ab2+pcapr−cb+cr2a111=6hkrb−r2ab3−6hdb−r3a2b3+6hb−r4apr−cb+cr1a2b3+p2capr−cb+cr3μ1,a112=2hb−r3a2b3,b1=hb−rd−b−rapr−cb+crμ1−akr,b2=1,b13=−hb−r2apr−cb+cr,b11=2hb−r2abakrb−r−d+2hb−r3abapr−cb+crμ1,b12=hb−r2ab,b113=2hb−r3abapr−cb+cr,b111=6hb−r3a2b2d−akrb−r−6hb−r4a2b2apr−cb+crμ1,b112=−2hb−r3a2b2.Next, we consider the following nonsingular matrix:(37)T=a2a2−1−a1λ2−a1,implementing the following translation:(38)z1z2=Tuv.Then, it is easy to see that system (24) can be written as follows:(39)uv⟶−100λ2uv+f1u,v,μ∗f2u,v,μ∗,where(40)f1u,v,μ∗=a11λ2−a1−a2b11a21+λ2z12+a12λ2−a1−a2b12a21+λ2z1z2+a13λ2−a1−a2b13a21+λ2z1μ∗+a113λ2−a1−a2b113a21+λ2z12μ∗+a111λ2−a1−a2b111a21+λ2z13+a112λ2−a1−a2b112a21+λ2z12z2+Oz1+z2+μ∗4f2u,v,μ∗=a111+a1+a2b11a21+λ2z12+a121+a1+a2b12a21+λ2z1z2+a131+a1+a2b13a21+λ2z1μ∗+a1131+a1+a2b113a21+λ2z12μ∗+a1111+a1+a2b111a21+λ2z13+a1121+a1+a2b112a21+λ2z12z2+Oz1+z2+μ∗4z1=a2u+v,z2=−1+a1u+λ2−a1v,z1z2=a2−1+a1u2+λ2−2a1−1uv+a2λ2−a1v2,z12=a22u2+uv+v2,z13=a23u3+3u2v+3uv2+v3,z12z2=a22−1+a1u3+λ2−3a1−2u2v+2λ2−3a1−1uv2+λ2−1v3.In order to determine the center manifoldWc0,0,0 of (39) about the equilibrium point 0,0 in a small neighborhood of μ∗, we implement an application of the center manifold theorem [29], and it is easy to see that there exists a center manifold of the following form:(41)Wc0,0,0=u,v,μ∗∈R3,v=hu,μ∗,h0,0=0,Dh0,0=0,foru and μ∗ sufficiently small. Furthermore, we assume that(42)hu,μ∗=c0u2+c1uμ∗+c2μ∗2+Ou+μ∗3.Then, it is easy to see that the center manifold satisfies the following:(43)h−u+f1u,hu,μ∗,μ∗,μ∗=λ2hu,μ∗+f2u,hu,μ∗,μ∗.Taking into account (42) and (43) and then comparing coefficients for (43), one can easily obtain that(44)c0=11−λ221+a1a2a11−a121+a1−a2b12+a22b11,c1=a131+a1+a2b131+λ22,c2=0.Furthermore, taking into account (39), it is restricted to the center manifold Wc0,0,0 as follows:(45)F:u⟶−u+d11u2+d12uμ∗+d112u2μ∗+d122uμ∗2+d111u3+Ou+μ∗4,where(46)d11=1λ2+1λ2−a1a11a2−a121+a1+a113a2−b11a22+1+a1b12a2−b113a22,d12=1λ2+1a13λ2−1−b13a2,d112=1λ2+1λ2−a1a11a2c1+a12c1λ2−2a1−1+a13c0−b11a22c1−b12a2c1λ2−2a1−1−b13a2c0,d122=1λ2+1a13c1λ2−1−b13a2c1,d111=1λ2+1λ2−a1a11a2c0+a111a22+a12c0λ2−2a1−1−a112a21+a1−b11a22c0−b12a2c0λ2−2a1−1−b111a23+b112a221+a1.In order for system (45) to undergo a flip bifurcation, we require that two discriminatory quantities α1and α2 are not zero, where(47)α1=2∂2F∂μ∗∂u+∂F∂μ∗∂F∂u0,0=2d12,α2=12∂2F∂u22+13∂3F∂u30,0=2d111+d112.Keeping in view the aforementioned computation and bifurcation theory given in [36, 37], we have the following theorem.Theorem 2.
Assume thata,b,c,d,k,p,r,μ1∈FB. If α2≠0, then system (17) undergoes a flip bifurcation at the fixed point x0,y0,e0 when the parameter μ varies in a small neighborhood of μ1. Moreover, if α2>0 (resp., α2<0), then the period-2 orbits that bifurcate from x0,y0,e0 are stable (resp., unstable).In the last Section5 related to numerical simulation, we choose some parametric values for system (4) such that it undergoes period-doubling bifurcation about positive equilibrium as μ varies in the suitable interval (see Figure 4).Figure 4
Bifurcation diagrams and LLE for system (3) with a=5.8, b=8.8, c=2.1, d=3.6, k=2.5, r=1.6, p=2.2, h=0.98, μ∈0.19,0.22, and x0,y0,e0=1.28,0.6428,0.287. (a) Bifurcation diagram for xn. (b) Bifurcation diagram for yn. (c) Bifurcation diagram for en. (d) Largest Lyapunov exponents.
(a)(b)(c)(d)


## 3.2. Neimark–Sacker Bifurcation

Finally, we discuss the Neimark–Sacker bifurcation about equilibriumx0,y0,e0 with variation of parameters a,b,c,d,k,p,r,μ2,h in a small neighborhood of HB. For this, we choose parameters a,b,c,d,k,p,r,μ2,h arbitrarily from HB, taking into account system (4) with a,b,c,d,k,p,r,μ2,h∈HB. In this case, system (4) is described by the following 2-dimensional map:(48)x⟶x+xd−kx−ya+x−e,y⟶y+y−r+bxa+x,μ2=epx−c.Then, it is easy to observe that map (48) has a unique positive fixed point x0,y0,e0. Next, we take μ∗ as a bifurcation parameter and consider a perturbation corresponding to map (48) given as follows:(49)x⟶x+xd−kx−ya+x−e,y⟶y+y−r+bxa+x,μ2+μ∗=epx−c,where μ∗≪1 is taken as a small perturbation parameter.Furthermore, from (24), we have(50)z1⟶a1z1+a2z2+a11z12+a12z1z2+a111z13+a112z12z2+Oz1+z24,z2⟶b1z1+b2z2+b11z12+b12z1z2+b111z13+b112z12z2+Oz1+z24,where coefficients a1,a2,a11,a12,a111,a112 and b1,b2,b11,b12,b111,b112 given in (24) are similar to the coefficients in (36) by replacing μ1 with μ2+μ∗. Moreover, it must be noted that the characteristic equation corresponding to the linearization of (50) about z1,z2=0,0 is given as follows:(51)λ2+pμ∗λ+qμ∗=0,where(52)pμ∗=hμ∗−qμ∗−1,qμ∗=1arp/b−r−c2R2μ−S2−1arpb−r−c2,hμ∗=abx0y0a+x03>0.Suppose thata,b,c,d,k,p,r,μ2,h∈HB; then eigenvalues about 0,0 are conjugate complex numbers denoted by λ and λ¯ with modulus 1 such that(53)λ,λ¯=−pμ∗2±i24qμ∗−p2μ∗,and, consequently, one has(54)λμ∗=0=q0=1,l=dλdμ∗|μ∗=0=R22r1−c2≠0.On the other hand, it is necessary that, atμ∗=0, one must have λm,λ¯m≠1,m=1,2,3,4, or equivalently one has p0≠−2,0,1,2. Next, it must be noted that p0=h0−q0−1=Hμ2−2 and Hμ2<4 because a,b,c,d,k,p,r,μ2,h∈HB. As a result, one has p0≠−2,2. Consequently, it is required that p0≠0,1, which is equivalently written as follows:(55)Hμ2≠2,3.Consequently, the eigenvalues of system (50) which are given by λ,λ¯ about fixed point 0,0 do not lie in the intersection of the unit circle with the coordinate axes at μ∗=0 whenever conditions (55) hold true. Moreover, in order to discuss the normal form for system (50) about μ∗=0, we choose μ∗=0, σ=1−Hμ2/2, ω=1/2Hμ24−Hμ2,(56)T=a20σ−a1−ω,taking into account the following translation:(57)z1z2=Tuv.Then, it is easy to see that system (50) takes the following form:(58)uv⟶σ−ωωσuv+f1¯u,vf2¯u,v,where(59)f1¯u,v=a11a2z12+a12a2z1z2+a111a2z13+a112a2z12z2+Oz1+z24f2¯u,v=σ−a1a11−a2b11a2ωz12+σ−a1a12−a2b12a2ωz1z2+σ−a1a111−a2b111a2ωz13+σ−a1a112−a2b112a2ωz12z2+Oz1+z24z1=a2u,z2=σ−a1u−ωv,z1z2=a2σ−a1u2+a2ωuv,z12=a22u2,z13=a23u3,z12z2=a22σ−a1u3−ωu2v.Therefore, at point0,0, we have(60)f¯1uu=2a11+2a12σ−a1,f¯1uv=a12ω,f¯1uuu=6a111a22,f¯1uuv=−a112a2ω,f¯1vv=f¯1uvv=f¯1vvv=0,f¯2uu=2ωσ−a1a11a2+a12σ−a1−a2b12−a22b11,f¯2uv=σ−a1a12−a2b12,f¯2uuu=6a2ωσ−a1a2a111+σ−a1a112−a2b112−a2b111,f¯2uuv=2a2a2b112−σ−a1a112,f¯2vv=f¯2uvv=f¯2vvv=0.Then, system (58) undergoes a Neimark–Sacker bifurcation, if we have the following nonzero discriminatory quantity [30, 32]:(61)θ=−Re1−2λλ¯21−λL11L12−12L112−L212+Reλ¯L22μ∗=0,where(62)L11=14f¯1uu+f¯1vv+if¯2uu+f¯2vv,L12=18f¯1uu−f¯1vv+2f¯2uv+if¯2uu−f¯2vv−2f¯1uv,L21=18f¯1uu−f¯1vv−2f¯2uv+if¯2uu−f¯2vv+2f¯1uv,L22=116f¯1uuu+f¯1uvv+f¯2uuv+f¯2vvv+if¯2uuu+f¯2uvv−f¯2uuv−f¯2vvv,Keeping in mind the aforementioned computation, one has the following result.Theorem 3.
Assume that condition (55) is satisfied and θ≠0, then the Neimark–Sacker bifurcation exists at the fixed point x0,y0,e0 in a small neighborhood of μ2. Further, if θ<0 (resp., θ>0 ), then an attracting (resp., repelling) closed invariant curve bifurcates from the fixed point for μ>μ2 (resp., μ<μ2 ).


## 4. Chaos Control

Bifurcating behavior, chaos, and unstable fluctuations have always been considered as adverse criteria in biology; therefore, these are damaging for the reproduction of the biological population. Certainly, we require to take action to stabilize the biological population. For further details and applications of chaos control methods, we refer to [37–50] and the references therein.In this section, we propose the following state delayed feedback control method for system (4):(63)xn+1=xn+hxnd−kxn−yna+xn−en+δxn−xn−1,yn+1=yn+hynbxna+xn−r,μ=enpxn−c,where δ is the feedback gain for the controlled system (63).Next, introducingun≔xn−xn−1, we obtain the following controlled system equivalent to system (63):(64)xn+1=xn+hxnd−kxn−yna+xn−en+δun,yn+1=yn+hynbxna+xn−r,un+1=hxnd−kxn−yna+xn−en+δun,μ=enpxn−c.The generalized variational matrixJx0,y0,u0,e0 of the controlled system (64) about its fixed point x0,y0,u0,e0=ar/b−r,a+x0d−kx0−e0,0,μ/px0−c is given as follows:(65)Jx0,y0,u0,e0=m11−hx0a+x0δ−hx0abhy0a+x021001−m11−hx0a+x0δ−hx0pe000px0−c,where(66)m11≔1+hx0y0a+x02−hkx0.Then, generalized characteristic polynomial forJx0,y0,u0,e0 is defined as follows:(67)Fλ=detm11−λ−hx0a+x0δ−hx0abhy0a+x021−λ001−m11−hx0a+x0δ−λ−hx0pe000px0−c=0,which, on simplification, yields(68)Fλ=λ3+Aλ2+Bλ+C,where(69)A≔px0δ+e0h+m11+1−cδ+m11+1c−px0,B≔cabh2x0y0/a+x03+2δm11+m11+px0−abh2x0y0/a+x03−e0h−2δm11−m11c−px0,C≔δ−2δm11.Keeping in view the controllability for system (63), we have the following lemma.Lemma 2.
Assume thatb>r,bc<rap+c,bd−rak+d/b−r2−μ/rap+c−bc>0; then the fixed point x0,y0,u0,e0=ar/b−r,a+x0d−kx0−e0,0,μ/px0−c of system (63) is a sink if the following conditions are satisfied:(70)A+C<1+B,A−3C<3−B,C2+B−AC<1.


## 5. Numerical Simulation and Discussion

In this section, our main purpose is to validate theoretical findings with numerical simulation. For this, first of all, the existence of the period-doubling bifurcation is illustrated through particular choice of biologically feasible parametric values. Choosea=5.8, b=8.8, c=2.1, d=3.6, k=2.5, r=1.6, p=2.2, and h=0.98 with the variation of bifurcation parameter μ in 0.19,0.22. Then, around μ≈0.211175, system (3) undergoes period-doubling bifurcation. On the other hand, at a=5.8, b=8.8, c=2.1, d=3.6, k=2.5, r=1.6, p=2.2, h=0.98, and μ≡μ1=0.211175, system (3) has a unique positive fixed point 1.28889,0.642835,0.287096, and the characteristic polynomial for the Jacobian matrix of singular system (3) is given as follows:(71)Pλ=λ2+0.057005λ−0.942995.Moreover, the roots ofPλ are −1 and 0.942995. Consequently, a,b,c,d,k,p,r,μ1,h∈FB, and it follows the correctness of Theorem 2. On the other hand, bifurcation diagrams for singular system (3) and corresponding largest Lyapunov exponents (LLE) are depicted in Figure 4.Secondly, choosea=41.5, b=85, c=2.7, d=3.4, k=1.4, r=3.5, p=2.3, and h=0.75 with a variation of bifurcation parameter μ in 0.4,1.1 for validity of Neimark–Sacker bifurcation. Then, around μ≈0.678396, system (3) undergoes a Neimark–Sacker bifurcation. On the other hand, at a=41.5, b=85, c=2.7, d=3.4, k=1.4, r=3.5, p=2.3, h=0.75, and μ≡μ2=0.678396, system (3) has a unique positive fixed point 1.78221,18.1794,0.484887, and the characteristic polynomial for the Jacobian matrix of singular system (3) is given as follows:(72)Pλ=λ2−1.20713λ+1.Moreover, the roots ofPλ are 0.603567−0.797313i and 0.603567+0.797313i with 0.603567±0.797313i=1. Consequently, a,b,c,d,k,p,r,μ2,h∈HB, and it follows the correctness of Theorem 3. On the other hand, the bifurcation diagrams for singular system (3) and corresponding largest Lyapunov exponents (LLE) are depicted in Figure 5. Moreover, in the chaotic region, that is, for μ∈0.678396,1.1, some phase portraits of system (4) are depicted in Figure 6.Figure 5
Bifurcation diagrams and LLE for system (3) with a=41.5, b=85, c=2.7, d=3.4, k=1.4, r=3.5, p=2.3, h=0.75, μ∈0.4,1.1, and x0,y0,e0=1.5,5.58,0.775. (a) Bifurcation diagram for xn. (b) Bifurcation diagram for yn. (c) Bifurcation diagram for en. (d) Largest Lyapunov exponents.
(a)(b)(c)(d)Figure 6
Phase portraits of system (4) for a=41.5, b=85, c=2.7, d=3.4, k=1.4, r=3.5, p=2.3, and h=0.75 and with different values of μ∈0.678396,1.1. (a) Phase portrait for μ = 0 : 678396. (b) Phase portrait for μ = 0 : 7. (c) Phase portrait for μ = 0 : 75. (d) Phase portrait for μ = 0 : 8. (e) Phase portrait for μ = 0 : 85. (f) Phase portrait for μ = 0 : 95. (g) Phase portrait for μ = 1. (h) Phase portrait for μ = 1 : 08.
(a)(b)(c)(d)(e)(f)(g)(h)Finally, in order to see the efficaciousness of delayed feedback control strategy for system (4), we choose a=41.5, b=85, c=2.7, d=3.4, k=1.4, r=3.5, p=2.3, h=0.75, and μ=0.95 in the chaotic region. At these parametric values, system (4) has unique positive fixed point x0,y0,e0=1.78221,9.77703,0.679018, and the complex conjugate roots of the characteristic equation for the Jacobian matrix are 0.813859±0.625908i with 0.813859±0.625908i=1.02671>1. Therefore, 1.78221,9.77703,0.679018 is a source for system (4). Furthermore, for these parametric values, system (64) can be written as follows:(73)xn+1=xn+0.75xn3.4−1.4xn−yn41.5+xn−en+δun,yn+1=yn+0.75yn85xn41.5+xn−3.5,un+1=0.75xn3.4−1.4xn−yn41.5+xn−en+δun,0.95=en2.3xn−2.7.On the other hand, the characteristic polynomial of the Jacobian matrix of system (73) is given as follows:(74)Pλ=λ3−δ+1.62772λ2+2δ+1.05413λ−δ.Taking into account the conditions of Lemma2, we have that the positive fixed point of system (73) is a sink if −0.920461<δ<−0.14539. Consequently, system (73) is controllable in δ∈−0.920461,−0.14539. More generally, if we take a=41.5, b=85, c=2.7, d=3.4, k=1.4, r=3.5, p=2.3, h=0.75, and μ∈0.678396,1.1, then the controllable region for system (64) is depicted in Figure 7.Figure 7
Controllability region of system (64) for a=41.5, b=85, c=2.7, d=3.4, k=1.4, r=3.5, p=2.3, h=0.75, and μ∈0.678396,1.1 in μδ-plane.


## 6. Conclusion

We discuss the dynamical behavior of a discrete-time singular bioeconomic model. Moreover, Euler’s approximation is implemented to a bioeconomic model governed by the differential-algebraic system proposed in [11]. The Neimark–Sacker bifurcation and period-doubling bifurcation are studied for the discrete bioeconomic model with the implementation of normal form theory, bifurcation theory, and the center manifold theorem. On the other hand, we select μ (the economic profit parameter) as a bifurcation parameter. Our investigations show richer dynamical behaviors for the discrete-time bioeconomic model compared with its continuous counterpart studied in [11]. Numerical computation for maximum Lyapunov exponents ensures further dynamical behavior and complexity of the model. Such type of complex phenomena might result from economic profit [25]. With the variation of the bifurcation parameter μ, the biologically feasible fixed point resembles the stability of economic profit, and later on, the system may sacrifice its stability by undergoing period-doubling or Neimark–Sacker bifurcation, and consequently trajectories tend to a period-doubling cascade or an invariant circle. Our theoretical discussion reveals that if the economic revenue μ increases beyond a certain threshold value μ1 (respectively, μ2), a phenomenon of period-doubling bifurcation (respectively, Neimark–Sacker bifurcation) occurs. From Theorem 2 (respectively, Theorem 3), if the economic revenue μ is equal to or larger than the bifurcation value μ1 (respectively, μ2), the predator population, the prey population, and the harvesting effort will not stay at steady states, which will result in serious biological economic environmental imbalance. Thus, it is sensible for fishermen to keep the economic revenue within a certain range for the purpose of maintaining the sustainable development of biological resources. On the other hand, the suggested state feedback control method of the delayed type vanishes such fluctuating phenomena and derives the bifurcating behavior of singular discrete-time economic prey-predator system towards a stable situation. Computation reveals that this delayed-type control strategy can be applied only by modifying the effort of harvesting for the prey population density by taking into account the present prey population density and the previous one. Consequently, economists can formulate some strategies to restrain or encourage the harvesting attempts in practical applications, for example, adjusting market price and revenue, abating pollution, making an allowance to fishermen, so that the biological populations can remain at their stable states.



---

*Source: 6679161-2021-07-13.xml*