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

This study examines the complexity of a discrete-time predator-prey systemwith 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.


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][4][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][7][8][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][11][12][13][14][15][16][17][18][19][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][22][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: where  and  stand for densities of prey and predator, respectively,  > 0 is predator's death rate, and () is prey's specific growth rate.(/) is ratio-dependent response function of predator and ,  > 0 are rate of capturing and conversion, respectively.If () = (1 − /) and () = /( + ), then system (1) takes in the form [6]: where , , , and  are all positive constants.By scaling the variables and parameters as  → /,  → /, and  → , we write system (2) as ẋ =  (1 − ) −   +  , ẏ =  (− +   +  ) with  = /,  = /, and  = / being positive constants.Forward Euler scheme is applied to system (3) to obtain the following system: ) , 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.Section 2 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.

Fixed Points: Existence and Their Stability
The fixed points of (4) are the solution of the following equations: A simple algebraic computation yields the following result.
To show the region in the space (, , ) for which positive fixed point of system (4) exists, we plot  = {(, , ) :  <  < /( − 1) with  > 1} in Figure 1.We see that  2 lies inside C but not in regions (I)-(II).Next, we investigate stability of (4) at fixed points.The Jacobian matrix of system (4) at fixed point (, ) is where Then the characteristic equation of ( 6) is where It is noted that the magnitude of eigenvalues of Jacobian matrix evaluated at fixed point (, ) determines the local stability of that fixed point.The following propositions represent the stability conditions of fixed points by Jury's criterion [24].where From Proposition 3, it follows that if term (iii.1)holds, two eigenvalues of ( 2 ) are  1 = −1 and  2 ̸ = −1, 1.We rewrite the term (iii.1) as follows:

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][26][27][28].The parameter  is chosen as a bifurcation parameter.
We next discuss a Neimark-Sacker bifurcation of system (4).
Since the parameters belong to NSB  2 , the complex eigenvalues of ( 6) are where the condition (tr ) 2 − 4 det  < 0 leads to  Then by algebraic calculation, we obtain We set  =  2 (1 +  NS  2 − ,  NS  1 )  , where Then it is obvious that ⟨, ⟩ = 1, where ⟨, ⟩ =  1  2 +  2  1 for ,  ∈ C 2 .We can decompose vector  ∈ R 2 as  =  + , for  close to  NS and  ∈ C. Obviously,  = ⟨, ⟩.Thus, we obtain the following transformed form of system (30) for || near  NS : where () = (1 + ()) () with ( NS ) = 0 and (, , ) is a smooth complex-valued function.Taylor expression of  with respect to (, ) yields According to multilinear symmetric vector functions, we can express the coefficients   as follows: We calculate the coefficient ( NS ) by the following formula to determine the direction in which invariant closed curve appears.

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.
Figure 2(a) displays the bifurcation diagram of system (4).This illustrates that stability of fixed point  1 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.
For case (ii): by taking parameters as in Table 1, we see that, at the fixed point (0.332, 0.221333), a NS bifurcation emerges at  =  NS = 2.14859 and it shows ( NS ) = −0.630982< 0. This confirms Theorem 6.
Figures 3(a) and 3(b) display the bifurcation diagrams of system (4).This illustrates that stability of fixed point  2 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 The phase portraits of bifurcation diagrams related to Figures 3(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).
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],  ∈ [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  is plotted in Figure 5(c) and its 2D projection onto (, ) 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  = 1.67 and nonchaotic dynamics for  = 2.7 and  = 1.67 (see Figure 4), which are compatible with the signs of Lyapunov exponents in Figure 5(c).This exhibits that as the parameter  increases, the map (4) shows more chaotic dynamics.

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 variable  of system (4) against the number of iterations with initial points ( 0 ,  0 ) and ( 0 + 0.001,  0 ) are generated and plotted in Figure 6(a).The difference (error) between two trajectories is computed by Δ =  1 − 2 and presented in Figure 6(b), where  1 and  2 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.[29,30] which characterizes strange attractors of a map is defined by

Fractal Dimension. The fractal dimension (FD)
where  1 ,  2 , . . .,   are Lyapunov exponents and  is the largest integer such that ∑  =1   ≥ 0 and ∑ +1 =1   < 0. For our two-dimensional map (4), the fractal dimension takes the form 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) 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.

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   to system (4), the controlled form of (4) becomes where the feedback gains are denoted by  1 and  2 and ( * ,  * ) represent positive fixed point of (4).The Jacobian matrix   of the controlled system (53) is given by where  11 = 1 +  ).Let  1 and  2 be the roots of (56), that is, eigenvalues of (55).Then The solution of the equations  1 = ±1 and  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   = 1; that is,      (61) Then the lines  1 ,  2 , and  3 (see Figure 8(a)) in the ( 1 ,  2 ) plane determine a triangular region which keeps stable eigenvalues.
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 ( 0 ,  0 ) = (0.324, 0.21633), and the feedback gains are  1 = 2.38 and  2 = −1.23.Figures 8(b) and 8(c) illustrate that, at the fixed point (0.332, 0.221333), the chaotic trajectory is stabilized.

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

Figure 1 :
Figure 1: Graphical depiction for positive fixed point.

Figure 8 :
Figure 8: Control of chaotic trajectories of controlled system (53).(a) Bounded region for stable eigenvalues and (b-c) the time series for states  and .

Table 1 :
Set of parameter values.