# The Bifurcation of Two Invariant Closed Curves in a Discrete Model

Authors: Yingying Zhang; Yicang Zhou
Journal: Discrete Dynamics in Nature and Society
 (2018)

---

## Abstract

A discrete population model integrated using the forward Euler method is investigated. The qualitative bifurcation analysis indicates that the model exhibits rich dynamical behaviors including the existence of the equilibrium state, the flip bifurcation, the Neimark-Sacker bifurcation, and two invariant closed curves. The conditions for existence of these bifurcations are derived by using the center manifold and bifurcation theory. Numerical simulations and bifurcation diagrams exhibit the complex dynamical behaviors, especially the occurrence of two invariant closed curves.


---

## Body

## 1. Introduction

Differential equations and difference equations are generally applied in mathematical biology modeling. Compared with differential equations, there are less different equations in ecological or epidemiological modeling; the main reason is the incompleteness of theories and methodologies for dynamical studies of discrete models. However, there are increasing attentions and research results on discrete models recently. One major reason is that discrete models as computational models are more efficient for numerical simulations and research results exhibit more richer dynamics than continuous ones. Another reason is that difference equations are more realistic and intelligible to characterize the population size of a single species when it has a fixed interval between generations or measurements.Up to now, the researches on discrete models are still focused on the dynamical behaviors (including stability, periodic solutions, bifurcations, chaos, and chaotic control; see [1–6]), and most of the scholars have studied the Neimark-Sacker bifurcation, a codimension-1 bifurcation, which has shown one invariant closed curve bifurcating from an equilibrium point. However, they do not consider further, the same as the phenomenon of limit cycles bifurcation in continuous models, the fact that the discrete models can exhibit the interesting dynamic behavior of two invariant closed curves bifurcating from an equilibrium point. One reason is the heavy computation, and another reason is the lack of proper theories and methods for the discrete models. In this paper we intend to research on this dynamic behavior in a discrete population model by using the bifurcation theory and normal form method; it is a fresh attempt different from previous related works. Meanwhile, we also propose a new research aspect by this work.We consider a continuous of Kolmogorov model [8–10] given by(1)dxdt=r1-xKx-mx-qx2x2+ay,dydt=px2x2+a-cy,where x, y represent the sizes of the prey and the predator population, respectively. It is assumed that the prey admits the Logistic growth if no predator exists. r is the intrinsic growth rate; K is the carrying capacity of the environment; m is the minimum of viable population; q is the maximal per capita consumption rate; a is the number of prey necessary to achieve one-half of the maximum rate q; p is the conversion efficiency of the consumed prey into new predator; c is the natural death rate of predators when predator y≡0 [11]. From the biological meaning all parameters are assumed to be positive and m≪K.The change of variables and the time rescaling given by the functionφ [12, 13](2)φu,v,τ=Ku,rK2qv,u2+a/K2rKτ=x,y,tturns model (1) into a polynomial system(3)dudτ=1-uu-MA+u2-uvu,dvdτ=Bu2-C2v,where A=a/K2, B=p-c/rK, C=ac/K2(p-c) for p-c>0, and M=m/K. The new parameters satisfy the following conditions:(4)0<A<1,B>0,C>0,0<M≪1.It is clear that det⁡Dφ(u,v,τ)=1/q(a+K2u2)>0, and φ is a diffeomorphism [14].We apply the forward Euler discrete scheme to system (3) and obtain the following discrete model:(5)ut+1=ut+δ1-utut-MA+ut2-utvtut,vt+1=vt+δBut2-C2vt.


## 2. The Existence and Stability of Equilibria

