# Bifurcation Analysis and Chaos Control in a Discrete-Time Predator-Prey System of Leslie Type with Simplified Holling Type IV Functional Response

Authors: S. M. Sohel Rana; Umme Kulsum
Journal: Discrete Dynamics in Nature and Society
 (2017)

---

## Abstract

The dynamic behavior of a discrete-time predator-prey system of Leslie type with simplified Holling type IV functional response is examined. We algebraically show that the system undergoes a bifurcation (flip or Neimark-Sacker) in the interior ofR+2. Numerical simulations are presented not only to validate analytical results but also to show chaotic behaviors which include bifurcations, phase portraits, period 2, 4, 6, 8, 10, and 20 orbits, invariant closed cycle, and attracting chaotic sets. Furthermore, we compute numerically maximum Lyapunov exponents and fractal dimension to justify the chaotic behaviors of the system. Finally, a strategy of feedback control is applied to stabilize chaos existing in the system.


---

## Body

## 1. Introduction

In ecology and mathematical biology, the dynamics of predator-prey interaction when the growth of predator depends on the ratio of predator and prey have been extensively investigated by ecologist and mathematician [1, 2] and the reference therein. Qualitative analyses of these works found many rich dynamics which include limit cycle, states of stability, codimension 1 subcritical Hopf bifurcation, codimension 2 Bogdanov-Takens bifurcation, and codimension 3 degenerate focus type Bogdanov-Takens bifurcation around positive equilibrium. But lots of exploratory works have suggested that the discretization of predator-prey models is more suitable compared to continuous ones when size of populations is small [3–10]. These researches have mainly focused on Gauss-type predator-prey interaction with monotonic functional responses. They obtained many complex properties in discrete-time models including the possibility of bifurcations (flip and Neimark-Sacker) and chaos phenomenon which had been derived either by using numerical simulations or by using center manifold theory.In recent times, a few number of articles in literature discussed the dynamics of discrete-time predator-prey systems of Leslie type [11, 12]. For example, a discrete-time predator-prey system of Holling and Leslie type with constant-yield prey harvesting was investigated in [11], and a discrete predator-prey model with Holling-Tanner functional response was studied in [12]. In their studies, the authors paid their attention to drive the existence of flip and Neimark-Sacker bifurcations by using center manifold theory.In this paper, we consider the following predator-prey system of Leslie type [1]:(1)x˙=rx1-xK-mxyx2+b,y˙=sy1-yhx,where x and y stand for densities of prey and predator, respectively; r,K,m,b,s, and h are all positive constants, and mxy/x2+b denotes the functional response of simplified Holling type IV. In [1], it is shown that model (1) has a degenerate Bogdanov-Takens bifurcation of codimension 3 at positive equilibrium. For simplicity, we scale the variables and parameters as (2)xK⟶x,myrK2⟶y,rt⟶t,a=bK2,α=sr,β=sKhm.Then we write system (1) in the form(3)x˙=x1-x-xyx2+a,y˙=yα-βyx.Forward Euler scheme is applied to system (3) to get the following discrete system:(4)xy⟼x+δx1-x-yx2+ay+δyα-βyx,where δ is the step size. The objective is to study systematically the existence condition of a bifurcation (flip or NS bifurcation) in the interior of R+2 using bifurcation theory and center manifold theory [13–16].This paper is organized as follows. Section2 deals with the existence condition for fixed points of system (4) and their stability criterion. In Section 3, we prove that under certain parametric condition system (4) admits a bifurcation. In Section 4, we implement numerical simulations of the system for one or more control parameters which include diagrams of bifurcation, phase portraits, maximum Lyapunov exponents, and fractal dimensions. In Section 5, we use the method of feedback control to stabilize chaos at unstable trajectories. Finally we carry out a short discussion in Section 6.


## 2. Fixed Points: Existence and Their Stability

