Codimension-Two Bifurcations of Fixed Points in a Class of Discrete Prey-Predator Systems

Abstract

The dynamic behaviour of a Lotka-Volterra system, described by a planar map, is analytically and numerically investigated. We derive analytical conditions for stability and bifurcation of the fixed points of the system and compute analytically the normal form coefficients for the codimension 1 bifurcation points (flip and Neimark-Sacker), and so establish sub- or supercriticality of these bifurcation points. Furthermore, by using numerical continuation methods, we compute bifurcation curves of fixed points and cycles with periods up to 16 under variation of one and two parameters, and compute all codimension 1 and codimension 2 bifurcations on the corresponding curves. For the bifurcation points, we compute the corresponding normal form coefficients. These quantities enable us to compute curves of codimension 1 bifurcations that branch off from the detected codimension 2 bifurcation points. These curves form stability boundaries of various types of cycles which emerge around codimension 1 and 2 bifurcation points. Numerical simulations confirm our results and reveal further complex dynamical behaviours.

1. Introduction

The dynamic relationship between predators and their prey is one of the dominant themes in both ecology and mathematical ecology due to its universal existence and importance; see [1–6]. The prototype Lotka-Volterra predator-prey system has received a lot of attention from theoretical and mathematical biologists; see [2–6]. This model is described by the following system of ordinary differential equations: 𝑥̇(𝑡)=𝑥𝑟₁−𝑎₁₁𝑥−𝑎₁₂𝑦̇(𝑡)𝑥𝑦,=𝑦−𝑟₂−𝑎₂₂𝑦+𝑎₂₁𝑥𝑦, (1.1) where 𝑥(𝑡) and 𝑦(𝑡) represent the densities of the prey and the predator, 𝑟₁, 𝑎₁₂, 𝑟₂, and 𝑎₂₁ are the intrinsic growth rate of the prey, the capture rate, the death rate of the predator, and the conversion rate, respectively, 𝑎₁₁ and 𝑎₂₂ denote the intraspecific competition coefficients of the prey and the predator, and 𝑟₁/𝑎₁₁ is the carrying capacity of the prey.

This and related models have been studied intensively in the previous decades and it has been noted, in particular, that they present a wealth of bifurcations of various codimensions. Organizing centers of codimension 3 have been investigated in great detail in [7] and references therein. For background on bifurcation theory, we refer to [8, 9].

Far-reaching generalizations of the model have been studied both analytically and numerically. We refer in particular to [10] where a five-parameter family of planar vector fields is studied that takes into account group defense strategy, competition between predators and a nonmonotonic response function. Though mathematical analysis is the prime tool also in this case, numerical methods are indispensable for a detailed study. In [10] not only the general-purpose languages Mathematica and Matlab were used, but also the dedicated bifurcation packages AUTO [11] and MatCont [12].

However, in recent years, many authors [4–6, 13] have suggested that discrete-time models are more appropriate than continuous ones, especially when the populations have nonoverlapping generations. Furthermore, discrete-time models often provide very effective approximations to continuous models which cannot be solved explicitly.

Though discrete models are by no means less complex than continuous ones, they have the computational (numerical) advantage that no differential equations have to be solved. In fact, the state-of-the-art in computational bifurcation study is in some respects more advanced for discrete systems than for continuous ones. So far, normal form coefficients for codimension 2 bifurcations of limit cycles are not computed by any available software. For discrete systems, such coefficients can be computed and they can even be used to start the computation of codimension 1 cycles rooted in such points; see [12, 14]. We will exploit this possibility.

In this paper we investigate a map in the plane for a pair of populations, first studied in [15, 16]. Danca et al. [15] investigated Neimark-Sacker bifurcations numerically. Murakami [16] discussed its branch points, flip bifurcations, and Neimark-Sacker bifurcations, establishing the sub- or supercritical character of the flip and Neimark-Sacker bifurcations by an explicit reduction to the center manifolds, obtaining a prediction for the invariant curve that branches off at the NS point. We undertake an analysis of the dynamics of the iterates of the map. This paper is organized as follows. In Section 2, we first discuss the map and the stability and bifurcations of its fixed points, proving the results in [16] by direct computations of the normal form coefficients. Next, in Section 3, we numerically compute curves of fixed points and bifurcation curves of the map and its iterates up to order 16 under variation of one and two parameters. We compute the critical normal form coefficients of all computed codim-1 and codim-2 bifurcations. These coefficients are powerful tools to compute stability boundaries of the map and its iterates and to switch to other bifurcation curves. In Section 4, we study more details on the bifurcation scenario of the system around a Neimark-Sacker bifurcation point. We conclude our work in Section 5 with a discussion of the obtained results.

2. The Model, the Fixed Points, Their Stability, and Bifurcations

In this paper, we study the map𝑥𝑦⟼()𝐹∶𝑎𝑥1−𝑥−𝑏𝑥𝑦𝑑𝑥𝑦, (2.1) which is analogous to (1.1) for the case of predators and prey with nonoverlapping generations; see [16, 17]. It can also be seen as an approximate discretization of the continuous-time Lotka-Volterra model𝑥̇(𝑡)=𝛼₀𝑥(𝑡(𝑡(𝑡)𝑦(𝑡),𝑦̇(𝑡)(𝑡)𝑦(𝑡)(𝑡),)(1−𝑥))−𝛼𝑚𝑥=𝑚𝑥−𝛽𝑦 (2.2) which is a simplification of (1.1) in which the intraspecific competition of the predator is ignored and the carrying capacity of the prey is 1. 𝑥 and 𝑦 represent the densities of the prey and predator, and 𝛼₀, 𝛼, 𝑚, and 𝛽 are nonnegative parameters. Applying the forward Euler scheme to system (2.2) with the stepsize 1/𝛽 and assuming 𝛼₀/𝛽≫1, we obtain the map (2.1) with nonnegative parameters 𝑎=𝛼₀/𝛽, b=𝛼𝑚/𝛽, and 𝑑=𝑚/𝛽.

We note that in (2.1) one of the two parameters 𝑏,𝑑 can be removed by a rescaling of 𝑦. So, the system is in reality a two-parameter system and fully equivalent to the system studied in [16]. We will generally choose 𝑑, 𝑎 as the unfolding parameters in the bifurcation study. In a way, it is the simplest possible discrete predator-prey model and, therefore, allows a reasonably complete analytical treatment as far as the fixed points of the map are concerned. However, we will see that even in this case the behaviour of cycles is very complicated and can only be studied by numerical methods.