The equilibrium point(u∗,v∗) of (5) satisfies following equations:(6)u∗=u∗+δ1-u∗u∗-MA+u∗2-u∗v∗u∗,v∗=v∗+δBu∗2-C2v∗.There are four equilibrium points; it is easy to get that they are (0,0),(1,0),(M,0), and (C,L), where(7)L=1C1-CC-MA+C2.The existence of the unique positive equilibrium point (C,L) requires that 0<M<C<1,0<A<1,B>0,δ>0. The Jacobian matrix of model (5) at equilibrium point (u∗,v∗) is(8)J=1+δG∗-δu∗22δBu∗v∗1+δBu∗2-C2,where G∗=-AM+2(A+AM)u∗-3(A+M)u∗2+4(1+M)u∗3-5u∗4-2u∗v∗.Letλ1 and λ2 be the two eigenvalues of matrix J; we use the following definitions and conclusions [15–17]:( 1 ) If λ1<1 and λ2<1, then (u∗,v∗) is called a sink, and it is locally asymptotical stable.( 2 ) If λ1>1 and λ2>1, then (u∗,v∗) is called a source, and it is unstable.( 3 ) If λ1>1 and λ2<1 (or λ1<1 and λ2>1), then (u∗,v∗) is called a saddle, and it is unstable.( 4 ) If λ1=1 and λ2≠1 (or λ1≠1 and λ2=1), then (u∗,v∗) is called nonhyperbolic, and the corresponding bifurcation may occur.Lemma 1.
The eigenvalue and the stability of the trivial fixed points are as follows [1]:(1)
The eigenvalues of(0,0) are λ1=1-δAM<1,λ2=1-δBC2<1.
(1.1)
(0,0) is a sink if δ<2/AM and δ<2/BC2.(1.2)
(0,0) is a source if δ>2/AM and δ>2/BC2.(1.3)
(0,0) is a saddle if δ>2/AM and δ<2/BC2(orδ<2/AMandδ>2/BC2).(1.4)
(0,0) is nonhyperbolic, and there may be a flip bifurcation, if δ=2/AM and δ≠2/BC2(or δ≠2/AM and δ=2/BC2).(1.5)
(0,0) is nonhyperbolic, and there may be a strong1 :2 resonance bifurcationifδ=2/AM and δ=2/BC2.(2)
The eigenvalues of(1,0) are λ1=1+δ(A+1)(M-1)<1,λ2=1+δB(1-C2)>1.
(2.1)
(1,0) is a source if δ>2/(A+1)(1-M).(2.2)
(1,0) is a saddle if δ<2/(A+1)(1-M).(2.3)
(1,0) is nonhyperbolic, and there may be a flip bifurcation if δ=2/(A+1)(1-M).(3)
The eigenvalues of(M,0) are λ1=1+δM(A+M2)(1-M)>1,λ2=1+δB(M2-C2)<1.
(3.1)
(M,0) is a source if δ>2/B(C2-M2).(3.2)
(M,0) is a saddle if δ<2/B(C2-M2).(3.3)
(M,0) is nonhyperbolic, and there may be a flip bifurcation if δ=2/B(C2-M2).The Jacobian matrix of model (5) at the positive fixed point (C,L) is(9)J∗=1+δG-δC22δBCL1.So then, the characteristic polynomial of the Jacobian matrixJ∗ is(10)Fλ=λ2-2+Gδλ+1+Gδ+Hδ2,where(11)G=-AM+2A+AMC-3A+MC2+41+MC3-5C4-2CL,H=2BC3L.It is easy to see that F(1)=Hδ2>0 and F(-1)=Hδ2+2Gδ+4.Lemma 2.
The stability of the positive equilibrium point(C,L) is as follows [1]:(1)
(C,L) is a sink if δ∈(0,-G-G2-4H/H) and G<-2H (or 0<δ<-G/H and -2H<G<0).(2)
(C,L) is a source if one of the following conditions holds:
(2.1)
δ>-G/H and -2H≤G<0.(2.2)
δ>-G+G2-4H/H.(2.3)
G≥0.(3)
(C,L) is a saddle if -G-G2-4H/H<δ<-G+G2-4H/H and G<-2H.(4)
(C,L) is nonhyperbolic, and there may be a flip bifurcation, if δ=-G±G2-4H/H,G<-2H,δ≠-2/G,δ≠-4/G.(5)
(C,L) is nonhyperbolic, and there may be a N-S bifurcation if δ=-G/H,-2H<G<0.ForA=0.4762, B=0.0689, and C=0.9112, the conclusions of Lemmas 1 and 2 are shown in δ-M plane (Figure 1).Figure 1
The stability of equilibrium point inδ-M plane. (a) The stability of (0,0) and the two curves denote the flip bifurcation manifolds and divide the δ-M plane into four regions with different stability of (0,0). The intersections R2 denote the strong 1 : 2 resonance. (b) The stability of (1,0). (c) The stability of (M,0). (d) The stability of (C,L). (C,L) exists when M values are underneath M=0.9112 curve. N-S curve denotes the Neimark-Sacker manifold. The N-S curve and flip curve divide the δ-M plane into three regions with different stability of (C,L).
(a) (b) (c) (d)


## 3. Bifurcations Analysis

### 3.1. Flip Bifurcation

