# Chaotic Dynamics and Control of Discrete Ratio-Dependent Predator-Prey System

Authors: Sarker Md. Sohel Rana
Journal: Discrete Dynamics in Nature and Society
 (2017)

---

## Abstract

This study examines the complexity of a discrete-time predator-prey system with ratio-dependent functional response. We establish algebraically the conditions for existence of fixed points and their stability. We show that under some parametric conditions the system passes through a bifurcation (flip or Neimark-Sacker). Numerical simulations are presented not only to justify theoretical results but also to exhibit new complex behaviors which include phase portraits, orbits of periods 9, 19, and 26, invariant closed circle, and attracting chaotic sets. Moreover, we measure numerically the Lyapunov exponents and fractal dimension to confirm the chaotic dynamics of the system. Finally, a state feedback control method is applied to control chaos which exists in the system.


---

## Body

## 1. Introduction

The interaction between predator and prey is one of the most studied topics in ecology and mathematical biology. The Lotka-Volterra model [1, 2] has received more attention from mathematician and ecologist. After them, many authors qualitatively analysed other types of predator-prey models in which response function of predator depends on prey densities only. But a number of respectable researchers [3–5] have claimed that in some environments (when predator needs to search and share foods), the response function of predator may depend on the ratio of prey to predator abundance. For detailed results in predator-prey systems with ratio-dependent response function, we refer to [6–9]. However, when the size of populations is small, the discretization of predator-prey systems is more suitable compared to continuous ones to understand unpredictable dynamic behaviors which exist in the system. Lots of empirical and theoretical works [10–20] have discussed the bifurcations and chaos phenomenon by using numerical simulations or center manifold theory and bifurcation theory.In recent times, there are a few number of articles discussing the dynamics of ratio-dependent discrete-time predator-prey systems [21–23]. For discrete ratio-dependent predator-prey systems, it is shown that positive equilibrium is globally asymptotically stable [21], the strong and the weak Allee effects are investigated in [22], and it is discussed that periodic solutions exist in [23].A ratio-dependent predator-prey system takes the form:(1)x˙=xgx-cfxyy,y˙=y-D+efxy,where x and y stand for densities of prey and predator, respectively, D>0 is predator’s death rate, and gx is prey’s specific growth rate. fx/y is ratio-dependent response function of predator and c,e>0 are rate of capturing and conversion, respectively.Ifgx=r1-x/K and fx=αx/(β+x), then system (1) takes in the form [6]:(2)x˙=rx1-xK-cαxyβy+x,y˙=y-D+eαxβy+x,where r, K, α, and β are all positive constants. By scaling the variables and parameters as x→x/K, y→βy/K, and t→rt, we write system (2) as (3)x˙=x1-x-axyy+x,y˙=y-d+bxy+xwith a=cα/βr, b=eα/r, and d=D/r being positive constants. Forward Euler scheme is applied to system (3) to obtain the following system:(4)H:xy⟶x+δx1-x-axyy+xy+δ-dy+bxyy+x,where δ is the step size. The motivation of this work is to study systematically the dynamical properties of (4) in detail which includes phase portraits, bifurcation in orbits of periods 9, 19, and 26, chaotic sets, Lyapunov exponents, and fractal dimension.This paper is organized as follows. Section2 deals with the existence condition for fixed points of (4) and their stability criterion. In Section 3, we prove that system (4) undergoes a flip bifurcation and a NS bifurcation in the interior of R+2 when δ changes its value in a small neighborhood of a specific parametric space. In Section 4, we provide numerical simulations for one or more control parameters to validate analytical results. In Section 5, we use the method of feedback control to stabilize chaos at the state of unstable trajectories. Finally, a brief discussion is carried out in Section 6.


## 2. Fixed Points: Existence and Their Stability

