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 of . 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.
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]:where and stand for densities of prey and predator, respectively; , and are all positive constants, and 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
Then we write system (1) in the formForward 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 using bifurcation theory and center manifold theory [13–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.
2. 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 for all parameter values. Next, it is of great interest to find the positive fixed point of system (4). Suppose that is a positive fixed point of system (4). Then, and are solutions of From (6), we can see that is the root of the following cubic equation: with coefficients Now, the transformation converts (7) to , where , . Using Cardano’s method, we get following result.
Lemma 1. If , then a unique positive fixed point 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 bywhereThe characteristic equation of matrix iswhere and
Using Jury’s criterion [17], we state the following stability conditions of fixed points.
Proposition 2. The predator free fixed point is a saddle if , source if , and nonhyperbolic if .
It is obvious that when , the two eigenvalues of are and . Therefore, a flip bifurcation can occur if parameters change in small vicinity of : At , we write (12) as where Then and
We state the following Proposition about stability criterion of .
Proposition 3. Suppose that fixed point of system (4) exists. Then it is (i)sink if one of the following conditions holds:(i.1) and ;(i.2) and ;(ii)source if one of the following conditions holds:(ii.1) and ;(ii.2) and ;(iii)nonhyperbolic if one of the following conditions holds:(iii.1) and ;(iii.2) and ;(iv)saddle if otherwise.
From Proposition 3, we see that two eigenvalues of are and if condition (iii.1) holds. We rewrite term (iii.1) as follows: or Therefore, a flip bifurcation can appear at if parameters vary around the set or .
Also we rewrite term (iii.2) as follows: and if the parameters change in small vicinity of , two eigenvalues of are complex having magnitude one and then NS bifurcation can emerge from fixed point .
3. Bifurcation Analysis
In this section, we will give attention to recapitulate bifurcations (flip and Neimark-Sacker) of system (4) around 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 with arbitrary parameter . One can consider the case of in a similar fashion. Since the parameters lie in , let , and then the eigenvalues of positive fixed point areThe condition leads to
Using the transformation , and writing , we shift the fixed point of system (4) to the origin. After Taylor expansion, system (4) reduces towhere andIt follows thatand .
Therefore, we obtain the following symmetric multilinear vector functions of :
Let be two eigenvectors of for eigenvalue such that and . Then by direct calculation we getWe set , where
Then by the standard scalar product in defined by , we show that . The direction of the flip bifurcation is obtained by sign , the coefficient of critical normal form [13], and is given byWe state the following result on flip bifurcation according to the above analysis.
Theorem 4. If (20) holds, , and the parameter changes its value around , then system (4) undergoes a flip bifurcation at positive fixed point . Moreover, the period 2 orbits that bifurcate from are stable (resp., unstable) if (resp., ).
3.2. Neimark-Sacker Bifurcation
We consider system (4) at fixed point with arbitrary . From (12), the eigenvalues are given by
Since the parameters belong to , the eigenvalues will be complex and LetTherefore, we have
Moreover, if , thenwhich obviously satisfies
Suppose that are two eigenvectors of and for eigenvalues and such that
Then by direct computation we obtainWe set , where
Then it is clear that where for . Now, we decompose vector as , for close to and . Obviously, . Thus, we obtain the following transformed form of system (21) for near : where with and is a smooth complex-valued function. After Taylor expression of with respect to , we obtain According to multilinear symmetric vector functions, the coefficients areThe invariant closed curve appear in the direction which is determined by the coefficient and calculated viawhere .
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, , and the parameter changes its value in small vicinity of , then system (4) passes through a Neimark-Sacker bifurcation at positive fixed point . Moreover, if (resp., ), then there exists a unique attracting (resp., repelling) invariant closed curve which bifurcates from .
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 , and fix .
Case 2. Vary in range , and fix .
Case 3. Vary in range , in range , and fix
For Case 1, taking parameters , and and covering , we find that, at fixed point , a flip bifurcation occurs at with multipliers , , and . This verifies Theorem 4.
According to bifurcation diagrams shown in Figures 1(a) and 1(b), we see that stability of fixed point happens for and loses its stability at and period doubling phenomena lead to chaos for . 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 , chaotic set occurs for , and the period 6 orbit occurs at within the window of chaotic region . The status of stable, periodic, or chaotic dynamics is compatible with sign in Figures 1(c) and 1(d).
(a)
(b)
(c)
(d)
For Case 2, taking parameters , , and and covering , we observe that a Neimark-Sacker (NS) bifurcation appears at fixed point for . Also, we have , , , , , , and . This verifies Theorem 5.
The bifurcation diagrams shown in Figures 2(a) and 2(b) demonstrate that stability of happens for and loses its stability at and an attracting invariant curve appears if . 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 , the sign of maximum Lyapunov exponent is confirming presence of chaos. Figure 2(d) is local amplification of Figure 2(a) for .
(a)
(b)
(c)
(d)
(e)
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 and , respectively. We also see that a fully developed chaos in system (4) occurs at .
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 and 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 and and the nonchaotic dynamics for and (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.
(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 where are Lyapunov exponents and is the largest integer such that and . 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.
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 to system (4), the controlled form of system (4) becomeswhere the feedback gains are denoted by and and represent positive fixed point of system (4).
The Jacobian matrix of the controlled system (43) is given bywhere , , , and are evaluated at . The characteristic equation of (44) iswhere and . Let and be the roots of (45). Then
The solution of the equations and determines the lines of marginal stability. These conditions confirm that . Suppose that , then from (47) we have
Assume that , then from (46) and (47) we get
Next, assume that , then from (46) and (47) we obtain
Then the lines , and (see Figure 5(a)) in the plane determine a triangular region which keeps eigenvalues with magnitude less than .
(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 and rest as in Case 2. The initial value is , and the feedback gains are and . Figures 5(b) and 5(c) show that, at the fixed point , 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 . 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 or and . 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.
Conflicts of Interest
The authors declare that there are no conflicts of interest.