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

Authors: R. Khoshsiar Ghaziani; W. Govaerts; C. Sonck
Journal: Discrete Dynamics in Nature and Society
 (2011)

---

## 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.


---

## Body

## 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)ẋ(t)=x(r1-a11x)-a12xy,ẏ(t)=y(-r2-a22y)+a21xy,
where x(t) and y(t) represent the densities of the prey and the predator, r1, a12, r2, and a21 are the intrinsic growth rate of the prey, the capture rate, the death rate of the predator, and the conversion rate, respectively, a11 and a22 denote the intraspecific competition coefficients of the prey and the predator, and r1/a11 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(2.1)F:(xy)⟼(ax(1-x)-bxydxy),
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(2.2)ẋ(t)=α0x(t)(1-x(t))-αmx(t)y(t),ẏ(t)=mx(t)y(t)-βy(t),
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. x and y represent the densities of the prey and predator, and α0, α, m, and β are nonnegative parameters. Applying the forward Euler scheme to system (2.2) with the stepsize 1/β and assuming α0/β≫1, we obtain the map (2.1) with nonnegative parameters a=α0/β, b=αm/β, and d=m/β.We note that in (2.1) one of the two parameters b,d can be removed by a rescaling of y. So, the system is in reality a two-parameter system and fully equivalent to the system studied in [16]. We will generally choose d, a 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 (x*,y*) to (2.3)ax*(1-x*)-bx*y*=x*,dx*y*=y*.
The origin E1=(0,0) is a fixed point of (2.1) but is of little interest. Two further nontrivial fixed points are E2=((a-1)/a,0) and E3=(1/d,1/b(a(1-(1/d))-1)). We note that E2 is biologically possible only if its coordinates are nonnegative, that is, a≥1. E3 is biologically possible only if the following conditions hold:  (i)
a>1,(ii)
d≥a/(a-1).We start the local bifurcation analysis of the map (2.1) by linearization of F around each of its fixed points. The Jacobian matrix J(x,y) is given by(2.4)J(x,y)=(a-2ax-by-bxdydx).The characteristic equation ofJ(x,y) is given by(2.5)λ2-tr⁡(J)+det⁡(J)=0,
where tr(J)=a-2ax-by+dx and det(J)=dxa-2dx2a.


### 2.1. Stability and Bifurcation ofE1

Proposition 2.1.
The fixed pointE1 is asymptotically stable for 0≤a<1. It loses stability via branching for a=1 and there bifurcates to E2.Proof.
Eigenvalues of the Jacobian atE1 are a and 0. So, E1 is stable if a<1 and loses stability at a=1. In (a,x)-space E1 forms the curve (a,0) with tangent vector (1,0). E2 is represented in (a,x)-space by the curve (a,(a-1)/a). When a=1, these curves intersect at (1,0) and the tangent vector in (a,x)-space is (1,1), so it is clear that E1 branches to E2 for a=1.


### 2.2. Stability and Bifurcation ofE2