The fixed points of (4) are the solution of the following equations: (5)x+δx1-x-axyy+x=x,y+δ-dy+bxyy+x=y.A simple algebraic computation yields the following result.Lemma 1.
System (4) has two fixed points:
(i) the predator free fixed pointE11,0, which always exists;
(ii) the interior fixed pointE2x∗,y∗, where x∗=1+a-1+d/b and y∗=(b-d)/d1+a-1+d/b. This fixed point exists if d<b<ad/(a-1) with a>1.To show the region in the spaced,a,b for which positive fixed point of system (4) exists, we plot C=d,a,b:d<b<ad/(a-1)witha>1 in Figure 1. We see that E2 lies inside C but not in regions (I)-(II).Figure 1
Graphical depiction for positive fixed point.Next, we investigate stability of (4) at fixed points. The Jacobian matrix of system (4) at fixed point Ex,y is(6)Jx,y=1+δa1-δb1δa21+δb2,where(7)a1=1-ayx+y+x-2+ayx+y2,b1=ax2x+y2,a2=by2x+y2,b2=bx2x+y2-d.Then the characteristic equation of (6) is(8)Fλ≔λ2-tr⁡Jλ+det⁡J=0,where(9)tr⁡J=2+δa1+b2,det⁡J=1+δa1+b2+δ2a1b2+a2b1.It is noted that the magnitude of eigenvalues of Jacobian matrix evaluated at fixed point Ex,y determines the local stability of that fixed point. The following propositions represent the stability conditions of fixed points by Jury’s criterion [24].Proposition 2.
For all permissible parameters values, there are four different topological types ofE1 and it is a(i)
sink if0<δ<2 and d-2/δ<b<d;(ii)
source ifδ>2 and (b<d-2/δ and b>d);(iii)
nonhyperbolic ifδ=2 or (b=d-2/δ and b≠d);(iv)
saddle if otherwise.Condition (iii) can be written as(10)FBE11=a,b,d,δ∈0,+∞:δ=2,b≠d-2δ,b≠d,or (11)FBE12=a,b,d,δ∈0,+∞:b=d-2δ,δ≠2,b≠d.It follows that if parametersa,b,d,δ belong to FBE11orFBE12, then two eigenvalues of JE1 are λ1=-1 and λ2≠-1,1. Thus a flip bifurcation may occur if parameters belong to FBE11orFBE12. This observation illustrates that predator goes to extinct and prey passes through flip bifurcation leading to chaos as bifurcation parameter δ varies.Proposition 3.
Ifd<b<ad/(a-1) with a>1, then positive fixed point E2 of (4) exists and it is a (i)
sink if either of the following conditions holds:
(i.1)
P∗≥0 and δ<(N∗-P∗)/M∗;(i.2)
P∗<0 and δ<N∗/M∗;(ii)
source if either of the following conditions holds:
(ii.1)
P∗≥0 and δ>(N∗+P∗)/M∗;(ii.2)
P∗<0 and δ>N∗/M∗;(iii)
nonhyperbolic if either of the following conditions holds:
(iii.1)
P∗≥0 and δ=(N∗±P∗)/M∗;(iii.2)
P∗<0 and δ=N∗/M∗;(iv)
saddle if otherwise,
where(12)M∗=a1b2+a2b1,N∗=-a1+b2,P∗=N∗2-4M∗.From Proposition3, it follows that if term (iii.1) holds, two eigenvalues of JE2 are λ1=-1 and λ2≠-1,1. We rewrite the term (iii.1) as follows: (13)FBE21=a,b,d,δ∈0,+∞:δ=N∗-P∗M∗,P∗≥0, or (14)FBE22=a,b,d,δ∈0,+∞:δ=N∗+P∗M∗,P∗≥0.Therefore, a flip bifurcation may emerge from fixed point E2 if parameters belong to FBE21 or FBE22.Also two eigenvaluesλ1,2 of JE2 are complex conjugate having magnitude one if (iii.2) holds. We rewrite the term (iii.2) as follows: (15)NSBE2=a,b,d,δ∈0,+∞:δ=N∗M∗,P∗<0,and if the parameters lie inNSBE2, there exists a NS bifurcation emerging from E2.


## 3. Analysis of Bifurcation

In this section, attention is paid to recapitulate bifurcations (flip and/or Neimark-Sacker) of system (4) around fixed points using the theories in center manifold and of bifurcation [25–28]. The parameter δ is chosen as a bifurcation parameter.


### 3.1. Flip Bifurcation at Fixed PointE11,0

We consider system (4) at fixed point E1 and take parameter a,b,d,δ arbitrarily from FBE11 (similarly for the case of FBE12).Letδ=δ1=2, or δ=δ1=-2/(b-d), for the case of FBE12. Then the eigenvalues of positive fixed point E11,0 are (16)λ1δ1=-1,λ2δ1=1+2b-d.The condition λ2δ1≠1 holds if(17)b-d≠0,-1.Using the transformationx~=x-1, y~=y and writing Aδ=J1,0, we shift the fixed point 1,0 of system (4) to the origin. After Taylor expansion, system (4) reduces to(18)x~y~⟶Aδx~y~+F1x~,y~,δF2x~,y~,δ,where X=x~,y~T and(19)F1x~,y~,δ=-ax~y~2-ay~3-x~2+ay~2δ+OX4,F2x~,y~,δ=bx~y~2+by~3-by~2δ+OX4.Therefore, we get the following symmetric multilinear vector functions ofx,y,u∈R2: (20)Bx,y=B1x,yB2x,y=2ax2y2-2x1y1δ-2bx2y2δ,(21)Cx,y,u=C1x,y,uC2x,y,u=-2x2y1u2-2x1y2u2-2x2y2u1-6x2y2u2aδ2x2y1u2+2x1y2u2+2x2y2u1+6x2y2u2bδ,where (22)Bix,y=∑j,k=12δ2Fiξ,δδξjδξkξ=0xjyk,Cix,y,u=∑j,k,l=12δ2Fiξ,δδξjδξkδξlξ=0xjykul,i=1,2,δ=δ1.Letp,q∈R2 be two eigenvectors (left and right) of A associated with eigenvalue λ1δ1=-1, respectively. Then Aδ1q=-q and ATδ1p=-p. After some algebraic calculation, we get(23)q~2+δ1b-d,0T,p~2+δ1b-d,2aT.We setp=γ12+δ1b-d,2aT, where γ1=1/2+δ1b-d2. Then by the standard scalar product in R2 defined by p,q=p1q1+p2q2, it is obvious to see that p,q=1. The direction of the flip bifurcation can be obtained by the sign of cδ1, the coefficient of critical normal form [25], and is computed via(24)cδ1=16p,Cq,q,q-12p,Bq,A-I-1Bq,q=2δ12+δ1b-d2≠0.The above discussion leads to the following result.Theorem 4.
If (17) holds and δ varies around the set FBE11, then system (4) passes through a flip bifurcation at fixed point E11,0. Moreover, the period-2 orbits that bifurcate from E11,0 are stable (resp., unstable) if cδ1>0 (resp., cδ1<0).


