# Complex Dynamics on the Routes to Chaos in a Discrete Predator-Prey System with Crowley-Martin Type Functional Response

Authors: Huayong Zhang; Shengnan Ma; Tousheng Huang; Xuebing Cong; Zichun Gao; Feifan Zhang
Journal: Discrete Dynamics in Nature and Society
 (2018)

---

## Abstract

We present in this paper an investigation on a discrete predator-prey system with Crowley-Martin type functional response to know its complex dynamics on the routes to chaos which are induced by bifurcations. Via application of the center manifold theorem and bifurcation theorems, occurrence conditions for flip bifurcation and Neimark-Sacker bifurcation are determined, respectively. Numerical simulations are performed, on the one hand, verifying the theoretical results and, on the other hand, revealing new interesting dynamical behaviors of the discrete predator-prey system, including period-doubling cascades, period-2, period-3, period-4, period-5, period-6, period-7, period-8, period-9, period-11, period-13, period-15, period-16, period-20, period-22, period-24, period-30, and period-34 orbits, invariant cycles, chaotic attractors, sub-flip bifurcation, sub-(inverse) Neimark-Sacker bifurcation, chaotic interior crisis, chaotic band, sudden disappearance of chaotic dynamics and abrupt emergence of chaos, and intermittent periodic behaviors. Moreover, three-dimensional bifurcation diagrams are utilized to study the transition between flip bifurcation and Neimark-Sacker bifurcation, and a critical case between the two bifurcations is found. This critical bifurcation case is a combination of flip bifurcation and Neimark-Sacker bifurcation, showing the nonlinear characteristics of both bifurcations.


---

## Body

## 1. Introduction

Predator-prey interaction shows widespread existence in nature and can take many forms, such as resource-consumer, plant-herbivore, and phytoplankton-zooplankton forms [1]. Due to the ubiquity and importance of the predator-prey interaction between populations, the dynamics of predator-prey systems have attracted the attention of many scholars, and the research on predator-prey systems has become one of the dominant themes in ecology [2–4]. In recent decades, mathematical models have been established to analyze various complex dynamics of the predator-prey systems in various circumstances [4–6].The Lotka-Volterra predator-prey model, proposed in the pioneering works of Lotka and Volterra [7, 8], is probably the simplest and most basic predator-prey model in the field. The development of subsequent predator-prey models is mostly based on the Lotka-Volterra model [3, 9, 10]. The development may be a change of functional response (e.g., Holling type I, II, and III classifications) or numerical response to describe different ecological processes in the predator-prey interaction [11, 12]. Several researchers argued that the nonautonomous model is closer to the realistic predator-prey system than the autonomous model [3]. Therefore, predator-prey models with a time delay were also studied because the predator-prey interaction may have time lag [13]. In order to incorporate the influence of seasonal variation into the predator-prey interaction, predator-prey models with periodically varying parameters were also proposed and investigated [14–16].In earlier approaches, most predator-prey models were time-continuous. However, the continuous models hardly described the dynamics when the populations have nonoverlapping generations or the number of populations is small [17, 18]. For these cases, applying discrete-time models should be more reasonable and appropriate [17]. In recent years, more and more scholars have payed attention to discrete-time predator-prey models. Compared with the continuous models, the discrete dynamic models often have advantages in describing richer nonlinear characteristics and complexity, such as various periodic orbits, quasiperiodic behaviors, and chaotic attractors [19–21]. Moreover, a few research works found that the discrete model can lead to more accurate results than the corresponding continuous model in describing predator-prey dynamics [6, 22, 23].Through the research of He and Lai [2], Liu and Xiao [19], Jing and Yang [20], Agiza et al. [5], Hu et al. [6], Huang and Zhang [21], and others, they found out that flip bifurcation and Neimark-Sacker bifurcation are two basic bifurcations occurring in discrete predator-prey systems. The two bifurcations often start the routes to chaos, revealing the dynamic transition from the nonchaotic state to the chaotic state. Moreover, complex dynamics on the routes to chaos can be observed, including cascades of period doubling, periodic windows, intermittent chaos, and chaotic crisis [19–21]. For example, via calculating Lyapunov exponents and fractal dimension, Agiza et al. determined diverse strange attractors of a predator-prey system [5].When we study the nonlinear characteristics of predator-prey systems, the functional response of the predator to the prey is an essential factor which affects the dynamical properties. From an ecological point of view, functional responses may be determined by the prey escape capability, prey habitat property, and predator hunting capability [24]. Skalski and Gilliam [12] suggested that three classical prey-predator-dependent functional responses—Beddington-DeAngelis, Crowley-Martin, and Hassell-Varley functional responses—can provide a better description of predator feeding over a range of predator-prey abundance presence.Particularly, Crowley-Martin functional response is much more suitable for the case where the predator feeding rate is decreased by a higher predator density even when the prey density is high [12]. The Crowley-Martin functional response depends on both the predator and the prey and takes into account the interaction between predators, regardless of whether a predator is looking for prey. It is similar to the classical Beddington-DeAngelis functional response; in particular, it shows the interspecies interference with each other. Usually, it is employed to demonstrate the predator-prey stability. And, therefore, it shows ecological significance for researching the predator-prey models with Crowley-Martin functional response. A classical predator-prey system with Crowley-Martin functional response can be described by the following equations [25]: (1a)dudt=ua-u-bv1+αu1+βv=fu,v,(1b)dvdt=v-c+du1+αu1+βv=gu,v,where u(t) and v(t) represent the densities of the prey and the predator, respectively; parameters a, b, c, d, α, and β are positive constants; the parameter a is the intrinsic growth rate of the prey and the parameter c is the mortality rate of the predator; b and d/b stand for the effect of capture rate and conversion factor denoting the newly born predators for each captured prey; the parameters α and β are the saturating parameters of Crowley-Martin functional response, α measures the magnitude of interference among prey, and β expresses the interference among the predators.In recent years, the research on the predator-prey systems with Crowley-Martin functional response has attracted a lot of attention. For example, Yang studied the stability and instability of positive equilibrium in a predator-prey model with Crowley-Martin functional response and time delay [24]. Asymptotic properties of a stochastic predator-prey model with Crowley-Martin functional response were studied by Liu et al., such as global existence and stochastic boundedness of the positive solution [26]. A stage-structured predator-prey system with Crowley-Martin functional response was investigated, and the persistence of the system and global asymptotic stability of the positive equilibrium were assessed [27]. Sivakumar et. al. analyzed a diffusive density-dependent predator-prey model with Crowley-Martin functional response and revealed Neimark-Sacker bifurcation, Turing bifurcation, and spatiotemporal patterns [28]. Dong et al. investigated the positive solutions of a spatiotemporal predator-prey system with Crowley-Martin functional response and obtained a complete understanding of the uniqueness and nonuniqueness of positive solutions [29, 30]. Li and Wu discussed the extinction and persistence results of time-dependent positive solutions to the system [31]. From previous studies, it can be seen that the predator-prey model with Crowley-Martin functional response often exhibits rich dynamics and attracts the attention of many researchers. However, the complex dynamics of the discrete-time predator-prey model with Crowley-Martin functional response are still in lack of research.The discretized form of system (1a) and (1b) was never investigated in previous studies. Since the growth, death, feeding, and migration of the predator and prey individuals always occur periodically, we can observe the dynamics of the predator-prey system by a particular time scale. This time scale can be defined by the generation span of the predator and prey populations and it measures the regeneration time of both populations. Its value is mainly determined by the population types; on the other hand, it is influenced by the change of population size or environmental conditions. In this research, we denote the time scale on which predator-prey dynamics are described and observed as parameter τ. Applying the forward Euler scheme with time interval τ, system (1a) and (1b) is now transformed to a discrete predator-prey system: (2a)uk+1=fuk,vk=uk+τuka-uk-bvk1+αuk1+βvk,(2b)vk+1=guk,vk=vk+τvk-c+duk1+αuk1+βvk,where k represents the sequence number of iterations. If we give the initial time as t0, the kth iteration represents the time at t0+kτ.In recent decades, the research on discrete predator-prey systems has become an important topic. Via the research of Li and Wu [31], the dynamical properties of system (1a) and (1b) have been revealed a lot. However, nonlinear characteristics of discrete system (2a) and (2b) are still scarcely known. In the literature, Ren et al. have actually investigated a Leslie-Gower type discrete predator-prey model with Crowley-Martin functional response and found the existence of Marotto chaos [32]. They also controlled chaotic orbits to be a fixed point by a feedback control method. Different from the study of Ren et al. [32], this research focuses on a detailed exploration of the complex dynamics, such as periodic windows, subbifurcations, and chaos crisis, on the routes to chaos induced by the bifurcations of system (2a) and (2b).In order to facilitate the study of system (2a) and (2b), we rewrite the equations of system (2a) and (2b) as a map; that is, (3)uv⟶u+τua-u-bv1+αu1+βvv+τv-c+du1+αu1+βv.In this research, the bifurcations and complex dynamics of the discrete predator-prey system with Crowley-Martin functional response are explored through map (3). The exploration is arranged as follows. Section 2 will give the stability analysis on the fixed points of map (3). Section 3 will provide bifurcation analysis and determine occurrence conditions for flip bifurcation and Neimark-Sacker bifurcation. Section 4 will present numerical simulations for complex dynamic behaviors of the discrete predator-prey system. Section 5 will describe the conclusions.