In this subsection,δ is chosen as a bifurcation parameter. In order to simplify calculation, we use the transformation u~(t)=u(t)-u∗,v~(t)=v(t)-v∗, then the equilibrium point (u∗,v∗) of model (5) is transformed into the origin spontaneously. The parameter perturbation δ=δ~+δ1 is used to write model (5)(12)u~t+1=u~t+δ~+δ11-u~t-x∗u~t+x∗-MA+u~t+x∗2-u~t+x∗v~t+y∗u~t+x∗,v~t+1=v~t+δ~+δ1Bu~t+x∗2-C2v~t+y∗.Model (12) is expanded to(13)u~t+1=a11u~t+a12v~t+a13u~2t+a14u~tv~t+b1u~tδ1+b2v~tδ1+a15u~3t+a16u~2tv~t+b3u~2tδ1+b4u~tv~tδ1+Ou~t+v~t+δ14,v~t+1=a21u~t+a22v~t+a23u~2t+a24u~tv~t+e1u~tδ1+a25u~2tv~t+e2u~2tδ1+e3u~tv~tδ1+Ou~t+v~t+δ14,where(14)a11=1+-5u∗4+4+4Mu∗3+-3M-3Au∗2+2+2MA-2v∗u∗-MAδ~,a12=-δ~u∗2,a13=-10u∗3+6M+6u∗2+-3M-3Au∗+1+MA-v∗δ~,a14=-2δ~u∗,b1=-5u∗4+4+4Mu∗3+-3M-3Au∗2+2+2MA-2v∗u∗-MA,b2=-u∗2,a15=-10u∗2+-4M-4u∗+A+Mδ~,a16=-δ~,b3=-10u∗3+6M+6u∗2+-3M-3Au∗+1+MA-v∗,b4=-2u∗,a21=2δ~Bu∗v∗L,a22=1+δ~Bu∗2-C2,a23=δ~Bu∗,a24=2δ~Bu∗,e1=2Bu∗v∗,a25=δ~B,e2=Bv∗,e3=2Bu∗.The transformation u~(t)=a12x~(t)+a12y~(t),v~(t)=-(1+a11)x~(t)+(λ2-a11)y~(t) changes (13) into the form(15)x~t+1y~t+1=-100λ2x~ty~t+fx~t,y~t,δ1gx~t,y~t,δ1,where(16)fx~t,y~t,δ1=λ2-a11a13a12λ2+1-a23λ2+1u~2t+λ2-a11a14a12λ2+1-a24λ2+1u~tv~t+λ2-a11b1a12λ2+1-e1λ2+1δ1u~t+λ2-a11b2a12λ2+1v~tδ1+λ2-a11a15a12λ2+1u~3t+λ2-a11a16a12λ2+1-a25λ2+1v~tu~2t+λ2-a11b3a12λ2+1-e2λ2+1δ1u~2t+λ2-a11b4a12λ2+1-e3λ2+1u~tv~tδ1,(17)gx~,y~,δ1=1+a11a13a12λ2+1+a23λ2+1u~2t+1+a11a14a12λ2+1+a24λ2+1u~tv~t+1+a11b1a12λ2+1+e1λ2+1u~tδ1+1+a11b2a12λ2+1v~tδ1+1+a11a15a12λ2+1u~3t+1+a11a16a12λ2+1+a25λ2+1u~2tv~t+1+a11b3a12λ2+1+e2λ2+1u~2tδ1+1+a11b4a12λ2+1+e3λ2+1u~tv~tδ1,and u~(t)=a12x~(t)+a12y~(t),v~(t)=-(1+a11)x~(t)+(λ2-a11)y~(t).Based on the center manifold theorem [18–20], for x~ and δ1 sufficiently small, we assume(18)y~t=hx~t,δ1=a1x~2t+a2x~tδ1+a3δ12+Ox~t+δ13,then we have(19)h-x~t+fx~t,hx~t,δ1,δ1-λ2hx~t,δ1-gx~t,hx~t,δ1,δ1=0.From the invariant property and the power series method [21], we have the expression(20)a1=1+a11a14-1-a11+a13a122+a24a12-1-a11+a23a1221-λ22,a2=-e1a122+1+a11b2-1-a11+b1a12a12λ2+12,a3=0.Restricted to the center manifold, the reduced model of (15) is(21)x~t+1=-x~t+h1x~2t+h2x~tδ1+h3x~2tδ1+h4x~tδ12+h5x~3t+Ox~t+δ14,where(22)h1=1λ2+1a11-λ2a14a11+1+a24a11+1a12-a11-λ2a13a12-a23a122,h2=1a12λ2+1λ2-a11b21+a11+a12b1λ2-a11-e1a122,h3=1λ2+12a2λ2-1a13a12-2a2a23a122+a11-λ21+a11b4+e31+a11a12+λ2-a11b1a1-e1a12a1+λ2-a11b3a12-e2a122,h4=1a12λ2+1λ2-a112b2a2+λ2-a11b1a12a2-e1a122a2,h5=1λ2+1λ2-a11a15a122+λ2-a112a14a1-a24λ2-a11a1a12+2λ2-a11a13a12a1-2a23a122a1+a11-λ21+a11a12a16+a251+a11a122+a11-λ21+a11a14a1+a241+a11a12a1.If we define(23)Fx~,δ1=-x~t+h1x~2t+h2x~tδ1+h3x~2tδ1+h4x~tδ12+h5x~3t+Ox~t+δ14,then model (21) will undergo a flip bifurcation, if the following nondegeneracy conditions are satisfied:(24)F.1∂2F∂x~∂δ10,0=h2≠0,F.216∂3F∂x~3+12∂2F∂x~220,0=h5+h12≠0.Theorem 3.
Assume that(F.1) and (F.2) hold. Model (5) will undergo a flip bifurcation if one of the following conditions holds [7]:(1)
δ~=2/AM (Lemma 1 (1.4)).(2)
δ~=2/(A+1)(1-M) (Lemma 1 (2.3)).(3)
δ~=2/B(C2-M2) (Lemma 1 (3.3)).(4)
δ~=-G±G2-4H/H,G<-2H,δ~≠-2/G,δ~≠-4/G (Lemma 2(4)). Moreover, if h5+h12>0(resp.,h5+h12<0), then the flip bifurcation is supercritical(resp.,subcritical)


### 3.2. N-S Bifurcation

