Bifurcations of a two-dimensional discrete-time predator–prey model

Abstract

We study the local dynamics and bifurcations of a two-dimensional discrete-time predator–prey model in the closed first quadrant \(\mathbb{R}_{+}^{2}\). It is proved that the model has two boundary equilibria: \(O(0,0)\), \(A (\frac{\alpha _{1}-1}{\alpha _{1}},0 )\) and a unique positive equilibrium \(B (\frac{1}{\alpha _{2}},\frac{ \alpha _{1}\alpha _{2}-\alpha _{1}-\alpha _{2}}{\alpha _{2}} )\) under some restriction to the parameter. We study the local dynamics along their topological types by imposing the method of linearization. It is proved that a fold bifurcation occurs about the boundary equilibria: \(O(0,0)\), \(A (\frac{\alpha _{1}-1}{\alpha _{1}},0 )\) and a period-doubling bifurcation in a small neighborhood of the unique positive equilibrium \(B (\frac{1}{\alpha _{2}},\frac{\alpha _{1} \alpha _{2}-\alpha _{1}-\alpha _{2}}{\alpha _{2}} )\). It is also proved that the model undergoes a Neimark–Sacker bifurcation in a small neighborhood of the unique positive equilibrium \(B (\frac{1}{ \alpha _{2}},\frac{\alpha _{1}\alpha _{2}-\alpha _{1}-\alpha _{2}}{\alpha _{2}} )\) and meanwhile a stable invariant closed curve appears. From the viewpoint of biology, the stable closed curve corresponds to the periodic or quasi-periodic oscillations between predator and prey populations. Numerical simulations are presented to verify not only the theoretical results but also to exhibit the complex dynamical behavior such as the period-2, -4, -11, -13, -15 and -22 orbits. Further, we compute the maximum Lyapunov exponents and the fractal dimension numerically to justify the chaotic behaviors of the discrete-time model. Finally, the feedback control method is applied to stabilize chaos existing in the discrete-time model.

Introduction

Different models have been invoked to understand the mechanism of competition between populations of two-species. In 1931, Volterra proposed a famous prey–predator model which is represented by the following system of ordinary differential equations [1]:

where X denotes the number of prey and Y denotes the number of predator. Moreover, a, b, c, d are positive parameters. It has been shown that the number of prey grows exponentially in the absence of predators, while the number of predators decreases exponentially in the absence of a prey population. The terms \(bXY\) and \(dXY\) explain the prey–predators encounters which are conducive to predators and lethal to prey. It is noted that model (1) then takes the following form, if one consider some harvesting effect [2]:

and as a result the reasonable harvesting effect is favorable to prey population. There are also some other prey–predator models which are more fascinating and effective for a number of interacting species greater than two or which assume a parasitic infection of the populations [3, 4].

It is a well-known fact that discrete-time models described by difference equations are more beneficial and reliable than continuous-time models whenever there are non-overlapping generations in the populations. Moreover, these models also provide efficient computational results for numerical simulations and provide a rich dynamics as compared to the continuous ones [5,6,7,8,9,10]. In the last few years, many interesting papers have appeared in the literature that discuss the stability, bifurcation and chaos phenomena in discrete-time models (see [11,12,13,14,15,16,17,18,19,20] and the references cited therein).

This paper deals with the study of stability, bifurcations and chaos control of the following discrete-time predator–prey model [21]:

The rest of the paper is organized as follows: Sect. 2 deals with the study of the existence of equilibria and local stability along their different topological types of the discrete-time model (4). In Sect. 3, we study the existence of bifurcations about equilibria of the model (4). Section 4 deals with a bifurcation analysis about the unique positive equilibrium of the model (4). In Sect. 5, numerical simulations are presented to verify the theoretical results. This also includes the study of fractal dimensions which characterize the strange attractors of the model (4). In Sect. 6, we study the chaos control by the feedback control method to stabilize chaos at unstable trajectories. A conclusion is given in Sect. 7.

Existence of equilibria and local stability of the discrete-time model (4)

Lemma 1

System (4) has at least two boundary equilibria and one unique positive equilibrium in\(\mathbb{R}_{+}^{2}\). More precisely,

