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

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 ofR2+. Numerical simulations are presented not only to validate analytical results but also to show chaotic behaviorswhich 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.


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][4][5][6][7][8][9][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 discretetime 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]: where  and  stand for densities of prey and predator, respectively; , , , , , and ℎ are all positive constants, and /( 2 + ) denotes the functional response of simplified Holling type IV.In [1], it is shown that model (1)  ( Then we write system (1) in the form Forward Euler scheme is applied to system (3) to get the following discrete system: 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][14][15][16].
This paper is organized as follows.Section 2 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.

Fixed Points: Existence and Their Stability
The fixed points of system (4) are solution of A simple algebraic computation shows that system (4) has a predator free fixed point  1 (1, 0) for all parameter values.Next, it is of great interest to find the positive fixed point of system (4).Suppose that  2 ( * ,  * ) is a positive fixed point of system (4).Then,  * and  * are solutions of From ( 6), we can see that  * ∈ (0, 1) is the root of the following cubic equation: with coefficients Now, the transformation  =  0  +  1 converts (7) to . Using Cardano's method, we get following result.

Lemma 1.
If  2 + 4 3 > 0, then a unique positive fixed point  2 ( * ,  * ) of system (4) exists, where and  denotes one of the three values of Next, we investigate stability of system (4) at fixed points.Note that the magnitude of eigenvalues of Jacobian matrix evaluated at fixed point (, ) determines the local stability of that fixed point.The Jacobian matrix of system (4) evaluated at fixed point (, ) is given by where The characteristic equation of matrix  is where (, ) = −tr  = −( It is obvious that when  = 2, the two eigenvalues of ( 1 ) are  1 = 1− and  2 = 1+ > 1.Therefore, a flip bifurcation can occur if parameters change in small vicinity of FB  1 : At  2 ( * ,  * ), we write (12) as where Then (1) = Ω 2 > 0 and (−1) = 4 + 2Δ + Ω 2 .We state the following Proposition about stability criterion of  2 .Proposition 3. Suppose that fixed point  2 ( * ,  * ) of system (4) exists.Then it is (i) sink if one of the following conditions holds: (ii) source if one of the following conditions holds: (iii) nonhyperbolic if one of the following conditions holds: From Proposition 3, we see that two eigenvalues of ( 2 ) are  1 = −1 and  2 ̸ = −1, 1 if condition (iii.1)holds.We rewrite term (iii.1) as follows: or Therefore, a flip bifurcation can appear at  2 if parameters vary around the set FB 1  2 or FB 2  2 .Also we rewrite term (iii.2) as follows: and if the parameters change in small vicinity of NSB  2 , two eigenvalues  1,2 of ( 2 ) are complex having magnitude one and then NS bifurcation can emerge from fixed point  2 .

Bifurcation Analysis
In this section, we will give attention to recapitulate bifurcations (flip and Neimark-Sacker) of system (4) around  2 by using the theory of bifurcation [13][14][15][16].We set  as a bifurcation parameter.
Then by direct calculation we get We set  =  1 (2 +   − 2   * / * ,    * /( * 2 + ))  , where Then by the standard scalar product in R 2 defined by ⟨, ⟩ =  1  1 +  2  2 , we show that ⟨, ⟩ = 1.The direction of the flip bifurcation is obtained by sign (  ), the coefficient of critical normal form [13], and is given by We state the following result on flip bifurcation according to the above analysis.
The invariant closed curve appear in the direction which is determined by the coefficient ( NS ) and calculated via where  ( 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.

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  = 5.5,  = 2.5,  = 4.0.
For Case 2, taking parameters  = 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 The bifurcation diagrams shown in Figures 2(a) and 2(b) demonstrate that stability of  2 happens for  < 1.94431 and loses its stability at  = 1.94431 and an attracting invariant curve appears if  > 1.94431.We dispose the maximum   The phase portraits of bifurcation diagrams in Figures 2(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.
For Case 3, 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.

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 where ℎ 1 , ℎ 2 , . . ., ℎ  are Lyapunov exponents and  is the largest integer such that ∑  =1 ℎ  ≥ 0 and ∑ +1 =1 ℎ  < 0. For our two-dimensional system (4), the fractal dimension takes the form With parameter values as in Case 2, 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.

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   to system (4), the controlled form of system (4) becomes where the feedback gains are denoted by  1 and  2 and ( * ,  * ) represent positive fixed point of system (4).The Jacobian matrix   of the controlled system (43) is given by where Then the lines  1 ,  2 , and  3 (see Figure 5(a)) in the ( 1 ,  2 ) plane determine a triangular region which keeps eigenvalues with magnitude less than 1.
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 ( 0 ,  0 ) = (0.65, 0.95), and the feedback gains are  1 = −1.3 and  2 = 0.16.Figures 5(b) and 5(c) show that, at the fixed point (0.662032, 0.993048), the chaotic trajectory is stabilized.

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 FB 1  2 or FB 2  2 and NSB  2 .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.

Figure 3 :
Figure 3: Phase portraits of bifurcation diagrams Figures 2(a) and 2(b) for different values of .