According to Lemma2, model (5) may undergo a Neimark-Sacker bifurcation if the condition (5) of Lemma 2 holds. We choose (C,L) as the bifurcation point and δ as the bifurcation parameter. In order to simplify calculation, we use the transformation u~(t)=u(t)-C,v~(t)=v(t)-L,δ=δ^+δ2. Model (5) becomes(25)u~t+1=u~t+δ^+δ21-u~t-Cu~t+C-MA+u~t+C2-u~t+Cv~t+Lu~t+C,v~t+1=v~t+δ^+δ2Bu~t+C2-C2v~t+L.We have the expansion(26)u~t+1=a11u~t+a12v~t+a13u~2t+a14u~tv~t+a15u~3t+a16u~2tv~t+a17u~4t+a18u~5t,v~t+1=a21u~t+a22v~t+a23u~2t+a24u~tv~t+a25u~2tv~t.where a11,a12,a13,a14,a15,a16,a21,a22,a23,a24,a25 are given in model (13) by substituting δ~ for δ^+δ2, u∗ for C, and v∗ for L, a17=(δ^+δ2)(1+M-5C),a18=-δ^-δ2.Whenδ2=0, the eigenvalues of the linearized system of model (26) are a pair of complex conjugate numbers λ and λ- with modulus 1, λ=μ+ωi, where(27)μ=1+Gδ^+δ22,ω=iδ^+δ24H-G22.Then we have(28)δ^=-GH,-2H<G<0,(29)λ=1+δ^+δ2G+Hδ^+δ22,dλdδ2δ2=0=-G2>0.When δ2=0, we require λm≠1 and λ-m≠1(m=1,2,3,4), which leads to(30)G2≠2H,3H,4H.The transformation u~(t)=a12x~(t),v~(t)=(μ-a11)x~(t)-ωy~(t) takes model (26) into(31)x~t+1y~t+1=μ-ωωμx~ty~t+f1x~t,y~tg1x~t,y~t,where(32)f1x~t,y~t=c13x~2t+c14x~ty~t+c15x~3t+c16x~2ty~t+c17x~4t+c18x~5t,g1x~t,y~t=c23x~2t+c24x~ty~t+c25x~3t+c26x~2ty~t+c27x~4t+c28x~5t,c13=a13a12+a14μ-a14a11,c14=-a14ω,c15=a12a15a12+a16μ-a16a11,c16=-a16a12ω,c17=a17a123,c18=a18a124,c23=1ω-a23a122+-a24+a13μ-a11a12+a14μ-a112,c24=-μ-a11a14+a24a12,c25=a12ωμ-a11-a25a12+a15a12+a16μ-a16a11,c26=a25a122-μ-a11a12a16,c27=μ-a11a17a123ω,c28=μ-a11a18a124ω.By the introduction of a complex variable, model (31) can be transformed into the form(33)z~t+1=μ+iωzt+∑2≤i+j≤51i!j!gijzitz¯jt,where(34)g20=12c24+c13-ic14-c23,g11=12c13+ic23,g02=12c13-c24+ic23+c14,g30=34c15+c26-ic16-c25,g21=143c15+c26-ic16-3c25,g12=143c15-c26+ic16+3c25,g03=34c15-c26+ic16+c25,g40=g31=g22=g13=g04=32c17+ic27.g50=g41=g32=g23=g14=g05=154c18+ic28,The following nondegeneracy conditiona(0)≠0 is necessary for Neimark-Sacker bifurcation, where(35)a0=Reλ-k10=-Re1-2λλ-221-λg20g11-12g112-g0224+Reλ-g212δ2=0,witha0=aδ2=0,k10=k1δ2=0,k10=g20g11λ¯-3+2λ2λ2-λλ¯-1+g1121-λ¯+g0222λ2-λ¯+g212δ2=0.Theorem 4.
If condition (28),(30) hold and a(0)≠0, then model (3) will undergo a Neimark-Sacker bifurcation when δ2 varies in a small neighborhood of origin. Moreover, if a(0)<0(a(0)>0), then a stable (unstable) invariant closed curve will bifurcate from the equilibrium point (C,L) for δ2>0(δ2<0).


### 3.3. The Two Invariant Closed Curves

In this subsection, we assume that the nondegeneracy condition of Neimark-Sacker bifurcation does not hold, that is,a(0)=0, then model (5) will undergo Chenciner bifurcation at the equilibrium point (C,L). The bifurcation is characterized by two parameters, δ and M. In order to simplify calculation, we use the change of variables u~(t)=u(t)-C,v~(t)=v(t)-L,δ=δˇ+δ3,M=Mˇ+M1. Plugging into model (5) gives(36)u~t+1=u~t+δˇ+δ31-u~t-Cu~t+C-Mˇ-M1A+u~t+C2-u~t+Cv~t+L~u~t+C,v~t+1=v~t+δˇ+δ3Bu~t+C2-C2v~t+L~,where L~=1/C(1-C)(C-Mˇ-M1)(A+C2).The expanded form of model (36) is(37)u~t+1=a11u~t+a12v~t+a13u~2t+a14u~tv~t+a15u~3t+a16u~2tv~t+a17u~4t+a18u~5t,v~t+1=a21u~t+a22v~t+a23u~2t+a24u~tv~t+a25u~2tv~t.where a11,a12,a13,a14,a15,a16,a17,a18,a21,a22,a23,a24,a25 are given in model (13) by substituting δ~ for δˇ+δ3, M for Mˇ+M1, u∗ for C, and v∗ for L, a17=(δˇ+δ3)(1+Mˇ+M1-5C),a18=-δˇ-δ3.Whenδ3=0,M1=0, the eigenvalues of (0,0) are a pair of complex conjugate numbers λ and λ- with modulus 1.Letλ=μ+ωi,H~=2BC3L~, and(38)G~=-AMˇ+M1+2A+AMˇ+M1C-3A+Mˇ+M1C2+41+Mˇ+M1C3-5C4-2CL~,where(39)μ=1+G~δˇ+δ32,ω=iδˇ+δ38BC3L~-G~22.We have(40)δˇ=-G~H~,-2H~<G~<0,(41)λ=1+δˇ+δ3G~+2BC3L~δˇ+δ32;when δ3=0,M1=0, we require λm≠1 and λ-m≠1(m=1,2,3,4,5,6) which leads to(42)G~2≠H~,2H~,3H~,4H~,5±52H~. As seen in Section 3.2, model (37) can be written by using a complex variable as(43)z~t+1=μ+iωzt+∑2≤i+j≤51i!j!gijzitz¯jt,where gij(2≤i+j≤5) are given in (33) by substituting δ^+δ2 for δˇ+δ3 and M for M˘+M1.We setβ1=|λ|-1,β2=a(δ3,M1),α=(δ3,M1),β=(β1,β2), the first inequality of (44) is required that it has the regularity at α=0, and moreover, the second one is an additional nondegeneracy condition for the bifurcation [7] of two invariant closed curves,(44)det⁡∂β∂αα=0≠0,L20=Reλ-k20+Imλ-k1022≠0,where(45)k10=g20g11λ¯-3+2λ2λ2-λλ¯-1+g1121-λ¯+g0222λ2-λ¯+g212δ3=0,M1=0,k20=g202g124λ-λ¯-2λ¯24λ2-λ2λ¯2-1+g20g22λ¯+λ-14λ2-λλ¯-1+g20g02g303λ3-6λ2+3λ+1-λ¯12λ-1λ¯2-λλ2-1+⋯δ3=0,M1=0.We do not provide the complete formula fork2(0) due to its length (the interested reader can refer to [22]).Theorem 5.
If the conditions (40),(42),(44) hold, then two invariant closed curves bifurcate from the equilibrium point (C,L) of model (5). The outside invariant closed curve is stable and the inner one is unstable (the outside one is unstable; the inner one is stable) if L2(0)<0 (L2(0)>0).