(i)

for all parametric values\(\alpha _{1}\)and\(\alpha _{2}\), system (4) has the boundary equilibrium\(O(0,0)\);

Now will study the local dynamics of (4) about \(O(0,0)\), \(A (\frac{\alpha _{1}-1}{\alpha _{1}},0 )\) and \(B (\frac{1}{ \alpha _{2}},\frac{\alpha _{1}\alpha _{2}-\alpha _{1}-\alpha _{2}}{\alpha _{2}} )\). The Jacobian matrix \(J_{(\widehat{x},\widehat{y})}\) of the discrete-time model (4) about equilibrium \((\widehat{x}, \widehat{y})\) is given by

Existence of bifurcations about equilibria of the discrete-time model (4)

In this section based on theoretical studies in Sect. 2, we will study the existence of bifurcations about equilibria. The existence of corresponding bifurcations about the equilibria \(O(0,0)\), \(A (\frac{ \alpha _{1}-1}{\alpha _{1}},0 )\) and \(B (\frac{1}{\alpha _{2}},\frac{ \alpha _{1}\alpha _{2}-\alpha _{1}-\alpha _{2}}{\alpha _{2}} )\) can be summarized as follows:

(i)

From Lemma 2, we can see that when \(\alpha _{1}=1\), one of the eigenvalues about the equilibrium \(O(0,0)\) is 1. So a fold bifurcation may occur when the parameter varies in the small neighborhood of \(\alpha _{1}=1\).

(ii)

From Lemma 3, we can easily see that if \(\alpha _{2}=\frac{\alpha _{1}}{\alpha _{1}-1}\) holds then one of the eigenvalues about \(A (\frac{\alpha _{1}-1}{\alpha _{1}},0 )\) is 1. So a fold bifurcation occurs when the parameter varies in a small neighborhood of \(\alpha _{2}=\frac{\alpha _{1}}{\alpha _{1}-1}\). And we denote the parameters satisfying \(\alpha _{2}=\frac{\alpha _{1}}{ \alpha _{1}-1}\) as

From Lemma 4, we see that if \(\alpha _{1}=\frac{ \alpha _{2}}{\alpha _{2}-2}\) holds then the eigenvalues of \(J_{B (\frac{1}{ \alpha _{2}},\frac{\alpha _{1}\alpha _{2}-\alpha _{1}-\alpha _{2}}{\alpha _{2}} )}\) about \(B (\frac{1}{\alpha _{2}},\frac{\alpha _{1} \alpha _{2}-\alpha _{1}-\alpha _{2}}{\alpha _{2}} )\) are a pair of complex conjugate with modulus 1. So a Neimark–Sacker bifurcation exists by the variation of parameter in a small neighborhood of \(\alpha _{1}=\frac{\alpha _{2}}{\alpha _{2}-2}\). Precisely we represent the parameters satisfying \(\alpha _{1}=\frac{\alpha _{2}}{\alpha _{2}-2}\) as

From Lemma 5, we see that if \(\alpha _{1}=\frac{3 \alpha _{2}}{3-\alpha _{2}}\) holds we see that one of the eigenvalues of \(J_{B (\frac{1}{\alpha _{2}},\frac{\alpha _{1}\alpha _{2}-\alpha _{1}-\alpha _{2}}{\alpha _{2}} )}\) about \(B (\frac{1}{\alpha _{2}},\frac{\alpha _{1}\alpha _{2}-\alpha _{1}-\alpha _{2}}{\alpha _{2}} )\) is −1 and other is neither 1 nor −1. So a period-doubling bifurcation exists by the variation of parameter in a small neighborhood of \(\alpha _{1}=\frac{3\alpha _{2}}{3-\alpha _{2}}\). More precisely we can also represent the parameters satisfying \(\alpha _{1}=\frac{3\alpha _{2}}{3-\alpha _{2}}\) as

where \(x^{*}=\frac{1}{\alpha _{2}}\), \(y^{*}=\frac{\alpha _{1}\alpha _{2}- \alpha _{1}-\alpha _{2}}{\alpha _{2}}\). Hereafter when \(\epsilon =0\), the normal form of system (8) is studied. Expanding (8) about \((u_{n},v_{n})=(0,0)\) by Taylor series, we get