### 3.2. Flip and Neimark-Sacker Bifurcation at Fixed PointE2x∗,y∗

We first discuss flip bifurcation of system (4) at fixed point E2. We take parameter a,b,d,δ arbitrarily from FBE21 (similarly for the case of FBE22); then it follows that P∗>0; that is,(25)a1-b22-4a2b1>0.Let(26)δ=δF=N∗-P∗M∗,or(27)δ=δF=N∗+P∗M∗,forthe  case  ofFBE22.Then the eigenvalues of positive fixed pointx∗,y∗ are(28)λ1δF=-1,λ2δF=3-δFN∗.The condition λ2δF≠1 leads to(29)δFN∗≠2,4.Using the transformationx~=x-x∗, y~=y-y∗ and writing Aδ=Jx∗,y∗, we shift the fixed point x∗,y∗ of system (4) to the origin. After Taylor expansion, system (4) reduces to(30)x~y~⟶Aδx~y~+F1x~,y~,δF2x~,y~,δ,whereX=x~,y~T and(31)F1x~,y~,δ=16-6aδx∗2x∗+y∗4y~3-6aδy∗2x∗+y∗4x~3+2aδ2x∗-y∗y∗x∗+y∗43x~2y~-2aδx∗-2y∗x∗x∗+y∗43x~y~2+122aδx∗2x∗+y∗3y~2+-2δ+2aδy∗2x∗+y∗3x~2-2aδx∗y∗x∗+y∗32x~y~+OX4,F2x~,y~,δ=166bδx∗2x∗+y∗4y~3+6bδy∗2x∗+y∗4x~3-2bδ2x∗-y∗y∗x∗+y∗43x~2y~+2bδx∗-2y∗x∗x∗+y∗43x~y~2+12-2bδx∗2x∗+y∗3y~2-2bδy∗2x∗+y∗3x~2+2bδx∗y∗x∗+y∗32x~y~+OX4.It follows that(32)B1x,y=∑j,k=12δ2F1ξ,δδξjδξkξ=0xjyk=2aδx∗2x∗+y∗3x2y2+-2δ+2aδy∗2x∗+y∗3x1y1-2aδx∗y∗x∗+y∗3x1y2+x2y1,(33)B2x,y=∑j,k=12δ2F2ξ,δδξjδξkξ=0xjyk=-2bδx∗2x∗+y∗3x2y2-2bδy∗2x∗+y∗3x1y1+2bδx∗y∗x∗+y∗3x1y2+x2y1,(34)C1x,y,u=∑j,k,l=12δ2F1ξ,δδξjδξkδξlξ=0xjykul=-6aδx∗2x∗+y∗4x2y2u2-6aδy∗2x∗+y∗4x1y1u1+2aδ2x∗-y∗y∗x∗+y∗4x1y1u2+x1y2u1+x2y1u1+x2y1u2+x1y2u2+x2y2u1,(35)C2x,y,u=∑j,k,l=12δ2F2ξ,δδξjδξkδξlξ=0xjykul=6bδx∗2x∗+y∗4x2y2u2+6bδy∗2x∗+y∗4x1y1u1-2bδx∗-2y∗x∗x∗+y∗4x1y1u2+x1y2u1+x2y1u1+x2y1u2+x1y2u2+x2y2u1,and δ=δF.Therefore, we get the following symmetric multilinear vector functions ofx,y,u∈R2: Bx,y=B1x,yB2x,y and Cx,y,u=C1x,y,uC2x,y,u.Letp,q∈R2 be two eigenvectors (left and right) of A associated with eigenvalue λ1δF=-1, respectively. Then AδFq=-q and ATδFp=-p. After some algebraic calculation, we get(36)q~-2-δFb2,δFa2T,p~-2-δFb2,-δFb1T.We setp=γ1-2-δFb2,-δFb1T, where (37)γ1=1-2-δFb22-δF2a2b1.Then by the standard scalar product inR2 defined by p,q=p1q1+p2q2, it is obvious to see that p,q=1. The direction of the flip bifurcation can be obtained by the sign of cδF, the coefficient of critical normal form [25], and is computed via(38)cδF=16p,Cq,q,q-12p,Bq,A-I-1Bq,q.The above discussion leads to the following result.Theorem 5.
If (29) holds, cδF≠0, and δ varies around the set FBE21, then system (4) passes through a flip bifurcation at fixed point E2x∗,y∗. Moreover, the period-2 orbits that bifurcate from E2x∗,y∗ are stable (resp., unstable) if cδF>0 (resp., cδF<0).We next discuss a Neimark-Sacker bifurcation of system (4).We first take parametera,b,d,δ arbitrarily from NSBE2 and consider system (4) at fixed point E2x∗,y∗.Since the parameters belong toNSBE2, the complex eigenvalues of (6) are(39)λ,λ¯=tr⁡J2±i24det⁡J-tr⁡J2,where the condition tr⁡J2-4det⁡J<0 leads to a1-b22-4a2b1<0. Let(40)δ=δNS=N∗M∗;then det⁡JδNS=1 and(41)λ=det⁡Jδ,dλδdδδ=δNS=N∗2≠0.Moreover, iftr⁡JδNS≠0,-1, then δNSN∗≠2,3 which leads to(42)λkδNS≠1fork=1,2,3,4.Letq,p∈C2 be eigenvectors of AδNS and ATδNS corresponding to the eigenvalues λδNS and λ¯δNS such that (43)AδNSq=λδNSq,AδNSq¯=λ¯δNSq¯,ATδNSp=λ¯δNSp,ATδNSp¯=λδNSp¯.Then by algebraic calculation, we obtain(44)q~1+δNSb2-λ,-δNSa2T,p~1+δNSb2-λ¯,δNSb1T.We setp=γ21+δNSb2-λ¯,δNSb1T, where (45)γ2=11+δNSb2-λ¯2-δNS2a2b1.Then it is obvious thatp,q=1, where p,q=p1¯q2+p2¯q1 for p,q∈C2. We can decompose vector X∈R2 as X=zq+z¯q¯, for δ close to δNS and z∈C. Obviously, z=p,X. Thus, we obtain the following transformed form of system (30) for δ near δNS: (46)z⟼λδz+gz,z¯,δ,where λδ=1+φδeiθδ with φδNS=0 and gz,z¯,δ is a smooth complex-valued function. Taylor expression of g with respect to z,z¯ yields (47)gz,z¯,δ=∑k+l≥21k!l!gklδzkz¯l,withgkl∈C,k,l=0,1,….According to multilinear symmetric vector functions, we can express the coefficientsgkl as follows:(48)g20δNS=p,Bq,q,g11δNS=p,Bq,q¯,g02δNS=p,Bq¯,q¯,g21δNS=p,Cq,q,q¯.We calculate the coefficientaδNS by the following formula to determine the direction in which invariant closed curve appears.(49)aδNS=Ree-iθδNSg212-Re1-2eiθδNSe-2iθδNS21-eiθδNSg20g11-12g112-14g022, where eiθδNS=λδNS.It is clear that two conditions (41) and (42) known as transversal and nondegenerate hold well for system (4). Therefore, we have the following result.Theorem 6.
IfaδNS≠0 and δ changes its value around NSBE2, then system (4) passes through a Neimark-Sacker bifurcation at fixed point E2. Moreover, if aδNS<0 (resp., >0), then there exists an attracting (resp., repelling) invariant closed curve which bifurcates from E2.