## 3.1. Flip Bifurcation

In this subsection,δ is chosen as a bifurcation parameter. In order to simplify calculation, we use the transformation u~(t)=u(t)-u∗,v~(t)=v(t)-v∗, then the equilibrium point (u∗,v∗) of model (5) is transformed into the origin spontaneously. The parameter perturbation δ=δ~+δ1 is used to write model (5)(12)u~t+1=u~t+δ~+δ11-u~t-x∗u~t+x∗-MA+u~t+x∗2-u~t+x∗v~t+y∗u~t+x∗,v~t+1=v~t+δ~+δ1Bu~t+x∗2-C2v~t+y∗.Model (12) is expanded to(13)u~t+1=a11u~t+a12v~t+a13u~2t+a14u~tv~t+b1u~tδ1+b2v~tδ1+a15u~3t+a16u~2tv~t+b3u~2tδ1+b4u~tv~tδ1+Ou~t+v~t+δ14,v~t+1=a21u~t+a22v~t+a23u~2t+a24u~tv~t+e1u~tδ1+a25u~2tv~t+e2u~2tδ1+e3u~tv~tδ1+Ou~t+v~t+δ14,where(14)a11=1+-5u∗4+4+4Mu∗3+-3M-3Au∗2+2+2MA-2v∗u∗-MAδ~,a12=-δ~u∗2,a13=-10u∗3+6M+6u∗2+-3M-3Au∗+1+MA-v∗δ~,a14=-2δ~u∗,b1=-5u∗4+4+4Mu∗3+-3M-3Au∗2+2+2MA-2v∗u∗-MA,b2=-u∗2,a15=-10u∗2+-4M-4u∗+A+Mδ~,a16=-δ~,b3=-10u∗3+6M+6u∗2+-3M-3Au∗+1+MA-v∗,b4=-2u∗,a21=2δ~Bu∗v∗L,a22=1+δ~Bu∗2-C2,a23=δ~Bu∗,a24=2δ~Bu∗,e1=2Bu∗v∗,a25=δ~B,e2=Bv∗,e3=2Bu∗.The transformation u~(t)=a12x~(t)+a12y~(t),v~(t)=-(1+a11)x~(t)+(λ2-a11)y~(t) changes (13) into the form(15)x~t+1y~t+1=-100λ2x~ty~t+fx~t,y~t,δ1gx~t,y~t,δ1,where(16)fx~t,y~t,δ1=λ2-a11a13a12λ2+1-a23λ2+1u~2t+λ2-a11a14a12λ2+1-a24λ2+1u~tv~t+λ2-a11b1a12λ2+1-e1λ2+1δ1u~t+λ2-a11b2a12λ2+1v~tδ1+λ2-a11a15a12λ2+1u~3t+λ2-a11a16a12λ2+1-a25λ2+1v~tu~2t+λ2-a11b3a12λ2+1-e2λ2+1δ1u~2t+λ2-a11b4a12λ2+1-e3λ2+1u~tv~tδ1,(17)gx~,y~,δ1=1+a11a13a12λ2+1+a23λ2+1u~2t+1+a11a14a12λ2+1+a24λ2+1u~tv~t+1+a11b1a12λ2+1+e1λ2+1u~tδ1+1+a11b2a12λ2+1v~tδ1+1+a11a15a12λ2+1u~3t+1+a11a16a12λ2+1+a25λ2+1u~2tv~t+1+a11b3a12λ2+1+e2λ2+1u~2tδ1+1+a11b4a12λ2+1+e3λ2+1u~tv~tδ1,and u~(t)=a12x~(t)+a12y~(t),v~(t)=-(1+a11)x~(t)+(λ2-a11)y~(t).Based on the center manifold theorem [18–20], for x~ and δ1 sufficiently small, we assume(18)y~t=hx~t,δ1=a1x~2t+a2x~tδ1+a3δ12+Ox~t+δ13,then we have(19)h-x~t+fx~t,hx~t,δ1,δ1-λ2hx~t,δ1-gx~t,hx~t,δ1,δ1=0.From the invariant property and the power series method [21], we have the expression(20)a1=1+a11a14-1-a11+a13a122+a24a12-1-a11+a23a1221-λ22,a2=-e1a122+1+a11b2-1-a11+b1a12a12λ2+12,a3=0.Restricted to the center manifold, the reduced model of (15) is(21)x~t+1=-x~t+h1x~2t+h2x~tδ1+h3x~2tδ1+h4x~tδ12+h5x~3t+Ox~t+δ14,where(22)h1=1λ2+1a11-λ2a14a11+1+a24a11+1a12-a11-λ2a13a12-a23a122,h2=1a12λ2+1λ2-a11b21+a11+a12b1λ2-a11-e1a122,h3=1λ2+12a2λ2-1a13a12-2a2a23a122+a11-λ21+a11b4+e31+a11a12+λ2-a11b1a1-e1a12a1+λ2-a11b3a12-e2a122,h4=1a12λ2+1λ2-a112b2a2+λ2-a11b1a12a2-e1a122a2,h5=1λ2+1λ2-a11a15a122+λ2-a112a14a1-a24λ2-a11a1a12+2λ2-a11a13a12a1-2a23a122a1+a11-λ21+a11a12a16+a251+a11a122+a11-λ21+a11a14a1+a241+a11a12a1.If we define(23)Fx~,δ1=-x~t+h1x~2t+h2x~tδ1+h3x~2tδ1+h4x~tδ12+h5x~3t+Ox~t+δ14,then model (21) will undergo a flip bifurcation, if the following nondegeneracy conditions are satisfied:(24)F.1∂2F∂x~∂δ10,0=h2≠0,F.216∂3F∂x~3+12∂2F∂x~220,0=h5+h12≠0.Theorem 3.
Assume that(F.1) and (F.2) hold. Model (5) will undergo a flip bifurcation if one of the following conditions holds [7]:(1)
δ~=2/AM (Lemma 1 (1.4)).(2)
δ~=2/(A+1)(1-M) (Lemma 1 (2.3)).(3)
δ~=2/B(C2-M2) (Lemma 1 (3.3)).(4)
δ~=-G±G2-4H/H,G<-2H,δ~≠-2/G,δ~≠-4/G (Lemma 2(4)). Moreover, if h5+h12>0(resp.,h5+h12<0), then the flip bifurcation is supercritical(resp.,subcritical)