Remark

According to bifurcation theory discussed in [30, 31], the bifurcation is called a supercritical Neimark–Sacker bifurcation if the discriminatory quantity \(\chi <0\). In the following section, numerical simulations guarantee that a supercritical Neimark–Sacker bifurcation occurs for the discrete-time model (4).

Hereafter we determine the center manifold \(W^{c}(0,0)\) of (17) about \((0,0 )\) in a small neighborhood of \(\alpha _{1}^{*}\). By center manifold theorem, there exists a center manifold \(W^{c}(0,0)\) that can be represented as follows:

and a non-degenerate condition for the existence of Neimark–Sacker bifurcation of the discrete-time model (4), i.e., \(\frac{ \alpha _{2}-2}{2\alpha _{2}}=0.5341880341880343>0\) hold. Moreover, after some manipulations from Mathematica one gets

In view of (22) and (23) the value of the discriminatory quantity from (12) is \(\chi =-0.11221718590370466<0\). Therefore if \(\alpha _{1}=2.34>2.33333\) then the discrete-time model (4) undergoes a supercritical Neimark–Sacker bifurcation and in fact a stable invariant close curve appears, which is presented in Fig. 2(a). Similarly for other choices of bifurcation parameter the value of discriminatory quantity is less than 0 (see Table 1) and their corresponding attracting close invariant curves are depicted in Figs. 2(b)–2(o). In the context of biology, attracting closed invariant curve bifurcations from the supercritical Neimark–Sacker bifurcation imply that host and parasitoid populations will coexist under periodic or quasi-periodic oscillations with long time.

Hereafter we will provide the numerical simulation in order to verify the theoretical results obtained in Sect. 4.2 by fixing \(\alpha _{2}=1.53\) and varying \(1.5\le \alpha _{1}\le 18.5\). Fixing \(\alpha _{2}=1.53\), then from the non-hyperbolic condition (iii) of Lemma 5 one gets \(\alpha _{1}=3.1224489795918364\). From a theoretical point of view the unique positive equilibrium point \((0.6535947712418301, 0.08163300653594788)\) of (4) is stable if \(\alpha _{1}< 3.1224489795918364\); bifurcation occurs if \(\alpha _{1}=3.1224489795918364\), and there is a period-doubling bifurcation if \(\alpha _{1}>3.1224489795918364\).

From Figs. 3(a)–3(b), we see that the equilibrium point is stable if \(\alpha _{1}< 3.1224489795918364\), and loses its stability at the period-doubling bifurcation parameter value \(\alpha _{1}=3.1224489795918364\). The maximum Lyapunov exponents corresponding to Figs. 3(a)–3(b) are plotted in Fig. 3(c). Moreover, 3D bifurcation diagrams are also plotted in Figs. 4(a)–4(c). The phase portraits which are associated with Figs. 3(a)–3(b) are depicted in Figs. 5(a)–5(f), which indicates that the discrete-time model (4) exhibits a complex dynamics such as period-2, -4, -11, -13, -15 and -22 orbits.

Fractal dimension

where \(\kappa _{1}, \kappa _{2},\ldots , \kappa _{n}\) are Lyapunov exponents and j is the largest integer such that \(\sum_{i=1}^{j} \kappa _{j}\ge 0\) and \(\sum_{i=1}^{j+1}\kappa _{j}<0\). For our under consideration discrete-time model (4), the fractal dimension takes the following form:

If \(\alpha _{2}=1.53\) then after some manipulation two Lyapunov exponents are \(\kappa _{1}= 1.3153221370247266\) (resp. \(\kappa _{1}=1.2928270741698256\)) and \(\kappa _{2}=-0.3806816141489094\) (resp. \(\kappa _{2}=-0.40393818528093656\)) for \(\alpha _{1}=1.63\) (resp. \(\alpha _{1}=1.7\)). So the fractal dimension for the discrete-time model (4) is