## 3.1. Flip Bifurcation at Fixed PointE11,0

We consider system (4) at fixed point E1 and take parameter a,b,d,δ arbitrarily from FBE11 (similarly for the case of FBE12).Letδ=δ1=2, or δ=δ1=-2/(b-d), for the case of FBE12. Then the eigenvalues of positive fixed point E11,0 are (16)λ1δ1=-1,λ2δ1=1+2b-d.The condition λ2δ1≠1 holds if(17)b-d≠0,-1.Using the transformationx~=x-1, y~=y and writing Aδ=J1,0, we shift the fixed point 1,0 of system (4) to the origin. After Taylor expansion, system (4) reduces to(18)x~y~⟶Aδx~y~+F1x~,y~,δF2x~,y~,δ,where X=x~,y~T and(19)F1x~,y~,δ=-ax~y~2-ay~3-x~2+ay~2δ+OX4,F2x~,y~,δ=bx~y~2+by~3-by~2δ+OX4.Therefore, we get the following symmetric multilinear vector functions ofx,y,u∈R2: (20)Bx,y=B1x,yB2x,y=2ax2y2-2x1y1δ-2bx2y2δ,(21)Cx,y,u=C1x,y,uC2x,y,u=-2x2y1u2-2x1y2u2-2x2y2u1-6x2y2u2aδ2x2y1u2+2x1y2u2+2x2y2u1+6x2y2u2bδ,where (22)Bix,y=∑j,k=12δ2Fiξ,δδξjδξkξ=0xjyk,Cix,y,u=∑j,k,l=12δ2Fiξ,δδξjδξkδξlξ=0xjykul,i=1,2,δ=δ1.Letp,q∈R2 be two eigenvectors (left and right) of A associated with eigenvalue λ1δ1=-1, respectively. Then Aδ1q=-q and ATδ1p=-p. After some algebraic calculation, we get(23)q~2+δ1b-d,0T,p~2+δ1b-d,2aT.We setp=γ12+δ1b-d,2aT, where γ1=1/2+δ1b-d2. Then by the standard scalar product in R2 defined by p,q=p1q1+p2q2, it is obvious to see that p,q=1. The direction of the flip bifurcation can be obtained by the sign of cδ1, the coefficient of critical normal form [25], and is computed via(24)cδ1=16p,Cq,q,q-12p,Bq,A-I-1Bq,q=2δ12+δ1b-d2≠0.The above discussion leads to the following result.Theorem 4.
If (17) holds and δ varies around the set FBE11, then system (4) passes through a flip bifurcation at fixed point E11,0. Moreover, the period-2 orbits that bifurcate from E11,0 are stable (resp., unstable) if cδ1>0 (resp., cδ1<0).