## 3.2. N-S Bifurcation

According to Lemma2, model (5) may undergo a Neimark-Sacker bifurcation if the condition (5) of Lemma 2 holds. We choose (C,L) as the bifurcation point and δ as the bifurcation parameter. In order to simplify calculation, we use the transformation u~(t)=u(t)-C,v~(t)=v(t)-L,δ=δ^+δ2. Model (5) becomes(25)u~t+1=u~t+δ^+δ21-u~t-Cu~t+C-MA+u~t+C2-u~t+Cv~t+Lu~t+C,v~t+1=v~t+δ^+δ2Bu~t+C2-C2v~t+L.We have the expansion(26)u~t+1=a11u~t+a12v~t+a13u~2t+a14u~tv~t+a15u~3t+a16u~2tv~t+a17u~4t+a18u~5t,v~t+1=a21u~t+a22v~t+a23u~2t+a24u~tv~t+a25u~2tv~t.where a11,a12,a13,a14,a15,a16,a21,a22,a23,a24,a25 are given in model (13) by substituting δ~ for δ^+δ2, u∗ for C, and v∗ for L, a17=(δ^+δ2)(1+M-5C),a18=-δ^-δ2.Whenδ2=0, the eigenvalues of the linearized system of model (26) are a pair of complex conjugate numbers λ and λ- with modulus 1, λ=μ+ωi, where(27)μ=1+Gδ^+δ22,ω=iδ^+δ24H-G22.Then we have(28)δ^=-GH,-2H<G<0,(29)λ=1+δ^+δ2G+Hδ^+δ22,dλdδ2δ2=0=-G2>0.When δ2=0, we require λm≠1 and λ-m≠1(m=1,2,3,4), which leads to(30)G2≠2H,3H,4H.The transformation u~(t)=a12x~(t),v~(t)=(μ-a11)x~(t)-ωy~(t) takes model (26) into(31)x~t+1y~t+1=μ-ωωμx~ty~t+f1x~t,y~tg1x~t,y~t,where(32)f1x~t,y~t=c13x~2t+c14x~ty~t+c15x~3t+c16x~2ty~t+c17x~4t+c18x~5t,g1x~t,y~t=c23x~2t+c24x~ty~t+c25x~3t+c26x~2ty~t+c27x~4t+c28x~5t,c13=a13a12+a14μ-a14a11,c14=-a14ω,c15=a12a15a12+a16μ-a16a11,c16=-a16a12ω,c17=a17a123,c18=a18a124,c23=1ω-a23a122+-a24+a13μ-a11a12+a14μ-a112,c24=-μ-a11a14+a24a12,c25=a12ωμ-a11-a25a12+a15a12+a16μ-a16a11,c26=a25a122-μ-a11a12a16,c27=μ-a11a17a123ω,c28=μ-a11a18a124ω.By the introduction of a complex variable, model (31) can be transformed into the form(33)z~t+1=μ+iωzt+∑2≤i+j≤51i!j!gijzitz¯jt,where(34)g20=12c24+c13-ic14-c23,g11=12c13+ic23,g02=12c13-c24+ic23+c14,g30=34c15+c26-ic16-c25,g21=143c15+c26-ic16-3c25,g12=143c15-c26+ic16+3c25,g03=34c15-c26+ic16+c25,g40=g31=g22=g13=g04=32c17+ic27.g50=g41=g32=g23=g14=g05=154c18+ic28,The following nondegeneracy conditiona(0)≠0 is necessary for Neimark-Sacker bifurcation, where(35)a0=Reλ-k10=-Re1-2λλ-221-λg20g11-12g112-g0224+Reλ-g212δ2=0,witha0=aδ2=0,k10=k1δ2=0,k10=g20g11λ¯-3+2λ2λ2-λλ¯-1+g1121-λ¯+g0222λ2-λ¯+g212δ2=0.Theorem 4.
If condition (28),(30) hold and a(0)≠0, then model (3) will undergo a Neimark-Sacker bifurcation when δ2 varies in a small neighborhood of origin. Moreover, if a(0)<0(a(0)>0), then a stable (unstable) invariant closed curve will bifurcate from the equilibrium point (C,L) for δ2>0(δ2<0).