## 2. Stability Analysis

The fixed points of map (3) can be calculated through the following equations: (4a)u=u+τua-u-bv1+αu1+βv,(4b)v=v+τv-c+du1+αu1+βv,which directly yield (5)u0,v0:0,0;u1,v1:a,0;u2,v2:u∗,v∗,where u∗ is a positive root of the following cubic equations: (6)αβdu3+βd1-αau2+bd-βad-αcbu-cb=0,v∗=da-u∗u∗cb.According to previous results in the literature [21], the sufficient conditions for map (3) possessing a unique positive fixed point are determined as follows: (7)1-αa2βd1-αa2-18αβd1-αa-9αbd-βad-αcd-27α2cd<0.To determine the stability of the fixed points of map (3), the Jacobian matrix associated with the map is applied and can be described as follows: (8)Ju,v=1+τa-2u-bv1+αu21+βv-τbu1+αu1+βv2τdv1+αu21+βv1+τ-c+du1+αu1+βv2.For the Jacobian matrix of (8), we substitute the values of the three fixed points into it and calculate the two eigenvalues, λ1 and λ2. If λ1<1 and λ2<1, the corresponding fixed point is considered to be stable; if λ1>1 or λ2>1, an unstable fixed point emerges. For the fixed point (u0,v0), the corresponding Jacobian matrix is (9)Ju0,v0=1+aτ001-cτ,the two eigenvalues of which are 1+aτ and 1-cτ. Since 1+aτ is larger than 1, then (u0,v0) is unstable. For the fixed point (u1,v1), the Jacobian matrix is obtained as (10)Ju1,v1=1-aτ-τab1+αa01+τ-c+da1+αa.The two eigenvalues of Ju1,v1 are 1-aτ and 1-cτ+daτ/(1+αa). When the following conditions are established,(11)c+αac-ad>0,0<τ<min2a,2c1+αa-da1+αa,the two eigenvalues are both less than one. Therefore, under conditions (11), (u1,v1) is stable. The Jacobian matrix of (u2,v2) is presented as the following matrix: (12)Ju2,v2=a11a12a21a22=1+τa-2u∗-τbda-u∗u∗1+αu∗2cb+βdu∗a-u∗-τc2b3u∗1+αu∗cb+βdu∗a-u∗2τd2u∗a-u∗1+αu∗2cb+βdu∗a-u∗1-τc+τb2c2du∗1+αu∗cb+βdu∗a-u∗2.Likewise, the two eigenvalues of matrix (3) are(13)λ1,2=12-p±p2-4q,where (14)p=-2-τb∗,q=1+τb∗+τ2a∗,a∗=a-2u∗-bda-u∗u∗1+αu∗2cb+βda-u∗u∗-c+b2c2du∗1+αu∗cb+βda-u∗u∗2+b3c2d2a-u∗u∗21+αu∗3cb+βda-u∗u∗3,b∗=a-2u∗-c-bda-u∗u∗1+αu∗2cb+βda-u∗u∗+b2c2du∗1+αu∗cb+βda-u∗u∗2.We calculate the stability conditions of (u2,v2), that is, λ1<1 and λ2<1, as follows: (15)q<1,-1+q<p<1+q.


## 3. Bifurcation Analysis

### 3.1. Flip Bifurcation