The fixed points of system (4) are solution of (5)x+δx1-x-yx2+a=x,y+δyα-βyx=y.A simple algebraic computation shows that system (4) has a predator free fixed point E11,0 for all parameter values. Next, it is of great interest to find the positive fixed point of system (4). Suppose that E2x∗,y∗ is a positive fixed point of system (4). Then, x∗ and y∗ are solutions of (6)1-x∗=y∗x∗2+a,αx∗=βy∗.From (6), we can see that x∗∈0,1 is the root of the following cubic equation: (7)p0w3+3p1w2+3p2w+p3=0,p0≠0,with coefficients (8)p0=1a,3p1=-1a,3p2=1+αaβ,p3=-1.Now, the transformation z=p0w+p1 converts (7) to z3+3Hz+G=0, where G=p02p3-3p0p1p2+2p13, H=p0p2-p12. Using Cardano’s method, we get following result.Lemma 1.
IfG2+4H3>0, then a unique positive fixed point E2x∗,y∗ of system (4) exists, where (9)x∗=1p0q-Hq-p1,y∗=αx∗β,and q denotes one of the three values of -G+G2+4H3/21/3.Next, we investigate stability of system (4) at fixed points. Note that the magnitude of eigenvalues of Jacobian matrix evaluated at fixed point Ex,y determines the local stability of that fixed point. The Jacobian matrix of system (4) evaluated at fixed point Ex,y is given by(10)Jx,y=a11a12a21a22,where(11)a11=1+δ-2δx+δyx2-ax2+a2,a12=-δxx2+a,a21=βδy2x2,a22=1+αδ-2βδyx.The characteristic equation of matrix J is(12)λ2+px,yλ+qx,y=0,where px,y=-trJ=-a11+a22 and det⁡J=a11a22-a12a21.Using Jury’s criterion [17], we state the following stability conditions of fixed points.Proposition 2.
The predator free fixed pointE11,0 is a saddle if 0<δ<2, source if δ>2, and nonhyperbolic if δ=2.It is obvious that whenδ=2, the two eigenvalues of JE1 are λ1=1-δ and λ2=1+αδ>1. Therefore, a flip bifurcation can occur if parameters change in small vicinity of FBE1: (13)FBE1=a,α,β,δ∈0,+∞:δ=2.At E2x∗,y∗, we write (12) as (14)Fλ≔λ2-2+Δδλ+1+Δδ+Ωδ2=0,where (15)Δ=1-2x∗+y∗x∗2-ax∗2+a2+α-2βy∗x∗,Ω=1-2x∗+y∗x∗2-ax∗2+a2α+2βy∗2x∗-1x∗-βy∗2x∗2-3ax∗x∗2+a2.Then F1=Ωδ2>0 and F-1=4+2Δδ+Ωδ2.We state the following Proposition about stability criterion ofE2.Proposition 3.
Suppose that fixed pointE2x∗,y∗ of system (4) exists. Then it is (i)
sink if one of the following conditions holds:
(i.1)
Δ2-4Ω≥0 and δ<-Δ-Δ2-4Ω/Ω;(i.2)
Δ2-4Ω<0 and δ<-Δ/Ω;(ii)
source if one of the following conditions holds:
(ii.1)
Δ2-4Ω≥0 and δ>-Δ+Δ2-4Ω/Ω;(ii.2)
Δ2-4Ω<0 and δ>-Δ/Ω;(iii)
nonhyperbolic if one of the following conditions holds:
(iii.1)
Δ2-4Ω≥0 and δ=-Δ±Δ2-4Ω/Ω;(iii.2)
Δ2-4Ω<0 and δ=-Δ/Ω;(iv)
saddle if otherwise.From Proposition3, we see that two eigenvalues of JE2 are λ1=-1 and λ2≠-1,1 if condition (iii.1) holds. We rewrite term (iii.1) as follows: (16)FBE21=a,α,β,δ∈0,+∞:δ=-Δ-Δ2-4ΩΩ,Δ2-4Ω≥0,or (17)FBE22=a,α,β,δ∈0,+∞:δ=-Δ+Δ2-4ΩΩ,Δ2-4Ω≥0.Therefore, a flip bifurcation can appear at E2 if parameters vary around the set FBE21 or FBE22.Also we rewrite term (iii.2) as follows:(18)NSBE2=a,α,β,δ∈0,+∞:δ=-ΔΩ,Δ2-4Ω<0,and if the parameters change in small vicinity of NSBE2, two eigenvalues λ1,2 of JE2 are complex having magnitude one and then NS bifurcation can emerge from fixed point E2.