The strange attractors for the above fixed parametric values are also plotted and presented in Figs. 6(a)–6(b), which illustrate that the discrete model (4) has a complex dynamical behavior as the parameter \(\alpha _{1}\) increases.

Chaos control

In this section, we will study the chaos control by applying the state feedback control method [35, 36]. By adding a feedback control law as the control force \(u_{n}\) to the discrete-time model (4), the controlled model (4) takes the following form:

where the feedback gains are denoted by \(k_{1}\) and \(k_{2}\) and \((x^{*},y^{*})\) is the unique positive equilibrium point, i.e., \((x^{*},y^{*})=B (\frac{1}{\alpha _{2}},\frac{\alpha _{1}\alpha _{2}-\alpha _{1}-\alpha _{2}}{\alpha _{2}} )\) of the discrete-time model (4). The Jacobian matrix \(J_{c}\) of the controlled system (27) is

The solutions of the equations \(\kappa _{1}=\pm 1\) and \(\kappa _{1} \kappa _{2}=1\) determine the lines of marginal stability. These conditions confirm that \(|\kappa _{1,2}|<1\). Suppose that \(\kappa _{1} \kappa _{2}=1\), then from (32) one gets

Then the lines \(l_{1}\), \(l_{2}\) and \(l_{3}\) in the \((k_{1}, k_{2})\) plane determine a triangular region which gives \(|\kappa _{1,2}|<1\) (see Fig. 7(a)).

Figure 7

Control of chaotic trajectories of the controlled discrete-time model (27) for \(\alpha _{1}=1.53\), \(\alpha _{2}=3.75\) with initial values \((0.2,0.25)\) (a) stability region in \((k_{1}, k_{2})\)-plan. (b)–(c) Time series for states \(x_{n}\) and \(y_{n}\), respectively

In order to check how the implementation of feedback control method works and how to control chaos at an unstable state, we have performed numerical simulations. Figures 7(b)–7(c) show that about the unique positive equilibrium the chaotic trajectories are stabilized.

Conclusion

This work deals with the study of local dynamics, bifurcations and chaos control of a discrete-time predator–prey model (4) in \(\mathbb{R}^{2}_{+}\). It is proved that the model has the boundary equilibria \(O(0,0)\), \(A (\frac{\alpha _{1}-1}{\alpha _{1}},0 )\) and a unique positive equilibrium \(B (\frac{1}{\alpha _{2}},\frac{\alpha _{1} \alpha _{2}-\alpha _{1}-\alpha _{2}}{\alpha _{2}} )\). We have investigated local stability along their topological types about \(O(0,0)\), \(A (\frac{\alpha _{1}-1}{\alpha _{1}},0 )\), \(B (\frac{1}{\alpha _{2}},\frac{\alpha _{1}\alpha _{2}-\alpha _{1}- \alpha _{2}}{\alpha _{2}} )\) by imposing the method of linearization and conclusions are presented in Table 2. We proved that about \(A (\frac{\alpha _{1}-1}{\alpha _{1}},0 )\) there exists a fold bifurcation when the parameters of the discrete model (4) are located in the following set:

We have also shown that about \(B (\frac{1}{\alpha _{2}},\frac{ \alpha _{1}\alpha _{2}-\alpha _{1}-\alpha _{2}}{\alpha _{2}} )\) the discrete-time model (4) undergoes both a Neimark–Sacker bifurcation and a period-doubling bifurcation when the parameters, respectively, are located in the following sets:

Numerical simulations have verified not only the theoretical results but also exhibited a complex dynamical behavior such as the period-2, -4, -11, -13, -15 and -22 orbits. Further, we have computed maximum Lyapunov exponents numerically. Finally, the feedback control method is applied to stabilize chaos existing in the discrete-time model (4).

Table 2 Number of equilibria along their qualitative behavior of the discrete-time model (4)

Acknowledgements

The author thanks the main editor and referees for their valuable comments and suggestions leading to improvement of this paper. This work was supported by the Higher Education Commission (HEC) of Pakistan.

Funding

The author declares that he got no funding on any part of this research.

Contributions

Corresponding author

Ethics declarations

Competing interests

The author declares that he has no conflict of interests regarding the publication of this paper.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.