Around the stable fixed point(u2,v2), flip bifurcation of map (3) is studied via the flip bifurcation theorem and center manifold theorem [27]. According to the flip bifurcation theorem, the predator-prey system undergoes flip bifurcation when one of the two eigenvalues of (u2,v2) is equal to −1. Therefore, it is explicitly said that q=p-1; that is, (16)1+τb∗+τ2a∗=-2-τb∗-1,which directly yields (17)τ∗=4-b∗∓b∗2-4a∗.At such critical point τ∗, the two eigenvalues of (u2,v2) are λ1=1 and λ2=1-p. Moreover, the value of λ2 should satisfy λ2≠1; that is,(18)τ∗c+2u∗-a+bda-u∗u∗1+αu∗2cb+βda-u∗u∗-b2c2du∗1+αu∗cb+βda-u∗u∗2≠2,4.Consider the bifurcation parameterτ also as a dependent variable. We then translate the fixed point (u2,v2,τ∗) to the origin via the following translation: (19)w=u-u∗,z=v-v∗,τ~=τ-τ∗.Consequently, via Taylor expansion near the fixed point, map (3) can be transformed as (20)wzτ~⟶a11w+a12z+12a13w2+12a14z2+a15wz+a16wτ~+a17zτ~+16a18w3+12a19wz2+12a110w2z+16a111z3+12a112w2τ~+12a113z2τ~+a114wzτ~+ow+z+τ~4a21w+a22z+12a23w2+12a24z2+a25wz+a26wτ~+a27zτ~+16a28w3+12a29wz2+12a210w2z+16a211z3+12a212w2τ~+12a213z2τ~+a214wzτ~+ow+z+τ~4τ~,where ow+z+τ~4 is a polynomial with at least four orders with variables w, z, and τ~. The coefficients aij can be calculated by (21)a1j=∂s1f∂us1∂s2f∂vs2∂s3f∂τs3u,v,τ=u∗,v∗,τ∗,a2j=∂s1g∂us1∂s2g∂vs2∂s3g∂τs3u,v,τ=u∗,v∗,τ∗,where s1, s2, and s3 are the orders of the variables w, z, and τ~, respectively, in the corresponding term of aij. Here, we define ∂0g/∂u0=1.Applying a reversible transformation as(22)wz⟶a12a12-1-a11λ2-a11w~z~,map (20) becomes (23)w~z~τ~⟶-1000λ20001w~z~τ~+1a121+λ2F1w~,z~,τ~F2w~,z~,τ~0,where (24)
F
1
w
~
,
z
~
,
τ
~
=
a
12
2
a
13
λ
2
-
a
11
2
-
a
12
a
23
2
w
~
+
z
~
2
+
a
14
λ
2
-
a
11
2
-
a
12
a
24
2
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
2
+
a
12
a
15
λ
2
-
a
11
-
a
12
a
25
w
~
+
z
~
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
+
a
12
a
16
λ
2
-
a
11
-
a
12
a
26
w
~
+
z
~
τ
~
+
a
17
λ
2
-
a
11
-
a
12
a
27
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
τ
~
+
a
12
3
a
18
λ
2
-
a
11
6
-
a
12
a
28
6
w
~
+
z
~
3
+
a
12
a
19
λ
2
-
a
11
2
-
a
12
a
29
2
w
~
+
z
~
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
2
+
a
12
2
a
110
λ
2
-
a
11
2
-
a
12
a
210
2
w
~
+
z
~
2
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
+
a
111
λ
2
-
a
11
6
-
a
12
a
211
6
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
3
+
a
12
2
a
112
λ
2
-
a
11
2
-
a
12
a
212
2
w
~
+
z
~
2
τ
~
+
a
113
λ
2
-
a
11
2
-
a
12
a
213
2
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
2
τ
~
+
a
12
a
114
λ
2
-
a
11
-
a
12
a
214
w
~
+
z
~
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
τ
~
+
o
w
~
+
z
~
+
τ
~
4
,
(25)
F
2
w
~
,
z
~
,
τ
~
=
a
12
2
a
13
1
+
a
11
2
+
a
12
a
23
2
w
~
+
z
~
2
+
a
14
1
+
a
11
2
+
a
12
a
24
2
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
2
+
a
12
a
15
1
+
a
11
+
a
12
a
25
w
~
+
z
~
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
+
a
12
a
16
;
21
+
a
11
+
a
12
a
26
w
~
+
z
~
τ
~
+
a
17
1
+
a
11
+
a
12
a
27
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
τ
~
+
a
12
3
a
18
1
+
a
11
6
+
a
12
a
28
6
w
~
+
z
~
3
+
a
12
a
19
1
+
a
11
2
+
a
12
a
29
2
w
~
+
z
~
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
2
+
a
12
2
a
110
1
+
a
11
2
+
a
12
a
210
2
w
~
+
z
~
2
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
+
a
111
1
+
a
11
6
+
a
12
a
211
6
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
3
+
a
12
2
a
112
1
+
a
11
2
+
a
12
a
212
2
w
~
+
z
~
2
τ
~
+
a
113
1
+
a
11
2
+
a
12
a
213
2
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
2
τ
~
+
a
12
a
114
1
+
a
11
+
a
12
a
214
w
~
+
z
~
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
τ
~
+
o
w
~
+
z
~
+
τ
~
4
.The center manifold of map (23) at the fixed point (0,0,0) is then determined. According to the center manifold theorem, a central manifold at the fixed point (0,0,0) can be approximated by the following: (26)WC0,0,0=w~,z~,τ~∈R3∣z~=h∗w~,τ~,h∗0,0=0,Dh∗0,0=0,where h∗w~,τ~ is assumed to be(27)h∗w~,τ~=e0τ~+e1w~2+e2w~τ~+e3τ~2+ow~+τ~3.Utilizing map (23) on both sides of z~=h∗w~,τ~ leads to (28)e0τ~+e1-w~+F1w~,h∗w~,τ~,τ~a121+λ22+e2τ~-w~+F1w~,h∗w~,τ~,τ~a121+λ2+e3τ~2-λ2e0τ~+e1w~2+e2w~τ~+e3τ~2-F2w~,h∗w~,τ~,τ~a121+λ2=ow~+τ~3,simultaneously balancing the variables in (28); then, we get v e0=e3=0, and the expressions of e1 and e2 are (29)e1=12a121-λ22a122a131+a11+a12a23+a141+a11+a12a241+a112-2a12a151+a11+a12a251+a11,e2=1a121+λ22a171+a112+a12a27-a161+a11-a122a26.Considering the dynamics of map (23) restricted in the center manifold WC(0,0,0), hence we obtain a one-dimensional mapping as (30)F:w~⟶-w~+μ1w~2+μ2w~τ~+μ3w~2τ~+μ4w~τ~2+μ5w~3+ow~+τ~4.The coefficients in map (30) are described as follows: (31a)μ1=1a121+λ2a1222a13λ2-a11-a12a23+12a14λ2-a11-a12a241+a112-a12a15λ2-a11-a12a251+a11,(31b)μ2=11+λ2a16λ2-a11-a12a26-1+a11a121+λ2a17λ2-a11-a12a27,(31c)μ3=1a121+λ2a122a13λ2-a11-a12a23e2-a14λ2-a11-a12a24λ2-a111+a11e2+a12a15λ2-a11-a12a25λ2-2a11-1e2+a12a16λ2-a11-a12a26e1+a17λ2-a11-a12a27λ2-a11e1+a122a112λ2-a112-a12a2122+a113λ2-a112-a12a21321+a112-a12a114λ2-a11-a12a2141+a11,(31d)μ4=e21+λ2a16λ2-a11-a12a26+e2λ2-a11a121+λ2a17λ2-a11-a12a27,(31e)μ5=1a121+λ2a122a13λ2-a11-a12a23e1-a14λ2-a11-a12a24λ2-a111+a11e1+a12a15λ2-a11-a12a25λ2-2a11-1e1+a123a18λ2-a116-a12a286+a12a19λ2-a112-a12a2921+a112-a122a110λ2-a112-a12a21021+a11-a111λ2-a116-a12a21161+a113.On the basis of map (30), the predator-prey system experiencing flip bifurcation needs the requirement of the following two conditions:(32)η1=∂2F∂w~∂τ~+12∂F∂τ~∂2F∂w~2≠0,η2=16∂3F∂w~3+12∂2F∂w~22≠0,w~=0,τ~=0,which are equivalent to (33)η1=μ2≠0,η2=μ5+μ12≠0.Hence, the occurrence of flip bifurcation at the fixed point(u2,v2) with parameter τ varying in a small neighborhood of τ∗ needs the satisfaction of conditions (18) and (33). Furthermore, if η2>0, a period-2 orbit bifurcates from (u2,v2) and is stable; if η2<0, the period-2 orbit bifurcating from (u2,v2) is unstable.


### 3.2. Neimark-Sacker Bifurcation