The Jacobian matrix of (2.1) at E2 is given by(2.6)J(E2)=(-a+2-b(a-1)a0d(a-1)a).Proposition 2.2.
The fixed pointE2 is linearly asymptotically stable if and only if a∈]1,3[ and d<a/(a-1). Moreover, it loses stability:  (i)
via branching fora=1 and there bifurcates to E1,(ii)
via branching ford=a/(a-1) and there bifurcates to E3 if 1<a<3,(iii)
via a supercritical flip fora=3 if d<3/2.Proof.
The multipliers ofJ(E2) are λ1=-a+2 and λ2=d(a-1)/a. The fixed point E2 is asymptotically stable if and only if |λ1|<1 and |λ2|<1, that is, if and only if a∈]1,3[ and d<a/(a-1). Boundary points of the stability region must satisfy one of three conditions: a=1, d=a/(a-1), or a=3.
In the first case, the conditionsd<a/(a-1) and a<3 are satisfied for nearby values a>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 E1 points.
In the second case, this is a stability boundary only if1<a<3. The Jacobian (2.6) then has an eigenvalue +1 and it is checked easily that these boundary points are also E3 points.
In the third case, this is a stability boundary only ifd<3/2. In this case, λ1=-1 which means that E2 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 b1,
(2.7)b1=16〈p,C(q,q,q)+3B(q,(I-A)-1B(q,q))〉,
derived by center manifold reduction is positive; see [9], Ch. 8 and [14]. Here, A=J(E2), and B(·,·),C(·,·,·) are the second- and third- order multilinear forms, respectively, and p and q are the left and right eigenvectors of A for the eigenvalue -1, respectively. These vectors are normalized by 〈p,q〉=1, 〈q,q〉=1, where 〈··〉 is the standard scalar product in ℝ2. We obtain
(2.8)q=(q1q2)=(10),p=(p1p2)=(12b3+2d).
The components of the multilinear form B(q,q) are given by
(2.9)[B(q,q)]1=∑j,k=12∂2(ax(1-x)-bxy)∂xj∂xkqjqk=-2aq1q1=-6,[B(q,q)]2=∑j,k=12∂2(dxy)∂xj∂xkqjqk=2dq1q2=0,
where the state variable vector is for ease of notation generically denoted by (x1,x2)T instead of (x,y)T.
Letζ=(I-A)-1B(q,q), then we have ζ=(-30) and find
(2.10)[B(q,ζ)]1=-2aq1ζ1=-2a(-3)=18,[B(q,ζ)]2=dq1ζ2=0.
The third-order multilinear form C(q,q,q) is identically zero. The critical normal form coefficient b1 is given by
(2.11)b1=16pT(540)=9,
which is clearly positive. This completes the proof of supercriticality of the flip point at E2.The stability regionΩE2S of E2, as obtained in Proposition 2.2, is shown in Figure 1.Figure 1
Stability regions in(d,a)-space. ΩE1S, ΩE2S are the stability regions of E1 and E2, respectively. The stability region ΩE3S of E3 is the union of ΩE3S1 and ΩE3S2 (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 ofE3

To study the stability ofE3, we use the Jury's criteria; see [6, Section  A2.1]. Let F(λ)=λ2-tr(J(E3))λ+det(J(E3)) be the characteristic polynomial of J(E3).According to the Jury’s criteria,E3 is asymptotically stable if the following conditions hold:   (2.12)F(-1)=1+tr⁡(J(E3))+det⁡(J(E3))>0,F(1)=1-tr⁡(J(E3))+det⁡(J(E3))>0,1-det⁡(J(E3))>0.AtE3, we have:(2.13)J(E3)=(d-ad-bdda-a-db1).
We note that tr(J(E3))=(2d-a)/d and det(J(E3))=a(d-2)/d are independent of b.Proposition 2.3.
E3 is linearly asymptotically stable if and only if one of the following mutually exclusive conditions holds:  (i)
3/2<d<9/4 and d/(d-1)<a<3d/(3-d),(ii)
d=9/4 and 1.8=d/(d-1)<a<d/(d-2)=3d/(3-d)=9,(iii)
d>9/4 and d/(d-1)<a<d/(d-2).Proof.
The criterionF(1)>0 is easily seen to be equivalent to the condition
(2.14)a>dd-1,d>1,
or equivalently,
(2.15)d>aa-1,a>1.
Next, the criteriondet(J(F3))<1 is easily seen to be equivalent to
(2.16)a<dd-2ifd>2.
The criterionF(-1)>0 translates as
(2.17)a<3d3-difd<3.In Figure1, the parts of the 3 curves a=d/(d-1), a=d/(d-2), and a=3d/(3-d) where a and d are positive are depicted (resp., for d>1, d>2, and 0<d<3). The curves a=d/(d-1) and a=3d/(3-d) intersect solely at (d=3/2,a=3) and the curves a=d/(d-2) and a=3d/(3-d) at (d=9/4,a=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 ΩE3S of E3 is thus the union of ΩE3S1 and ΩE3S2 (which correspond to (i) and (iii) in this proposition, resp.), and the open interval that separates them (and corresponds to (ii)).Proposition 2.4.
E3 loses stability:  (i)
via a flip point when3/2<d<9/4 and a=3d/(3-d),(ii)
via a Neimark-Sacker point when9/4<d and a=d/(d-2),(iii)
via a branch point whend>3/2 and a=d/(d-1) where it bifurcates to E2,(iv)
via a branch-flip (BPPD) point whend=3/2 and a=3,(v)
via a resonance 1 : 2 point whend=9/4 and a=9.Proof.
By Proposition2.3 the stability boundary of E3 consists of parts of three curves, namely,  (1)
Curve 1:a=3d/(3-d).(2)
Curve 2:a=d/(d-2).(3)
Curve 3:a=d/(d-1).The points of Curve 1 which are on the stability boundary of E3 satisfy F(-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(J(E3))=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 F(1)=0, that is, they have an eigenvalue 1. It can be checked easily that E3 then branches to E2.
Combining this with Proposition2.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 (whend=3/2 and a=3), the point is a branch-flip point (BPPD) with eigenvalues −1 and 1. The intersection point of Curves 1 and 2 (when d=9/4 and a=9) is a resonance 1 : 2 point with double eigenvalue −1.Proposition 2.5.
The flip point in Proposition2.4, part (i), is subcritical.Proof.
For the subcriticality of the flip bifurcation, it is sufficient to show that the normal form coefficientb1 in (2.7) is always negative. By using the same procedure as we used in the proof of Proposition 2.2, we obtain
(2.18)q=(q1q2)=(1d(2d-3)b(d-3)),p=(p1p2)=(1b2d).
To simplify computations, we do not normalize p and q, since rescaling does not change the sign of b1 provided that 〈p,q〉 is positive (it can be proved easily that this is the case). B(q,q) is computed as
(2.19)B(q,q)=(-4d2d2(2d-3)b(d-3)).
Letζ=(I-A)-1B(q,q), then we have ζ=(d-d2(4+(3/(3-d)))/b) and find
(2.20)B(q,ζ)=(2d23-d(6-d)-2ad-2d3(6-d)b(3-d)).
The third-order multilinear form C(q,q,q) is zero. The critical normal form coefficient b1=(-1/2)d3/(3-d), which is clearly negative for 3/2<d<9/4. This completes the proof of subcriticality of the flip point.Proposition 2.6.
The Neimark-Sacker point in Proposition2.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 coefficientd1,
(2.21)d1=12e-iθ0〈p,C(q,q,q¯)+2B(q,h1100)+B(q¯,h2000)〉,h1100=(In-A)-1B(q,q¯),h2000=(e2iθ0In-A)-1B(q,q),
is negative; see [14]. Here, θ0(0<θ0<π) is the argument of the critical multiplier, A=J(E3), and B(·,·), C(·,·,·) are the second- and third-order multilinear forms, respectively. p and q are the left and right eigenvectors of A(2.22)Aq=eiθ0q,ATp=e-iθ0p.
These vectors are normalized by 〈p,q〉=1,〈q,q〉=1, where 〈p,q〉=pHq is the Hermitian inner product in ℂ2. In the Neimark-Sacker point, a=d/(d-2) (d>9/4, so d-2>0). We get
(2.23)A=(d-3d-2-bddb(d-2)1),
with characteristic polynomial λ2-((2d-5)/(d-2))λ+1. It follows that ℜ(λ)=(d-(5/2))/(d-2) and ℑ(λ)=d-(9/4)/(d-2), so eiθ0=1/(d-2)(d-(5/2)+id-(9/4)). The eigenvectors are
(2.24)q=((d-2)bd-12-id-94),p=d(2d-(9/2)-id-(9/4))b(4d2-17d+18)(1bd(12-id-94)),
where q is not normalized since rescaling of q does not change the sign of d1 (this can be proved easily) and we can simplify the computations by not normalizing. The calculation of the components of B(q,q¯) gives us
(2.25)B(q,q¯)=-b(d-2)(bd1),
and it follows that
(2.26)h1100=b(d-2)(b(d-2)d-2).
The remaining calculations are done with the help of Maple. The calculation of h2000 gives us
(2.27)h2000=-2(d-2)3b-50d2+135d+6i4d-9d2-29i4d-9d-119+35i4d-9+6d3⋅(bd(17d-7i4d-9d-23+11i4d-9-3d2+i4d-9d2)i4d-9d+22-2i4d-9+4d2-19d),
where
(2.28)e2iθ0=1(d-2)2(d2-6d+172+i(2d-5)d-94).
For B(q¯,h2000) and B(q,h1100), we obtain
(2.29)B(q¯,h2000)=2(d-2)3b2-50d2+135d+6i4d-9d2-29i4d-9d-119+35i4d-9+6d3⋅(bd(32d-6i4d-9d-29+9i4d-9-13d2+i4d-9d2+2d3)i4d-9d2+2d-8i4d-9d-2d3+7d2-17+13i4d-9),B(q,h1100)=(d-2)2b22(b(1+i4d-9)-5-i4d-9).
The third-order multilinear form C(q,q,q) is identically zero, so we get
(2.30)d1=12e-iθ0pH(2B(q,h1100)+B(q¯,h2000))=14(-2d+5+i4d-9)(4d-9+i4d-9)b2d×(31d2+i4d-9d3+44-5d3-64d-12i4d-9-7i4d-9d2+16i4d-9d)×(4d-9)-1(-50d2+135d+6i4d-9d2-29i4d-9d-119+35i4d-9+6d3)-1.
The real part of d1 equals -1/2b2d (exact up to a positive factor due to the nonnormalization of q), 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) ofE3 obtained in Proposition 2.3 are depicted as ΩE3S1 and ΩE3S2, respectively, in Figure 1. The region (ii) is the open interval on the common boundary of ΩE3S1 and ΩE3S2. We have a complete description of the stability region of E3 for all parameter combinations.The biological interpretation of Figure1 is as follows. The situation with no prey and no predators exists for all parameter values but is stable only in ΩE1S, that is, for a<1. The situation with a fixed number of prey but no predators exists for all a>1 but is stable only in ΩE2S. Coexistence of fixed nonzero numbers of prey and predators is possible whenever a>1 and d>a/(a-1) but is stable only in the union of ΩE3S1 and ΩE3S2.


## 2.1. Stability and Bifurcation ofE1

Proposition 2.1.
The fixed pointE1 is asymptotically stable for 0≤a<1. It loses stability via branching for a=1 and there bifurcates to E2.Proof.
Eigenvalues of the Jacobian atE1 are a and 0. So, E1 is stable if a<1 and loses stability at a=1. In (a,x)-space E1 forms the curve (a,0) with tangent vector (1,0). E2 is represented in (a,x)-space by the curve (a,(a-1)/a). When a=1, these curves intersect at (1,0) and the tangent vector in (a,x)-space is (1,1), so it is clear that E1 branches to E2 for a=1.


## 2.2. Stability and Bifurcation ofE2

The Jacobian matrix of (2.1) at E2 is given by(2.6)J(E2)=(-a+2-b(a-1)a0d(a-1)a).Proposition 2.2.
The fixed pointE2 is linearly asymptotically stable if and only if a∈]1,3[ and d<a/(a-1). Moreover, it loses stability:  (i)
via branching fora=1 and there bifurcates to E1,(ii)
via branching ford=a/(a-1) and there bifurcates to E3 if 1<a<3,(iii)
via a supercritical flip fora=3 if d<3/2.Proof.
The multipliers ofJ(E2) are λ1=-a+2 and λ2=d(a-1)/a. The fixed point E2 is asymptotically stable if and only if |λ1|<1 and |λ2|<1, that is, if and only if a∈]1,3[ and d<a/(a-1). Boundary points of the stability region must satisfy one of three conditions: a=1, d=a/(a-1), or a=3.
In the first case, the conditionsd<a/(a-1) and a<3 are satisfied for nearby values a>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 E1 points.
In the second case, this is a stability boundary only if1<a<3. The Jacobian (2.6) then has an eigenvalue +1 and it is checked easily that these boundary points are also E3 points.
In the third case, this is a stability boundary only ifd<3/2. In this case, λ1=-1 which means that E2 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 b1,
(2.7)b1=16〈p,C(q,q,q)+3B(q,(I-A)-1B(q,q))〉,
derived by center manifold reduction is positive; see [9], Ch. 8 and [14]. Here, A=J(E2), and B(·,·),C(·,·,·) are the second- and third- order multilinear forms, respectively, and p and q are the left and right eigenvectors of A for the eigenvalue -1, respectively. These vectors are normalized by 〈p,q〉=1, 〈q,q〉=1, where 〈··〉 is the standard scalar product in ℝ2. We obtain
(2.8)q=(q1q2)=(10),p=(p1p2)=(12b3+2d).
The components of the multilinear form B(q,q) are given by
(2.9)[B(q,q)]1=∑j,k=12∂2(ax(1-x)-bxy)∂xj∂xkqjqk=-2aq1q1=-6,[B(q,q)]2=∑j,k=12∂2(dxy)∂xj∂xkqjqk=2dq1q2=0,
where the state variable vector is for ease of notation generically denoted by (x1,x2)T instead of (x,y)T.
Letζ=(I-A)-1B(q,q), then we have ζ=(-30) and find
(2.10)[B(q,ζ)]1=-2aq1ζ1=-2a(-3)=18,[B(q,ζ)]2=dq1ζ2=0.
The third-order multilinear form C(q,q,q) is identically zero. The critical normal form coefficient b1 is given by
(2.11)b1=16pT(540)=9,
which is clearly positive. This completes the proof of supercriticality of the flip point at E2.The stability regionΩE2S of E2, as obtained in Proposition 2.2, is shown in Figure 1.Figure 1
Stability regions in(d,a)-space. ΩE1S, ΩE2S are the stability regions of E1 and E2, respectively. The stability region ΩE3S of E3 is the union of ΩE3S1 and ΩE3S2 (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 ofE3

To study the stability ofE3, we use the Jury's criteria; see [6, Section  A2.1]. Let F(λ)=λ2-tr(J(E3))λ+det(J(E3)) be the characteristic polynomial of J(E3).According to the Jury’s criteria,E3 is asymptotically stable if the following conditions hold:   (2.12)F(-1)=1+tr⁡(J(E3))+det⁡(J(E3))>0,F(1)=1-tr⁡(J(E3))+det⁡(J(E3))>0,1-det⁡(J(E3))>0.AtE3, we have:(2.13)J(E3)=(d-ad-bdda-a-db1).
We note that tr(J(E3))=(2d-a)/d and det(J(E3))=a(d-2)/d are independent of b.Proposition 2.3.
E3 is linearly asymptotically stable if and only if one of the following mutually exclusive conditions holds:  (i)
3/2<d<9/4 and d/(d-1)<a<3d/(3-d),(ii)
d=9/4 and 1.8=d/(d-1)<a<d/(d-2)=3d/(3-d)=9,(iii)
d>9/4 and d/(d-1)<a<d/(d-2).Proof.
The criterionF(1)>0 is easily seen to be equivalent to the condition
(2.14)a>dd-1,d>1,
or equivalently,
(2.15)d>aa-1,a>1.
Next, the criteriondet(J(F3))<1 is easily seen to be equivalent to
(2.16)a<dd-2ifd>2.
The criterionF(-1)>0 translates as
(2.17)a<3d3-difd<3.In Figure1, the parts of the 3 curves a=d/(d-1), a=d/(d-2), and a=3d/(3-d) where a and d are positive are depicted (resp., for d>1, d>2, and 0<d<3). The curves a=d/(d-1) and a=3d/(3-d) intersect solely at (d=3/2,a=3) and the curves a=d/(d-2) and a=3d/(3-d) at (d=9/4,a=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 ΩE3S of E3 is thus the union of ΩE3S1 and ΩE3S2 (which correspond to (i) and (iii) in this proposition, resp.), and the open interval that separates them (and corresponds to (ii)).Proposition 2.4.
E3 loses stability:  (i)
via a flip point when3/2<d<9/4 and a=3d/(3-d),(ii)
via a Neimark-Sacker point when9/4<d and a=d/(d-2),(iii)
via a branch point whend>3/2 and a=d/(d-1) where it bifurcates to E2,(iv)
via a branch-flip (BPPD) point whend=3/2 and a=3,(v)
via a resonance 1 : 2 point whend=9/4 and a=9.Proof.
By Proposition2.3 the stability boundary of E3 consists of parts of three curves, namely,  (1)
Curve 1:a=3d/(3-d).(2)
Curve 2:a=d/(d-2).(3)
Curve 3:a=d/(d-1).The points of Curve 1 which are on the stability boundary of E3 satisfy F(-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(J(E3))=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 F(1)=0, that is, they have an eigenvalue 1. It can be checked easily that E3 then branches to E2.
Combining this with Proposition2.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 (whend=3/2 and a=3), the point is a branch-flip point (BPPD) with eigenvalues −1 and 1. The intersection point of Curves 1 and 2 (when d=9/4 and a=9) is a resonance 1 : 2 point with double eigenvalue −1.Proposition 2.5.
The flip point in Proposition2.4, part (i), is subcritical.Proof.
For the subcriticality of the flip bifurcation, it is sufficient to show that the normal form coefficientb1 in (2.7) is always negative. By using the same procedure as we used in the proof of Proposition 2.2, we obtain
(2.18)q=(q1q2)=(1d(2d-3)b(d-3)),p=(p1p2)=(1b2d).
To simplify computations, we do not normalize p and q, since rescaling does not change the sign of b1 provided that 〈p,q〉 is positive (it can be proved easily that this is the case). B(q,q) is computed as
(2.19)B(q,q)=(-4d2d2(2d-3)b(d-3)).
Letζ=(I-A)-1B(q,q), then we have ζ=(d-d2(4+(3/(3-d)))/b) and find
(2.20)B(q,ζ)=(2d23-d(6-d)-2ad-2d3(6-d)b(3-d)).
The third-order multilinear form C(q,q,q) is zero. The critical normal form coefficient b1=(-1/2)d3/(3-d), which is clearly negative for 3/2<d<9/4. This completes the proof of subcriticality of the flip point.Proposition 2.6.
The Neimark-Sacker point in Proposition2.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 coefficientd1,
(2.21)d1=12e-iθ0〈p,C(q,q,q¯)+2B(q,h1100)+B(q¯,h2000)〉,h1100=(In-A)-1B(q,q¯),h2000=(e2iθ0In-A)-1B(q,q),
is negative; see [14]. Here, θ0(0<θ0<π) is the argument of the critical multiplier, A=J(E3), and B(·,·), C(·,·,·) are the second- and third-order multilinear forms, respectively. p and q are the left and right eigenvectors of A(2.22)Aq=eiθ0q,ATp=e-iθ0p.
These vectors are normalized by 〈p,q〉=1,〈q,q〉=1, where 〈p,q〉=pHq is the Hermitian inner product in ℂ2. In the Neimark-Sacker point, a=d/(d-2) (d>9/4, so d-2>0). We get
(2.23)A=(d-3d-2-bddb(d-2)1),
with characteristic polynomial λ2-((2d-5)/(d-2))λ+1. It follows that ℜ(λ)=(d-(5/2))/(d-2) and ℑ(λ)=d-(9/4)/(d-2), so eiθ0=1/(d-2)(d-(5/2)+id-(9/4)). The eigenvectors are
(2.24)q=((d-2)bd-12-id-94),p=d(2d-(9/2)-id-(9/4))b(4d2-17d+18)(1bd(12-id-94)),
where q is not normalized since rescaling of q does not change the sign of d1 (this can be proved easily) and we can simplify the computations by not normalizing. The calculation of the components of B(q,q¯) gives us
(2.25)B(q,q¯)=-b(d-2)(bd1),
and it follows that
(2.26)h1100=b(d-2)(b(d-2)d-2).
The remaining calculations are done with the help of Maple. The calculation of h2000 gives us
(2.27)h2000=-2(d-2)3b-50d2+135d+6i4d-9d2-29i4d-9d-119+35i4d-9+6d3⋅(bd(17d-7i4d-9d-23+11i4d-9-3d2+i4d-9d2)i4d-9d+22-2i4d-9+4d2-19d),
where
(2.28)e2iθ0=1(d-2)2(d2-6d+172+i(2d-5)d-94).
For B(q¯,h2000) and B(q,h1100), we obtain
(2.29)B(q¯,h2000)=2(d-2)3b2-50d2+135d+6i4d-9d2-29i4d-9d-119+35i4d-9+6d3⋅(bd(32d-6i4d-9d-29+9i4d-9-13d2+i4d-9d2+2d3)i4d-9d2+2d-8i4d-9d-2d3+7d2-17+13i4d-9),B(q,h1100)=(d-2)2b22(b(1+i4d-9)-5-i4d-9).
The third-order multilinear form C(q,q,q) is identically zero, so we get
(2.30)d1=12e-iθ0pH(2B(q,h1100)+B(q¯,h2000))=14(-2d+5+i4d-9)(4d-9+i4d-9)b2d×(31d2+i4d-9d3+44-5d3-64d-12i4d-9-7i4d-9d2+16i4d-9d)×(4d-9)-1(-50d2+135d+6i4d-9d2-29i4d-9d-119+35i4d-9+6d3)-1.
The real part of d1 equals -1/2b2d (exact up to a positive factor due to the nonnormalization of q), 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) ofE3 obtained in Proposition 2.3 are depicted as ΩE3S1 and ΩE3S2, respectively, in Figure 1. The region (ii) is the open interval on the common boundary of ΩE3S1 and ΩE3S2. We have a complete description of the stability region of E3 for all parameter combinations.The biological interpretation of Figure1 is as follows. The situation with no prey and no predators exists for all parameter values but is stable only in ΩE1S, that is, for a<1. The situation with a fixed number of prey but no predators exists for all a>1 but is stable only in ΩE2S. Coexistence of fixed nonzero numbers of prey and predators is possible whenever a>1 and d>a/(a-1) but is stable only in the union of ΩE3S1 and ΩE3S2.


## 3. Numerical Bifurcation Analysis ofE2 and E3

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 Figures2, 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 a=3.571920967580968, b=0.2, d=1.4.


### 3.1. Numerical Bifurcation ofE2

By continuation ofE2=(0.5;0) starting from a=2, b=0.2, d=1.4, in the stable region of E2 with a free, we see that E2 is stable when 1<a<3. It loses stability via a supercritical period doubling point (PD, the corresponding normal form coefficient is 9>0) when a=3, and via a branch point (BP) when a 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+000The first two entries ofx are the coordinate values of the fixed point E2, and the last entry of x is the value of the free parameter a 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+002A stable 4-cycle is born whena>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+004Thus, whena>3.544090 a stable 8-cycle emerges. A stable 8-cycle is given by C8={X18,…,X88}, where X18=(0.880450408535962,0) where a=3.571920967580968, b=0.2, and d=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 fora=3.5, b=0.2, and d=2.83.We note that forE2 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 a; see [19].Continuation ofE2 starting from the same parameter values as in Run 1, with d as free parameter, leads to: label = BP, x = (0.500000 0.000000 2.000000)The appearance of a branch point is consistent with Proposition2.2 part (ii) which states that E2 bifurcates to E3 when b=d=a/(a-1)=2.


### 3.2. Numerical Bifurcation ofE3

We now considerE3=(0.454545;4.545454) which is in the stable region for the parameter values a=3.5, b=0.2, and d=2.2 (stability follows from Proposition 2.3 part (i)). We do a numerical continuation of E3 with control parameter d, 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−001label = NS, x = (0.357143 6.250000 2.800000)normal form coefficient of NS = −6.971545e−002E3 is stable when 1.615385<d<2.8. It loses stability via a supercritical Neimark-Sacker (NS) point when d=2.8, which is consistent with Proposition 2.4 part (ii) (d/(d-2)=3.5=a). It also loses stability through a subcritical PD point when d=1.615385, which is consistent with Propositions 2.5 and 2.4 part (i) since 3d/(3-d)=3.5=a.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 d=2.83, a stable closed invariant curve is created around the unstable fixed point E3 (see Figure 4). Now, we compute the period doubling curve, with a and d 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−006label = R2, x = (0.444444 20.000000 9.000000 2.250000)Normal form coefficient for R2:[c, d] = 1.197634e−001, −7.185806e−001Two codim-2 bifurcation points are detected on the flip curve, namely, a fold-flip LPPD and a resonance 2 bifurcationR2. 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 a and d 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 ilabel = R3, x = (0.428571 15.000000 7.000000 2.333333 −0.500000)Normal form coefficient of R3: Re(c_1) = −5.000000e−001label = R2, x = (0.444444 20.000000 9.000000 2.250000 −1.000000)Normal form coefficient of R2: [c, d] = 1.197646e−001, −7.185574e−001Figure 5
Flip and Neimark-Sacker bifurcation curves starting from points inRun 2. We note that the LPPD point is in reality a BPPD point.Figure 6
A stable 4-cycle near anR4 point for a=4.99, b=0.2, and d=2.56.The computed curve of NS points is also shown in Figure5 (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 a=d/(d-2) at a resonance 1 : 2 (R2) point, which is consistent with Proposition 2.4 part (v). It also hits the BP curve a=d/(d-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 coefficientA of the R4 point in Run 4 satisfies |A|>1, hence two cycles of period 4 of the map are born. A stable 4-cycle for a=4.99, b=0.2, and d=2.56 is given by: C4={X1,X2,X3,X4}, where X1=(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 |A|>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+001label = LPPD, x = (0.540029 8.419070 5.511767 2.615029)Normal form coefficient for LPPD:[a/e, be] = 7.735689e−001, −8.261238e−002Figure 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 LP4 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+001label = LPPD, x = (0.540029 8.419070 5.511767 2.615029)Normal form coefficient for LPPD:[a/e, be] = 7.735689e−001, −8.261238e−002label = GPD, x = (0.541277 8.238085 5.478361 2.609243)Normal form coefficient of GPD = 6.483721e−001label = R2, x = (0.582936 7.330070 5.625162 2.664169)Normal form coefficient for R2:[c, d] = 6.130829e−002, −3.509615e+000By superposing the flip curve on Figure7, we obtain Figure 8. We further compute a curve of fixed points of the fourth iterate starting from the 4-cycle C4 presented in Figure 6, with control parameter d. 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−001label = PD, x = (0.554614 8.335785 2.864961)normal form coefficient of PD = 1.909178e−001label = NS, x = (0.546431 9.030507 2.608544)normal form coefficient of NS = −2.177981e+000Figure 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-cycleC4.The 4-cycle remains stable when2.502419<d<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 ilabel = R3, x = (0.529139 9.928228 5.218594 2.614545 −0.500000)Normal form coefficient of R3: Re(c_1) = −6.065820e+000label = R2, x = (0.513005 11.383897 5.625162 2.664169 −1.000000)Normal form coefficient of R2: [c, d] = 1.004097e−001, −5.657220e+000Figure 10
Stability regionΩS4 of the 4-cycle C4. The boundaries consist of NS, LP4, PD4, and NS4 curves.Figure 11
A stable 16-cycle around theR4 point of Run 9 for a=4.68, b=0.2, and d=2.66.The stability regionΩS4 of the 4-cycle C4 is bounded by the LP4, PD4, and NS4 curves (see Figure 10).Now, we consider theR4 point computed in Run 9. Since the corresponding normal form coefficient A satisfies |A|>1, two cycles of period 16 of the map are born. A stable 16-cycle for a=4.68, b=0.2, and d=2.66 is given by C16={X1,X2,…,X16}, where X1=(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 theR4 point. These curves exist since |A|>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−001First Lyapunov coefficient for second iterate = 3.799853e−001label = LPPD, x = (0.519225 8.927562 4.765661 2.658381)Normal form coefficient for LPPD:[a/e, be] = 9.129886e+000, 5.166118e+000First Lyapunov coefficient for second iterate = 5.166118e+000label = LPPD, x = (0.469475 10.678780 5.153847 2.636682)Normal form coefficient for LPPD:[a/e, be] = 1.298255e−001, −4.388341e+000label = LPPD, x = (0.580415 9.268192 5.370935 2.641181)Normal form coefficient for LPPD:[a/e, be] = 1.528369e+000, 6.631508e+001First Lyapunov coefficient for second iterate = 6.631508e+001label = LPPD, x = (0.578148 9.853334 5.628124 2.676611)Normal form coefficient for LPPD:[a/e, be] = −9.656910e−001, −2.522009e+004label = R1, x = (0.577592 9.941486 5.660840 2.682153)normal form coefficient of R1 = −1Figure 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-cycleC16 presented in Figure 11, with control parameter a. 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+000label = PD, x = (0.595604 7.701225 4.700640)normal form coefficient of PD = 2.829176e+002Figure 13
Curve of fixed points of the sixteenth iterate starting from the 16-cycleC16.Figure 14
A stable 32-cycle fora=4.7071, b=0.2, and d=2.66.The 16-cycle remains stable when4.646130<d<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 fora=4.7071, b=0.2, and d=2.66 is given by: C32={X1,X2,…,X32}, where X1=(0.619778764655944,7.243979650146042). We present this cycle in Figure 14.Now, we compute a PD-curve starting from the computed PD point inRun 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 ΩS16 of the 16-cycle consist of the LP16 and PD16 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+003First Lyapunov coefficient for second iterate = 1.579190e+003label = LPPD, x = (0.587235 7.998201 4.765661 2.658381)Normal form coefficient for LPPD:[a/e, be] = 9.129867e+000, 9.328278e+002First Lyapunov coefficient for second iterate = 9.328278e+002label = GPD, x = (0.529363 9.922888 5.130266 2.635413)Normal form coefficient of GPD = 1.292549e+007label = LPPD, x = (0.526233 10.046755 5.153847 2.636682)Normal form coefficient for LPPD:[a/e, be] = 1.298254e−001, −7.143714e+002label = LPPD, x = (0.521191 10.840009 5.370935 2.641181)Normal form coefficient for LPPD:[a/e, be] = 1.528369e+000, 1.217953e+002First Lyapunov coefficient for second iterate = 1.217953e+002label = GPD, x = (0.521992 10.863772 5.384512 2.642770)Normal form coefficient of GPD = −6.188850e+003Figure 15
Stability regionΩS16 of the 16-cycle C16. 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 forb=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 Figure7. 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 Figure17, we present an Arnold Tongue rooted in a weak 2 : 7 resonant Neimark-Sacker point. Its computation started from a stable 7-cycle with x=0.401985064603439, y=12.930999111981340, a=5.844728289174310, b=0.2, and d=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.


## 3.1. Numerical Bifurcation ofE2

By continuation ofE2=(0.5;0) starting from a=2, b=0.2, d=1.4, in the stable region of E2 with a free, we see that E2 is stable when 1<a<3. It loses stability via a supercritical period doubling point (PD, the corresponding normal form coefficient is 9>0) when a=3, and via a branch point (BP) when a 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+000The first two entries ofx are the coordinate values of the fixed point E2, and the last entry of x is the value of the free parameter a 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+002A stable 4-cycle is born whena>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+004Thus, whena>3.544090 a stable 8-cycle emerges. A stable 8-cycle is given by C8={X18,…,X88}, where X18=(0.880450408535962,0) where a=3.571920967580968, b=0.2, and d=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 fora=3.5, b=0.2, and d=2.83.We note that forE2 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 a; see [19].Continuation ofE2 starting from the same parameter values as in Run 1, with d as free parameter, leads to: label = BP, x = (0.500000 0.000000 2.000000)The appearance of a branch point is consistent with Proposition2.2 part (ii) which states that E2 bifurcates to E3 when b=d=a/(a-1)=2.


## 3.2. Numerical Bifurcation ofE3

We now considerE3=(0.454545;4.545454) which is in the stable region for the parameter values a=3.5, b=0.2, and d=2.2 (stability follows from Proposition 2.3 part (i)). We do a numerical continuation of E3 with control parameter d, 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−001label = NS, x = (0.357143 6.250000 2.800000)normal form coefficient of NS = −6.971545e−002E3 is stable when 1.615385<d<2.8. It loses stability via a supercritical Neimark-Sacker (NS) point when d=2.8, which is consistent with Proposition 2.4 part (ii) (d/(d-2)=3.5=a). It also loses stability through a subcritical PD point when d=1.615385, which is consistent with Propositions 2.5 and 2.4 part (i) since 3d/(3-d)=3.5=a.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 d=2.83, a stable closed invariant curve is created around the unstable fixed point E3 (see Figure 4). Now, we compute the period doubling curve, with a and d 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−006label = R2, x = (0.444444 20.000000 9.000000 2.250000)Normal form coefficient for R2:[c, d] = 1.197634e−001, −7.185806e−001Two codim-2 bifurcation points are detected on the flip curve, namely, a fold-flip LPPD and a resonance 2 bifurcationR2. 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 a and d 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 ilabel = R3, x = (0.428571 15.000000 7.000000 2.333333 −0.500000)Normal form coefficient of R3: Re(c_1) = −5.000000e−001label = R2, x = (0.444444 20.000000 9.000000 2.250000 −1.000000)Normal form coefficient of R2: [c, d] = 1.197646e−001, −7.185574e−001Figure 5
Flip and Neimark-Sacker bifurcation curves starting from points inRun 2. We note that the LPPD point is in reality a BPPD point.Figure 6
A stable 4-cycle near anR4 point for a=4.99, b=0.2, and d=2.56.The computed curve of NS points is also shown in Figure5 (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 a=d/(d-2) at a resonance 1 : 2 (R2) point, which is consistent with Proposition 2.4 part (v). It also hits the BP curve a=d/(d-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 coefficientA of the R4 point in Run 4 satisfies |A|>1, hence two cycles of period 4 of the map are born. A stable 4-cycle for a=4.99, b=0.2, and d=2.56 is given by: C4={X1,X2,X3,X4}, where X1=(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 |A|>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+001label = LPPD, x = (0.540029 8.419070 5.511767 2.615029)Normal form coefficient for LPPD:[a/e, be] = 7.735689e−001, −8.261238e−002Figure 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 LP4 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+001label = LPPD, x = (0.540029 8.419070 5.511767 2.615029)Normal form coefficient for LPPD:[a/e, be] = 7.735689e−001, −8.261238e−002label = GPD, x = (0.541277 8.238085 5.478361 2.609243)Normal form coefficient of GPD = 6.483721e−001label = R2, x = (0.582936 7.330070 5.625162 2.664169)Normal form coefficient for R2:[c, d] = 6.130829e−002, −3.509615e+000By superposing the flip curve on Figure7, we obtain Figure 8. We further compute a curve of fixed points of the fourth iterate starting from the 4-cycle C4 presented in Figure 6, with control parameter d. 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−001label = PD, x = (0.554614 8.335785 2.864961)normal form coefficient of PD = 1.909178e−001label = NS, x = (0.546431 9.030507 2.608544)normal form coefficient of NS = −2.177981e+000Figure 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-cycleC4.The 4-cycle remains stable when2.502419<d<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 ilabel = R3, x = (0.529139 9.928228 5.218594 2.614545 −0.500000)Normal form coefficient of R3: Re(c_1) = −6.065820e+000label = R2, x = (0.513005 11.383897 5.625162 2.664169 −1.000000)Normal form coefficient of R2: [c, d] = 1.004097e−001, −5.657220e+000Figure 10
Stability regionΩS4 of the 4-cycle C4. The boundaries consist of NS, LP4, PD4, and NS4 curves.Figure 11
A stable 16-cycle around theR4 point of Run 9 for a=4.68, b=0.2, and d=2.66.The stability regionΩS4 of the 4-cycle C4 is bounded by the LP4, PD4, and NS4 curves (see Figure 10).Now, we consider theR4 point computed in Run 9. Since the corresponding normal form coefficient A satisfies |A|>1, two cycles of period 16 of the map are born. A stable 16-cycle for a=4.68, b=0.2, and d=2.66 is given by C16={X1,X2,…,X16}, where X1=(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 theR4 point. These curves exist since |A|>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−001First Lyapunov coefficient for second iterate = 3.799853e−001label = LPPD, x = (0.519225 8.927562 4.765661 2.658381)Normal form coefficient for LPPD:[a/e, be] = 9.129886e+000, 5.166118e+000First Lyapunov coefficient for second iterate = 5.166118e+000label = LPPD, x = (0.469475 10.678780 5.153847 2.636682)Normal form coefficient for LPPD:[a/e, be] = 1.298255e−001, −4.388341e+000label = LPPD, x = (0.580415 9.268192 5.370935 2.641181)Normal form coefficient for LPPD:[a/e, be] = 1.528369e+000, 6.631508e+001First Lyapunov coefficient for second iterate = 6.631508e+001label = LPPD, x = (0.578148 9.853334 5.628124 2.676611)Normal form coefficient for LPPD:[a/e, be] = −9.656910e−001, −2.522009e+004label = R1, x = (0.577592 9.941486 5.660840 2.682153)normal form coefficient of R1 = −1Figure 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-cycleC16 presented in Figure 11, with control parameter a. 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+000label = PD, x = (0.595604 7.701225 4.700640)normal form coefficient of PD = 2.829176e+002Figure 13
Curve of fixed points of the sixteenth iterate starting from the 16-cycleC16.Figure 14
A stable 32-cycle fora=4.7071, b=0.2, and d=2.66.The 16-cycle remains stable when4.646130<d<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 fora=4.7071, b=0.2, and d=2.66 is given by: C32={X1,X2,…,X32}, where X1=(0.619778764655944,7.243979650146042). We present this cycle in Figure 14.Now, we compute a PD-curve starting from the computed PD point inRun 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 ΩS16 of the 16-cycle consist of the LP16 and PD16 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+003First Lyapunov coefficient for second iterate = 1.579190e+003label = LPPD, x = (0.587235 7.998201 4.765661 2.658381)Normal form coefficient for LPPD:[a/e, be] = 9.129867e+000, 9.328278e+002First Lyapunov coefficient for second iterate = 9.328278e+002label = GPD, x = (0.529363 9.922888 5.130266 2.635413)Normal form coefficient of GPD = 1.292549e+007label = LPPD, x = (0.526233 10.046755 5.153847 2.636682)Normal form coefficient for LPPD:[a/e, be] = 1.298254e−001, −7.143714e+002label = LPPD, x = (0.521191 10.840009 5.370935 2.641181)Normal form coefficient for LPPD:[a/e, be] = 1.528369e+000, 1.217953e+002First Lyapunov coefficient for second iterate = 1.217953e+002label = GPD, x = (0.521992 10.863772 5.384512 2.642770)Normal form coefficient of GPD = −6.188850e+003Figure 15
Stability regionΩS16 of the 16-cycle C16. 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 forb=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 Figure7. 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 Figure17, we present an Arnold Tongue rooted in a weak 2 : 7 resonant Neimark-Sacker point. Its computation started from a stable 7-cycle with x=0.401985064603439, y=12.930999111981340, a=5.844728289174310, b=0.2, and d=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 d. We fix the parameters a=3.5, b=0.2 and assume that d is free.Figure18 shows that E3 is a stable attractor for d=2.67. We note that for the given parameters value, E3 is a stable fixed point consistent with Proposition 2.3 part (iii). The behaviour of (2.2) at d=2.78, so before the NS point at d=2.8, is depicted in Figure 19.Figure 18
An attracting fixed point for system (2.2) for a=3.5, b=0.2, and d=2.67.Figure 19
Phase portrait for the system (2.2) before the NS point (a=3.5, b=0.2, and d=2.78).Figure20 demonstrates the behaviour of the model after the NS bifurcation, for d=2.81. From Figures 19 and 20 it turns out that the fixed point E3 loses its stability through an NS bifurcation, when d 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 E3.Figure 20
Phase portrait for the system (2.2) for a=3.5, b=0.2, and d=2.81.Asd 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 d=3.13.Figure 21
The breakdown of the closed invariant curve of the system (2.2) for a=3.5, b=0.2, and d=3.13.For increasingd, 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 d. 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 a=3.5, b=0.2, and d=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 d=3.25 and d=3.8, respectively, which exhibit a fractal structure.Figure 23
Chaotic attractor for the system (2.2) for a=3.5, b=0.2, and d=3.25.Figure 24
Chaotic attractor for the system (2.2) for a=3.5, b=0.2, and d=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,E1, E2, and E3. 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 E2 and E3 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.



---

*Source: 862494-2011-06-09.xml*