Abstract
The complex dynamics of a nonlinear discretized predator-prey model with the nonlinear Allee effect in prey and both populations are investigated. First, the rigorous results are derived from the existence and stability of the fixed points of the model. Second, we establish a model with the Allee effect in prey undergoing codimension-one bifurcations (flip bifurcation and Neimark–Sacker bifurcation) and codimension-two bifurcation associated with 1 : 2 strong resonance by using center manifold theorem and bifurcation theory, and the direction of bifurcations is also evaluated. In particular, chaos in the sense of Marotto is proved at some certain conditions. Third, numerical simulations are performed to illustrate the effectiveness of the theoretical results and other complex dynamical behaviors, such as the period-3, 4, 6, 8, 9, 30, and 43 orbits, attracting invariant cycles, coexisting chaotic sets, and so forth. Of most interest is the finding of coexisting attractors and multistability. Moreover, a moderate Allee effect in predators can stabilize the dynamical behavior. Finally, the hybrid feedback control strategy is implemented to stabilize chaotic orbits existing in the model.
1. Introduction
The studies on biological mathematics indicate that discrete-time models described by different equations are more pertinent mathematical tools to formulate and simulate the dynamical properties than continuous ones, particularly when populations have nonoverlapping generations [1]. Moreover, using the discrete case is also crucial because it may reveal richer dynamical behaviors and provide more efficient computations. For instance, chaos can be observed in single-species discrete models, whereas for continuous-time models, minimum three-species models are needed [1]. The models of discrete-time take the following form:
As early as the 1970s, May [2, 3] had clearly illustrated the complex dynamics possible in simple discrete-time models, whether it is from the standpoint of ecology or a mathematical point of view, to investigate the complex dynamics of discrete models are very meaningful [4–9].
The most important and fascinating topic of current research in biological mathematics is the impact of the Allee effect on population dynamics. Ecologist Allee [10] first described that the Allee effect shows a positive relationship between population density and the growth rate per capita, and the existence of it indicates that there is a minimum population size necessary for a population to sustain itself in the natural world [11, 12]. There are various kinds of biological and physical factors behind this Allee mechanism, such as complications in finding mates, lessened defenses against predators, reproductive facilitation, inbreeding depressions, social thermoregulation, environment conditioning, and so on [13, 14]. This effect has already been discovered in numerous organisms in the real world, including insects [15], birds and mammals [16], plants [17], bacteria [18], marine invertebrates [19], and plant-herbivore interactions [20]. However, the Allee effects of different species, because of their different mechanisms [21], may affect the dynamics of populations differently. It is now well understood that the Allee effect plays an important role in the conservation of endangered and exploited species. Therefore, the study of the Allee effect is important to conservation biology.
In recent decades, the Allee effect in the discrete-time predator-prey model [22–28] has developed quickly. Celik and Duman [29] studied the stability of the positive equilibrium point for the discrete predator-prey model without the Allee effect:and the model with the Allee effect in the prey population:where and represent the densities of prey and predator populations at time , respectively. is the intrinsic growth rate parameter of prey and is the conversion efficiency parameter of prey to predator. and are positive constants. in model (3) denotes the Allee effect function [12, 13, 30], where is the Allee constant which means the population size at which fitness is half its maximum value [12, 14] and . Thereafter, Yuan and Yang [31] derived the existence and bifurcation directions of flip bifurcation, Neimark–Sacker bifurcation, and codimension-two bifurcations for model (2) with 1 : 2 resonance by using center manifold theorem, bifurcation theory, and the normal form method. Chen et al. [32] found some similar dynamic characteristics for models (2) and (3) and the important role played by the Allee constant in model (3). On the basis of model (2), Wu and Zhao [33] took as the Allee effect function on prey population and analyzed the detailed bifurcation behaviors of codimension-one and codimension-two with 1 : 2 resonance of the positive fixed point. Wang et al. [34] took as the mate-finding Allee effect (i.e., it was hard to find a compatible and receptive mate at low population size or density [14, 19], is the Allee effect constant) on prey population and considered the following model:
Further, they considered the Allee effects on both prey and predator populations as the modelwhere incorporates the Allee effect on predator population and is the Allee effect constant. Bigger reflects that the stronger the Allee effect and the predator density with respect to the prey population will be decreased. For these two models, they obtained the asymptotically stable conditions of the equilibrium points and showed the bifurcation or chaotic phenomena by numerical simulations but did not theoretically prove the existence of them. Enlightened by the abovementioned works, in this paper, we are committed to analyze the complex dynamics (such as codimension-one bifurcations, codimension-two bifurcation, and Marotto’s chaos) of models (4) and (5) in detail by both analytical and numerical methods and obtain some new points (such as the coexisting attractors and multistability). The results obtained demonstrate that the Allee effect may promote persistence for both prey and predator and remain its stable dynamical behavior.
The paper is organized as follows. In Section 2, the rigorous results are derived from the existence of various fixed points of models (4) and (5) and their local dynamics. In Section 3, the details about the bifurcation analysis of codimension-one and codimension-two are given. In Section 4, chaos in the sense of Marotto is proved when certain conditions are satisfied. In Section 5, the extensive results of numerical simulations are used to verify the analytical results and to display rich dynamical behaviors. In Section 6, the hybrid feedback control strategy is utilized to control the chaos caused by the bifurcations of the model (4). In Section 7, our findings are summarized.
2. Analysis of Fixed Points
2.1. Existence and Stability of Fixed Points of Model (4)
Fixed points of the model (4) are determined by solving the following nonlinear equations:
For all parameter values , and , the trivial fixed point and the axial fixed point always exist, they are two boundary fixed points.
In the following theorem, we give the existence and the number of the positive (coexistence) fixed point in .
Theorem 1. Assume , and are positive real numbers; if , the model (4) has a unique positive (coexistence) fixed point , where andIf , the model (4) has no positive fixed point.
Proof. LetFor convenience, we assumeIt is easy to know that , andfor , so is a concave function in (0, 1). is the straight line crossing with slope . Thus, there will be only two cases for the roots of as follows:(i)If , then there is no root in (0, 1) of , in other words, no positive (coexistence) fixed point exists, as shown in Figure 1(a).(ii)If , then there exist exactly one root in (0, 1) of , in other words, exactly one positive (coexistence) fixed point exists, where and , as shown in Figure 1(b).The Jacobian matrix of the model (4) evaluated at any fixed point is as follows:whereIn equation (12), the symbol denotes a partitioned matrix with row blocks and column blocks, and the (i, j)-th block is , a smaller matrix, all other blocks are zero matrices (see more details in [35]).
The characteristic equation of the Jacobian matrix is as follows:whereIn order to obtain a clear insight into the local stability analysis of the fixed point , we shall introduce the following definition and lemma from references [8, 36].
(a)
(b)
Definition 1. The fixed pointis called(i)sink if and , and it is locally asymptotic stable;(ii)source if and , and it is locally unstable;(iii)saddle if and or vice versa ( and ), and it is locally unstable;(iv)nonhyperbolic if either or .
Lemma 1. Let, suppose that, andare two roots of equation, then(i)andif and only ifand;(ii)andif and only if and;(iii)and (orand ) if and only if;(iv)andif and only ifand;(v)andare complex conjugates withif and only ifand.
Using Definition 1 and Lemma 1, we have the following results.
Theorem 2. Model (4) admits:(1)The trivial fixed point is always nonhyperbolic.(2)The axial fixed point is(i)Source if and only if ;(ii)Saddle if and only if ;(iii)Nonhyperbolic if and only if .
Proof. The Jacobian matrix of model (4) evaluated at and is given by equation (15):respectively. Obviously, the trivial fixed point is always nonhyperbolic. It can be easily observed that the eigenvalues of are and , thus is saddle point if , i.e., ; is source if , i.e., and ; is nonhyperbolic if , i.e., and .
Theorem 3. Suppose that , then the unique positive fixed point is(1)sink if and only if(2)source if and only if(3)saddle if and only if where .
Proof. We first obtain that holds if . (1) holds if and only if . On the other hand, holds if and only if . Hence, is sink if and only if(2)The conclusion is the same as the case (1). We obtain holds if and only if . Hence, is source if and only if(3)We obtain holds if and only if . Hence, is saddle if and only if It is well known that a sink which is locally asymptotic stable; thus, the condition (16) is also the locally asymptotically stable condition for the positive fixed point .
2.2. Existence and Stability of Fixed Points of Model (5)
Fixed points of model (5) are determined by solving the following nonlinear equations.
Using the same method in the previous subsection, the existence result about fixed points of the model (5) in is summarized as follows.
Proposition 1. Supposeandare positive parameter values, model (5) has two boundary fixed points , and unique positive (coexistence) fixed point providing that , where satisfys
The Jacobian matrix of model (5) evaluated at any fixed point is as follows:where
The characteristic equation of the Jacobian matrix is as follows:where
It follows from Lemma 1 that the boundary fixed point and are always nonhyperbolic. Using the same method of Theorem 2, some simple manipulation yields the asymptotically stable condition of the positive fixed point in the following theorem.
Theorem 4. The unique positive fixed pointis locally asymptotically stable if and only ifwhere .
3. Analysis of Bifurcations
In this section, the qualitative behavior of model (4) is analyzed by studying codimension-one bifurcations (flip bifurcation and Neimark–Sacker bifurcation) and codimension-two bifurcation at the unique positive (coexistence) fixed point . The phenomena of bifurcations are verified by using the center manifold theorem and bifurcation theory [37–39].
3.1. Flip Bifurcation
In this subsection, we show that there exist some values of parameters such that model (4) undergoes flip bifurcation by choosing as a bifurcation parameter. The characteristic equation of model (4) at a positive fixed point is as follows:where
Let the critical parameter value is , if satisfiesthen, the roots of characteristic (4) at the fixed point are as follows:
The condition implies that , i.e.,
We define the following set:
The variation of the parameters in the small neighborhood of the set gives emergence of flip bifurcation of model (4). Suppose that in the following analysis.
Let and . We transform the fixed point of model (4) into the origin and expand the Taylor series around to the third order, take the as a new dependent variable, we have the following:where the coefficients of model (35) are given in Appendix A.
We construct an invertible matrixand use the translation
then, model (35) can be changed intowherewhere the coefficients of are given in Appendix A.
Now, we use the center manifold theory to determine the center manifold of model (38) at the fixed point (0, 0) in a small neighborhood of . The center manifold can be presented as follows:
Assume that the center manifold takes the following form:then, apply the model (38) on two sides of simultaneously, we have
By comparing coefficients of (42), we obtain
Therefore, model (38) restricted to the center manifold as follows:where
In order for map (44) to undergo a flip bifurcation, we require that the following two discriminatory quantities and to be nonzero, whereand
From the previously mentioned analysis, we obtain the following theorem on flip bifurcation.
Theorem 5. If the bifurcation parameter varies in the small neighborhood of, andthen, model (4) undergoes a flip bifurcation at positive (coexistence) fixed point . What’s more, if () in equation (47), the period-2 orbits that bifurcate from this fixed point are stable (respectively unstable).
3.2. Neimark–Sacker Bifurcation
Next, we give the existence criteria and the direction for Neimark–Sacker bifurcation by choosing as a bifurcation parameter. The characteristic equation of the model (4) at positive fixed point is as follows:and the roots of characteristic (49) have the following form:where
The eigenvalues will be a pair of complex conjugate for , that is,
Let the critical parameter value be , suppose satisfies
We can obtain that and .
We define the following set
The variation of the parameters in the small neighborhood of the set gives emergence of the Neimark–Sacker bifurcation of model (4). Suppose that in the following analysis.
Under the conditions (52) and (53), we have
In addition, if , which implies
Correspondingly, we have .
At , we shift the fixed point of model (4) into the origin by letting . After the substitution, we obtainwhere are given in Appendix B.
Let the invertible matrixand under the transformation
model (58) becomeswherewhere the coefficients are given in Appendix B.
In order to guarantee the existence of Neimark–Sacker bifurcation of model (4), the following discriminatory quantity will be used to determine the direction of the bifurcation should not be zero:where
In which,
After some calculation, we obtainwhere
From the aforementioned discussions, we have the following theorem related to Neimark–Sacker bifurcation.
Theorem 6. If the bifurcation parameter varies in the small neighborhood ofsatisfying the conditions (56), (57) andthen, model (4) undergoes a Neimark–Sacker bifurcation at positive (coexistence) fixed point . More concretely, if (respectively, ), then Neimark–Sacker bifurcation is supercritical (respectively, subcritical) and an attracting (respectively, repelling) invariant closed curve bifurcates from the fixed point for (respectively, ).
3.3. Codimension-Two Bifurcation
Now, we discuss the conditions for the occurrence of codimension-two bifurcation associated with 1 : 2 strong resonance of the model (4) by selecting and as bifurcation parameters. The characteristic equation of model (4) at a positive fixed point is as follows:where
By the relation between roots and coefficients, has a root −1 with multiplicity 2 if and only if and . We define the following set:
It is clear see that if , then the two eigenvalues of are . Thus, model (4) undergoes codimension-two bifurcation with 1 : 2 strong resonance at .
Taking parameters arbitrarily from , we consider model (4) with as followswhere are small perturbations of parameters .
Let , then we can transform the fixed point to the origin and expand the Taylor series at to the third-order then, model (72) takes the following form:where the coefficients and are given in Appendix C.
Make an invertible matrixand use the transformationthen, model (73) can be rewritten as follows:wherewith , , and ; the coefficients of are given in Appendix C.
We denote the invertible matrixand take the nonsingular linear coordinate transformationthen, model (76) becomeswhere , , andthe coefficients of are given in Appendix C.
Lastly, we take the following transformationwhere coefficients and will be confirmed in the following discussion. The inverse transformation of the variable from transformation (82) isin which
Using transformation (82) and (83) for model (80), we can obtain as follows:whereand the coefficients of are given in Appendix C.
In order to eliminate all quadratic terms, we takeand can obtain and for all nonnegative integers and . Further, in order to remove all cubic but those resonant terms, we takeand can obtain and for all nonnegative integers and . After the previously mentioned transformations, model (85) can finally be reduced into the following formwhere and . When , we have
By making use of the results on the existence of 1 : 2 strong resonance bifurcation theorems in [37, 39] and the analysis performed, we have the following theorem.
Theorem 7. If the bifurcation parametersandvary in the small neighborhood of, andthen, model (4) undergoes a codimension-two bifurcation with 1 : 2 strong resonance at positive (coexistence) fixed point .
4. Analysis of Marotto’s Chaos
In the present section, Marotto’s theorem in references [40, 41] is used to prove that the unique positive (coexistence) fixed point of model (4) is a snap-back repeller and model (4) possesses a chaotic behavior in the sense of Marotto.
Before giving the conditions for the existence of chaotic phenomena, we first present the definitions and theorem of Marotto’s chaos given in references [32, 33]. For any map , and any positive integer , let represent the composition of with itself times. For a differentiable function , let denote the Jacobian matrix of evaluated at the point , and its determinant. Let denote the closed ball in of radius centered at the point and its interior. Also, let be the usual Euclidean norm of in .
Definition 2. Letbe differentiable in. The pointis an expanding fixed point ofin, ifand all eigenvalues ofexceed 1 in norm for all.
Definition 3. Suppose thatis an expanding fixed point of a given mapinwith, thenis called a snap-back repeller of if there exists a point satisfies , , andfor some positive integer.
Theorem 8. If has a snap-back repeller, then the mapis chaotic.
The characteristic equation of the Jacobian matrix of model (4) has a pair of conjugate complex roots and , and if and only ifwhere .
In order to give the conditions under which the fixed point is a snap-back repeller, a neighborhood of is found to satisfy (92) for all . For convenience, the following equations are proposed:
By some calculations, we haveandwhere
Then, the equation has two real rootsonly if . Since , we denote . Thus, if and .
On the other hand,withwhere
Then, the equation has two real rootswhenwhere . Denote , thus, if one of the following conditions is satisfied.(1) and ;(2) and .
According to the previous analysis, we obtain the following lemma.
Lemma 2. Let, and, denote the following conditions.
If the fixed point of model (4) satisfiesthen, is called an expanding fixed point of in , is defined by the map (4).
By Definition 3 of a snap-back repeller, we need to find one point which satisfies the conditions and , where is the times iterations, is an integer and .
Note thatand
If equations (105) and (106) have solutions which are all different from , then a map is constructed to map the point to the fixed point after two iterations.
The solutions of (104) satisfy the following equations:
Putting and in (105), and solving , we obtain
By a simple calculation, we obtainwhere
According to the previous analysis, we have the following theorem.
Theorem 9. Assume that the conditions in Lemma 2 hold. If the solutions of equations (107) and (108) satisfy thatthen is a snap-back repeller of model (4) in , and in this case the model (4) is chaotic in the sense of Marotto.
Next, we give an example to illustrate the existence of conditions in Theorem 9.
Example 1. For , and , the map (4) has a positive fixed point and the eigenvalues associated with is . According to the conditions of Lemma 2 and Theorem 9, we find that a region of (see Figure 2) isand there exists a point satisfying that andwherewhere , , and . We can easily observe that . Therefore, is a snap-back repeller and the model (4) is chaotic in the sense of Marotto.
5. Numerical Simulations
Diagrams for bifurcation, phase-plane, maximum Lyapunov exponents, the basin of attractor, two-parameter dynamical map, and time sequence of the model (4) and model (5) are drawn to validate the theoretical results in this section. Furthermore, we find some new, complex, and surprising dynamical behaviors as the parameters vary.
5.1. Numerical Simulations for Bifurcations and Marotto’s Chaos at
5.1.1. Flip Bifurcation
The case of flip bifurcation is given in Figure 3. Suppose that parameters , and , based on Theorem 1, we know model (4) has only one positive fixed point . Figure 3(a) shows the bifurcation diagram of the model with respect to the bifurcating parameter in plane with initial values (0.9, 0.9), the flip bifurcation (labeled “FB”) takes place for the fixed point with and . The chaotic attractor for is given in Figure 3(b) and the corresponding time sequence diagram is given in Figure 3(c) which shows the chaotic nature. Figure 3 confirms the correctness of the theoretical results in Theorem 5.
(a)
(b)
(c)
5.1.2. Neimark–Sacker Bifurcation
The case of Neimark–Sacker bifurcation is given in Figure 4. Suppose that parameters and , based on Theorem 1, we know model (4) has only one positive fixed point. Figure 4(a) shows that the bifurcation diagram of model (4) with respect to the bifurcating parameter in plane with initial values (0.5, 0.5), the Neimark–Sacker bifurcation (labeled “NSB”) takes place for fixed point and with and , an attracting invariant cycles appears. The chaotic attractor for is given in Figure 4(b) and the corresponding time sequence diagram is given in Figure 4(c) which shows the chaotic nature. Figure 4 confirms the correctness of the theoretical results in Theorem 6.
(a)
(b)
(c)
5.1.3. Codimension-Two Bifurcation
The case of codimension-two bifurcation is given in Figure 5. Suppose that parameters , and , based on Theorem 1, we know the model (4) has only one positive fixed point. After calculating, we obtain the positive fixed point , and the eigenvalues of the Jacobian matrix are . Furthermore, we have and , thus , and . Figure 5(a) shows that the three-dimensional bifurcation of the model (4) in plane for and vary in a neighborhood of with initial values (0.6, 0.6). Therefore, it follows from Theorem 7 that a fixed point is a 1 : 2 strong resonance bifurcation point. Figure 5(b) shows the two-dimensional bifurcation diagram in plane when . Figure 5 confirms the correctness of the theoretical results in 3.3.
(a)
(b)
5.1.4. Marotto’s Chaos
Using Example 1, the case of Marotto’s chaos is given in Figure 6. Suppose that parameters , and , based on Theorem 1, we know the model (4) has only one positive fixed point . Figure 6(a) shows the bifurcation diagram in plane with initial values (0.3, 0.3). The maximum Lyapunov exponents corresponding to (a) are calculated and plotted in Figure 6(b). It can be seen from Figure 6(b) that all the maximal Lyapunov exponents are larger than 0, so the model (4) with the given initial values are all chaotic. From Example 1 and Theorem 9, one can know that the chaotic attractor is located in the chaotic region , which is in line with the maximal Lyapunov exponents in Figure 6(b). The Marotto’s chaotic attractor is exhibited in Figure 6(c) for . Figure 6 verifies the correctness of the theoretical results in Theorem 9.
(a)
(b)
(c)
5.2. Numerical Simulations for Model (4) in Different Bifurcation Parameters
In this part, we investigate the influence of the parameters on the dynamics of the system (4), and some new and more complex dynamic behaviors are presented as the parameters vary. Bifurcation parameters are considered in the following four cases, respectively:(1)Vary in range , and fix ;(2)Vary in range , and fix ;(3)Vary in range , and fix ;(4)Vary in range , and fix .
5.2.1. Case (1)
Taking parameters and covering [0, 1.9] with initial values (0.66, 0.66), the bifurcation diagram of model (4) with the Allee effect on prey population in plane is given in Figure 7(a) and the local amplification is given in Figure 7(b) for . The maximum Lyapunov exponents corresponding to (a) are plotted in Figure 7(c) which confirms the existence of the chaos and period orbits in the parametric space. In Figure 7(a), there is a stable fixed point for more large regions , the stable fixed point loses its stability at and closed invariant circle appears when the parameter exceeds 1.57 and eventually leads to chaos at . There are two different routes (invariant circle and inverse period-doubling bifurcation) to the chaos which occurs simultaneously. The inverse period-doubling bifurcation phenomena occurs at , and there are some cascades of period-doubling bifurcations lead to chaos. The chaotic behavior suddenly appears at . Figure 7(b) shows that there is an unending sequence of bifurcations until chaos. The phase-plane diagrams corresponding to Figure 7(a) are shown in Figures 7(d)–7(f) for the two-coexisting chaotic sets at , four-coexisting chaotic sets at , and the invariant circle at , respectively.
(a)
(b)
(c)
(d)
(e)
(f)
5.2.2. Case (2)
Taking parameters and covering [0, 1.9] with initial values (0.08, 0.08), the bifurcation diagram of model (4) with the Allee effect on prey population in plane is given in Figure 8(a) and the local amplifications are plotted in Figures 8(b) and 8(c) for and , respectively. The maximum Lyapunov exponents corresponding to (a) are plotted in Figure 8(d) which confirms the existence of the chaos and period orbits in the parametric space. Figure 8(a) shows that there is a stable fixed point for and the fixed point loses its stability as increases. The period-3 orbits are observed for , and there happen two different routes to chaos at the same time; they are inverse period-doubling bifurcation and period-doubling bifurcation. When increases at , a period-6 orbit appears, and some cascades of period-doubling bifurcations lead to chaos. The chaotic behaviors suddenly appear at . One can also observe that there is an inverse period-doubling bifurcation which occurs at changes into chaos with complex period windows and the chaotic behavior suddenly appears at . Figure 8(c) depicts the alternates of model (4) between periodic and chaotic motions. The phase-plane diagrams corresponding to Figure 8(a) are shown in Figures 8(e)–8(i) for the two-coexisting chaotic sets at , period-3 orbits at , period-6 orbits at , three-coexisting chaotic sets at , and chaotic attractor is at , respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
5.2.3. Case (3)
Taking parameters and covering [1, 3.4] with initial values (0.4, 0.4), the bifurcation diagram of the model (4) with the Allee effect on prey population in plane is given in Figure 9(a). Figure 9(b) is the local amplifications for . The maximum Lyapunov exponents corresponding to (a) are plotted in Figure 9(c) which confirms the existence of the chaos and period orbits in the parametric space. Figure 9(a) shows that the stability of fixed point happens for , loses its stability as increases and through a Neimark–Sacker bifurcation for . For , the fixed point becomes unstable and a closed invariant circle appears, the radius of the circle becomes bigger and bigger, and eventually leads to chaos. In particular, for , the invariant circle disappears, then it becomes period-4 orbits at suddenly. Model (4) undergoes period-doubling phenomena, which leads to chaos with complex period windows. The phase-plane diagrams for various values of are shown in Figures 9(d)–9(l), which clearly illustrate the act of how an invariant circle bifurcates from the stable fixed point and increases its radius, the invariant circle leads to chaotic attractors and the chaotic behavior suddenly disappears at . Figure 9 shows that there are period-4, period-8, period-30, four-coexisting chaotic sets, and chaotic attractor.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
5.2.4. Case (4)
Taking parameters and covering [4, 6.5] with initial values (0.4, 0.4), the bifurcation diagram of the model (4) with the Allee effect on prey population in plane is given in Figure 10(a) and the local amplifications for and are given in Figures 10(b) and 10(c), respectively. The maximum Lyapunov exponents corresponding to (a) are plotted in Figure 10(d) which confirms the existence of the chaos and period orbits in the parametric space. Figure 10(a) shows that the stability areas of a fixed point is , the fixed point becomes unstable as increases. This would suggest that a population increases the carrying capacity at a certain growth rate and then lose its stability. The phase-plane diagrams in Figures 10(e)–10(m) clearly depicts the process of an invariant curve becoming increasingly complex as parameter growth and chaotic attractors appear. When increases at certain values, for example, at , a closed curve appears, indicating that model (4) undergoes a Neimark–Sacker bifurcation, and the quasiperiodic orbit arise. At , the circle disappears and a period-9 orbit appears, then evolves into chaos. Figure 10 shows the rich types of dynamics for the model (4); there are period-9, period-43, quasiperiodic orbit, nonattracting chaotic set, and chaotic attractor.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
(m)
5.3. Coexisting Attractors and Multistability of Model (4)
As mentioned in Subsection 5.2, model (4) can exhibit a variety of rich dynamical behaviors when varying the parameters , and , respectively. In this subsection, we will present some coexisting attractors and multistability to obtain more insightful information on the complex dynamics of the model (4).
With a view of exploring more detailed dynamics of model (4), we first plot a two-parameter dynamical diagram [42] as illustrated in Figure 11 to present the impact of the parameters , , and . It can be clearly seen from Figure 11 that there are two color areas in which the yellow area denotes the periodic state and the dark-blue region represents the chaotic state. Figures 11(a) and 11(b) show the dynamics regions of the model (4) in the plane with the initial conditions (0.08, 0.08) and initial conditions (0.7, 0.7), respectively. Figures 11(c) and 11(d) show the dynamics regions of the model (4) in the plane with the initial conditions (0.08, 0.08) and initial conditions (0.7, 0.7), respectively. Observed from Figures 11(a) and 11(b), there are some finely banded yellow regions occurring in the dark-blue area, indicating that the orbits of the model (4) alter many times in these regions between period and chaos. This also shows that model (4) is extremely sensitive in these regions, meaning the emergence of the complicated state transition behavior. Furthermore, when we vary both the parameters and within a certain range, one can observe that the color area is always relatively concentrated and complete in the yellow region (chaos) or dark-blue region (period), which means model (4) has strong structural stability. Especially, viewing the color of the dynamical map horizontally along the direction of , the dynamical behaviors are in good agreement with the bifurcation diagram of Figure 12(a). Figures 11(c) and 11(d) show that when the parameter is small, model (4) is in a chaotic state. Then it converts into the periodic state, next drop into chaos again, and finally go out of chaos into periodic orbits with gradually increasing. We also see that period is mainly concentrated in the left half of the graph while chaos has mainly appeared in the right half. Remark that seeing the color of the dynamical map horizontally along the direction of , the dynamical behaviors also coincide well with the bifurcation diagram of Figure 12(a). A comparison of Figures 11(a) and 11(b) (or Figures 11(c) and 11(d)) shows that different initial conditions may be able to have different dynamical behaviors when taking the same parameters and . In conclusion, the two-parameter dynamical map, as a useful tool, can offer a better parameter selection for population evolution in the real world.
(a)
(b)
(c)
(d)
(a)
(b)
Based on the previously discussed analysis, our attention is further focused on the study of coexisting attractors and multistability. The coexisting attractors is an important nonlinear dynamics which means two or more attractors generate simultaneously under different initial conditions but with the same parameter set [43]. Particularly, owing to only one positive stable fixed point, the resulting attractors are all hidden [43, 44]. The trajectories of the model with multistability can asymptotically approach different states when varying initial conditions and fixing parametric values. In order to observe this interesting phenomenon in the model (4), suppose that , and select as the adjustable parameter. The coexisting bifurcation diagram with varying parameter is displayed in Figure 12(a) and the corresponding maximum Lyapunov exponents to (a) are plotted in Figure 12(b). Note that the cyan trajectory and blue trajectory in Figure 12 begin with the initial conditions (0.08, 0.08) and (0.7, 0.7), respectively. The parametric region in which Figure 12(a) presents different behaviors is called the coexisting zone. Obviously, from Figure 12(a), one can observe the coexistence of period and chaos in an extremely narrow area between the two vertical red dashed lines with . With the parameter varying in the coexisting zone, there are chaotic attractors coexist with the periodic attractors (e.g., period-3 orbits, period-6 orbits, and period-12 orbits). When setting the parameter and keeping other parameters unchanged, the phase portrait of the coexistence of the period-3 attractor and chaotic attractor are illustrated in Figure 13. It is worth noting that the bifurcation diagram matches well with its corresponding maximum Lyapunov exponents.
It is well known that the coexisting attractors can be characterized by the basins of attraction, and the multistability depends on the different basins of attraction [42]. An attractor’s attraction basin is the region of the initial condition plane, which tend to the attractor as the time approaches infinity. That is when we start with points that are within the basin of attraction, they will finally be iterated into the attractor [45, 46]. Fix the parameters and . The basin of attraction of the coexisting period-3 attractor and chaotic attractor of Figure 13 is shown in Figure 14. Different colors represent different attraction regions, where the orange-red and the royal blue indicate period-3 attractor and chaotic attractor, respectively. Figure 14 depicts such model is crucially dependent on the initial states, and selecting the initial value in different colors, the model (4) will generate chaotic or periodic attractors. What is more, it’s easy to observe that the model is delicate with respect to the initial conditions. Especially, the orange-red region and royal blue region are interweaved with each other in the bottom half of the graph, implying that a tiny perturbation to the initial state is likely to lead to a state transition behavior.
5.4. Impact of the Allee Effect on the Stability of Model (5)
For better visualization of the impact of the Allee parameter on the stability of model (5), we draw the bifurcation diagrams of model (5) in plane, plane and plane, respectively. Bifurcation parameters are considered in four cases which is the same as subsection 5.2.
In Figure 15, the four different values of are selected. Figures 15(a)–15(d) show that if we gradually increase the value of for the model (5) with the Allee effect on both populations, then model (5) becomes stable, i.e., the Allee effect in predators has a stabilizing effect on the model (5) with the set of parameter values, the chaotic dynamics is observed at large values of than the model without Allee constant (see Figures 7(a) and 8(a)). However, changing the values of has no significant effect on the inverse period-doubling bifurcation (see Figures 15(a)–15(h)). In fact, we can use the feedback control technique to change the chaotic dynamics of inverse period-doubling bifurcation (see Section 6). A similar behavior is observed in the remaining subfigures ((d)-(h)), ((i)-(l)), and ((m)-(p)) also. The results we obtained mean that the Allee constant plays a very important role in managing the population.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
(m)
(n)
(o)
(p)
6. Chaos Control
Controlling chaos is one of the most vital and developed areas for population models in the current research. Chaos due to the emergence of bifurcations is considered to be unfavorable phenomenon for the survival of populations. Hence, it is necessary to execute chaos control techniques to avoid the population from unpredictable situations. In this section, the hybrid feedback control strategy (state feedback and parameter perturbation) in references [31, 47–49] is applied to control the chaotic behavior which develops because of the appearance of bifurcation of model (4). This strategy has significant characteristics in population models, especially in discrete-time population models associated with biological species, which can also make the chaotic behavior predictable and stable.
Assuming that model (4) substantiates bifurcation at the unique positive (coexistence) fixed point , the controlled form of model (4) is given as follows:where is the control parameter.
The Jacobian matrix of controllable model (115) estimated at is given as follows:then, its characteristic equation is as follows:where
Next, we give the following theorem on the stability of the controlled model (113).
Theorem 10. The unique positive (coexistence) fixed point of the controlled model (115) is locally asymptotically stable when the following condition holds:
Proof. According to Lemma 1 (i), the point is locally asymptotically stable if and only if and . Thus, holds.
The following are two examples to demonstrate how this strategy controls the unstable fixed point.
Example 2. As seen in subsection 5.1, the model (4) undergoes flip bifurcation when . We apply the hybrid feedback control strategy to control chaos which appears because of the flip bifurcation. The controlled model (115) takes the following form:Table 1 shows that the control interval of stability is . In addition, Figure 16(a) presents the bifurcation diagram for the controlled model (120), which also indicates that the flip bifurcation of the fixed point of the controlled model (120) can be advanced.
(a)
(b)
Example 3. As seen in subsection 5.1, the model (4) undergoes Neimark–Sacker bifurcation when . We apply the hybrid feedback control strategy to control the chaos which appears because of the Neimark–Sacker bifurcation. For this, the controlled model (113) takes the following form:Table 2 shows that there exists a control interval of stability for . Moreover, Figure 16(b) illustrates the bifurcation diagram for the controlled model (121), which also indicates that the Neimark–Sacker bifurcation of the fixed point of the controlled system (121) can be delayed.
By choosing the appropriate value of the control parameter , the chaotic behavior is changed into different periodic orbits. For example, when we choose , which is drawn in Figure 17 that chaotic trajectories of flip bifurcation and Neimark–Sacker bifurcation are controlled in period-1 orbit. Figures 17(a) and 17(b) show a sharp contrast to Figures 3(c) and 4(c), respectively.
(a)
(b)
7. Summary
For many species, population growth is seasonal rather than a continuous process. Meanwhile, studying population dynamics involving Allee effects becomes an important subject of contemporary research.
In this study, a nonlinear discrete-time predator-prey model with a nonlinear Allee effect has been established and analyzed in detail. By combing mathematical analysis and numerical simulations, we have shown that the model with the Allee effect in prey can undergo flip bifurcation, Neimark–Sacker bifurcation, codimension-two bifurcation associated with 1 : 2 strong resonance, and Marotto’s chaos under certain parametric conditions. With the help of a bifurcation diagram, phase portraits, time series, maximum Lyapunov exponents, dual-parameter dynamical map, the basin of attraction, the model reveals copious dynamics behaviors, including orbits of period-3, 4, 6, 8, 9, 30, and 43, invariant cycles, quasiperiodic orbits, types of two, three and four coexisting chaotic sets, and nonattracting chaotic set. More interestingly, diverse coexisting attractors and multistability are also observed when selecting appropriate parameter sets and different initial conditions. Particularly, we can discover that the Allee effect has a stabilizing effect. More precisely, an appropriate increase in the Allee effect strength on predators can stabilize population dynamics, thus reducing the probability of extinction. These results could be useful for understanding the dynamics of the discrete-time predator-prey model. Furthermore, these highly complicated dynamics of the model may provide the possibility of real-world applications [43] for the model under study and provide reasonable guidance in population dynamics.
From an ecological perspective, attracting invariant closed curves, the cascade of period-doubling, and the chaotic sets imply that the prey and predator will coexist within a certain time. However, the predator-prey model is unstable if chaotic behavior happens. Chaos can lead the population to run a higher risk of extinction; thus they are harmful to the breeding of biological populations [48]. Specifically, the hybrid feedback control strategy is applied to stabilize the chaotic orbits existing in the model (4), which will make predator and prey maintain stable coexistence. Remarkably, it should be stressed that such model is often affected by external disturbances or noise [45], so there is a richness of complicated dynamical properties of this model still to be uncovered. Thus, we will pay attention to that questions in our future research work.
Appendix
A. The Coefficients of Model (35) and Model (38) in Subsection 3.1
(1)The coefficients of model (35) are shown as follows:(2)The coefficients of in model (38) are shown as follows:
B. Appendix
(1)The coefficients of model (58) are shown as follows:(2)The coefficients of in model (61) are shown as follows:
C. The Coefficients of Model (71), (74), (78) and Model (83) in Subsection 3.3
(1)The coefficients of model (75) are shown as follows:(2)The coefficients of in model (76) are shown as follows:(3)The coefficients of in model (80) are shown as follows:(4)The coefficients of in model (85) are shown as follows:
Data Availability
No data were used to support this study.
Conflicts of Interest
The author declares that there are no conflicts of interest.