## 3. Bifurcation Analysis

In this section, we will give attention to recapitulate bifurcations (flip and Neimark-Sacker) of system (4) around E2 by using the theory of bifurcation [13–16]. We set δ as a bifurcation parameter.


### 3.1. Flip Bifurcation

Consider system (4) at the fixed point E2x∗,y∗ with arbitrary parameter a,α,β,δ∈FBE21. One can consider the case of FBE22 in a similar fashion. Since the parameters lie in FBE21, let δ=δF=-Δ-Δ2-4Ω/Ω, and then the eigenvalues of positive fixed point x∗,y∗ are(19)λ1δF=-1,λ2δF=3+ΔδF.The condition λ2δF≠1 leads to(20)ΔδF≠-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(21)x~y~⟶Aδx~y~+F1x~,y~,δF2x~,y~,δ,where X=x~,y~T and(22)F1x~,y~,δ=166δy∗a2-6ax∗2+x∗4x∗2+a4x~3+3δ6ax∗-2x∗3x∗2+a3x~2y~+12-2δ+δy∗6ax∗-2x∗3x∗2+a3x~2+2δx∗2-ax∗2+a2x~y~+OX4,F2x~,y~,δ=βδx∗4x~x∗y~-y∗x~2-βδx∗3x~x∗y~-y∗x~2+OX4.It follows that(23)B1x,y=∑j,k=12δ2F1ξ,δδξjδξkξ=0xjyk=δx∗2-ax∗2+a2x1y2+x2y1+-2δ+δy∗6ax∗-2x∗3x∗2+a3x1y1,B2x,y=∑j,k=12δ2F2ξ,δδξjδξkξ=0xjyk=-2βδx∗x2y2+2βδy∗x∗2x1y2+x2y1-2βδy∗2x∗3x1y1,C1x,y,u=∑j,k,l=12δ2F1ξ,δδξjδξkδξlξ=0xjykul=δ6ax∗-2x∗3x∗2+a3x1y2u1+x2y1u1+x1y1u2+6δy∗a2-6ax∗2+x∗4x∗2+a4x1y1u1,C2x,y,u=∑j,k,l=12δ2F2ξ,δδξjδξkδξlξ=0xjykul=2βδx∗2x1y2u2+x2y1u2+x2y2u1-4βδy∗x∗3x1y1u2+x1y2u1+x2y1u1+6βδy∗2x∗4x1y1u1,and δ=δF.Therefore, we obtain the following symmetric multilinear vector functions ofx,y,u∈R2:(24)Bx,y=B1x,yB2x,y,Cx,y,u=C1x,y,uC2x,y,u.Letp,q∈R2 be two eigenvectors of A for eigenvalue λ1δF=-1 such that AδFq=-q and ATδFp=-p. Then by direct calculation we get(25)q~2+αδF-2βδFy∗x∗,-βδFy∗2x∗2T,p~2+αδF-2βδFy∗x∗,δFx∗x∗2+aT.We set p=γ12+αδF-2βδFy∗/x∗,δFx∗/x∗2+aT, where (26)γ1=12+αδF-2βδFy∗/x∗2-βδF2y∗2/x∗x∗2+a.Then by the standard scalar product inR2 defined by p,q=p1q1+p2q2, we show that p,q=1. The direction of the flip bifurcation is obtained by sign cδF, the coefficient of critical normal form [13], and is given by(27)cδF=16p,Cq,q,q-12p,Bq,A-I-1Bq,q.We state the following result on flip bifurcation according to the above analysis.Theorem 4.
If (20) holds, cδF≠0, and the parameter δ changes its value around δF, then system (4) undergoes a flip bifurcation at positive 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).