## 3.2. Flip and Neimark-Sacker Bifurcation at Fixed PointE2x∗,y∗

We first discuss flip bifurcation of system (4) at fixed point E2. We take parameter a,b,d,δ arbitrarily from FBE21 (similarly for the case of FBE22); then it follows that P∗>0; that is,(25)a1-b22-4a2b1>0.Let(26)δ=δF=N∗-P∗M∗,or(27)δ=δF=N∗+P∗M∗,forthe  case  ofFBE22.Then the eigenvalues of positive fixed pointx∗,y∗ are(28)λ1δF=-1,λ2δF=3-δFN∗.The condition λ2δF≠1 leads to(29)δFN∗≠2,4.Using the transformationx~=x-x∗, y~=y-y∗ and writing Aδ=Jx∗,y∗, we shift the fixed point x∗,y∗ of system (4) to the origin. After Taylor expansion, system (4) reduces to(30)x~y~⟶Aδx~y~+F1x~,y~,δF2x~,y~,δ,whereX=x~,y~T and(31)F1x~,y~,δ=16-6aδx∗2x∗+y∗4y~3-6aδy∗2x∗+y∗4x~3+2aδ2x∗-y∗y∗x∗+y∗43x~2y~-2aδx∗-2y∗x∗x∗+y∗43x~y~2+122aδx∗2x∗+y∗3y~2+-2δ+2aδy∗2x∗+y∗3x~2-2aδx∗y∗x∗+y∗32x~y~+OX4,F2x~,y~,δ=166bδx∗2x∗+y∗4y~3+6bδy∗2x∗+y∗4x~3-2bδ2x∗-y∗y∗x∗+y∗43x~2y~+2bδx∗-2y∗x∗x∗+y∗43x~y~2+12-2bδx∗2x∗+y∗3y~2-2bδy∗2x∗+y∗3x~2+2bδx∗y∗x∗+y∗32x~y~+OX4.It follows that(32)B1x,y=∑j,k=12δ2F1ξ,δδξjδξkξ=0xjyk=2aδx∗2x∗+y∗3x2y2+-2δ+2aδy∗2x∗+y∗3x1y1-2aδx∗y∗x∗+y∗3x1y2+x2y1,(33)B2x,y=∑j,k=12δ2F2ξ,δδξjδξkξ=0xjyk=-2bδx∗2x∗+y∗3x2y2-2bδy∗2x∗+y∗3x1y1+2bδx∗y∗x∗+y∗3x1y2+x2y1,(34)C1x,y,u=∑j,k,l=12δ2F1ξ,δδξjδξkδξlξ=0xjykul=-6aδx∗2x∗+y∗4x2y2u2-6aδy∗2x∗+y∗4x1y1u1+2aδ2x∗-y∗y∗x∗+y∗4x1y1u2+x1y2u1+x2y1u1+x2y1u2+x1y2u2+x2y2u1,(35)C2x,y,u=∑j,k,l=12δ2F2ξ,δδξjδξkδξlξ=0xjykul=6bδx∗2x∗+y∗4x2y2u2+6bδy∗2x∗+y∗4x1y1u1-2bδx∗-2y∗x∗x∗+y∗4x1y1u2+x1y2u1+x2y1u1+x2y1u2+x1y2u2+x2y2u1,and δ=δF.Therefore, we get the following symmetric multilinear vector functions ofx,y,u∈R2: Bx,y=B1x,yB2x,y and Cx,y,u=C1x,y,uC2x,y,u.Letp,q∈R2 be two eigenvectors (left and right) of A associated with eigenvalue λ1δF=-1, respectively. Then AδFq=-q and ATδFp=-p. After some algebraic calculation, we get(36)q~-2-δFb2,δFa2T,p~-2-δFb2,-δFb1T.We setp=γ1-2-δFb2,-δFb1T, where (37)γ1=1-2-δFb22-δF2a2b1.Then by the standard scalar product inR2 defined by p,q=p1q1+p2q2, it is obvious to see that p,q=1. The direction of the flip bifurcation can be obtained by the sign of cδF, the coefficient of critical normal form [25], and is computed via(38)cδF=16p,Cq,q,q-12p,Bq,A-I-1Bq,q.The above discussion leads to the following result.Theorem 5.
If (29) holds, cδF≠0, and δ varies around the set FBE21, then system (4) passes through a flip bifurcation at fixed point E2x∗,y∗. Moreover, the period-2 orbits that bifurcate from E2x∗,y∗ are stable (resp., unstable) if cδF>0 (resp., cδF<0).We next discuss a Neimark-Sacker bifurcation of system (4).We first take parametera,b,d,δ arbitrarily from NSBE2 and consider system (4) at fixed point E2x∗,y∗.Since the parameters belong toNSBE2, the complex eigenvalues of (6) are(39)λ,λ¯=tr⁡J2±i24det⁡J-tr⁡J2,where the condition tr⁡J2-4det⁡J<0 leads to a1-b22-4a2b1<0. Let(40)δ=δNS=N∗M∗;then det⁡JδNS=1 and(41)λ=det⁡Jδ,dλδdδδ=δNS=N∗2≠0.Moreover, iftr⁡JδNS≠0,-1, then δNSN∗≠2,3 which leads to(42)λkδNS≠1fork=1,2,3,4.Letq,p∈C2 be eigenvectors of AδNS and ATδNS corresponding to the eigenvalues λδNS and λ¯δNS such that (43)AδNSq=λδNSq,AδNSq¯=λ¯δNSq¯,ATδNSp=λ¯δNSp,ATδNSp¯=λδNSp¯.Then by algebraic calculation, we obtain(44)q~1+δNSb2-λ,-δNSa2T,p~1+δNSb2-λ¯,δNSb1T.We setp=γ21+δNSb2-λ¯,δNSb1T, where (45)γ2=11+δNSb2-λ¯2-δNS2a2b1.Then it is obvious thatp,q=1, where p,q=p1¯q2+p2¯q1 for p,q∈C2. We can decompose vector X∈R2 as X=zq+z¯q¯, for δ close to δNS and z∈C. Obviously, z=p,X. Thus, we obtain the following transformed form of system (30) for δ near δNS: (46)z⟼λδz+gz,z¯,δ,where λδ=1+φδeiθδ with φδNS=0 and gz,z¯,δ is a smooth complex-valued function. Taylor expression of g with respect to z,z¯ yields (47)gz,z¯,δ=∑k+l≥21k!l!gklδzkz¯l,withgkl∈C,k,l=0,1,….According to multilinear symmetric vector functions, we can express the coefficientsgkl as follows:(48)g20δNS=p,Bq,q,g11δNS=p,Bq,q¯,g02δNS=p,Bq¯,q¯,g21δNS=p,Cq,q,q¯.We calculate the coefficientaδNS by the following formula to determine the direction in which invariant closed curve appears.(49)aδNS=Ree-iθδNSg212-Re1-2eiθδNSe-2iθδNS21-eiθδNSg20g11-12g112-14g022, where eiθδNS=λδNS.It is clear that two conditions (41) and (42) known as transversal and nondegenerate hold well for system (4). Therefore, we have the following result.Theorem 6.
IfaδNS≠0 and δ changes its value around NSBE2, then system (4) passes through a Neimark-Sacker bifurcation at fixed point E2. Moreover, if aδNS<0 (resp., >0), then there exists an attracting (resp., repelling) invariant closed curve which bifurcates from E2.