When map (3) undergoes Neimark-Sacker bifurcation, the two eigenvalues of the fixed point (u2,v2) are a pair of conjugate complex numbers whose module is one. Hence, we directly have p2-4q<0 and q=1, which are further determined as (34)b∗2-4a∗<0,τ0=-b∗a∗.Under conditions (34), we translate the fixed point (u2,v2) to the origin by the following translation: (35)w=u-u∗,z=v-v∗.With application of Taylor expansion, map (3) is transformed to (36)wz⟶a11w+a12z+12a13w2+12a14z2+a15wz+16a18w3+12a19wz2+12a110w2z+16a111z3+ow+z4a21w+a22z+12a23w2+12a24z2+a25wz+16a28w3+12a29wz2+12a210w2z+16a211z3+ow+z4.The coefficients in the map are also calculated by (21) but u,v,τ=u∗,v∗,τ∗. The two eigenvalues at the fixed point (0,0) of map (36) are also a pair of conjugate complexes with a module of one. These two eigenvalues are written as (37)λτ0,λ-τ0=-pτ02±i~24qτ0-p2τ0=θ±i~φ,where i~=-1 and λ=q(τ0)=1. The modulus of the two eigenvalues satisfies(38)γ=dλτdττ=τ0=121a∗τ2+b∗τ+12a∗τ+b∗τ=τ0=-b∗/a∗=-b∗2≠0. Simultaneously, the two eigenvalues should not be real numbers or imaginary numbers, which means λτ0n≠1, n=1,2,3,4. Since qτ0=1, we have pτ0≠-2,0,1,2. In the meantime, p2τ0<4qτ0=4. Therefore, pτ0≠0,1, which is equivalent to (39)τ0c+2u∗-a+bda-u∗u∗1+αu∗2cb+βda-u∗u∗-b2c2du∗1+αu∗cb+βda-u∗u∗2≠2,3.Through the transformation(40)wz⟶a110θ-a11-φw~z~,the normal form of map (36) is obtained as (41)w~z~⟶θ-φφθw~z~+1a12φG1w~,z~G2w~,z~,in which (42a)
G
1
w
~
,
z
~
=
a
12
2
a
13
2
φ
+
a
14
2
φ
θ
-
a
11
2
+
a
12
a
15
θ
-
a
11
φ
w
~
2
+
a
14
2
φ
3
z
2
-
a
14
θ
-
a
11
φ
2
+
a
12
a
15
φ
2
w
~
z
~
+
a
12
3
a
18
6
+
a
12
a
19
2
φ
θ
-
a
11
2
+
a
12
2
a
110
2
φ
θ
-
a
11
+
1
6
a
111
φ
θ
-
a
11
3
w
~
3
+
a
12
a
19
2
φ
3
+
a
111
2
φ
3
θ
-
a
11
w
~
z
~
2
-
a
12
a
19
θ
-
a
11
φ
2
+
a
12
2
a
110
2
φ
2
+
a
111
2
θ
-
a
11
2
φ
2
w
~
2
z
~
-
a
111
6
φ
4
z
~
3
+
o
w
~
+
z
~
4
,
(42b)
G
2
w
~
,
z
~
=
a
12
2
2
a
13
θ
-
a
11
-
a
12
a
23
+
1
2
a
14
θ
-
a
11
-
a
12
a
24
θ
-
a
11
2
+
a
12
θ
-
a
11
a
15
θ
-
a
11
-
a
12
a
25
w
~
2
+
a
14
θ
-
a
11
2
-
a
12
a
24
2
φ
2
z
~
2
-
a
14
θ
-
a
11
-
a
12
a
24
θ
-
a
11
φ
+
a
12
a
15
θ
-
a
11
-
a
12
a
25
φ
w
~
z
~
+
a
12
3
6
a
18
θ
-
a
11
-
a
12
a
28
+
a
12
2
a
19
θ
-
a
11
-
a
12
a
29
θ
-
a
11
2
+
a
12
2
a
110
θ
-
a
11
2
-
a
12
a
210
2
θ
-
a
11
+
a
111
θ
-
a
11
6
-
a
12
a
211
6
θ
-
a
11
3
w
~
3
+
a
12
2
a
19
θ
-
a
11
-
a
12
a
29
φ
2
+
a
111
θ
-
a
11
2
-
a
12
a
211
6
θ
-
a
11
φ
2
w
~
z
~
2
-
a
12
a
19
θ
-
a
11
-
a
12
a
29
θ
-
a
11
φ
+
a
12
2
a
110
θ
-
a
11
2
-
a
12
a
210
2
φ
+
a
111
θ
-
a
11
2
-
a
12
a
211
2
θ
-
a
11
2
φ
w
~
2
z
~
-
a
111
θ
-
a
11
6
-
a
12
a
211
6
φ
3
z
~
3
+
o
w
~
+
z
~
4
.Calculating the two- and three-order derivatives ofG1(w~,z~) and G2(w~,z~) at w~=0 and τ~=0, the following results are obtained: (43)G1w~w~=a122a132φ+a142φθ-a112+a12a15θ-a11φ,G1w~z~=-a14θ-a11φ2+a12a15φ2,G1z~z~=a142φ3,G1w~w~w~=a123a186+a12a192φθ-a112+a122a1102φθ-a11+16a111φθ-a113,G1w~w~z~=-a12a19θ-a11φ2+a122a1102φ2+a1112θ-a112φ2,G1w~z~z~=a12a192φ3+a1112φ3θ-a11,G1z~z~z~=-a1116φ4,G2w~w~=a1222a13θ-a11-a12a23+12a14θ-a11-a12a24θ-a112+a12θ-a11a15θ-a11-a12a25,G2w~z~=-a14θ-a11-a12a24θ-a11φ+a12a15θ-a11-a12a25φ,G2z~z~=a14θ-a112-a12a242φ2,G2w~w~w~=a123a18θ-a116-a12a286+a12a19θ-a112-a12a292θ-a112+a122a110θ-a112-a12a2102θ-a11+a111θ-a116-a12a2116θ-a113,G2w~w~z~=-a12a19θ-a11-a12a29θ-a11φ+a122a110θ-a112-a12a2102φ+a111θ-a112-a12a2112θ-a112φ,G2w~z~z~=a122a19θ-a11-a12a29φ2+a111θ-a112-a12a2112α-a11φ2,G2z~z~z~=-a111θ-a116-a12a2116φ3.On the basis of map (41), the predator-prey system experiencing Neimark-Sacker bifurcation needs the requirement of the following condition:(44)a0=-Re⁡1-2λ-λ-21-λξ11ξ20-12ξ112-ξ022+Re⁡λ-ξ21≠0,where (45)ξ20=18G1w~w~-G1z~z~+2G2w~z~+i~G2w~w~-G2z~z~-2G1w~z~,ξ11=14G1w~w~+G1z~z~+i~G2w~w~+G2z~z~,ξ02=18G1w~w~-G1z~z~-2G2w~z~+i~G2w~w~-G2z~z~+2G1w~z~,ξ21=116G1w~w~w~+G1w~z~z~+G2w~w~z~+G2z~z~z~+i~G2w~w~w~+G2w~z~z~-G1w~w~z~-G1z~z~z~.Using the above formula, we can get (46)a0=-1A0M0G1w~w~+G1z~z~G1w~w~-G1z~z~+2G2w~z~-G2w~w~+G2z~z~G2w~w~-G2z~z~-2G1w~z~+N0G2w~w~+G2z~z~G1w~w~-G1z~z~+2G2w~z~-G1w~w~+G1z~z~G2w~w~-G2z~z~-2G1w~z~-132a122φ2G1w~w~+G1z~z~2+G2w~w~+G2z~z~2-164a122φ2G1w~w~-G1z~z~-2G2w~z~2+G2w~w~-G2z~z~+2G1w~z~2+116a12φαG1w~w~w~+G1w~z~z~+G2w~w~z~+G2z~z~z~+φG2w~w~w~+G2w~z~z~-G1w~w~z~-G1z~z~z~≠0,in which (47)A0=32a122φ21-θ2+φ2,N0=2θφ1-θ1-2θ-2φ2-φ3-4θθ2-φ2,M0=1-3θ+2θ2-2φ2θ2-φ2+6θ-8θ2φ2.Taking the above calculations together, the occurrence of Neimark-Sacker bifurcation at the fixed point(u2,v2) in the discrete predator-prey system needs the satisfaction of conditions (34), (39), and (46). Moreover, if a0<0 and γ>0, an attracting invariant closed curve bifurcates from the fixed point for τ>τ0; if a0>0 and γ>0, a repelling invariant closed curve bifurcates from the fixed point for τ<τ0.


## 3.1. Flip Bifurcation