## 3.3. The Two Invariant Closed Curves

In this subsection, we assume that the nondegeneracy condition of Neimark-Sacker bifurcation does not hold, that is,a(0)=0, then model (5) will undergo Chenciner bifurcation at the equilibrium point (C,L). The bifurcation is characterized by two parameters, δ and M. In order to simplify calculation, we use the change of variables u~(t)=u(t)-C,v~(t)=v(t)-L,δ=δˇ+δ3,M=Mˇ+M1. Plugging into model (5) gives(36)u~t+1=u~t+δˇ+δ31-u~t-Cu~t+C-Mˇ-M1A+u~t+C2-u~t+Cv~t+L~u~t+C,v~t+1=v~t+δˇ+δ3Bu~t+C2-C2v~t+L~,where L~=1/C(1-C)(C-Mˇ-M1)(A+C2).The expanded form of model (36) is(37)u~t+1=a11u~t+a12v~t+a13u~2t+a14u~tv~t+a15u~3t+a16u~2tv~t+a17u~4t+a18u~5t,v~t+1=a21u~t+a22v~t+a23u~2t+a24u~tv~t+a25u~2tv~t.where a11,a12,a13,a14,a15,a16,a17,a18,a21,a22,a23,a24,a25 are given in model (13) by substituting δ~ for δˇ+δ3, M for Mˇ+M1, u∗ for C, and v∗ for L, a17=(δˇ+δ3)(1+Mˇ+M1-5C),a18=-δˇ-δ3.Whenδ3=0,M1=0, the eigenvalues of (0,0) are a pair of complex conjugate numbers λ and λ- with modulus 1.Letλ=μ+ωi,H~=2BC3L~, and(38)G~=-AMˇ+M1+2A+AMˇ+M1C-3A+Mˇ+M1C2+41+Mˇ+M1C3-5C4-2CL~,where(39)μ=1+G~δˇ+δ32,ω=iδˇ+δ38BC3L~-G~22.We have(40)δˇ=-G~H~,-2H~<G~<0,(41)λ=1+δˇ+δ3G~+2BC3L~δˇ+δ32;when δ3=0,M1=0, we require λm≠1 and λ-m≠1(m=1,2,3,4,5,6) which leads to(42)G~2≠H~,2H~,3H~,4H~,5±52H~. As seen in Section 3.2, model (37) can be written by using a complex variable as(43)z~t+1=μ+iωzt+∑2≤i+j≤51i!j!gijzitz¯jt,where gij(2≤i+j≤5) are given in (33) by substituting δ^+δ2 for δˇ+δ3 and M for M˘+M1.We setβ1=|λ|-1,β2=a(δ3,M1),α=(δ3,M1),β=(β1,β2), the first inequality of (44) is required that it has the regularity at α=0, and moreover, the second one is an additional nondegeneracy condition for the bifurcation [7] of two invariant closed curves,(44)det⁡∂β∂αα=0≠0,L20=Reλ-k20+Imλ-k1022≠0,where(45)k10=g20g11λ¯-3+2λ2λ2-λλ¯-1+g1121-λ¯+g0222λ2-λ¯+g212δ3=0,M1=0,k20=g202g124λ-λ¯-2λ¯24λ2-λ2λ¯2-1+g20g22λ¯+λ-14λ2-λλ¯-1+g20g02g303λ3-6λ2+3λ+1-λ¯12λ-1λ¯2-λλ2-1+⋯δ3=0,M1=0.We do not provide the complete formula fork2(0) due to its length (the interested reader can refer to [22]).Theorem 5.
If the conditions (40),(42),(44) hold, then two invariant closed curves bifurcate from the equilibrium point (C,L) of model (5). The outside invariant closed curve is stable and the inner one is unstable (the outside one is unstable; the inner one is stable) if L2(0)<0 (L2(0)>0).


## 4. Numerical Simulations