## 4. Numerical Simulations

Here, diagrams for bifurcation, phase portraits, Lyapunov exponents, and fractal dimension of system (4) will be presented by performing numerical simulation to validate our analytical results. The parameter δ is chosen as a bifurcation parameter to see the effect of ratio-dependent functional response in the system. We restrict our attention in the case of NS bifurcation only.Now we consider bifurcation parameters in the following cases.For case (i): by taking parameters as in Table1, we see that, at the fixed point 1,0, a flip bifurcation occurs at δ=δ1=2 and it shows cδ1=12.96>0. This confirms Theorem 4.Table 1
Set of parameter values.
Varying parameter
Fixed parameters
Case (i)
1.8
≤
δ
≤
3
a
=
1.0, b=0.5, d=0.6
Case (ii)
2.05
≤
δ
≤
3.2
a
=
1.67, b=1.0, d=0.6
Case (iii)
2.05
≤
δ
≤
3.2, 1.6≤a≤1.7
b
=
1.0, d=0.6Figure2(a) displays the bifurcation diagram of system (4). This illustrates that stability of fixed point E1 happens for δ<2, bifurcation occurs at δ=2, and there appears period doubling bifurcation if δ>2. Figure 2(b) results in the maximum Lyapunov exponents relating bifurcation in Figure 2(a) confirming the parametric space in which the behaviors of system (4) are chaotic or stable period window.Figure 2
Flip bifurcation and Lyapunov exponent of system (4). (a) Bifurcation for prey and (b) maximum Lyapunov exponents related to (a). Initial value: x0,y0=0.98,0.1.
(a)
(b)For case (ii): by taking parameters as in Table1, we see that, at the fixed point 0.332,0.221333, a NS bifurcation emerges at δ=δNS=2.14859 and it shows aδNS=-0.630982<0. This confirms Theorem 6.Figures3(a) and 3(b) display the bifurcation diagrams of system (4). This illustrates that stability of fixed point E2 happens for δ<2.14859, bifurcation occurs at δ=2.14859, and there appears attracting invariant circle if δ>2.14859. Figure 3(c) results in the maximum Lyapunov exponents relating bifurcation in Figures 3(a) and 3(b) confirming the parametric space in which the behaviors of system (4) are chaotic or stable period window. Figure 3(d) is local amplification of Figure 3(a) for δ∈2.6503,2.999.Figure 3
Neimark-Sacker bifurcation and Lyapunov exponent of system (4). (a) Bifurcation for prey, (b) bifurcation for predator, (c) maximum Lyapunov exponents related to (a-b), and (d) local amplification related to (a) for δ∈2.6503,2.999. Initial value: x0,y0=0.324,0.2163.
(a)
(b)
(c)
(d)The phase portraits of bifurcation diagrams related to Figures3(a) and 3(b) for various values of δ are plotted in Figure 4, which clearly illustrates the act of smooth invariant circle and how it bifurcates from the stable fixed point. We observe that as δ increases, the invariant circle suddenly disappears and periods 9, 26, and 19 orbits appear at δ=2.7,2.902,2.93, respectively. We also see that for δ=3.2, there exists a fully developed chaos in system (4).Figure 4
Phase portraits of bifurcation diagrams Figures3(a) and 3(b) for several values of δ.For case (iii): the dynamics of map (4) can change when more parameters vary. The diagrams of bifurcation of map (4) for control parameters δ∈2.05,3.2, a∈1.6,1.7, respectively, and the fixing remaining parameters as in case (ii) are disposed in Figures 5(a) and 5(b). The 3D maximum Lyapunov exponent for two control parameters δ and a is plotted in Figure 5(c) and its 2D projection onto δ,a is shown in Figure 5(d). It is easy to find values of control parameters for which the dynamics of map (4) is in chaotic, periodic, or nonchaotic status. Obviously, there are chaotic dynamics for δ=3.1492 and a=1.67 and nonchaotic dynamics for δ=2.7 and a=1.67 (see Figure 4), which are compatible with the signs of Lyapunov exponents in Figure 5(c). This exhibits that as the parameter a increases, the map (4) shows more chaotic dynamics.Figure 5
Numerical diagnostic of system (4) for control parameters δ and a. (a-b) Bifurcation for prey and predator covering δ∈2.05,3.2, a∈1.6,1.7, and b=1.0, d=0.6 (c) The 3D view of maximum Lyapunov exponents related to (a-b). (d) The 2D projection onto δ,a. Initial value: x0,y0=0.324,0.2163.
(a)
(b)
(c)
(d)