Around the stable fixed point(u2,v2), flip bifurcation of map (3) is studied via the flip bifurcation theorem and center manifold theorem [27]. According to the flip bifurcation theorem, the predator-prey system undergoes flip bifurcation when one of the two eigenvalues of (u2,v2) is equal to −1. Therefore, it is explicitly said that q=p-1; that is, (16)1+τb∗+τ2a∗=-2-τb∗-1,which directly yields (17)τ∗=4-b∗∓b∗2-4a∗.At such critical point τ∗, the two eigenvalues of (u2,v2) are λ1=1 and λ2=1-p. Moreover, the value of λ2 should satisfy λ2≠1; that is,(18)τ∗c+2u∗-a+bda-u∗u∗1+αu∗2cb+βda-u∗u∗-b2c2du∗1+αu∗cb+βda-u∗u∗2≠2,4.Consider the bifurcation parameterτ also as a dependent variable. We then translate the fixed point (u2,v2,τ∗) to the origin via the following translation: (19)w=u-u∗,z=v-v∗,τ~=τ-τ∗.Consequently, via Taylor expansion near the fixed point, map (3) can be transformed as (20)wzτ~⟶a11w+a12z+12a13w2+12a14z2+a15wz+a16wτ~+a17zτ~+16a18w3+12a19wz2+12a110w2z+16a111z3+12a112w2τ~+12a113z2τ~+a114wzτ~+ow+z+τ~4a21w+a22z+12a23w2+12a24z2+a25wz+a26wτ~+a27zτ~+16a28w3+12a29wz2+12a210w2z+16a211z3+12a212w2τ~+12a213z2τ~+a214wzτ~+ow+z+τ~4τ~,where ow+z+τ~4 is a polynomial with at least four orders with variables w, z, and τ~. The coefficients aij can be calculated by (21)a1j=∂s1f∂us1∂s2f∂vs2∂s3f∂τs3u,v,τ=u∗,v∗,τ∗,a2j=∂s1g∂us1∂s2g∂vs2∂s3g∂τs3u,v,τ=u∗,v∗,τ∗,where s1, s2, and s3 are the orders of the variables w, z, and τ~, respectively, in the corresponding term of aij. Here, we define ∂0g/∂u0=1.Applying a reversible transformation as(22)wz⟶a12a12-1-a11λ2-a11w~z~,map (20) becomes (23)w~z~τ~⟶-1000λ20001w~z~τ~+1a121+λ2F1w~,z~,τ~F2w~,z~,τ~0,where (24)
F
1
w
~
,
z
~
,
τ
~
=
a
12
2
a
13
λ
2
-
a
11
2
-
a
12
a
23
2
w
~
+
z
~
2
+
a
14
λ
2
-
a
11
2
-
a
12
a
24
2
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
2
+
a
12
a
15
λ
2
-
a
11
-
a
12
a
25
w
~
+
z
~
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
+
a
12
a
16
λ
2
-
a
11
-
a
12
a
26
w
~
+
z
~
τ
~
+
a
17
λ
2
-
a
11
-
a
12
a
27
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
τ
~
+
a
12
3
a
18
λ
2
-
a
11
6
-
a
12
a
28
6
w
~
+
z
~
3
+
a
12
a
19
λ
2
-
a
11
2
-
a
12
a
29
2
w
~
+
z
~
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
2
+
a
12
2
a
110
λ
2
-
a
11
2
-
a
12
a
210
2
w
~
+
z
~
2
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
+
a
111
λ
2
-
a
11
6
-
a
12
a
211
6
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
3
+
a
12
2
a
112
λ
2
-
a
11
2
-
a
12
a
212
2
w
~
+
z
~
2
τ
~
+
a
113
λ
2
-
a
11
2
-
a
12
a
213
2
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
2
τ
~
+
a
12
a
114
λ
2
-
a
11
-
a
12
a
214
w
~
+
z
~
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
τ
~
+
o
w
~
+
z
~
+
τ
~
4
,
(25)
F
2
w
~
,
z
~
,
τ
~
=
a
12
2
a
13
1
+
a
11
2
+
a
12
a
23
2
w
~
+
z
~
2
+
a
14
1
+
a
11
2
+
a
12
a
24
2
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
2
+
a
12
a
15
1
+
a
11
+
a
12
a
25
w
~
+
z
~
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
+
a
12
a
16
;
21
+
a
11
+
a
12
a
26
w
~
+
z
~
τ
~
+
a
17
1
+
a
11
+
a
12
a
27
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
τ
~
+
a
12
3
a
18
1
+
a
11
6
+
a
12
a
28
6
w
~
+
z
~
3
+
a
12
a
19
1
+
a
11
2
+
a
12
a
29
2
w
~
+
z
~
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
2
+
a
12
2
a
110
1
+
a
11
2
+
a
12
a
210
2
w
~
+
z
~
2
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
+
a
111
1
+
a
11
6
+
a
12
a
211
6
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
3
+
a
12
2
a
112
1
+
a
11
2
+
a
12
a
212
2
w
~
+
z
~
2
τ
~
+
a
113
1
+
a
11
2
+
a
12
a
213
2
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
2
τ
~
+
a
12
a
114
1
+
a
11
+
a
12
a
214
w
~
+
z
~
λ
2
-
a
11
z
~
-
1
+
a
11
w
~
τ
~
+
o
w
~
+
z
~
+
τ
~
4
.The center manifold of map (23) at the fixed point (0,0,0) is then determined. According to the center manifold theorem, a central manifold at the fixed point (0,0,0) can be approximated by the following: (26)WC0,0,0=w~,z~,τ~∈R3∣z~=h∗w~,τ~,h∗0,0=0,Dh∗0,0=0,where h∗w~,τ~ is assumed to be(27)h∗w~,τ~=e0τ~+e1w~2+e2w~τ~+e3τ~2+ow~+τ~3.Utilizing map (23) on both sides of z~=h∗w~,τ~ leads to (28)e0τ~+e1-w~+F1w~,h∗w~,τ~,τ~a121+λ22+e2τ~-w~+F1w~,h∗w~,τ~,τ~a121+λ2+e3τ~2-λ2e0τ~+e1w~2+e2w~τ~+e3τ~2-F2w~,h∗w~,τ~,τ~a121+λ2=ow~+τ~3,simultaneously balancing the variables in (28); then, we get v e0=e3=0, and the expressions of e1 and e2 are (29)e1=12a121-λ22a122a131+a11+a12a23+a141+a11+a12a241+a112-2a12a151+a11+a12a251+a11,e2=1a121+λ22a171+a112+a12a27-a161+a11-a122a26.Considering the dynamics of map (23) restricted in the center manifold WC(0,0,0), hence we obtain a one-dimensional mapping as (30)F:w~⟶-w~+μ1w~2+μ2w~τ~+μ3w~2τ~+μ4w~τ~2+μ5w~3+ow~+τ~4.The coefficients in map (30) are described as follows: (31a)μ1=1a121+λ2a1222a13λ2-a11-a12a23+12a14λ2-a11-a12a241+a112-a12a15λ2-a11-a12a251+a11,(31b)μ2=11+λ2a16λ2-a11-a12a26-1+a11a121+λ2a17λ2-a11-a12a27,(31c)μ3=1a121+λ2a122a13λ2-a11-a12a23e2-a14λ2-a11-a12a24λ2-a111+a11e2+a12a15λ2-a11-a12a25λ2-2a11-1e2+a12a16λ2-a11-a12a26e1+a17λ2-a11-a12a27λ2-a11e1+a122a112λ2-a112-a12a2122+a113λ2-a112-a12a21321+a112-a12a114λ2-a11-a12a2141+a11,(31d)μ4=e21+λ2a16λ2-a11-a12a26+e2λ2-a11a121+λ2a17λ2-a11-a12a27,(31e)μ5=1a121+λ2a122a13λ2-a11-a12a23e1-a14λ2-a11-a12a24λ2-a111+a11e1+a12a15λ2-a11-a12a25λ2-2a11-1e1+a123a18λ2-a116-a12a286+a12a19λ2-a112-a12a2921+a112-a122a110λ2-a112-a12a21021+a11-a111λ2-a116-a12a21161+a113.On the basis of map (30), the predator-prey system experiencing flip bifurcation needs the requirement of the following two conditions:(32)η1=∂2F∂w~∂τ~+12∂F∂τ~∂2F∂w~2≠0,η2=16∂3F∂w~3+12∂2F∂w~22≠0,w~=0,τ~=0,which are equivalent to (33)η1=μ2≠0,η2=μ5+μ12≠0.Hence, the occurrence of flip bifurcation at the fixed point(u2,v2) with parameter τ varying in a small neighborhood of τ∗ needs the satisfaction of conditions (18) and (33). Furthermore, if η2>0, a period-2 orbit bifurcates from (u2,v2) and is stable; if η2<0, the period-2 orbit bifurcating from (u2,v2) is unstable.


## 3.2. Neimark-Sacker Bifurcation