### 3.2. Neimark-Sacker Bifurcation

We consider system (4) at fixed point E2x∗,y∗ with arbitrary a,α,β,δ∈NSBE2. From (12), the eigenvalues are given by (28)λ,λ¯=-pδ±pδ2-4qδ2.Since the parameters belong toNSBE2, the eigenvalues will be complex and (29)λ,λ¯=1+Δδ2±iδ24Ω-Δ2.Let(30)δ=δNS=-ΔΩ.Therefore, we have(31)λ=qδ,qδNS=1,dλδdδδ=δNS=-Δ2≠0.Moreover, ifpδNS≠0,1, then(32)Δ2Ω≠2,3,which obviously satisfies(33)λkδNS≠1fork=1,2,3,4.Suppose thatq,p∈C2 are two eigenvectors of AδNS and ATδNS for eigenvalues λδNS and λ¯δNS such that (34)AδNSq=λδNSq,AδNSq¯=λ¯δNSq¯,ATδNSp=λ¯δNSp,ATδNSp¯=λδNSp¯.Then by direct computation we obtain(35)q~1+αδNS-2βδy∗x∗-λ,-βδy∗2x∗2T,p~1+αδNS-2βδy∗x∗-λ¯,δx∗x∗2+aT.We set p=γ21+αδNS-2βδy∗/x∗-λ¯,δx∗/x∗2+aT, where (36)γ2=11+αδNS-2βδy∗/x∗-λ¯2-βδNS2y∗2/x∗x∗2+a.Then it is clear thatp,q=1 where p,q=p1¯q2+p2¯q1 for p,q∈C2. Now, we 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 (21) for δ near δNS: (37)z⟼λδz+gz,z¯,δ,where λδ=1+φδeiθδ with φδNS=0 and gz,z¯,δ is a smooth complex-valued function. After Taylor expression of g with respect to z,z¯, we obtain (38)gz,z¯,δ=∑k+l≥21k!l!gklδzkz¯l,withgkl∈C,k,l=0,1,….According to multilinear symmetric vector functions, the coefficients gkl are(39)g20δNS=p,Bq,q,g11δNS=p,Bq,q¯,g02δNS=p,Bq¯,q¯,g21δNS=p,Cq,q,q¯.The invariant closed curve appear in the direction which is determined by the coefficient aδNS and calculated via(40)aδNS=Ree-iθδNSg212-Re1-2eiθδNSe-2iθδNS21-eiθδNSg20g11-12g112-14g022,where eiθδNS=λδNS.It is clear that conditions (31) and (33) known as transversal and nondegenerate for system (4) hold well. We obtain the following result.Theorem 5.
If (32) holds, aδNS≠0, and the parameter δ changes its value in small vicinity of NSBE2, then system (4) passes through a Neimark-Sacker bifurcation at positive fixed point E2. Moreover, if aδNS<0 (resp., >0), then there exists a unique attracting (resp., repelling) invariant closed curve which bifurcates from E2.


## 3.1. Flip Bifurcation