We naturally start the bifurcation analysis of (2.1) with the calculation of the fixed points. These are the solutions (𝑥^∗,𝑦^∗) to 𝑎𝑥^∗1−𝑥^∗−𝑏𝑥^∗𝑦^∗=𝑥^∗,𝑑𝑥^∗𝑦^∗=𝑦^∗. (2.3) The origin 𝐸₁=(0,0) is a fixed point of (2.1) but is of little interest. Two further nontrivial fixed points are 𝐸₂=((𝑎−1)/𝑎,0) and 𝐸₃=(1/𝑑,1/𝑏(𝑎(1−(1/𝑑))−1)). We note that 𝐸₂ is biologically possible only if its coordinates are nonnegative, that is, 𝑎≥1. 𝐸₃ is biologically possible only if the following conditions hold: (i)𝑎>1,(ii)𝑑≥𝑎/(𝑎−1).

We start the local bifurcation analysis of the map (2.1) by linearization of 𝐹 around each of its fixed points. The Jacobian matrix 𝐽(𝑥,𝑦) is given by𝐽𝑥,𝑦=𝑎−2𝑎𝑥−𝑏𝑦−𝑏𝑥𝑑𝑦𝑑𝑥. (2.4)

The characteristic equation of 𝐽(𝑥,𝑦) is given by𝜆²−tr(𝐽)+det(𝐽)=0, (2.5) where tr(𝐽)=𝑎−2𝑎𝑥−𝑏𝑦+𝑑𝑥 and det(𝐽)=𝑑𝑥𝑎−2𝑑𝑥²𝑎.

2.1. Stability and Bifurcation of 𝐸₁

Proposition 2.1. The fixed point 𝐸₁ is asymptotically stable for 0≤𝑎<1. It loses stability via branching for 𝑎=1 and there bifurcates to 𝐸₂.

Proof. Eigenvalues of the Jacobian at 𝐸₁ are 𝑎 and 0. So, 𝐸₁ is stable if 𝑎<1 and loses stability at 𝑎=1. In (𝑎,𝑥)-space 𝐸₁ forms the curve (𝑎,0) with tangent vector (1,0). 𝐸₂ is represented in (𝑎,𝑥)-space by the curve (𝑎,(𝑎−1)/𝑎). When 𝑎=1, these curves intersect at (1,0) and the tangent vector in (𝑎,𝑥)-space is (1,1), so it is clear that 𝐸₁ branches to 𝐸₂ for 𝑎=1.

2.2. Stability and Bifurcation of 𝐸₂

The Jacobian matrix of (2.1) at 𝐸₂ is given by𝐽𝐸₂=()−𝑎+2−𝑏𝑎−1𝑎0𝑑(𝑎−1)𝑎. (2.6)