When map (3) undergoes Neimark-Sacker bifurcation, the two eigenvalues of the fixed point (u2,v2) are a pair of conjugate complex numbers whose module is one. Hence, we directly have p2-4q<0 and q=1, which are further determined as (34)b∗2-4a∗<0,τ0=-b∗a∗.Under conditions (34), we translate the fixed point (u2,v2) to the origin by the following translation: (35)w=u-u∗,z=v-v∗.With application of Taylor expansion, map (3) is transformed to (36)wz⟶a11w+a12z+12a13w2+12a14z2+a15wz+16a18w3+12a19wz2+12a110w2z+16a111z3+ow+z4a21w+a22z+12a23w2+12a24z2+a25wz+16a28w3+12a29wz2+12a210w2z+16a211z3+ow+z4.The coefficients in the map are also calculated by (21) but u,v,τ=u∗,v∗,τ∗. The two eigenvalues at the fixed point (0,0) of map (36) are also a pair of conjugate complexes with a module of one. These two eigenvalues are written as (37)λτ0,λ-τ0=-pτ02±i~24qτ0-p2τ0=θ±i~φ,where i~=-1 and λ=q(τ0)=1. The modulus of the two eigenvalues satisfies(38)γ=dλτdττ=τ0=121a∗τ2+b∗τ+12a∗τ+b∗τ=τ0=-b∗/a∗=-b∗2≠0. Simultaneously, the two eigenvalues should not be real numbers or imaginary numbers, which means λτ0n≠1, n=1,2,3,4. Since qτ0=1, we have pτ0≠-2,0,1,2. In the meantime, p2τ0<4qτ0=4. Therefore, pτ0≠0,1, which is equivalent to (39)τ0c+2u∗-a+bda-u∗u∗1+αu∗2cb+βda-u∗u∗-b2c2du∗1+αu∗cb+βda-u∗u∗2≠2,3.Through the transformation(40)wz⟶a110θ-a11-φw~z~,the normal form of map (36) is obtained as (41)w~z~⟶θ-φφθw~z~+1a12φG1w~,z~G2w~,z~,in which (42a)
G
1
w
~
,
z
~
=
a
12
2
a
13
2
φ
+
a
14
2
φ
θ
-
a
11
2
+
a
12
a
15
θ
-
a
11
φ
w
~
2
+
a
14
2
φ
3
z
2
-
a
14
θ
-
a
11
φ
2
+
a
12
a
15
φ
2
w
~
z
~
+
a
12
3
a
18
6
+
a
12
a
19
2
φ
θ
-
a
11
2
+
a
12
2
a
110
2
φ
θ
-
a
11
+
1
6
a
111
φ
θ
-
a
11
3
w
~
3
+
a
12
a
19
2
φ
3
+
a
111
2
φ
3
θ
-
a
11
w
~
z
~
2
-
a
12
a
19
θ
-
a
11
φ
2
+
a
12
2
a
110
2
φ
2
+
a
111
2
θ
-
a
11
2
φ
2
w
~
2
z
~
-
a
111
6
φ
4
z
~
3
+
o
w
~
+
z
~
4
,
(42b)
G
2
w
~
,
z
~
=
a
12
2
2
a
13
θ
-
a
11
-
a
12
a
23
+
1
2
a
14
θ
-
a
11
-
a
12
a
24
θ
-
a
11
2
+
a
12
θ
-
a
11
a
15
θ
-
a
11
-
a
12
a
25
w
~
2
+
a
14
θ
-
a
11
2
-
a
12
a
24
2
φ
2
z
~
2
-
a
14
θ
-
a
11
-
a
12
a
24
θ
-
a
11
φ
+
a
12
a
15
θ
-
a
11
-
a
12
a
25
φ
w
~
z
~
+
a
12
3
6
a
18
θ
-
a
11
-
a
12
a
28
+
a
12
2
a
19
θ
-
a
11
-
a
12
a
29
θ
-
a
11
2
+
a
12
2
a
110
θ
-
a
11
2
-
a
12
a
210
2
θ
-
a
11
+
a
111
θ
-
a
11
6
-
a
12
a
211
6
θ
-
a
11
3
w
~
3
+
a
12
2
a
19
θ
-
a
11
-
a
12
a
29
φ
2
+
a
111
θ
-
a
11
2
-
a
12
a
211
6
θ
-
a
11
φ
2
w
~
z
~
2
-
a
12
a
19
θ
-
a
11
-
a
12
a
29
θ
-
a
11
φ
+
a
12
2
a
110
θ
-
a
11
2
-
a
12
a
210
2
φ
+
a
111
θ
-
a
11
2
-
a
12
a
211
2
θ
-
a
11
2
φ
w
~
2
z
~
-
a
111
θ
-
a
11
6
-
a
12
a
211
6
φ
3
z
~
3
+
o
w
~
+
z
~
4
.Calculating the two- and three-order derivatives ofG1(w~,z~) and G2(w~,z~) at w~=0 and τ~=0, the following results are obtained: (43)G1w~w~=a122a132φ+a142φθ-a112+a12a15θ-a11φ,G1w~z~=-a14θ-a11φ2+a12a15φ2,G1z~z~=a142φ3,G1w~w~w~=a123a186+a12a192φθ-a112+a122a1102φθ-a11+16a111φθ-a113,G1w~w~z~=-a12a19θ-a11φ2+a122a1102φ2+a1112θ-a112φ2,G1w~z~z~=a12a192φ3+a1112φ3θ-a11,G1z~z~z~=-a1116φ4,G2w~w~=a1222a13θ-a11-a12a23+12a14θ-a11-a12a24θ-a112+a12θ-a11a15θ-a11-a12a25,G2w~z~=-a14θ-a11-a12a24θ-a11φ+a12a15θ-a11-a12a25φ,G2z~z~=a14θ-a112-a12a242φ2,G2w~w~w~=a123a18θ-a116-a12a286+a12a19θ-a112-a12a292θ-a112+a122a110θ-a112-a12a2102θ-a11+a111θ-a116-a12a2116θ-a113,G2w~w~z~=-a12a19θ-a11-a12a29θ-a11φ+a122a110θ-a112-a12a2102φ+a111θ-a112-a12a2112θ-a112φ,G2w~z~z~=a122a19θ-a11-a12a29φ2+a111θ-a112-a12a2112α-a11φ2,G2z~z~z~=-a111θ-a116-a12a2116φ3.On the basis of map (41), the predator-prey system experiencing Neimark-Sacker bifurcation needs the requirement of the following condition:(44)a0=-Re⁡1-2λ-λ-21-λξ11ξ20-12ξ112-ξ022+Re⁡λ-ξ21≠0,where (45)ξ20=18G1w~w~-G1z~z~+2G2w~z~+i~G2w~w~-G2z~z~-2G1w~z~,ξ11=14G1w~w~+G1z~z~+i~G2w~w~+G2z~z~,ξ02=18G1w~w~-G1z~z~-2G2w~z~+i~G2w~w~-G2z~z~+2G1w~z~,ξ21=116G1w~w~w~+G1w~z~z~+G2w~w~z~+G2z~z~z~+i~G2w~w~w~+G2w~z~z~-G1w~w~z~-G1z~z~z~.Using the above formula, we can get (46)a0=-1A0M0G1w~w~+G1z~z~G1w~w~-G1z~z~+2G2w~z~-G2w~w~+G2z~z~G2w~w~-G2z~z~-2G1w~z~+N0G2w~w~+G2z~z~G1w~w~-G1z~z~+2G2w~z~-G1w~w~+G1z~z~G2w~w~-G2z~z~-2G1w~z~-132a122φ2G1w~w~+G1z~z~2+G2w~w~+G2z~z~2-164a122φ2G1w~w~-G1z~z~-2G2w~z~2+G2w~w~-G2z~z~+2G1w~z~2+116a12φαG1w~w~w~+G1w~z~z~+G2w~w~z~+G2z~z~z~+φG2w~w~w~+G2w~z~z~-G1w~w~z~-G1z~z~z~≠0,in which (47)A0=32a122φ21-θ2+φ2,N0=2θφ1-θ1-2θ-2φ2-φ3-4θθ2-φ2,M0=1-3θ+2θ2-2φ2θ2-φ2+6θ-8θ2φ2.Taking the above calculations together, the occurrence of Neimark-Sacker bifurcation at the fixed point(u2,v2) in the discrete predator-prey system needs the satisfaction of conditions (34), (39), and (46). Moreover, if a0<0 and γ>0, an attracting invariant closed curve bifurcates from the fixed point for τ>τ0; if a0>0 and γ>0, a repelling invariant closed curve bifurcates from the fixed point for τ<τ0.


## 4. Numerical Results