Consider system (4) at the fixed point E2x∗,y∗ with arbitrary parameter a,α,β,δ∈FBE21. One can consider the case of FBE22 in a similar fashion. Since the parameters lie in FBE21, let δ=δF=-Δ-Δ2-4Ω/Ω, and then the eigenvalues of positive fixed point x∗,y∗ are(19)λ1δF=-1,λ2δF=3+ΔδF.The condition λ2δF≠1 leads to(20)ΔδF≠-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(21)x~y~⟶Aδx~y~+F1x~,y~,δF2x~,y~,δ,where X=x~,y~T and(22)F1x~,y~,δ=166δy∗a2-6ax∗2+x∗4x∗2+a4x~3+3δ6ax∗-2x∗3x∗2+a3x~2y~+12-2δ+δy∗6ax∗-2x∗3x∗2+a3x~2+2δx∗2-ax∗2+a2x~y~+OX4,F2x~,y~,δ=βδx∗4x~x∗y~-y∗x~2-βδx∗3x~x∗y~-y∗x~2+OX4.It follows that(23)B1x,y=∑j,k=12δ2F1ξ,δδξjδξkξ=0xjyk=δx∗2-ax∗2+a2x1y2+x2y1+-2δ+δy∗6ax∗-2x∗3x∗2+a3x1y1,B2x,y=∑j,k=12δ2F2ξ,δδξjδξkξ=0xjyk=-2βδx∗x2y2+2βδy∗x∗2x1y2+x2y1-2βδy∗2x∗3x1y1,C1x,y,u=∑j,k,l=12δ2F1ξ,δδξjδξkδξlξ=0xjykul=δ6ax∗-2x∗3x∗2+a3x1y2u1+x2y1u1+x1y1u2+6δy∗a2-6ax∗2+x∗4x∗2+a4x1y1u1,C2x,y,u=∑j,k,l=12δ2F2ξ,δδξjδξkδξlξ=0xjykul=2βδx∗2x1y2u2+x2y1u2+x2y2u1-4βδy∗x∗3x1y1u2+x1y2u1+x2y1u1+6βδy∗2x∗4x1y1u1,and δ=δF.Therefore, we obtain the following symmetric multilinear vector functions ofx,y,u∈R2:(24)Bx,y=B1x,yB2x,y,Cx,y,u=C1x,y,uC2x,y,u.Letp,q∈R2 be two eigenvectors of A for eigenvalue λ1δF=-1 such that AδFq=-q and ATδFp=-p. Then by direct calculation we get(25)q~2+αδF-2βδFy∗x∗,-βδFy∗2x∗2T,p~2+αδF-2βδFy∗x∗,δFx∗x∗2+aT.We set p=γ12+αδF-2βδFy∗/x∗,δFx∗/x∗2+aT, where (26)γ1=12+αδF-2βδFy∗/x∗2-βδF2y∗2/x∗x∗2+a.Then by the standard scalar product inR2 defined by p,q=p1q1+p2q2, we show that p,q=1. The direction of the flip bifurcation is obtained by sign cδF, the coefficient of critical normal form [13], and is given by(27)cδF=16p,Cq,q,q-12p,Bq,A-I-1Bq,q.We state the following result on flip bifurcation according to the above analysis.Theorem 4.
If (20) holds, cδF≠0, and the parameter δ changes its value around δF, then system (4) undergoes a flip bifurcation at positive 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).


## 3.2. Neimark-Sacker Bifurcation