### 4.1. Initial Perturbation

Sensitivity to initial values of a system demonstrates that at the beginning the two trajectories are arbitrarily closely overlapped but a significant difference builds up for future trajectories. For numerical simulation, two perturbed trajectories for state variablex of system (4) against the number of iterations with initial points x0,y0 and x0+0.001,y0 are generated and plotted in Figure 6(a). The difference (error) between two trajectories is computed by Δx=x1-x2 and presented in Figure 6(b), where x1 and x2 are solution trajectories of system (4) without and with initial perturbation. We take 1000 iterations for each initial value. From Figure 6, it is clear that the error is negligible initially but it increases for trajectories in future time. Therefore, a slightly initial perturbation may lead to complex dynamic behavior of the system.Figure 6
Initial perturbation for system (4). (a) Trajectories for state variable x, blue for x0,y0; red for x0+0.001,y0. (b) Error diagnostic of the two trajectories. Parameters: a=1.67,b=1.0,d=0.6,δ=2.98. Initial value: x0,y0=0.324,0.21633.
(a)
(b)


### 4.2. Fractal Dimension

The fractal dimension (FD) [29, 30] which characterizes strange attractors of a map is defined by (50)dL=j+∑i=1jλiλj,where λ1,λ2,…,λn are Lyapunov exponents and j is the largest integer such that ∑i=1jλi‍≥0 and ∑i=1j+1λi‍<0. For our two-dimensional map (4), the fractal dimension takes the form (51)dL=1+λ1λ2,λ1>0>λ2.Taking parameter values given as in case (i) and by computer simulation, two Lyapunov exponents areλ1≈0.0281 and λ2≈-0.0741 for δ=2.9815 and λ1≈0.0832 and λ2≈-0.0950 for δ=3.1494. So the fractal dimensions of map (4) are(52)dL≈1+0.02810.0741=1.3793forδ=2.9815,dL≈1+0.08320.0950=1.8759forδ=3.1494. The strange attractors given in Figure 4 are consistent with the delta values’ corresponding fractal dimension and this illustrates that the ratio-dependent discrete predator-prey system (4) has a very complex dynamic behavior as the parameter δ increases.  We refer to Figure 7 to understand the system dynamics using fractal dimension.Figure 7
Fractal dimension of system (4). (a) FD related to Figure 2(a) and (b) FD related to Figure 3(a).
(a)
(b)


## 4.1. Initial Perturbation

Sensitivity to initial values of a system demonstrates that at the beginning the two trajectories are arbitrarily closely overlapped but a significant difference builds up for future trajectories. For numerical simulation, two perturbed trajectories for state variablex of system (4) against the number of iterations with initial points x0,y0 and x0+0.001,y0 are generated and plotted in Figure 6(a). The difference (error) between two trajectories is computed by Δx=x1-x2 and presented in Figure 6(b), where x1 and x2 are solution trajectories of system (4) without and with initial perturbation. We take 1000 iterations for each initial value. From Figure 6, it is clear that the error is negligible initially but it increases for trajectories in future time. Therefore, a slightly initial perturbation may lead to complex dynamic behavior of the system.Figure 6
Initial perturbation for system (4). (a) Trajectories for state variable x, blue for x0,y0; red for x0+0.001,y0. (b) Error diagnostic of the two trajectories. Parameters: a=1.67,b=1.0,d=0.6,δ=2.98. Initial value: x0,y0=0.324,0.21633.
(a)
(b)