To show the bifurcations and complex dynamics of the discrete predator-prey system, numerical simulations are carried out. We fix the parameters asa=3, b=2, c=0.2, d=1, α=1, and β=2 and assume that τ varies. Under such conditions, the fixed point is (2.8086,1.3436) and the critical point for flip bifurcation is τ∗=0.7501. At the critical bifurcation point, the two eigenvalues are λ1=-1 and λ2=0.8899. Considering the dynamics of map (23) restricted to the center manifold, we then obtain the values of two quantities as η1=-2.6665 and η2=-0.0036. Therefore, the predator-prey system indeed undergoes flip bifurcation at (u2,v2) in a small neighborhood of τ=τ∗, and the system dynamics converge to a period-2 orbit. Taking τ=0.8 as an example, the two stable periodic points are (3.1929,1.1228) and (2.2583,1.094).Figure1(a) displays the flip bifurcation diagram of the discrete predator-prey system withτvarying in [0.6,1.08]. It shows a period-doubling cascade in orbits of periods 2, 4, 8, and 16 (Figures 2(a)–2(d)). As determined by the maximum Lyapunov exponent in Figure 1(b), we find that the system firstly experiences chaos at τ=0.9698. Figures 1(c)–1(e) are local amplifications of the bifurcation diagram, displaying three periodic windows occurring in the chaotic region. In these periodic windows, period-6, period-5, and period-3 orbits appear (Figures 2(g), 2(j), and 2(m)). Moreover, a sub-period-doubling cascade emerges in each periodic window, leading to interior chaotic crisis and dynamic transition from periodic behaviors to chaos (Figures 2(m)–2(p)). Figure 2 shows various attractors in the flip bifurcation diagram with the increase of τ value.Figure 1
(a) Flip bifurcation diagram withτ covering [0.6,1.08]; (b) maximum Lyapunov exponent corresponding to (a); (c)–(e) local amplifications of (a) showing periodic windows. a=3, b=2, c=0.2, d=1, α=1, and β=2.
(a) (b) (c) (d) (e)Figure 2
Phase portraits for various values ofτ corresponding to Figure 1(a).
(a)
τ = 0.7 (b)
τ = 0.8 (c)
τ = 0.94 (d)
τ = 0.96 (e)
τ = 0.964 (f)
τ = 0.98 (g)
τ = 0.985 (h)
τ = 0.99 (i)
τ = 1.02 (j)
τ = 1.025 (k)
τ = 1.03 (l)
τ = 1.05 (m)
τ = 1.054 (n)
τ = 1.055 (o)
τ = 1.058 (p)
τ = 1.061Figure3 displays Neimark-Sacker bifurcation of the discrete predator-prey system with parametric conditions given as a=2, b=2, c=2, d=1.75, α=0.1, and β=0.1. In this case, the Neimark-Sacker bifurcation emerges at the fixed point (1.3467, 0.3849) in a small neighborhood of τ∗=1.1166. We can verify that, at the critical bifurcation point, λ1,2=0.2500±0.9682i~ and λ1,2=1. The values of γ and a0 are determined as γ=0.6716 and a0=-1.4078. Since γ>0 and a<0, we know that an attracting invariant closed curve bifurcates from the fixed point for τ>τ0.Figure 3
(a) Neimark-Sacker bifurcation diagram withτ covering [1.05,1.505]; (b) maximum Lyapunov exponent corresponding to (a); (c)–(h) local amplifications of (a). a=2, b=2, c=2, d=1.75, α=0.1, and β=0.1.
(a) (b) (c) (d) (e) (f) (g) (h)The bifurcation diagram in Figure3(a) and the maximum Lyapunov exponent diagram in Figure 3(b) suggest that the Neimark-Sacker bifurcation induces a route to chaos, leading to a dynamic transition from a fixed point, to invariant curves, with periodic windows occurring in between, to chaotic dynamics. The maximum Lyapunov exponent diagram says that the predator-prey system firstly shows chaotic dynamics at τ=1.4072. Local amplifications in Figures 3(c)–3(h) illustrate that the invariant curves may show many isolated cycles and periodic windows can appear in the zones of both invariant curves and chaos. In the invariant curve zone, the periodic windows often emerge intermittently, whereas in the chaotic zone, the periodic windows often exhibit sub-period-doubling cascade. In the periodic windows of Figures 3(c)–3(h), we find period-9, period-13, period-20, period-22 orbits. Figure 4 explicitly displays the alternation between periodic behaviors and invariant curves or chaotic dynamics with the increase of τ value.Figure 4
Phase portraits for various values ofτ corresponding to Figure 3(a).
(a)
τ = 1.12 (b)
τ = 1.225 (c)
τ = 1.27 (d)
τ = 1.28 (e)
τ = 1.29 (f)
τ = 1.35 (g)
τ = 1.410 (h)
τ = 1.425 (i)
τ = 1.431 (j)
τ = 1.44 (k)
τ = 1.4548 (l)
τ = 1.48Phase portraits for various values ofτ corresponding to Figure 3(a) are shown in Figure 4. It is worth noting that Figure 4(l) is a chaotic attractor in the sense of Marotto. According to the definition of snap-back repeller and Marotto chaos [32], we can prove that the fixed point (u2,v2) is a snap-back repeller by an iterative method and the emergence of Marotto chaos. In the case of a=2, b=2, c=2, d=1.75, α=0.1, β=0.1, and τ=1.48, we have λ1,2=0.006±1.2833i~ and λ1,2=1.2833>1. Notice that the fixed point of f with the norm of all eigenvalues of Df(u,v) exceeding 1 in magnitude is equivalent to 1+p+q>0, 1–p+q>0, and q>1. Let H1, H2, and H3 represent the sets of three inequalities in the front and Ur(u2,v2)=H1∩H2∩H3; then, we obtain the neighborhood Ur(u2,v2)={(u,v)∣1.04<u<1.53,0.21<v<0.55}. We need to find a point (u′,v′)∈Ur(u2,v2), satisfying u′,v′3=(u2,v2) and det(Df((u′,v′)k))≠0 for 1⩽k⩽3, where u′,v′k=fk(u′,v′). Through three iterative calculations and numerical simulations, we get the point (u′,v′)=(1.1004,0.4816) which meets f3(u′,v′)=(u2,v2) and det(Df(u′,v′k))=2.4284,1.3198,-2.5154≠0 for 1⩽k⩽3. Through these calculations, we proved that (u2,v2) is a snap-back repeller, and map (3) exhibits a chaotic behavior in the sense of Marotto.With the variation of other parameter values (e.g., parameterd), the predator-prey system may exhibit richer dynamical behaviors in the Neimark-Sacker bifurcation diagram. When we set the parameter values as a=2, b=2, c=2, d=1.85, α=0.1, and β=0.1, a new Neimark-Sacker bifurcation diagram is obtained as shown in Figure 5(a). The system undergoes Neimark-Sacker bifurcation at τ=0.9492, and it firstly enters chaotic dynamics at about τ=1.1964 (Figures 5(a) and 6(b)). Similar nonlinear characteristics to Figures 3 and 4 are found in this case, such as route to chaos, invariant curves, chaotic attractors, and periodic windows, which here exhibit period-5, period-9, period-24, and period-34 orbits. However, in the periodic window of [1.3134,1.3189], different dynamical behaviors are observed (Figure 5(c)). At about τ=1.3134, the chaotic dynamics of the predator-prey system suddenly disappear and the system dynamics jump to a periodic behavior with period 7 (Figures 6(a) and 6(b)). On each branch of the period-7 orbit, sub-Neimark-Sacker bifurcation occurs (Figure 5(d)). In the phase portrait, this reflects as each periodic point changing to a cycle (Figure 6(c)). The Neimark-Sacker bifurcation in each branch of period-7 orbit also triggers a route to chaos like Figure 3(a), leading to 7 small chaotic subsets. Chaotic interior crisis then takes place and leads to an abrupt transition from chaotic subsets to a narrow chaotic band (Figures 5(c) and 6(d)). The chaotic band is in the range [1.317,1.3171], and after that, the system dynamics turn to be periodic again (see period-15 orbit as shown in Figure 6(e)). Similarly, on each branch of the periodic orbit, sub-flip bifurcation occurs and induces a period-doubling cascade (Figures 5(e) and 6(e)–6(g)). Chaotic subsets and interior crisis appear sequently, taking the system entering chaotic dynamics again (Figure 6(h)).Figure 5
(a) Neimark-Sacker bifurcation diagram with parametric conditions given asa=2, b=2, c=2, d=1.85, α=0.1, and β=0.1; (b) maximum Lyapunov exponent corresponding to (a); (c) local amplification of (a) in the range of [1.313,1.32]; (d)-(e) amplifications of sub-Neimark-Sacker bifurcation and sub-flip bifurcation diagrams on each branch of periodic orbits.
(a) (b) (c) (d) (e)Figure 6
Phase portraits for various values ofτ corresponding to Figure 5(a).
(a)
τ = 1.25 (b)
τ = 1.314 (c)
τ = 1.316 (d)
τ = 1.3171 (e)
τ = 1.318 (f)
τ = 1.3182 (g)
τ = 1.3188 (h)
τ = 1.32Figure7 shows the Neimark-Sacker bifurcation diagram when the parameter values are given as a=2, b=2, c=2, d=1.95, α=0.1, and β=0.1. The critical Neimark-Sacker bifurcation point is τ=0.9492, and the first chaotic point is at around τ=1.1012 (Figures 7(a) and 7(b)). On the route to chaos, periodic windows with period-5, period-6, period-7, period-11, and period-16 orbits and narrow chaotic band are found. When we amplify the range [1.152,1.196], new dynamical behaviors are observed (Figure 7(c)). On each branch of the period-6 orbit, which occurs after the chaotic band (Figures 8(a) and 8(b)), the predator-prey system sequently undergoes sub-Neimark-Sacker bifurcation and periodic window and inverse sub-Neimark-Sacker bifurcation with the increase ofτ value (Figures 7(d) and 7(e)). This process results in a dynamic transition from period-6 orbit, to invariant curves with 6 cycles, to period-24 orbit, and then the system experiences the above dynamical behaviors again but in an inverse way (Figures 8(b)–8(f)), until the emergence of chaotic behaviors (Figure 8(g)).Figure 7
(a) Neimark-Sacker bifurcation diagram with parameter values given asa=2, b=2, c=2, d=1.95, α=0.1, and β=0.1; (b) maximum Lyapunov exponent corresponding to (a); (c) local amplification of (a) in the range of [1.152,1.196]; (d)-(e) amplifications of dynamics on each branch of period-6 orbit in (c).
(a) (b) (c) (d) (e)Figure 8
Phase portraits for various values ofτ corresponding to Figure 7(a).
(a)
τ = 1.15 (b)
τ = 1.155 (c)
τ = 1.16 (d)
τ = 1.168 (e)
τ = 1.17 (f)
τ = 1.185 (g)
τ = 1.21Figure9 displays a case in which the bifurcation and route to chaos are restricted in the range of 0<τ<1. The parametric conditions for this case are given as a=2, b=4, c=1, d=2.2, α=0.05, and β=0.5. The fixed point (0.5727,0.4496) is a stable attractor at τ=0.55. With the change of τ value, the fixed point loses its stability through Neimark-Sacker bifurcation at τ=0.5728. The maximum Lyapunov exponent corresponding to Figure 9(a) is shown (Figure 9(b)); we find that the predator-prey system enters chaos at τ=0.9314. Before the chaotic region, the predator-prey system undergoes a periodic window with period-10 orbits, and as the value ofτincreases, the period orbit doubles (Figures 10(a) and 10(b)), causing the system to jump to quasiperiodic behaviors (Figure 9(d)), following which the system returns to another periodic window with period halving (Figures 9(e) and 10(d)–10(g)). And then, the system experiences intermittent chaos, which suggests the abrupt appearance of chaotic dynamics. The chaotic attractor is demonstrated by the cloud of points in Figure 10(h).Figure 9
(a) Neimark-Sacker bifurcation diagram with parameter values given asa=2, b=4, c=1, d=2.2, α=0.05, and β=0.5; (b) maximum Lyapunov exponent corresponding to (a); (c)–(e) local amplifications of (a) in the range of [0.52,0.96].
(a) (b) (c) (d) (e)Figure 10
Phase portraits for various values ofτ corresponding to Figure 9(a).
(a)
τ = 0.86 (b)
τ = 0.88 (c)
τ = 0.91 (d)
τ = 0.9175 (e)
τ = 0.919 (f)
τ = 0.925 (g)
τ = 0.93 (h)
τ = 0.96Figure11 shows three-dimensional bifurcation diagrams in (τ,d,u) space and (τ,β,u) space with varying parameter values. As shown in Figure 11(a), we find that, with the increase of parameter d, the Neimark-Sacker bifurcation diagram moves toward the negativeτ axis. Such movement takes the decrease of Neimark-Sacker bifurcation critical point and first chaotic point and simultaneously changes the periodic windows on the route to chaos. Figure 11(b) exhibits that, with the increase of parameter β, the predator-prey system presents a dynamic transition from Neimark-Sacker bifurcation to flip bifurcation. Moreover, we find a critical case between the two bifurcations, where the critical points of Neimark-Sacker bifurcation and flip bifurcation are the same and the two bifurcations combine together. With the given parametric conditions, we have τ∗=τ0≈1.3851 in such case. The predator-prey system firstly experiences flip bifurcation and then abruptly jumps to the latter part of the Neimark-Sacker bifurcation diagram.Figure 11
Bifurcation diagrams in(τ,d,u) space and (τ,β,u) space, describing (a) the change of Neimark-Sacker bifurcation when a=2, b=2, c=2, α=0.1, β=0.1, and d=1.75, 1.85, and 1.95 and (b) the transition between Neimark-Sacker bifurcation and flip bifurcation when a=1, b=2, c=0.2, d=1, α=1, and β=0.8, 1.1, 1.3851, 1.7, and 2.
(a) (b)