We consider system (4) at fixed point E2x∗,y∗ with arbitrary a,α,β,δ∈NSBE2. From (12), the eigenvalues are given by (28)λ,λ¯=-pδ±pδ2-4qδ2.Since the parameters belong toNSBE2, the eigenvalues will be complex and (29)λ,λ¯=1+Δδ2±iδ24Ω-Δ2.Let(30)δ=δNS=-ΔΩ.Therefore, we have(31)λ=qδ,qδNS=1,dλδdδδ=δNS=-Δ2≠0.Moreover, ifpδNS≠0,1, then(32)Δ2Ω≠2,3,which obviously satisfies(33)λkδNS≠1fork=1,2,3,4.Suppose thatq,p∈C2 are two eigenvectors of AδNS and ATδNS for eigenvalues λδNS and λ¯δNS such that (34)AδNSq=λδNSq,AδNSq¯=λ¯δNSq¯,ATδNSp=λ¯δNSp,ATδNSp¯=λδNSp¯.Then by direct computation we obtain(35)q~1+αδNS-2βδy∗x∗-λ,-βδy∗2x∗2T,p~1+αδNS-2βδy∗x∗-λ¯,δx∗x∗2+aT.We set p=γ21+αδNS-2βδy∗/x∗-λ¯,δx∗/x∗2+aT, where (36)γ2=11+αδNS-2βδy∗/x∗-λ¯2-βδNS2y∗2/x∗x∗2+a.Then it is clear thatp,q=1 where p,q=p1¯q2+p2¯q1 for p,q∈C2. Now, we 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 (21) for δ near δNS: (37)z⟼λδz+gz,z¯,δ,where λδ=1+φδeiθδ with φδNS=0 and gz,z¯,δ is a smooth complex-valued function. After Taylor expression of g with respect to z,z¯, we obtain (38)gz,z¯,δ=∑k+l≥21k!l!gklδzkz¯l,withgkl∈C,k,l=0,1,….According to multilinear symmetric vector functions, the coefficients gkl are(39)g20δNS=p,Bq,q,g11δNS=p,Bq,q¯,g02δNS=p,Bq¯,q¯,g21δNS=p,Cq,q,q¯.The invariant closed curve appear in the direction which is determined by the coefficient aδNS and calculated via(40)aδNS=Ree-iθδNSg212-Re1-2eiθδNSe-2iθδNS21-eiθδNSg20g11-12g112-14g022,where eiθδNS=λδNS.It is clear that conditions (31) and (33) known as transversal and nondegenerate for system (4) hold well. We obtain the following result.Theorem 5.
If (32) holds, aδNS≠0, and the parameter δ changes its value in small vicinity of NSBE2, then system (4) passes through a Neimark-Sacker bifurcation at positive fixed point E2. Moreover, if aδNS<0 (resp., >0), then there exists a unique attracting (resp., repelling) invariant closed curve which bifurcates from E2.


## 4. Numerical Simulations