In this section, the bifurcation diagrams and phase portraits of model (5) clearly demonstrate our theoretical results by numerical simulations. Parameters A,B,C are fixed at A=0.4762,B=0.0689,C=0.9112, respectively. The bifurcation parameters δ and M are considered in the following three cases.Case  a(δ=2.5698,M=0.1372). The equilibrium point of model (5) is (0.9112,0.0985), the conditions of (F.1), (F.2) are h2=-0.8050 and h5+h12=-259.2976, respectively. Based on Theorem 3, the flip bifurcation occurs when δ has a small perturbation in the interval (2.5698,2.78). The numerical simulation shows that there are orbits of periods 2, 4, 8, 16, and 32, chaotic behavior, and so on. Furthermore, Figures 2(a), 2(b), and 2(c) show the stable positive fixed point, period-2 orbit, and period-4 orbit, respectively, and Figure 2(d) exhibits the bifurcation diagram in (δ-u) plane.Figure 2
The initial value is(0.92,0.1), (a) δ=2.56, the middle line u=0.9112 denotes the equilibrium point and the equilibrium point is stable, (b) δ=2.5705, there is a period-2 orbit, (c) δ=2.7, there is a period-4 orbit, and (d) the bifurcation diagram in (δ-u) plane, showing the process of period doubling and chaos.
(a) (b) (c) (d)Caseb(δ=1,M=0.8189). The positive equilibrium point of model (5) is (0.9112,0.0117), and the eigenvalues of the equilibrium point are 0.9994±0.035i, a(0)=-1.125≠0. Based on Theorem 4, Neimark-Sacker bifurcation occurs when δ has a small perturbation. There is a unique and stable invariant closed curve bifurcating from the equilibrium point. Figure 3(a) exhibits the bifurcation diagram in (δ-u) plane. In addition, Figures 3(b) and 3(c) show the phase portraits of the model for different δ. Figure 3(d) shows the corresponding solution curves at δ=16.Figure 3
The initial value is(0.9112,0.01228), (a) the bifurcation diagram in δ-u plane, it exhibits the dynamical behaviors of model (5), (b) δ=1, the equilibrium point is stable, (c) δ=1.5, one invariant closed curve emerges from the equilibrium point, and (d) the solution curve when δ=1.5, which shows that the invariant closed curve is not orbit of model (5).
(a) (b) (c) (d)Casec(A=0.4762,B=0.0689,C=0.9112,M=0.81999,δ=9.3745×10-4). The positive equilibrium point of model (5) is (0.9112,0.0116), the eigenvalues of the equilibrium point are 0.9999±3×10-5i, and a=6.4315×10-21≈0; by using MatcontM [23], we get L2(0)=5.3263×10-5≠0. According to Theorem 5, two invariant closed curves will bifurcate from the equilibrium point. The bifurcation diagram in δ-M plane is shown in Figure 4(a), the green curve marked by N-S denotes the Neimark-Sacker bifurcation manifold [24], the blue curve denotes the sketchy border where two invariant closed curves collide and disappear, and similarly, C-H denotes the Chenciner bifurcation point. From Figure 4(b), we can see that when δ=0.3,M=0.821, the equilibrium point is unstable and there are no invariant closed curves bifurcating from the point. As M decreases, crossing the blue curve, at δ=0.3,M=0.81999, two invariant closed curves emerge from the equilibrium point. The larger closed curve is unstable and the smaller one is stable (Figure 4(c)). When δ=0.3,M=0.819, the smaller closed curve invariant closed curve enlarges gradually with the change of M and then breaks up; the larger one still retains unstable (Figure 4(d)).Figure 4
Parameters areA=0.4762,B=0.0689,C=0.9112, the initial values are (0.9112,0.01162), (0.9112,0.012247), (0.9112,0.012249),(0.95,0.011),(0.955,0.011), (a) the bifurcation diagram in (δ-M) plane [7], (b) δ=0.3,M=0.821, the equilibrium point is unstable, (c) δ=0.3,M=0.81999, the equilibrium point is unstable, two invariant closed curves emerge from the equilibrium point, the larger circle is unstable, and the smaller one is stable, and (d) δ=0.3,M=0.819, the equilibrium point is stable and the smaller circle disappears.
(a) (b) (c) (d)


## 5. Conclusion

We investigate the dynamic behaviors of the discrete model (5), obtained from the topologically equivalent model (3) of a differential equation system of Kolmogorov type (1). Initially, we analyze the local stability of the equilibrium point of model (5) and display the bifurcation diagrams. Furthermore, we present a one-parameter bifurcation analysis of the model (5) (for parameter δ) theoretically and show that model (5) can undergo flip bifurcation and Neimark-Sacker bifurcation. Moreover, we do a two-parameter bifurcation analysis of the model (5) (for parameter δ and M) theoretically and show that the model can undergo two invariant closed curves’ bifurcation for small perturbation of parameter δ and M. Finally, we carry out the numerical simulation and demonstrate our results. Compared with continuous model (3), if the step size δ value is small, the dynamical behavior of model (5) is similar to model (3). Otherwise, the dynamical behavior of model (5) is more complex than model (3), it is shown clearly in Figures 2(d) and 3(a).The novel result obtained in this paper is the existence of two invariant closed curves of model (5), it is a new result that other papers had not been mentioned, and it provided and extended the research direction on discrete models. Here we mention that we do not focus on chaos, because there are a lot of papers devoted to this topic. We not only verify the existence of two invariant closed curves theoretically and numerically, but also analyze and simulate the stability of the two curves. Although there are many examples [11, 25, 26] to elaborate the existence of multiple limit cycles for continuous models, the research on the existence of two or more invariant closed curves for discrete ecological models is few, even rare due to the complicated calculation. It is not only a tremendous challenge to verify these conditions to produce two invariant closed curves, but also not easy to find proper parameters to simulate the behavior; simultaneously, we need to keep the solution of the model positive. In our work, we provide a new idea to investigate the behavior, similar to continuous situation, by calculating normal form coefficients and finding the n-order focal values or Lyapunov coefficients [7, 19, 22, 27]; this is the first attempt to link limit cycles theory of continuous models with discrete models’ invariant closed curves bifurcations by normal form coefficients and Lyapunov coefficients. From above analysis, we can see that this method is feasible and intelligible.When the given conditions are satisfied, continuous models, as we know, can produce no more thann limit cycles if origin as the equilibrium point of the model is fine focus of order n. Analogously, discrete models can produce two or more invariant closed curves [22] if corresponding conditions are satisfied, but this is rarely studied in actual models; moreover, the invariant closed curves are not the periodic solutions of the discrete model generally. Therefore, there are many works to do.Future research will focus on the global stability of the equilibrium point and the existence of the discrete model with three or more invariant closed curves, and the conditions for the invariant closed curves are the periodic solutions of the discrete model.



---

*Source: 1613709-2018-05-30.xml*