## 5. Conclusions

In this research, the complex behaviors of discrete predator-prey systems with Crowley-Martin functional response are studied. Like many studies on the discrete predator-prey systems [5, 6, 20, 21, 32], we found that the discrete system can undergo flip bifurcation and Neimark-Sacker bifurcation at the stably coexistent fixed point. The two bifurcations can both trigger the route to chaos; that is, chaotic dynamics appear along with the emergence of bifurcations. On the routes to chaos, we also find a period-doubling cascade in orbits of periods 2, 4, and 8, periodic windows, invariant cycles, chaotic attractors, and Marotto chaos. Compared with the study of Ren et al. [32], who also investigated a discrete predator-prey system with Crowley-Martin functional response, many new interesting dynamic behaviors are revealed, including period-7, period-9, period-11, period-13, period-15, period-22, period-30, and period-34 orbits, sub-(inverse) Neimark-Sacker bifurcation, chaotic interior crisis, chaotic band, and intermittent periodic behaviors. Moreover, in the three-dimensional bifurcation diagram, we find the transition between flip bifurcation and Neimark-Sacker bifurcation and a critical bifurcation case between the two bifurcations. These results indicate that the discrete predator-prey system has more complex and richer dynamics than the corresponding continuous system.The generation of complex dynamics in the discrete predator-prey system is related to two important nonlinear mechanisms, flip bifurcation and Neimark-Sacker bifurcation, which cause the system to jump from steady state to oscillatory states, such as periodic and quasiperiodic states, and trigger routes to chaos. The occurrence of these behaviors induced by the two bifurcations can ecologically explain the complex dynamic states of populations. For example, the Australian entomologist Nicholson initiated a series of long-term experiments in which he kept populations of blowflies(Lucilia cuprina) in cages under a variety of resource renewal regimes. In one experiment, he found that the numbers of adult blowflies fluctuated as perturbed periodic state for the first 400 days [33]. Other experiments have investigated the long-term dynamics of freshwater organisms such as the cladoceranDaphnia in the laboratory. In one of the first such experiments, Slobodkin observed irregular fluctuations in time series that lasted between 200 and 400 days [34]. Application of the approaches for analyzing time series by computing Lyapunov exponents indicates that these time series may be quasiperiodic. Moreover, chaotic dynamics were also found in the natural world or laboratory experiments. Many authors have emphasized the presence of chaotic dynamics in a wide variety of communities with interacting species [34], such as outbreaks of planktons in the phytoplankton-zooplankton system, which is one of the important predator-prey systems. With the complex dynamics found in this research, many complicated states or oscillations of predator-prey populations in the real world may be possibly explained. Also, this research could be useful for researchers who will further explore the discrete-time predator-prey systems.



---

*Source: 2386954-2018-04-29.xml*