Here, diagrams for bifurcation, phase portraits, maximum Lyapunov exponents, and fractal dimension of system (4) will be drawn to validate our theoretical results using numerical simulation. We consider parameter values in the following cases for bifurcation analysis.Case 1.
Varyδ in range 0.75≤δ≤1.17, and fix a=5.5,α=2.5,β=4.0.Case 2.
Varyδ in range 1.9≤δ≤2.15, and fix a=2.5,α=0.75,β=0.5.Case 3.
Varyδ in range 1.9≤δ≤2.15, β in range 0.5≤β≤0.8, and fix a=2.5,α=0.75.For Case1, taking parameters a=5.5,α=2.5, and β=4.0 and δ covering 0.75,1.17, we find that, at fixed point E20.910115,0.568822, a flip bifurcation occurs at δ=δF~0.852505 with multipliers λ1=-1,λ2=0.112921, aδF=22.3878, and a,α,β,δ∈FBE21. This verifies Theorem 4.According to bifurcation diagrams shown in Figures1(a) and 1(b), we see that stability of fixed point E2 happens for δ<0.852505 and loses its stability at δ=0.852505 and period doubling phenomena lead to chaos for δ>0.852505. The maximum Lyapunov exponents and fractal dimension related to Figures 1(a) and 1(b) are computed and shown in Figures 1(c) and 1(d). We observe that the period 2, 4, and 8 orbits occur for δ∈0.75,1.104, chaotic set occurs for δ∈1.104,1.15, and the period 6 orbit occurs at δ=1.112 within the window of chaotic region δ∈1.104,1.15. The status of stable, periodic, or chaotic dynamics is compatible with sign in Figures 1(c) and 1(d).Figure 1
Flip 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) fractal dimension corresponding to (a). Initial value x0,y0=0.9,0.56.
(a)
(b)
(c)
(d)For Case2, taking parameters a=2.5, α=0.75, and β=0.5 and δ covering 1.9,2.15, we observe that a Neimark-Sacker (NS) bifurcation appears at fixed point 0.662032,0.993048 for δ=δNS~1.94431. Also, we have λ,λ¯=-0.274696±0.961531i, g20=3.36329-2.22294i, g11=3.85086-1.04183i, g02=-1.24765+5.18251i, g21=-4.31749+7.22596i, aδNS=-11.9921, and a,α,β,δ∈NSBE2. This verifies Theorem 5.The bifurcation diagrams shown in Figures2(a) and 2(b) demonstrate that stability of E2 happens for δ<1.94431 and loses its stability at δ=1.94431 and an attracting invariant curve appears if δ>1.94431. We dispose the maximum Lyapunov exponents in Figure 2(c) relating bifurcation in Figures 2(a) and 2(b), which confirm the existences of chaos and period window as parameter δ varying. When δ~2.126, the sign of maximum Lyapunov exponent is confirming presence of chaos. Figure 2(d) is local amplification of Figure 2(a) for δ∈2.05,2.11.Figure 2
NS bifurcation and Lyapunov exponent of system (4). (a) NS bifurcation for prey, (b) NS bifurcation for predator, (c) maximum Lyapunov exponents related to (a-b), (d) local amplification diagram in (a) for δ∈2.05,2.11, and (e) fractal dimension associated with (a). Initial value x0,y0=0.65,0.95.
(a)
(b)
(c)
(d)
(e)The phase portraits of bifurcation diagrams in Figures2(a) and 2(b) for different values of δ are displayed in Figure 3, which clearly illustrates the act of smooth invariant curve how it bifurcates from the stable fixed point and increases its radius. As δ grows, disappearance of closed curve occurs suddenly and a period 10 and 20 orbits appear at δ~2.06 and δ~2.093, respectively. We also see that a fully developed chaos in system (4) occurs at δ~2.15.Figure 3
Phase portraits of bifurcation diagrams Figures2(a) and 2(b) for different values of δ.For Case3, the dynamic complexity of system (4) can be observed when more parameters vary. The diagrams of bifurcation of system (4) for control parameters δ∈1.9,2.15 and β∈0.5,0.8 and fixing remaining parameters as in Case 2 are shown in Figures 4(a) and 4(b). The 3D maximum Lyapunov exponent for two control parameters is plotted in Figure 4(c) and its 2D projection onto δ,β plane is shown in Figure 4(d). It is easy to find values of control parameters for which the dynamics of system (4) are in status of nonchaotic, periodic, or chaotic. For instance, there are chaotic dynamics for δ=2.15 and β=0.5 and the nonchaotic dynamics for δ=1.93 and β=0.5 (see Figure 3), which are compatible with the signs of maximum Lyapunov exponents in Figure 4(c). This exhibits that the parameter β may play a role for chaotic dynamics in the system.Figure 4
Diagnostic of system (4) for control parameters
δ
and
β. (a-b) Bifurcation for prey and predator covering δ∈1.9,2.15, β=0.5,0.55,0.6,0.65,0.7,0.75,0.8∈0.5,0.8 in δ-β-x space and δ-β-y space. (c) The 3D view of maximum Lyapunov exponents related to (a-b). (d) The 2D projection onto δ,β plane. Initial value x0,y0=0.65,0.95.
(a)
(b)
(c)
(d)


### 4.1. Fractal Dimension of the Map

The measure of fractal dimensions characterizes the strange attractors of a system. By using Lyapunov exponents, the fractal dimension [18, 19] is defined by (41)dL=j+∑i=1jhihj,where h1,h2,…,hn are Lyapunov exponents and j is the largest integer such that ∑i=1jhi≥0 and ∑i=1j+1hi<0. For our two-dimensional system (4), the fractal dimension takes the form (42)dL=1+h1h2,h1>0>h2.With parameter values as in Case2, the fractal dimension of system (4) is plotted in Figure 2(e). The strange attractors given in Figure 3 and its corresponding fractal dimension illustrate that the Leslie type predator-prey system (4) has a chaotic dynamics as the parameter δ increases.


## 4.1. Fractal Dimension of the Map

The measure of fractal dimensions characterizes the strange attractors of a system. By using Lyapunov exponents, the fractal dimension [18, 19] is defined by (41)dL=j+∑i=1jhihj,where h1,h2,…,hn are Lyapunov exponents and j is the largest integer such that ∑i=1jhi≥0 and ∑i=1j+1hi<0. For our two-dimensional system (4), the fractal dimension takes the form (42)dL=1+h1h2,h1>0>h2.With parameter values as in Case2, the fractal dimension of system (4) is plotted in Figure 2(e). The strange attractors given in Figure 3 and its corresponding fractal dimension illustrate that the Leslie type predator-prey system (4) has a chaotic dynamics as the parameter δ increases.