## 4.2. Fractal Dimension

The fractal dimension (FD) [29, 30] which characterizes strange attractors of a map is defined by (50)dL=j+∑i=1jλiλj,where λ1,λ2,…,λn are Lyapunov exponents and j is the largest integer such that ∑i=1jλi‍≥0 and ∑i=1j+1λi‍<0. For our two-dimensional map (4), the fractal dimension takes the form (51)dL=1+λ1λ2,λ1>0>λ2.Taking parameter values given as in case (i) and by computer simulation, two Lyapunov exponents areλ1≈0.0281 and λ2≈-0.0741 for δ=2.9815 and λ1≈0.0832 and λ2≈-0.0950 for δ=3.1494. So the fractal dimensions of map (4) are(52)dL≈1+0.02810.0741=1.3793forδ=2.9815,dL≈1+0.08320.0950=1.8759forδ=3.1494. The strange attractors given in Figure 4 are consistent with the delta values’ corresponding fractal dimension and this illustrates that the ratio-dependent discrete predator-prey system (4) has a very complex dynamic behavior as the parameter δ increases.  We refer to Figure 7 to understand the system dynamics using fractal dimension.Figure 7
Fractal dimension of system (4). (a) FD related to Figure 2(a) and (b) FD related to Figure 3(a).
(a)
(b)


## 5. Chaos Control

We shall apply a state feedback control method [24, 31] in order to stabilize chaos of system (4) at the state of unstable trajectories. By adding a feedback control law as the control force un to system (4), the controlled form of (4) becomes(53)xn+1=xn+δxn1-xn-axnynyn+xn+un,yn+1=yn+δ-dyn+bxnynyn+xn,(54)un=-k1xn-x∗-k2yn-y∗,where the feedback gains are denoted by k1 and k2 and x∗,y∗ represent positive fixed point of (4).The Jacobian matrixJc of the controlled system (53) is given by(55)Jcx∗,y∗=a11-k1a12-k2a21a22,where a11=1+δa1, a12=-δb1, a21=δa2, a22=1+δb2, and a1, b1, a2, and b2 are determined by (7). The characteristic equation of matrix Jc is(56)λ2-tr⁡Jcλ+det⁡Jc=0,where trJc=a11+a22-k1 and detJc=a22a11-k1-a21a12-k2. Let λ1andλ2 be the roots of (56), that is, eigenvalues of (55). Then(57)λ1+λ2=tr⁡Jc,(58)λ1λ2=det⁡Jc.The solution of the equationsλ1=±1andλ1λ2=1 determines the lines of marginal stability. These conditions confirm that λ1,2<1. Assume that λ1λ2=1; then from (58) we have det⁡Jc=1; that is, (59)l1:a22k1-a21k2=a11a22-a12a21-1.Now, suppose thatλ1=1; then from (57) and (58), we get tr⁡Jc-det⁡Jc=1; that is, (60)l2:1-a22k1+a21k2=a11+a22-1-a11a22+a12a21.Next, assume thatλ1=-1; then from (57) and (58), we obtain tr⁡Jc+det⁡Jc=-1; that is, (61)l3:1+a22k1-a21k2=a11+a22+1+a11a22-a12a21.Then the lines l1, l2, and l3 (see Figure 8(a)) in the k1,k2 plane determine a triangular region which keeps stable eigenvalues.Figure 8
Control of chaotic trajectories of controlled system (53). (a) Bounded region for stable eigenvalues and (b-c) the time series for states x and y.
(a)
(b)
(c)In order to observe how feedback method works and controls chaos at unstable state, we have presented some numerical simulations. We setδ=3.15 and the rest as in case(i) as fixed parameters. The initial value is x0,y0=0.324,0.21633, and the feedback gains are k1=2.38andk2=-1.23. Figures 8(b) and 8(c) illustrate that, at the fixed point 0.332,0.221333, the chaotic trajectory is stabilized.


## 6. Discussions

We investigated the dynamic complexities of discrete predator-prey system with ratio-dependent functional response (4) in the closed first quadrant of R+2. We established the existence conditions for a flip bifurcation and a NS bifurcation of system (4) by using center manifold theorem and bifurcation theory. Specifically, we showed that system (4) displays the change of complex dynamical behaviors via flip bifurcation and NS bifurcation including an invariant cycle, orbits of periods 9, 19, and 26, and attracting chaotic sets when parameter δ changes its value around sets FBE1 and NSBE2, respectively. This indicates that, at the state of chaos, the system is unstable. In particular, it is observed that if the dynamics of prey is chaotic, then the dynamics of predator will ultimately tend to zero or tend to unstable periodic orbit. This results in the interaction between predator and prey which can approach oscillatory balance behavior. Moreover, the variation of two control parameters which displayed captivating chaotic dynamics exists in the system and one can directly observe the chaotic phenomenon from two-dimensional parameter-spaces. Finally, we applied a strategy of feedback control to stabilize chaos at the state of unstable trajectories resulting in predictable stable orbits.



---

*Source: 4537450-2017-07-10.xml*