Proposition 2.2. The fixed point 𝐸₂ is linearly asymptotically stable if and only if 𝑎∈]1,3[ and 𝑑<𝑎/(𝑎−1). Moreover, it loses stability: (i)via branching for 𝑎=1 and there bifurcates to 𝐸₁, (ii)via branching for 𝑑=𝑎/(𝑎−1) and there bifurcates to 𝐸₃ if 1<𝑎<3, (iii)via a supercritical flip for 𝑎=3 if 𝑑<3/2.

Proof. The multipliers of 𝐽(𝐸₂) are 𝜆₁=−𝑎+2 and 𝜆₂=𝑑(𝑎−1)/𝑎. The fixed point 𝐸₂ is asymptotically stable if and only if |𝜆₁|<1 and |𝜆₂|<1, that is, if and only if 𝑎∈]1,3[ and 𝑑<𝑎/(𝑎−1). Boundary points of the stability region must satisfy one of three conditions: 𝑎=1, 𝑑=𝑎/(𝑎−1), or 𝑎=3.
In the first case, the conditions 𝑑<𝑎/(𝑎−1) and 𝑎<3 are satisfied for nearby values 𝑎>1, hence this is a real stability boundary. In Proposition 2.1, we proved that this is a branch point and the new branch consists of 𝐸₁ points.
In the second case, this is a stability boundary only if 1<𝑎<3. The Jacobian (2.6) then has an eigenvalue +1 and it is checked easily that these boundary points are also 𝐸₃ points.
In the third case, this is a stability boundary only if 𝑑<3/2. In this case, 𝜆₁=−1 which means that 𝐸₂ loses stability via a period doubling point. For supercriticality of the period doubling point, it is sufficient to show that the corresponding critical normal form coefficient 𝑏₁, 𝑏₁=16𝑝,𝐶𝑞,𝑞,𝑞+3𝐵𝑞,(𝐼−𝐴)⁻¹𝐵𝑞,𝑞, (2.7) derived by center manifold reduction is positive; see [9], Ch. 8 and [14]. Here, 𝐴=𝐽(𝐸₂), and 𝐵(⋅,⋅),𝐶(⋅,⋅,⋅) are the second- and third- order multilinear forms, respectively, and 𝑝 and 𝑞 are the left and right eigenvectors of 𝐴 for the eigenvalue −1, respectively. These vectors are normalized by ⟨𝑝,𝑞⟩=1, ⟨𝑞,𝑞⟩=1, where ⟨⋅⋅⟩ is the standard scalar product in ℝ². We obtain 𝑞𝑞=¹𝑞₂=10,𝑝𝑝=¹𝑝₂=12𝑏.3+2𝑑 (2.8) The components of the multilinear form 𝐵(𝑞,𝑞) are given by 𝐵𝑞,𝑞₁=²_𝑗,𝑘=1𝜕²𝑎𝑥(1−𝑥)−𝑏𝑥𝑦𝜕𝑥_𝑗𝜕𝑥_𝑘𝑞_𝑗𝑞_𝑘=−2𝑎𝑞₁𝑞₁𝐵=−6,𝑞,𝑞₂=²_𝑗,𝑘=1𝜕²𝑑𝑥𝑦𝜕𝑥_𝑗𝜕𝑥_𝑘𝑞_𝑗𝑞_𝑘=2𝑑𝑞₁𝑞₂=0, (2.9) where the state variable vector is for ease of notation generically denoted by (𝑥₁,𝑥₂)^𝑇 instead of (𝑥,𝑦)^𝑇.
Let 𝜁=(𝐼−𝐴)⁻¹𝐵(𝑞,𝑞), then we have 𝜁=₀⁻³ and find 𝐵𝑞,𝜁₁=−2𝑎𝑞₁𝜁₁=−2𝑎(−3)=18,𝐵𝑞,𝜁₂=𝑑𝑞₁𝜁₂=0. (2.10) The third-order multilinear form 𝐶(𝑞,𝑞,𝑞) is identically zero. The critical normal form coefficient 𝑏₁ is given by 𝑏₁=16𝑝^𝑇054=9, (2.11) which is clearly positive. This completes the proof of supercriticality of the flip point at 𝐸₂.

The stability region Ω^𝑆_𝐸2 of 𝐸₂, as obtained in Proposition 2.2, is shown in Figure 1.

Figure 1 

Stability regions in (𝑑,𝑎)-space. Ω𝑆𝐸1, Ω𝑆𝐸2 are the stability regions of 𝐸1 and 𝐸2, respectively. The stability region Ω𝑆𝐸3 of 𝐸3 is the union of Ω𝑆1𝐸3 and Ω𝑆2𝐸3 (which correspond to (i) and (iii) in Proposition 2.3, resp.), and the open interval that separates them (and corresponds to (ii) in Proposition 2.3).

2.3. Stability and Bifurcation of 𝐸₃

To study the stability of 𝐸₃, we use the Jury's criteria; see [6, Section  A2.1]. Let 𝐹(𝜆)=𝜆²−tr(𝐽(𝐸₃))𝜆+det(𝐽(𝐸₃)) be the characteristic polynomial of 𝐽(𝐸₃).

According to the Jury’s criteria, 𝐸₃ is asymptotically stable if the following conditions hold: 𝐹(−1)=1+tr𝐽𝐸₃+det𝐽𝐸₃𝐹(1)𝐽𝐸>0,=1−tr₃+det𝐽𝐸₃𝐽𝐸>0,1−det₃>0. (2.12)

At 𝐸₃, we have:𝐽𝐸₃=𝑑−𝑎𝑑−𝑏𝑑𝑑𝑎−𝑎−𝑑𝑏1. (2.13) We note that tr(𝐽(𝐸₃))=(2𝑑−𝑎)/𝑑 and det(𝐽(𝐸₃))=𝑎(𝑑−2)/𝑑 are independent of 𝑏.

Proposition 2.3. 𝐸₃ is linearly asymptotically stable if and only if one of the following mutually exclusive conditions holds: (i)3/2<𝑑<9/4 and 𝑑/(𝑑−1)<𝑎<3𝑑/(3−𝑑), (ii)𝑑=9/4 and 1.8=𝑑/(𝑑−1)<𝑎<𝑑/(𝑑−2)=3𝑑/(3−𝑑)=9, (iii)𝑑>9/4 and 𝑑/(𝑑−1)<𝑎<𝑑/(𝑑−2).

Proof. The criterion 𝐹(1)>0 is easily seen to be equivalent to the condition 𝑑𝑎>𝑑−1,𝑑>1, (2.14) or equivalently, 𝑎𝑑>𝑎−1,𝑎>1. (2.15)
Next, the criterion det(𝐽(𝐹3))<1 is easily seen to be equivalent to 𝑑𝑎<𝑑−2if𝑑>2. (2.16)
The criterion 𝐹(−1)>0 translates as 𝑎<3𝑑3−𝑑if𝑑<3. (2.17)

In Figure 1, the parts of the 3 curves 𝑎=𝑑/(𝑑−1), 𝑎=𝑑/(𝑑−2), and 𝑎=3𝑑/(3−𝑑) where 𝑎 and 𝑑 are positive are depicted (resp., for 𝑑>1, 𝑑>2, and 0<𝑑<3). The curves 𝑎=𝑑/(𝑑−1) and 𝑎=3𝑑/(3−𝑑) intersect solely at (𝑑=3/2,𝑎=3) and the curves 𝑎=𝑑/(𝑑−2) and 𝑎=3𝑑/(3−𝑑) at (𝑑=9/4,𝑎=9). It follows from the figure that the 3 inequalities are fullfilled when (i), (ii), or (iii) holds. This can also easily be shown algebraic. In Figure 1, the stability region Ω^𝑆_𝐸3 of 𝐸₃ is thus the union of Ω^𝑆1_𝐸3 and Ω^𝑆2_𝐸3 (which correspond to (i) and (iii) in this proposition, resp.), and the open interval that separates them (and corresponds to (ii)).

Proposition 2.4. 𝐸₃ loses stability: (i)via a flip point when 3/2<𝑑<9/4 and 𝑎=3𝑑/(3−𝑑), (ii)via a Neimark-Sacker point when 9/4<𝑑 and 𝑎=𝑑/(𝑑−2), (iii)via a branch point when 𝑑>3/2 and 𝑎=𝑑/(𝑑−1) where it bifurcates to 𝐸₂, (iv)via a branch-flip (BPPD) point when 𝑑=3/2 and 𝑎=3, (v)via a resonance 1 : 2 point when 𝑑=9/4 and 𝑎=9.

Proof. By Proposition 2.3 the stability boundary of 𝐸₃ consists of parts of three curves, namely, (1)Curve 1: 𝑎=3𝑑/(3−𝑑). (2)Curve 2: 𝑎=𝑑/(𝑑−2). (3)Curve 3: 𝑎=𝑑/(𝑑−1). The points of Curve 1 which are on the stability boundary of 𝐸₃ satisfy 𝐹(−1)=0, that is, they have an eigenvalue −1 and, thus, are period doubling points. The points of Curve 2 which are on the stability boundary satisfy det(𝐽(𝐸₃))=1, that is, they have two eigenvalues with product 1 and, thus, are Neimark-Sacker points. The points of Curve 3 which are on the stability boundary satisfy 𝐹(1)=0, that is, they have an eigenvalue 1. It can be checked easily that 𝐸₃ then branches to 𝐸₂.
Combining this with Proposition 2.2 we find that the interior points of the boundary parts of Curves 1, 2, and 3 form the sets described in parts (i), (ii), and (iii) of Proposition 2.3, respectively.
At the intersection of Curves 1 and 3 (when 𝑑=3/2 and 𝑎=3), the point is a branch-flip point (BPPD) with eigenvalues −1 and 1. The intersection point of Curves 1 and 2 (when 𝑑=9/4 and 𝑎=9) is a resonance 1 : 2 point with double eigenvalue −1.

Proposition 2.5. The flip point in Proposition 2.4, part (i), is subcritical.

Proof. For the subcriticality of the flip bifurcation, it is sufficient to show that the normal form coefficient 𝑏₁ in (2.7) is always negative. By using the same procedure as we used in the proof of Proposition 2.2, we obtain 𝑞𝑞=¹𝑞₂=1𝑑(2𝑑−3)𝑏()𝑝𝑑−3,𝑝=¹𝑝₂=1𝑏2𝑑. (2.18) To simplify computations, we do not normalize 𝑝 and 𝑞, since rescaling does not change the sign of 𝑏₁ provided that ⟨𝑝,𝑞⟩ is positive (it can be proved easily that this is the case). 𝐵(𝑞,𝑞) is computed as 𝐵𝑞,𝑞=−4𝑑2𝑑²(2𝑑−3)𝑏(𝑑−3). (2.19)
Let 𝜁=(𝐼−𝐴)⁻¹𝐵(𝑞,𝑞), then we have 𝜁=^𝑑_−𝑑²_{(4+(3/(3−𝑑)))/𝑏} and find 𝐵𝑞,𝜁=2𝑑²()3−𝑑6−𝑑−2𝑎𝑑−2𝑑³(6−𝑑)𝑏(3−𝑑). (2.20) The third-order multilinear form 𝐶(𝑞,𝑞,𝑞) is zero. The critical normal form coefficient 𝑏₁=(−1/2)𝑑³/(3−𝑑), which is clearly negative for 3/2<𝑑<9/4. This completes the proof of subcriticality of the flip point.

Proposition 2.6. The Neimark-Sacker point in Proposition 2.4, part (ii), is supercritical.

Proof. To prove the supercriticality of the Neimark-Sacker point, it is sufficient to show that the real part of the corresponding critical normal form coefficient 𝑑₁, 𝑑₁=12𝑒^−𝑖𝜃0𝑝,𝐶𝑞,𝑞,𝑞+2𝐵𝑞,ℎ₁₁₀₀+𝐵𝑞,ℎ₂₀₀₀,ℎ₁₁₀₀=𝐼_𝑛−𝐴⁻¹𝐵𝑞,𝑞,ℎ₂₀₀₀=𝑒^2𝑖𝜃0𝐼_𝑛−𝐴⁻¹𝐵𝑞,𝑞, (2.21) is negative; see [14]. Here, 𝜃₀(0<𝜃₀<𝜋) is the argument of the critical multiplier, 𝐴=𝐽(𝐸₃), and 𝐵(⋅,⋅), 𝐶(⋅,⋅,⋅) are the second- and third-order multilinear forms, respectively. 𝑝 and 𝑞 are the left and right eigenvectors of 𝐴𝐴𝑞=𝑒^𝑖𝜃0𝐴𝑞,_𝑇𝑝=𝑒^−𝑖𝜃0𝑝. (2.22) These vectors are normalized by ⟨𝑝,𝑞⟩=1,⟨𝑞,𝑞⟩=1, where ⟨𝑝,𝑞⟩=𝑝^𝐻𝑞 is the Hermitian inner product in ℂ². In the Neimark-Sacker point, 𝑎=𝑑/(𝑑−2) (𝑑>9/4, so 𝑑−2>0). We get 𝐴=𝑑−3−𝑏𝑑−2𝑑𝑑𝑏()1𝑑−2, (2.23) with characteristic polynomial 𝜆²−((2𝑑−5)/(𝑑−2))𝜆+1. It follows that ℜ(𝜆)=(𝑑−(5/2))/(𝑑−2) and ℑ(𝜆)=𝑑−(9/4)/(𝑑−2), so 𝑒^𝑖𝜃0=1/(𝑑−2)(𝑑−(5/2)+𝑖𝑑−(9/4)). The eigenvectors are ()𝑏𝑞=𝑑−2𝑑−12−𝑖9𝑑−4,𝑑()𝑝=2𝑑−9/2−𝑖𝑑−(9/4)𝑏4𝑑²1−17𝑑+18𝑏𝑑12−𝑖9𝑑−4, (2.24) where 𝑞 is not normalized since rescaling of 𝑞 does not change the sign of 𝑑₁ (this can be proved easily) and we can simplify the computations by not normalizing. The calculation of the components of 𝐵(𝑞,𝑞) gives us 𝐵𝑞,𝑞()𝑏=−𝑏𝑑−2𝑑1, (2.25) and it follows that ℎ₁₁₀₀()𝑏()=𝑏𝑑−2𝑑−2𝑑−2. (2.26) The remaining calculations are done with the help of Maple. The calculation of ℎ₂₀₀₀ gives us ℎ₂₀₀₀=−2(𝑑−2)³𝑏−50𝑑²+135𝑑+6𝑖√4𝑑−9𝑑²−29𝑖√4𝑑−9𝑑−119+35𝑖√4𝑑−9+6𝑑³⋅𝑏𝑑17𝑑−7𝑖√4𝑑−9𝑑−23+11𝑖√4𝑑−9−3𝑑²+𝑖√4𝑑−9𝑑²𝑖√4𝑑−9𝑑+22−2𝑖√4𝑑−9+4𝑑²,−19𝑑 (2.27) where 𝑒^2𝑖𝜃0=1(𝑑−2)²𝑑²−6𝑑+172+𝑖(2𝑑−5)9𝑑−4. (2.28) For 𝐵(𝑞,ℎ₂₀₀₀) and 𝐵(𝑞,ℎ₁₁₀₀), we obtain 𝐵𝑞,ℎ₂₀₀₀=2(𝑑−2)³𝑏²−50𝑑²+135𝑑+6𝑖√4𝑑−9𝑑²−29𝑖√4𝑑−9𝑑−119+35𝑖√4𝑑−9+6𝑑³⋅𝑏𝑑32𝑑−6𝑖√4𝑑−9𝑑−29+9𝑖√4𝑑−9−13𝑑²+𝑖√4𝑑−9𝑑²+2𝑑³𝑖√4𝑑−9𝑑²+2𝑑−8𝑖√4𝑑−9𝑑−2𝑑³+7𝑑²−17+13𝑖√,𝐵4𝑑−9𝑞,ℎ₁₁₀₀=(𝑑−2)²𝑏²2𝑏1+𝑖√√4𝑑−9−5−𝑖.4𝑑−9 (2.29) The third-order multilinear form 𝐶(𝑞,𝑞,𝑞) is identically zero, so we get 𝑑₁=12𝑒^−𝑖𝜃0𝑝^𝐻2𝐵𝑞,ℎ₁₁₀₀+𝐵𝑞,ℎ₂₀₀₀=14−2𝑑+5+𝑖√4𝑑−94𝑑−9+𝑖√4𝑑−9𝑏²𝑑×31𝑑²+𝑖√4𝑑−9𝑑³+44−5𝑑³−64𝑑−12𝑖√4𝑑−9−7𝑖√4𝑑−9𝑑²+16𝑖√×()4𝑑−9𝑑4𝑑−9⁻¹−50𝑑²+135𝑑+6𝑖√4𝑑−9𝑑²−29𝑖√4𝑑−9𝑑−119+35𝑖√4𝑑−9+6𝑑³⁻¹. (2.30) The real part of 𝑑₁ equals −1/2𝑏²𝑑 (exact up to a positive factor due to the nonnormalization of 𝑞), which is clearly negative. This completes the proof of supercriticality of the Neimark-Sacker point in Proposition 2.4, part (ii).

The stability regions (i) and (iii) of 𝐸₃ obtained in Proposition 2.3 are depicted as Ω^𝑆1_𝐸3 and Ω^𝑆2_𝐸3, respectively, in Figure 1. The region (ii) is the open interval on the common boundary of Ω^𝑆1_𝐸3 and Ω^𝑆2_𝐸3. We have a complete description of the stability region of 𝐸₃ for all parameter combinations.

The biological interpretation of Figure 1 is as follows. The situation with no prey and no predators exists for all parameter values but is stable only in Ω^𝑆_𝐸1, that is, for 𝑎<1. The situation with a fixed number of prey but no predators exists for all 𝑎>1 but is stable only in Ω^𝑆_𝐸2. Coexistence of fixed nonzero numbers of prey and predators is possible whenever 𝑎>1 and 𝑑>𝑎/(𝑎−1) but is stable only in the union of Ω^𝑆1_𝐸3 and Ω^𝑆2_𝐸3.

3. Numerical Bifurcation Analysis of 𝐸₂ and 𝐸₃

In this section, we perform a numerical bifurcation analysis by using the MATLAB package Cl_MatContM; see [12, 14]. The bifurcation analysis is based on continuation methods, whereby we trace solution manifolds of fixed points while some of the parameters of the map vary; see [18]. We note that in a two-parameter setting the boundaries of stability domains of cycles are usually curves of codimension 1 bifurcations, which necessarily have to be computed by numerical continuation methods.

To validate the model, we compute a number of stable cycles for parameter values in the regions where we claim they exist; these stable cycles were found by orbit convergence, independently of the continuation methods. They are represented in the Figures 2, 6, 11, and 14. Similarly, a stable closed invariant curve was computed by orbit convergence and is represented in Figure 4. The computation of an Arnold tongue in Section 3.5 by a continuation method also validates our computations.

Figure 2 

A stable 8-cycle of (2.2) for 𝑎=3.571920967580968, 𝑏=0.2, 𝑑=1.4.

3.1. Numerical Bifurcation of 𝐸₂

By continuation of 𝐸₂=(0.5;0) starting from 𝑎=2, 𝑏=0.2, 𝑑=1.4, in the stable region of 𝐸₂ with 𝑎 free, we see that 𝐸₂ is stable when 1<𝑎<3. It loses stability via a supercritical period doubling point (PD, the corresponding normal form coefficient is 9>0) when 𝑎=3, and via a branch point (BP) when 𝑎 crosses 1.4. The output of Run 1 is given by:label = BP, x = (0.000000 0.000000 1.000000) label = PD, x = (0.666667 0.000000 3.000000) normal form coefficient of PD = 9.000000e+000

The first two entries of x are the coordinate values of the fixed point 𝐸₂, and the last entry of x is the value of the free parameter 𝑎 at the corresponding bifurcation point. We note that the normal form coefficient of the PD point is 9, confirming (2.11). We note furthermore that the detected bifurcation points in this run are in accordance with the statement of Proposition 2.2.

Beyond the PD point the dynamics of (2.1) is a stable 2-cycle. Cl_MatContM allows to switch to the continuation of this 2-cycle. It loses stability at a supercritical PD point: label = PD, x = (0.849938 0.000000 3.449490) normal form coefficient of PD = 4.062566e+002

A stable 4-cycle is born when 𝑎>3.449490. It loses stability at a supercritical PD point: label = PD, x = (0.884050 0.000000 3.544090) normal form coefficient of PD = 1.617268e+004

Thus, when 𝑎>3.544090 a stable 8-cycle emerges. A stable 8-cycle is given by 𝐶₈={𝑋⁸₁,…,𝑋⁸₈}, where 𝑋⁸₁=(0.880450408535962,0) where 𝑎=3.571920967580968, 𝑏=0.2, and 𝑑=1.4. This 8-cycle is represented in Figure 2. In this situation there are no predators and the number of prey repeats itself every 8 time spans. Switchings at PD points of the second and fourth iterates are given in Figure 3.

Figure 3 

Curves of fixed points of the first, second and fourth iterates.

Figure 4 

A stable closed invariant curve for 𝑎=3.5, 𝑏=0.2, and 𝑑=2.83.

We note that for 𝐸₂ the map (2.1) is a logistic map. In fact, in this case the predator becomes extinct and the prey undergoes the period-doubling bifurcation to chaos when further increasing the parameter 𝑎; see [19].

Continuation of 𝐸₂ starting from the same parameter values as in Run 1, with 𝑑 as free parameter, leads to: label = BP, x = (0.500000 0.000000 2.000000)

The appearance of a branch point is consistent with Proposition 2.2 part (ii) which states that 𝐸₂ bifurcates to 𝐸₃ when 𝑏=𝑑=𝑎/(𝑎−1)=2.

3.2. Numerical Bifurcation of 𝐸₃

We now consider 𝐸₃=(0.454545;4.545454) which is in the stable region for the parameter values 𝑎=3.5, 𝑏=0.2, and 𝑑=2.2 (stability follows from Proposition 2.3 part (i)). We do a numerical continuation of 𝐸₃ with control parameter 𝑑, and call this Run 2. The output of Run 2 is given by: label = PD, x = (0.619048 1.666667 1.615385) normal form coefficient of PD = −5.905026e−001 label = NS, x = (0.357143 6.250000 2.800000) normal form coefficient of NS = −6.971545e−002

𝐸₃ is stable when 1.615385<𝑑<2.8. It loses stability via a supercritical Neimark-Sacker (NS) point when 𝑑=2.8, which is consistent with Proposition 2.4 part (ii) (𝑑/(𝑑−2)=3.5=𝑎). It also loses stability through a subcritical PD point when 𝑑=1.615385, which is consistent with Propositions 2.5 and 2.4 part (i) since 3𝑑/(3−𝑑)=3.5=𝑎.

The dynamics of the system prior to the PD point consists of an unstable 2-cycle that coexists with a stable fixed point. Beyond the NS point the dynamics of the system consists of a stable closed invariant which coexists with unstable fixed points of (2.1). For 𝑑=2.83, a stable closed invariant curve is created around the unstable fixed point 𝐸₃ (see Figure 4). Now, we compute the period doubling curve, with 𝑎 and 𝑑 free, by starting from the PD point detected in Run 2. We call this Run 3. label = LPPD, x = (0.666667 0.00000 3 1.500000) Normal form coefficient for LPPD:[a/e, be] = −3.333334e−001, 4.541342e−006 label = R2, x = (0.444444 20.000000 9.000000 2.250000) Normal form coefficient for R2:[c, d] = 1.197634e−001, −7.185806e−001

Two codim-2 bifurcation points are detected on the flip curve, namely, a fold-flip LPPD and a resonance 2 bifurcation R2. We note that the LPPD point is in reality a branch-flip (BPPD) point, which by the software is detected as an LPPD point since BPPD points are ungeneric on a curve of PD points. This curve is shown in Figure 5 (left curve). Now, we compute the NS curve, with 𝑎 and 𝑑 free parameters, by starting from the NS point of Run 2. We call this Run 4. label = R4, x = (0.400000 10.000000 5.000000 2.500000 −0.000000) Normal form coefficient of R4: A = −1.240347e+000 + 6.201737e−001 i label = R3, x = (0.428571 15.000000 7.000000 2.333333 −0.500000) Normal form coefficient of R3: Re(c_1) = −5.000000e−001 label = R2, x = (0.444444 20.000000 9.000000 2.250000 −1.000000) Normal form coefficient of R2: [c, d] = 1.197646e−001, −7.185574e−001

Figure 5 

Flip and Neimark-Sacker bifurcation curves starting from points in Run 2. We note that the LPPD point is in reality a BPPD point.

Figure 6 

A stable 4-cycle near an R4 point for 𝑎=4.99, 𝑏=0.2, and 𝑑=2.56.

The computed curve of NS points is also shown in Figure 5 (right curve). The codim-2 bifurcations that are computed along the Neimark-Sacker curve are a resonance 1 : 2 (R2), resonance 1 : 3 (R3), and a resonance 1 : 4 (R4) point. In addition to the coordinates of the bifurcation point, parameter values and the real part of the Neimark-Sacker multiplier at the bifurcation point are output. We note that the PD curve crosses over the NS curve 𝑎=𝑑/(𝑑−2) at a resonance 1 : 2 (R2) point, which is consistent with Proposition 2.4 part (v). It also hits the BP curve 𝑎=𝑑/(𝑑−1) at a branch-flip point BPPD which is consistent with Proposition 2.4 part (iv).

3.3. Orbits of Period 4, 8, 16 and 32

The normal form coefficient 𝐴 of the R4 point in Run 4 satisfies |𝐴|>1, hence two cycles of period 4 of the map are born. A stable 4-cycle for 𝑎=4.99, 𝑏=0.2, and 𝑑=2.56 is given by: 𝐶₄={𝑋₁,𝑋₂,𝑋₃,𝑋₄}, where 𝑋₁=(0.506291269909196,9.483616960226117). We present this cycle in Figure 6. In order to compute the stability region of this 4-cycle, we compute the twofold curves of the fourth iterate rooted at the R4 point. These curves exist since |𝐴|>1, see [9], and switching from an R4 point to the fold curves of the fourth iterate is supported by Cl_MatContM. The stable fixed points of the fourth iterate exist in the wedge between the twofold curves. The output of this continuation, Run 5, is given below. The curves are shown in Figure 7. label = LPPD, x = (0.077546 7.084151 3.948550 3.637018) Normal form coefficient for LPPD:[a/e, be] = −1.638977e−001, −1.231773e+001 label = LPPD, x = (0.540029 8.419070 5.511767 2.615029) Normal form coefficient for LPPD:[a/e, be] = 7.735689e−001, −8.261238e−002

Figure 7 

Twofold curves (LP4) emanate from an R4 point on an NS curve.

We can further compute the stability boundaries of the 4-cycle. This region is bounded by the two just computed limit point curves and a period doubling curve of the fourth iterate rooted at the detected LPPD points on the branches of LP⁴ curves. Continuation of the flip curve of the fourth iterate emanated at the LPPD of Run 5 is given below. We call this Run 6.

label = LPPD, x = (0.862908 2.257241 3.948550 3.637018) Normal form coefficient for LPPD:[a/e, be] = −1.638977e−001, −9.218471e+001 label = LPPD, x = (0.540029 8.419070 5.511767 2.615029) Normal form coefficient for LPPD:[a/e, be] = 7.735689e−001, −8.261238e−002 label = GPD, x = (0.541277 8.238085 5.478361 2.609243) Normal form coefficient of GPD = 6.483721e−001 label = R2, x = (0.582936 7.330070 5.625162 2.664169) Normal form coefficient for R2:[c, d] = 6.130829e−002, −3.509615e+000

By superposing the flip curve on Figure 7, we obtain Figure 8. We further compute a curve of fixed points of the fourth iterate starting from the 4-cycle 𝐶₄ presented in Figure 6, with control parameter 𝑑. The output of this continuation, Run 7, is given below. The curve is presented in Figure 9. label = LP, x = (0.418135 9.911757 2.502419) normal form coefficient of LP = −2.095861e−001 label = PD, x = (0.554614 8.335785 2.864961) normal form coefficient of PD = 1.909178e−001 label = NS, x = (0.546431 9.030507 2.608544) normal form coefficient of NS = −2.177981e+000

Figure 8 

Twofold curves (LP4) of the fourth iterate emanate from an R4 point, and a flip curve of the fourth iterate (PD4) is rooted at an LPPD point.

Figure 9 

Curve of fixed points of the fourth iterate starting from the 4-cycle 𝐶4.

The 4-cycle remains stable when 2.502419<𝑑<2.608544. Now, we compute an NS-curve starting from the computed NS point in Run 8 and call this Run 9. This curve is superposed on Figure 8 and is depicted in Figure 10. label = R4, x = (0.594620 7.423051 4.615965 2.666797 0.000000)Normal form coefficient of R4: A = −4.453277e+000 + −1.491255e+000 i label = R3, x = (0.529139 9.928228 5.218594 2.614545 −0.500000)Normal form coefficient of R3: Re(c_1) = −6.065820e+000 label = R2, x = (0.513005 11.383897 5.625162 2.664169 −1.000000)Normal form coefficient of R2: [c, d] = 1.004097e−001, −5.657220e+000

Figure 10 

Stability region Ω4𝑆 of the 4-cycle 𝐶4. The boundaries consist of NS, LP4, PD4, and NS4 curves.

Figure 11 

A stable 16-cycle around the R4 point of Run 9 for 𝑎=4.68, 𝑏=0.2, and 𝑑=2.66.

The stability region Ω⁴_𝑆 of the 4-cycle 𝐶₄ is bounded by the LP⁴, PD⁴, and NS⁴ curves (see Figure 10).

Now, we consider the R4 point computed in Run 9. Since the corresponding normal form coefficient 𝐴 satisfies |𝐴|>1, two cycles of period 16 of the map are born. A stable 16-cycle for 𝑎=4.68, 𝑏=0.2, and 𝑑=2.66 is given by 𝐶₁₆={𝑋₁,𝑋₂,…,𝑋₁₆}, where 𝑋₁=(0.596607820551550,7.625016168527653). We present this cycle in Figure 11.

In order to compute the stability region of this 16-cycle, we compute twofold curves of the sixteenth iterate rooted at the R4 point. These curves exist since |𝐴|>1. The stable fixed points of the sixteenth iterate exist in the wedge between the twofold curves. The output of this continuation, Run 10, is given below, and the fold curves are shown in Figure 12. label = LPPD, x = (0.557836 8.068783 4.651131 2.670940) Normal form coefficient for LPPD:[a/e, be] = 1.235611e+001, 3.799853e−001 First Lyapunov coefficient for second iterate = 3.799853e−001 label = LPPD, x = (0.519225 8.927562 4.765661 2.658381) Normal form coefficient for LPPD:[a/e, be] = 9.129886e+000, 5.166118e+000 First Lyapunov coefficient for second iterate = 5.166118e+000 label = LPPD, x = (0.469475 10.678780 5.153847 2.636682) Normal form coefficient for LPPD:[a/e, be] = 1.298255e−001, −4.388341e+000 label = LPPD, x = (0.580415 9.268192 5.370935 2.641181) Normal form coefficient for LPPD:[a/e, be] = 1.528369e+000, 6.631508e+001 First Lyapunov coefficient for second iterate = 6.631508e+001 label = LPPD, x = (0.578148 9.853334 5.628124 2.676611) Normal form coefficient for LPPD:[a/e, be] = −9.656910e−001, −2.522009e+004 label = R1, x = (0.577592 9.941486 5.660840 2.682153) normal form coefficient of R1 = −1

Figure 12 

Twofold curves (LP16) emanate from an R4 point on the NS4 curve.

We further compute a curve of fixed points of the sixteenth iterate starting from the 16-cycle 𝐶₁₆ presented in Figure 11, with control parameter 𝑎. The output of this continuation, Run 11, is given below. The curve is presented in Figure 13. label = LP, x = (0.598114 7.473352 4.646130)normal form coefficient of LP = 2.965904e+000 label = PD, x = (0.595604 7.701225 4.700640) normal form coefficient of PD = 2.829176e+002

Figure 13 

Curve of fixed points of the sixteenth iterate starting from the 16-cycle 𝐶16.

Figure 14 

A stable 32-cycle for 𝑎=4.7071, 𝑏=0.2, and 𝑑=2.66.

The 16-cycle remains stable when 4.646130<𝑑<4.700640. The dynamics of the system beyond the supercritical PD is a stable 32-cycle which coexists with the unstable fixed points of the sixteenth iterate of the system.

A stable 32-cycle for 𝑎=4.7071, 𝑏=0.2, and 𝑑=2.66 is given by: 𝐶₃₂={𝑋₁,𝑋₂,…,𝑋₃₂}, where 𝑋₁=(0.619778764655944,7.243979650146042). We present this cycle in Figure 14.

Now, we compute a PD-curve starting from the computed PD point in Run 10, and call this Run 12. This curve is superposed on Figure 12, and is depicted in Figure 15. The boundaries of the stability region Ω_𝑆¹⁶ of the 16-cycle consist of the LP¹⁶ and PD¹⁶ curves, as illustrated in Figure 15.label = LPPD, x = (0.600292 7.504485 4.651131 2.670940) Normal form coefficient for LPPD:[a/e, be] = 1.235612e+001, 1.579190e+003 First Lyapunov coefficient for second iterate = 1.579190e+003 label = LPPD, x = (0.587235 7.998201 4.765661 2.658381) Normal form coefficient for LPPD:[a/e, be] = 9.129867e+000, 9.328278e+002 First Lyapunov coefficient for second iterate = 9.328278e+002 label = GPD, x = (0.529363 9.922888 5.130266 2.635413) Normal form coefficient of GPD = 1.292549e+007 label = LPPD, x = (0.526233 10.046755 5.153847 2.636682) Normal form coefficient for LPPD:[a/e, be] = 1.298254e−001, −7.143714e+002 label = LPPD, x = (0.521191 10.840009 5.370935 2.641181) Normal form coefficient for LPPD:[a/e, be] = 1.528369e+000, 1.217953e+002 First Lyapunov coefficient for second iterate = 1.217953e+002 label = GPD, x = (0.521992 10.863772 5.384512 2.642770) Normal form coefficient of GPD = −6.188850e+003

Figure 15 

Stability region Ω𝑆16 of the 16-cycle 𝐶16. The boundaries consist of LP16 and PD16 curves.

3.4. Orbits of Period 3

Next, we consider the resonance 1 : 3 (R3) point in Run 4. Since its normal form coefficient is negative, the bifurcation picture near the R3 point is qualitatively the same as presented in [9, Figure 9.12]. In particular, there is a region near the R3 point where a stable invariant closed curve coexists with an unstable fixed point. For parameter values close to the R3 point, the map has a saddle cycle of period three.

Furthermore, a Neutral Saddle bifurcation curve of fixed points of the third iterate emanates [9, Chapter  9]. We compute this curve by branch switching at the R3 point. This curve is presented in Figure 16.

Figure 16 

Curve of Neutral Saddle bifurcations of the third iterate for 𝑏=0.2.

3.5. Computation of Arnold Tongues

It is well known that near a Neimark-Sacker curve there exists a dense array of resonance tongues, generalizing the isolated tongue of period 4 in Figure 7. The tongues locally form an open and dense set of the parameter plane. There are also quasiperiodic invariant circles in between that correspond to a set of positive measure in the parameter plane.

So far, no numerical methods have been implemented to specifically compute the boundaries of the resonance tongues that are rooted in weakly resonant Neimark-Sacker points (unlike the strong resonant 1 : 4 case). However, since they are limit point curves of fixed points of cycles with known periods, they can be computed relatively easily if the cycles inside the tongue are globally stable (which depends on the criticality of the Neimark-Sacker curve and the noncritical multipliers as well). It is sufficient to find a fixed point of cycles inside the tongue by orbit convergence and to continue it in one free parameter to find a point on the boundary of the Arnold tongue as a limit point of cycles. From this, the boundary curves can be computed by a continuation in two free parameters.

In Figure 17, we present an Arnold Tongue rooted in a weak 2 : 7 resonant Neimark-Sacker point. Its computation started from a stable 7-cycle with 𝑥=0.401985064603439, 𝑦=12.930999111981340, 𝑎=5.844728289174310, 𝑏=0.2, and 𝑑=2.45. We note that the boundary curves contain further bifurcation points.

Figure 17 

An Arnold tongue rooted in a weak 2 : 7 resonant Neimark-Sacker point.

4. Numerical Simulation

To reveal the qualitative dynamical behaviours of (2.2) near the NS point in Run 1, we present a complete bifurcation sequence that is observed for different values of 𝑑. We fix the parameters 𝑎=3.5, 𝑏=0.2 and assume that 𝑑 is free.

Figure 18 shows that 𝐸₃ is a stable attractor for 𝑑=2.67. We note that for the given parameters value, 𝐸₃ is a stable fixed point consistent with Proposition 2.3 part (iii). The behaviour of (2.2) at 𝑑=2.78, so before the NS point at 𝑑=2.8, is depicted in Figure 19.

Figure 18 

An attracting fixed point for system (2.2) for 𝑎=3.5, 𝑏=0.2, and 𝑑=2.67.

Figure 19 

Phase portrait for the system (2.2) before the NS point (𝑎=3.5, 𝑏=0.2, and 𝑑=2.78).

Figure 20 demonstrates the behaviour of the model after the NS bifurcation, for 𝑑=2.81. From Figures 19 and 20 it turns out that the fixed point 𝐸₃ loses its stability through an NS bifurcation, when 𝑑 varies from 2.78 to 2.81. The dynamics of (2.2) beyond the NS point is a stable closed invariant curve which coexists with unstable fixed point 𝐸₃.

Figure 20 

Phase portrait for the system (2.2) for 𝑎=3.5, 𝑏=0.2, and 𝑑=2.81.

As 𝑑 is increased further, however, the phase portrait starts to fold. We see that the circle, after being stretched, shrunk, and folded, creates new phenomena due to the breakdown of the closed curve; see Figure 21 for 𝑑=3.13.

Figure 21 

The breakdown of the closed invariant curve of the system (2.2) for 𝑎=3.5, 𝑏=0.2, and 𝑑=3.13.

For increasing 𝑑, we obtain the multiple invariant closed curves brought about the NS bifurcation point of iterates of (2.2). In these cases, higher bifurcations of the torus occurs as the system moves out of quasiperiodic region by increasing 𝑑. The dynamics move from one closed curve to another periodically, but the dynamics in each closed curve, may be periodic or quasiperiodic. Figure 22 presents the set of closed curves around the NS bifurcation.

Figure 22 

The existence of multiple closed curves of the system (2.2) for 𝑎=3.5, 𝑏=0.2, and 𝑑=3.215.

Moreover, the closed curves may break and lead to multiple fractal tori on which the dynamics of (2.2) are chaotic. Figures 23 and 24 present strange attractors for (2.2) with 𝑑=3.25 and 𝑑=3.8, respectively, which exhibit a fractal structure.

Figure 23 

Chaotic attractor for the system (2.2) for 𝑎=3.5, 𝑏=0.2, and 𝑑=3.25.

Figure 24 

Chaotic attractor for the system (2.2) for 𝑎=3.5, 𝑏=0.2, and 𝑑=3.8.

5. Concluding Remarks

In this paper, we studied a planar map that models a predator-prey interaction with nonoverlapping generations. We derived analytically a complete description of the stability regions of the fixed points of the system, namely, 𝐸₁, 𝐸₂, and 𝐸₃. We showed that the system undergoes branching, period doubling, and Neimark-Sacker bifurcations. In particular, we determined criticality of the flip, and Neimark-Sacker bifurcations for 𝐸₂ and 𝐸₃ analytically by direct computation of the normal form coefficients without explicit reduction to the center manifold. To support the analytical results and reveal the further complex behaviour of the system, we employed numerical continuation methods to compute curves of codimension 1 and 2 bifurcation points. In particular, we computed curves of fixed points of different cycles as well as curves of fold, flip and Neimark-Sacker bifurcations of the fourth iterate, and fold and flip bifurcations of the sixteenth iterate. These curves form stability boundaries of different cycles of the system. These tasks were possible by means of the computation of normal form coefficients and branch-switching algorithms. We further used numerical simulation methods to reveal chaotic behaviours and a strange attractor near the Neimark-Sacker bifurcation.

Acknowledgments

This paper was greatly improved by the comments of four reviewers who read the paper from different viewpoints. The authors thank them all and thank the editorial staff of DDNS for the fast and efficient handling of the manuscript. R. K. Ghaziani also would like to thank Shahrekord University for the financial support through a research grant.