## 5. Chaos Control

To stabilize chaos at the state of unstable trajectories of system (4), a state feedback control method [17, 20] is applied. By adding a feedback control law as the control force un to system (4), the controlled form of system (4) becomes(43)xn+1=xn+δxn1-xn-ynxn2+a+un,yn+1=yn+δynα-βynxn,un=-k1xn-x∗-k2yn-y∗,where the feedback gains are denoted by k1 and k2 and x∗,y∗ represent positive fixed point of system (4).The Jacobian matrixJc of the controlled system (43) is given by(44)Jcx∗,y∗=a11-k1a12-k2a21a22,where a11=1+δ-2δx+δyx2-a/x2+a2, a12=-δx/x2+a, a21=βδy2/x2, and a22=1+αδ-2βδy/x are evaluated at x∗,y∗. The characteristic equation of (44) is(45)λ2-trJcλ+det⁡Jc=0,where trJc=a11+a22-k1 and det⁡Jc=a22a11-k1-a21a12-k2. Let λ1 and λ2 be the roots of (45). Then(46)λ1+λ2=a11+a22-k1,(47)λ1λ2=a22a11-k1-a21a12-k2.The solution of the equationsλ1=±1 and λ1λ2=1 determines the lines of marginal stability. These conditions confirm that λ1,2<1. Suppose that λ1λ2=1, then from (47) we have (48)l1:a22k1-a21k2=a11a22-a12a21-1.Assume thatλ1=1, then from (46) and (47) we get (49)l2:1-a22k1+a21k2=a11+a22-1-a11a22+a12a21.Next, assume thatλ1=-1, then from (46) and (47) we obtain (50)l3:1+a22k1-a21k2=a11+a22+1+a11a22-a12a21.Then the linesl1,l2, and l3 (see Figure 5(a)) in the k1,k2 plane determine a triangular region which keeps eigenvalues with magnitude less than 1.Figure 5
Control of chaotic trajectories of system (43). (a) Stability region in k1,k2 plane. (b-c) Time series for states x and y, respectively.
(a)
(b)
(c)In order to check how the implementation of feedback control method works and controls chaos at unstable state, we have performed numerical simulations. Parameter values are fixed asδ=2.126 and rest as in Case 2. The initial value is x0,y0=0.65,0.95, and the feedback gains are k1=-1.3 and k2=0.16. Figures 5(b) and 5(c) show that, at the fixed point 0.662032,0.993048, the chaotic trajectory is stabilized.


## 6. Discussions

We investigate complex behaviors of discretized Leslie type predator-prey system (4) with simplified Holling type IV functional response in the closed first quadrant R+2. By using center manifold theorem and bifurcation theory we establish that system (4) can undergo a bifurcation (flip or NS) at unique positive fixed point if δ varies around the sets FBE21 or FBE22 and NSBE2. Numerical simulations present unpredictable behaviors of the system through a flip bifurcation which includes orbits of period 2, period 4, period 6, and period 8 and through a NS bifurcation which includes an invariant cycle, orbits of period 10 and period 20, and chaotic sets, respectively. These indicate that, at the state of chaos, the system is unstable and particularly, the predator goes to extinct or goes to a stable fixed point when the dynamic of prey is chaotic. We confirm about the existence of chaos through the computation of maximum Lyapunov exponents and fractal dimension. Moreover, system (4) exhibits rich and chaotic dynamics by the variation of two control parameters and one can directly observe the chaotic phenomenon from this two-dimensional parameter-space. Finally, the chaotic trajectories at unstable state are controlled by implementing the strategy of feedback control. However, it is still a challenging problem to explore multiple parameter bifurcation in the system. We expect to obtain some more analytical results on this issue in the future.



---

*Source: 9705985-2017-07-25.xml*