# Equilibria and Bogdanov-Takens Bifurcation Analysis in the Bazykin’s Predator-Prey System

Authors: Shuangte Wang; Hengguo Yu
Journal: Discrete Dynamics in Nature and Society
 (2022)

---

## Abstract

In the paper, we proposed a Bazykin’s predator-prey system to explore the equilibrium point and Bogdanov-Takens bifurcation problems. Firstly, we derived some key parameter threshold conditions to ensure that the Bazykin’s predator-prey system had a multiple focus of multiplicity one, weak focus of order 2, cusps of codimension 2 and a degenerate Bogdanov-Takens singularity (focus or center case) of codimension 3. Furthermore, the distinction of two types of codimension 2 cusps was also discussed, which showed that the threshold of the two types of cusps could exhibit a cusp, which was a special case of the mentioned degenerate Bogdanov-Takens singularity (focus or center case) of codimension 3. Secondly, we systematically calculated that the Bazykin’s predator-prey system could undergo two types of Bogdanov-Takens bifurcations of codimension 2 and a degenerate focus type Bogdanov-Takens bifurcation of codimension 3. Finally, some numerical examples were implemented to verify the correctness and feasibility of mathematical theory derivation, which also directly showed all possible equilibrium points and Bogdanov-Takens bifurcations of Bazykin’s predator-prey system. In a word, all the research results could play an important theoretical support role in the study of controlling cyanobacteria bloom.


---

## Body

## 1. Introduction

At present and in the future, water eutrophication is still one of the major water environmental problems in the world, especially the phenomenon of cyanobacteria bloom [1]. Studying the formation mechanism of cyanobacteria bloom has important ecological and environmental significance for scientifically predicting the occurrence of cyanobacteria bloom in lakes and taking corresponding measures to reduce its impact [2]. In order to explore the formation mechanism of cyanobacteria blooms in eutrophic lakes, it is necessary to comprehensively understand the dynamic evolution process of cyanobacteria bloom and systematically analyze the synergistic influence mechanism of major factors (chemistry, physics and biology) on cyanobacteria bloom with the help of mathematical models, especially the Bazykin’s predator-prey system. This is because that Bazykin’s predator-prey system started from the Lotka-Volterra system, various regulating factors were considered, such as rates of birth and death, predation and competition, these different factors can have a stabilizing or a destabilizing effect on the community, and their interplay leads to increasingly complicated behaviors, which can describe the growth dynamic mechanism of some cyanobacteria populations [3].Bazykin [4] presented a variation of Volterra’s classical predator-prey model under the help of Michaelis-Menton equation, which can be called Bazykin’s predator-prey model according to Bazykin’s modification. Since then, the Bazykin’s predator-prey system has received the attention of a large number of scholars and obtained some excellent research results [5–10]. The paper [5] investigated bifurcations of equilibria in Bazykin’s predator-prey model. The paper [6] researched stability and bifurcation of Bazykin’s predator-prey model with memory effect, which can reveal the effect of memory based growth on global bifurcation threshold. The paper [7] considered Bazykin’s predator-prey model to address harvesting induced stability exchanges through bifurcation analysis, which can be useful to understand conservation policy and fishery management. The paper [8] proposed Bazykin predator-prey model with Beddington-DeAngelis response function to study persistence and global stability. The paper [9] inquired into the critical normal form coefficients and bifurcations of a discrete-time Bazykin’s predator-prey model. The paper [10] theoretically and numerically studied stability and bifurcation behaviors of the Bazykin’s predator-prey ecosystem with Holling type II functional response. It is obvious that the Bazykin’s predator-prey system can more realistically describe the dynamic relationship between predator population and prey population, and produce relatively rich bifurcation dynamics, such as transcritical bifurcation, saddle-node bifurcation, Hopf bifurcation and Bogdanov-Takens bifurcation, which can in fact be seen as bridging the gap between mathematical biology and bifurcation theory, Furthermore, the Bazykin’s predator-prey system can lead to sudden changes in ecosystem state and parameter space, which can analyze and predict qualitative changes in dynamics. Usually, these sudden changes are a very needed tool to study the outbreak of cyanobacteria bloom, because of the increasing adverse load on the cyanobacteria biosphere.In relation to Michaelis-Menton equation, we will consider a predator-prey ordinary differential equations (ODEs) system with prey-dependent Holling type II functional response and density-dependent death rate for the predator, which is a variation of Volterra’s classical predator-prey model and can be referred to Bazykin’s works.(1)x˙=r1x1−xK1−αxya+x−m1x,y˙=αexya+x−m2y−dy2,here the parameters r1, K1, α, a, m1, e, m2 and d are all positive constants with practically biological meanings, respectively. r1 denotes the intrinsic growth rate of the prey, K1 represents the carrying capacity of the environment, a is the half-saturation constant, α is the search efficiency of predator for prey, m1 and m2 are mortality rate of the prey and predator species respectively, e is the biomass conversion, d is the intra-specific competition coefficient. The functions x=xt and y=yt represent densities of the prey and predator at time t, respectively. The term dy2 represents interspecific density-restricted effect on the predator; the term αx/a+x (Holling type II functional response) is named after C. S. Holling, who proposed several functional responses for different kinds of species to model the phenomena of predation in 1965, since then the classical Lotka-Volterra predator-prey system in biomathematics was more advanced and realistic [11–14]. These functional responses described how predators transform harvested prey into the growth of themselves and were discussed by many researchers [15], especially the Rosenzweig-MacArthur model (R-M model) or a predator-prey model with Holling type II functional response [16].To investigate complex dynamical behaviors of the system (1) more precisely, in this paper, we concentrate on Bogdanov-Takens (BT) bifurcations at corresponding interior equilibria, respectively. It is clear to see that all solutions of the system (1) are positive and bounded with respect to the positive initial conditions x0>0, y0>0. Also, the x-axis and y-axis are invariant, namely, the system (1) is dissipative in the first quadrant ℝ+2 and well-defined on the domain ℝ+2=ℝ+2¯. The rest of this paper is organized as followed. In Section 2, we mainly consider stability of interior (positive) equilibria with respect to following special cases: (i) A1=0, A2>0; (ii) A1=A2=0; (iii) A2=0, A1≠0. The first case (i) ensures the potential Hopf bifurcation with the help of standard bifurcation theory, while the second case (ii) exhibits BT bifurcations of codimension 2 with two bifurcation parameters and a degenerate focus type BT bifurcation of codimension 3 with three bifurcation parameters. The statements of bifurcations are illustrated in Section 3. Finally, we also give some numerical simulations and discussions to support our theories. The last Section 4 is a short summary.


## 2. Equilibria Analysis

In this section, we mainly discuss the existence conditions and stability analysis of interior equilibria in some special cases. It is clear that the system (1) has bound equilibria E0:=0,0, E1:=0,−m2/d and E2:=K11−m1/r1,0. For biological consideration, we omit the point E1.


### 2.1. Preliminaries

In this subsection, we fully discuss preliminaries of the interior equilibrium point in the system (1). Here we denote it as E∗=x∗,y∗ for later use. An equilibrium point E∗ must satisfy following algebraic polynomial form equations(2)fx,y=r11−xK1a+x−αy−m1a+x=0,gx,y=αex−m2a+x−dya+x=0.In other words,x∗ and y∗ must be positive roots of following third-order polynomials (cubic equations) px=∑i=03aixi and py=∑i=03biyi, respectively, where coefficients are(3)a3=dr1,a2=dm1−r1K1+2adr1,a1=2am1−r1d+ααe−m2K1+a2dr1,a0=am1−r1d−m2αaK1,b3=K1d2,b2=−2K1dαe−m2,b1=αe−m22+ader1−m1K1+a2der1,b0=em1−r1αe−m2K1+am2r1a.We will mention an eliminant method in linear algebraic merely for completeness and shall not involve it in detail. If we sort above equations (2) as the form of fx,y=a0xy+a1x=0, gx,y=b0xy+b1x=0, the cubic equation px=0 can also be obtained again from the first Sylvester’s resultant Ryf,g=a0xb1x−a1xb0x. Similarly, if we sort the equations (2) as the form of fx,y=c0yx2+c1yx+c2y=0, gx,y=d0yx+d1y=0, the second resultant(4)Rxf,g=c0yc1yc2yd0yd1y00d0yd1y.Also deduces the cubic equationqy=0.Thanks to Scipione del Ferro, Niccolo Fontana(Tartaglia), Gerolamo Cardano and Shengjin Fan’s groundbreaking masterpieces, a third-order algebraic equation’s root(s) can be formulated by Cardano’s formula or Shengjin’s formula, and may qualitatively have one, two or three positive roots in the light of Rene Descartes’ rule of signs or Shengjin’s discriminant. Here we let these complicated expressions(5)px=a1a3−13a2a32,qx=227a2a33−13a2a3a1a3+a0a3,Δx=qx22+px33,and(6)py=b1b3−13b2b32,qy=227b2b33−13b2b3b1b3+b0b3,Δy=qy22+py33,be discriminants of above cubic equations px=0 and qy=0, respectively [5].The interior equilibrium pointE∗ does not exist when r1≤m1 or αe≤m2, thus we always assume r1>m1 and αe>m2 in the rest of this paper. The isoclines from the equations (2) and cobweb model show that an interior equilibrium point E∗ exists if condition(7)0<m2aαe−m2<K11−m1r1,holds. Noticing that the condition (7) deduces a0<0, b0<0, according to the zero theorem and(8)pm2aαe−m2=−αex2+m2a+x2e2r1a2α2dαe−m23<0,qαe−m2d=a2αe2r1>0.It is supposed that such condition (7) is well-done. In addition, the condition (7) and one of following conditions(9)ix2−a2≤m2aαe−m2,iia1,a2≥0,iiia22−3a1a3≤0,ivb22−3b1b3≤0,can lead to the uniqueness of the equilibrium point E∗, here x2=K11−m1/r1.


### 2.2. Multiple Focus with Multiplicity One

Due to the Routh-Hurwitz criterion and the Perron’s theorems, we merely need to consider the Jacobian matrix at a non-hyperbolic interior equilibrium pointE∗:(10)JE∗=Jij2×2=αx∗y∗a+x∗2−r1x∗K1−αx∗a+x∗αeay∗a+x∗2−dy∗.Furthermore, the trace, determinant and discriminant of matrixJE∗ are denoted as(11)A1:=A1E∗=trJE∗=−2r1x∗2+−αe−m1+m2+r1K1−r1ax∗+K1am2K1a+x∗,A2:=A2E∗=detJE∗ and Δ∗:=Δ∗E∗=A12−4A2, respectively. In order to obtain potential Hopf bifurcation, i.e. A1=0 and A2>0 (matrix JE∗ has a pair of pure imaginary eigenvalues), we take the threshold of parameter m1=r1−αe−m2, which can ensure that the numerator of A1 has a perfect square discriminant with respect of variable x∗. Thus the system (1) has an equilibrium point E3:=x3,y3=aλμ,λμ+1ea with parameters(12)m2=λαe,K1=μar1αe,d=−αλ−1μ+1λaλμ+12,and control variables λ∈0,1, μ>μm=1/1−λ. This special case is now denoted as (C1) for later use.Let us firstly discuss the equilibriumE3 more qualitatively, under the assumption of this case, the determinant is simplified as(13)A2E3=−λα2e2φA2λ,μλμ+12,φA2λ,μ=λ3μ2+−2μ2+2μλ2+μ−12λ−μ.The equationφA2λ,μ=0 has a quite clear positive root.(14)μ1=−2λ2+2λ+1+s2λ1−λ2>μm,s=−4λ2+4λ+1.That is to say: (i)A2E3>0 and E3 is a center or focus when μ∈μm,μ1; (ii) A2E3=0 and E3 is a potential cusp of codimension at least 2 when μ=μ1; (iii) A2E3<0 and E3 is just a saddle point when μ>μ1.This subsection will deal with above case (i). By using transformationx=u+x3,y=v+y3, the point E3 is translated to the origin O=0,0. Making a change of variables as u=−dy3X+βY, v=−αeay3/a+x32X, in which β=A2E3>0, then the system (1) becomes(15)X˙=F2X,Y=−βY+∑i+j=23aijXiYj+OX,Y4,Y˙=G2X,Y=βX+∑i+j=23bijXiYj+OX,Y4.Following [17, 18], the first Lyapunov number of the system (15) at the point E3, which is used to determine the stability of limit cycles as well, is presented as.(16)σ=3π2β3a30+b03+a12+b21−1β2a20b20−a02b02−a11a02+a20+b11b02+b20=3πα3e3λμφσλ,μ2a2λμ+14βφA2λ,μ,where(17)φσλ,μ=2+2λ4μ2+4μ−2μ2λ3+2−μλ2+μ+1λ.The transformationu=βX+dy3Y, v=αeay3/a+x32Y can also be used to deduce σ in (16).For any arbitraryλ, in view of a quadratic function φσλ,μ with respect to variable μ, it is a downward opening parabola with a symmetry axis μ=4λ2−λ+1/4λ21−λ>0 and has a positive maximum λ2+14λ+1/8λ1−λ. The equation φσλ,μ=0 straightly has a unique positive root(18)μσ=4λ2+λ2+14λ+1−λ+14λ21−λ>μm.What is more,μ1>μσ (or μσ exists) if λ>2/5; μ1<μσ if λ<2/5; μ1=μσ only if λ=2/5. Therefore, the equilibrium point E3 is respectively a stable(unstable if μσ exists) multiple focus with multiplicity one if μ<μσ>rbinμσ. When μ=μσ (if μσ exists), σ vanishes and one need to calculate the second Lyapunov coefficient l2. The generalized Hopf bifurcation of codimension 2 (Bautin bifurcation) may occur and two homocentric limit cycles may appear in some particular cases around the Bautin point E3 if l2≠0.


### 2.3. Weak Focus of Order 2

Whenμ=μσ, the first Lyapunov coefficient will be zero and we need to treat a weak (fine) focus E3 of order at least 2 in the so-called center-focus problem. This subsection starts with the system (15) or a standard form(19)x˙=−βy+φx,y,y˙=βx+ψx,y,where φx,y, ψx,y are both analytical functions and φ0,0=ψ0,0=0. Taking the polar coordinate transformation x=rcosθ, y=rsinθ, we derive a new system r˙=rRr,θ, θ.=β+Qr,θ, where Rr,θ=φcosθ+ψsinθ, Qr,θ=1/rψcosθ−φsinθ are also analytical functions. Now we rewrite the above system as(20)drdθ=rRβ+Q=∑k=2∞Rkθrk,here Rkθ are just polynomials of triangle functions cosθ and sinθ.We suppose that the special solution of the system (20) with initial condition rθ=0=c is rθ,c=∑k=1∞rkθck. Substituting it into the equation (20) and comparing all terms in ck, we have a series of coupled first-order differential equations(21)dr1dθ=0,dr2dθ=R2r12,dr3dθ=2R2r1r2+R3r13,⋯.These equations naturally yield solutionsr1=1, rk=gkθ+φkθk>1, where involved constants are(22)g2=12π∫02πR2τdτ,g3=12π∫02πR3τ+2R2τr2τdτ,and φkθ are periodic functions with a period 2π as well. For the system (19) with standard form, we obviously derive the first focal quantity g3 with g2=0, which repeats the formula (16) again, i.e. g3/σ≡constant.>0.Ifg3=0 or σ=0, from the integrals (22), we know g4=0 and the second focal quantity g5 (see the Appendix). From now on, with the system (15) at hand, the second focal quantity g5 reads(23)g5=−2097152λ8e4α42λ−16φg5λ3a4s1+3λ+1124λ2−λ+s1+141−λψg5λ7/2,in which auxiliary functions are ψg5λ=9λ3−λ2s1+3λ2+4λs1−3λ−s1−1 and(24)φg5λ=1+s1+54s1+61λ+1179s1+1533λ2+13252s1+20377λ3+79878s1+151954λ4+230064s1+604458λ5+57696s1+873264λ6+−1682544s1−2369616λ7+−6203262s1−14628750λ8+−11558132s1−36469478λ9+−6708162s1−43299230λ10+22132104s1+19510890λ11+88954156s1+206053636λ12+179705232s1+536172580λ13+268668984s1+899927400λ14+310805040s1+1188987000λ15+295112037s1+1226575797λ16+224630382s1+1060506033λ17+141939919s1+727205209λ18+69055348s1+415540245λ19+27134470s1+180678370λ20+6942224s1+61201034λ21+1283776s1+13517744λ22+19840s1+1703872λ23+−800s1−14240λ24+800λ25.Symbols1 is λ2+14λ+1. The resultant(25)Rs1φg5=10485765λ−22λ17λ−11412λ3+2λ2+2λ−1λ2+2λ+322λ2+λ+24≠0.Implies thatφg5λ>0 when λ∈2/5,1. While ψg5λ is also positive since Rs1ψg5=16λ2λ−15λ−2λ2+2λ+3≠0, then g5<0≠0. In other words, the equilibrium point E3 with σ=0 is a stable weak focus of order 2.


### 2.4. Cusp of Codimension 2

From the case (C1), undaunted by the scale of the distinction, this subsection reasonably set about finishing a puzzle when μ=μ1. Letting A2=0 in this case, we have threshold of parameters, a nilpotent E∗2=x∗2,y∗2 and an elementary equilibrium E∗1=x∗1,y∗1 as follows:(26)m1=r1−λ+1αe,d=λαμ1−λμ1−1aλμ1+12,x∗2=μ1aλ,y∗2=eaλμ1+1,x∗1=−λ2−λ+1s−λ2+λ+1aλs+1λ−1,y∗1=2aesλ+2λ2−2s−λ−2λ−1−2λ2+2λ+1+s.Firstly, we can make a transformation (I): x=X+x∗2, y=Y+y∗2 for the system (1) with above threshold of parameters:(27)X˙=F1X,Y,Y˙=G1X,Y.Secondly, we take the transformation(28)II:X=14s−2λ+1αeu+v,Y=αe2s−2λ+1s−1u8λ.Thus above system becomes Jordan’s standard form(29)u˙=F2u,v=v+a20u2+a11uv+a02v2+Ou,v3,v˙=G2u,v=b20u2+b11uv+b02v2+Ou,v3.By using the Lemma1 (or the Lemma 1 in [19]), this system (29) is equivalent to following system(30)x˙=y,y˙=d1λx2+d2λxy+Ox,y3,where discriminants manipulated by λ are(31)d1=d1λ=b20=α3e3λ−14φd1λ4a2λs+4λ−3s−5,d2=d2λ=b11+2a20=α2e2λ1−λ3s3+5−2λs2+19−12λs+23−18λas+1ψd2λ.Here temporarily auxiliary functions areφd1λ=2λs−2s2+8λ−s+1 and ψd2λ=s3+2λ−1s2+12λ−15s+14λ−17. We have added control variable s∈1,2.Lemma 1.
The following system(32)x˙=y+Ax2+Bxy+Cy2+Ox,y3,y˙=Dx2+Exy+Fy2+Ox,y3,is equivalent to the system(33)x˙=y,y˙=Dx2+E+2Axy+Ox,y3.After some nonsingular transformations in the neighbourhood of the originO.Letting functionφd1λ be zero and solving out the threshold of λ from the resultant (determinant)(34)Rsφd1=11−2λ2−1−8λ20011−2λ2−1−8λ210hλ0010hλ=4λ35λ−2,or equation φd1λ=0, it leads to a unique solution λ=2/5, where hλ=4λ2−4λ−1. And the derivative(35)dφd1λdλ=2ss2−1+8sλ+4λ1−λ>0.Implies that the functionφd1λ is a monotonic increasing function defined on interval 0,1. Combining the denominator of d1:(36)4a2λs+4λ−3s−5=4a2s+4λ−1−s+1<0.Here we know thatd1<>0 if λ><2/5. Notice that the numerator of d2 is obviously positive since the range of λ. In fact, we can assert that d2 is non zero. Indeed, for the function ψd2λ in the denominator of d2, its negativity will be found and well verified from an inequality(37)ψd2λ=s2s+λ−3s+1+λ−1s2+12s+14<s2−3s+1.That is to say,ψd2λ is always negative, namely d2 is negative here. All in all, d1d2≠0 when m2≠2/5αe, while d1d2=0 when m2=2/5αe.Next, we mainly concentrate on the stability of the equilibrium pointE∗1 whose trace, determinant and discriminant are parameterized as(38)A1E∗1=αe5−sφA1λ2s+12λs+4λ−3s−5,A2E∗1=αeαe−m2s−3φA2λ2s+12λs+4λ−3s−5,Δ∗E∗1=−m22φΔλs+14ψΔλ.In which some continuous auxiliary functions can be respectively expressed by(39)φA1λ=−2s3+2λ−7s2+−6λ+6s−16λ+19,φA2λ=s3+14λ−19s2+12λ−9s−34λ+43,φΔλ=∑j=19xjs9−j,ψΔλ=s4+4s3+8λ−6s2+28λ−34s+24λ−29,with(40)x1=1,x2=64λ−226,x3=124λ−354,x4=−1532λ+3156,x5=−2452λ+5010,x6=8248λ−11682,x7=13172λ−18786,x8=−10108λ+12208,x9=−16476λ+19889.The determinants(41)RsφA1=1−λ+723λ−38λ−192001−λ+723λ−38λ−19210hλ00010hλ00010hλ=165λ−2λ−15,RsφA2=114λ−1912λ−9−34λ+4300114λ−1912λ−9−34λ+4310hλ00010hλ00010hλ=1285λ−22λ−14,and(42)RsφΔ=1x2x3x4x5x6x7x8x9001x2x3x4x5x6x7x8x910hλ0000000010hλ0000000010hλ0000000010hλ0000000010hλ0000000010hλ0000000010hλ0000000010hλ=6553641λ2−54λ+81λ−1125λ−22,which are respectively from equations φA1λ=0, φA2λ=0 and φΔλ=0, yield the unique solution (zero point) λ=2/5 once more.For the sake of mathematical natures of the functionsφA1λ, φA2λ, φΔλ and ψΔλ defined on interval 0,1, we calculate their derivatives up to second order with respect to λ:(43)dφA1λdλ=φDφA1λs,d2φA1λdλ2=24s4+−48λ+72s3+−96λ2+168λ−72s2+96λ−1λ−1/22s3,dφA2λdλ=φDφA2λs,d2φA2λdλ2=−12s4+−336λ+264s3+48λ2−192λ+96s2−192λ3+336λ2−192λ+36s3,dφΔλdλ=φDφΔλs,d2φΔλdλ2=1s3−32s9+−5376λ+8120s8+768λ2−9696λ+11664s7+35840λ3−162400λ2+227440λ−125400s6+47616λ3−183552λ2+265536λ−153376s5+−367680λ3+1125120λ2−1146288λ+428520s4+−313856λ3+955136λ2−1035872λ+415984s3+395904λ3−956640λ2+781008λ−229448s2+161728λ−122λ−436361,dψΔλdλ=φDψΔλs,where new auxiliary continuous functions are(44)φDφA1λ=2s3+24λ−18s2+−16λ2+64λ−44s+24λ2−36λ+12,φDφA2λ=14s3+−12λ+18s2+−112λ2+208λ−110s−48λ2+60λ−18,φDφΔλ=64∑j=19yjs9−j,φDψΔλ=−16λ+16s3+−48λ+52s2+−64λ2+80λs−112λ2+192λ−68,in which(45)y1=1,y2=−32λ+14064,y3=−1792λ2+7224λ−469664,y4=−2976λ2+9984λ−670064,y5=30640λ2−78440λ+3980864,y6=39232λ2−99776λ+5325264,y7=−98976λ2+189672λ−8020064,y8=−105376λ2+202976λ−9162064,y9=40432λ2−69048λ+2441664.The following determinants(46)RsφDφA1=112λ−9−8λ2+32λ−2212λ2−18λ+600112λ−9−8λ2+32λ−2212λ2−18λ+610hλ00010hλ00010hλ=144λ−1420λ2−12λ−3,RsφDφA2=1−6λ+97−8λ2+104λ−557−24λ2+30λ−97001−6λ+97−8λ2+104λ−557−24λ2+30λ−9710hλ00010hλ00010hλ=1152495λ−25λ+1λ−14,RsφDφΔ=1y2y3y4y5y6y7y8y9001y2y3y4y5y6y7y8y910hλ0000000010hλ0000000010hλ0000000010hλ0000000010hλ0000000010hλ0000000010hλ0000000010hλ=1286560λ4−12808λ3+14305λ2−4395λ−15205λ−2λ−111,and(47)RsφDψΔ=1−48λ+52−16λ+16−64λ2+80λ−16λ+16−112λ2+192λ−68−16λ+16001−48λ+52−16λ+16−64λ2+80λ−16λ+16−112λ2+192λ−68−16λ+1610hλ00010hλ00010hλ=64λλ−15,which are respectively constructed from equations φDφA1λ=0, φDφA2λ=0, φDφΔλ=0 and φDψΔλ=0, yield that:(i)
The functionφA1λ has a (local) minimum φA13+26/10=227−1476/25 on interval 0,1 since d2/dλ2φA13+26/10=4386−1008>0, and, the detailed numerical simulation results are shown in Figure 1(a);(ii)
The functionφA2λ has a (local) minimum φA22/5=0 on interval 0,1 since d2/dλ2φA22/5=6480/49>0, and, the detailed numerical simulation results are shown in Figure 1(b);(iii)
The functionφΔλ has a local minimum φΔ2/5=0 and a local maximum φΔλ3≈534.507461 on interval 0,1 since d2/dλ2φΔ2/5=1156415616/30625>0 and d2/dλ2φΔλ3≈−13531.790941<0, where λ3 is another positive root of equation φDφΔλ=0 or the unique positive root of equation 6560λ4−12808λ3+14305λ2−4395λ−1520=0 in open interval 0,1, i.e.(48)λ3=19840M1/6N1/46−205NM2/3+−4040497M1/3+17198445275N+69931867473M1/2+4803M1/6N1/4−3N3/4≈0.764691,with(49)M=2783067866245+1917582022669769211,N=410M2/3−4040497M1/3−34396890550,and then, the detailed numerical simulation results are shown in Figure1(c); (iv) The function ψΔλ is a monotonic function on interval 0,1 since dψΔλ/dλ≠0, and, the detailed numerical simulation results are shown in Figure 1(d).Hence, maxima and minima of the functionsφA1λ, φA2λ, φΔλ and ψΔλ on unit interval 0,1 are listed as follows:(50)maxφA1λ=16,minφA1λ=227−147625<0,maxφA2λ=16,minφA2λ=0,maxφΔλ=9216,minφΔλ=0,maxψΔλ=−4,minψΔλ=−64.Combining Figures1(a)–1(d) and the determinant(51)RsψΔ=148λ−628λ−3424λ−2900148λ−628λ−3424λ−2910hλ000010hλ000010hhλλ000010hλ=256λ−18.Thus we have a rough summary:(i)
φA1λ is positive (negative) when λ<2/5λ>2/5, and it has a unique zero point 2/5;(ii)
φA2λ and φΔλ are nonnegative functions defined on interval 0,1, and they both have the unique zero point λ=2/5;(iii)
ψΔλ is a negative function defined on interval 0,1.Thus, it follows that: (i)A1E∗1 is negative (positive) when λ<2/5λ>2/5; (ii) A2E∗1 and Δ∗E∗1 are always positive when λ≠2/5. In a word, the equilibrium point E∗1 is a stable or unstable elementary node when λ≠2/5.On the whole, whenm2≠2/5αe, we have d1d2≠0, an interior equilibrium point E∗2 with A1=A2=0 from the case (C1) is a cusp of codimension 2, and an interior equilibrium point E∗1 is a stable (unstable) node when m2<2/5αem2>2/5αe. When m2=2/5αe, the unique interior equilibrium point E∗ becomes a cusp of codimension at least 3.Theorem 1.
From the case (C1), when m2≠2/5αe, an interior equilibrium point E∗2 with A1=A2=0 is a cusp of codimension 2 (Bogdanov-Takens bifurcation point) since d1d2≠0, and an interior equilibrium point E∗1 is an asymptotic stable (unstable) node when m2<2/5αem2>2/5αe. When m2=2/5αe, the unique interior equilibrium point E∗=4a,5ea is a cusp of codimension at least 3, indeed it is a codimension 3 BT singularity (focus or center case), see next subsection.Remark 1.
When we use the transformations (III):u=p+a02pq,v=q−a20p2, and (IV): p=w, q=z−c11wz, the system (29) firstly becomes(52)p˙=q+c11pq+Op,q3,q˙=d20p2+d11pq+d02q2+Op,q3,and then becomes(53)w˙=z+Ow,z3,z˙=f20w2+f11wz+f02z2+Ow,z3.
Finally, in order to obtain standard form, we construct a transformation (V):w=x1+1/2f02x12, z=f02x1y1+y1, above system becomes(54)x1.=y1+Ox1,y13,y1.=d1x12+d2x1y1+Ox1,y13.
This also ensures the Theorem1 and (31).Figure 1
(a) Figure of functionφA1λ on the interval 0,1; (b) figure of function φA2λ on the interval 0,1; (c) figure of function φΔλ on the interval 0,1; (d) figure of function ψΔλ on the interval 0,1.
(a)(b)(c)(d)


### 2.5. Bogdanov-Takens Singularity (Focus or Center) of Codimension 3

Lettingpx=Δx=0 or py=Δy=0, we obviously derive thresholds of parameters(55)K1=ar18αe+m2r1−m1αe−m2,d=8αe+m2αe−m2227αar1−m1e2,and a unique degenerate equilibrium point E6≔x6,y6=a2αe+m2/αe−m2,18r1−m1e2αa/8αe+m2αe−m2. In this case, it is quite obvious to see that A2E6=0, while the trace reduces to(56)A1E6=−16α2e2+6αem1−14αem2−6αer1+3m1m2−2m22−3m2r138αe+m2.The combinations in this case withA1E6≠0 is denoted as (C4). Furthermore, if A1E6=0 in the case (C4), thresholds of parameters r1, K1, d and a unique nilpotent equilibrium point E7 are(57)r1=2αe−m28αe+m2+3m12αe+m232αe+m2,K1=a2αe−m28αe+m2+3m12αe+m22αe−m22,d=2αe+m2αe−m218e2αa,E7≔x7,y7=a2αe+m2αe−m2,12e2αa2αe+m2.Similarly, we denote this special case as (C5). Note that the conditions (9) (iii) and (9) (iv) all hold in cases (C4) or (C5) since a22−3a1a3=b22−3b1b3=0.For the case (C4), we take transformations (I): x=X+x6, y=Y+y6 and(58)II:X=αx6ua+x6+v,Y=dy6u+3r1−m1e8αe+m2v,τ=tA1.The system (1) becomes a normal form(59)u˙=u+Φu,v,(60)v˙=Ψu,v,where Φu,v,Ψu,v=Ou,v2. The implicit function u=hv=h2v2+⋯ from right hand side of the equation (59) implies ψv=Ψhv,v=a3v3+⋯ with(61)a3=2r1−m1αe−m232αe+m2278αe+m2e2α2a2A12≠0.There is only one zero eigenvalue and a nonzero eigenvalue for the matrixJE6, by using the Theorem 7.1 in Zhifen Zhang’s book [20], it is easy to see that the equilibrium point E6 is a stable (unstable) node if A1<0A1>0.We will provide some explicitly smooth transformations to obtain a normal form with terms up to fourth order and determine the exact type of the equilibrium pointE7. Firstly, we take a transformation (I): x=X+x7, y=Y+y7 and translate the equilibrium point E7 to the origin O. Secondly, we take a linear transformation(62)II:u=−27αe−m22αe+m2Xaα8αe+m22−9αe−m22αe+m2Yaαe8αe+m22,v=−6αe−m22Xαa8αe+m2+3αe−m22αe+m2Yeaα8αe+m2,and transform matrix JE7 to its real Jordan’s canonical form, then above system becomes(63)u˙=F2u,v=v+∑i+j=24aijuivj+Ou,v5,v˙=G2u,v=∑i+j=24bijuivj+Ou,v5,where a11=b20=0 and a20a02b11b02≠0.Coming to the equation (63), we denote it in obvious notations:(64)u˙=v+Φu,v,(65)v˙=Ψu,v,where Φu,v,Ψu,v=Ou,v2. Substituting an implicit function v=hu=α2u2+α3u3+⋯ from right hand side of the first equation (64) into the second equation (65), where(66)α2=−28αe+m2αe−m281e,α3=4αe−m228αe+m22187e2.We immediately derive seriesψu=Ψu,hu=a3u3+⋯, Φuu,v+Ψvu,v∣v=hu=b1u+⋯ and a discriminant μ=b12+8a3, where(67)a3=b30−a20b11=−4αe−m228αe+m226561e2<0,b1=αe−m28αe+7m28αe+m281e2αe+m2≠0,μ=−αe−m228αe+m2264α2e2+16αem2−17m226561e22αe+m22<0.Thus the equilibrium pointE7 is a center or a focus.Thirdly, we take near identity transformations(68)III:u=p,v=q−a20p2−a02q2,IV:p=w+12d02w2,q=d02wz+z.To eliminate second order terms inu˙ and q2 term in q˙, thus above system firstly becomes(69)p˙=F3p,q=q+Op,q3,q˙=G3p,q=d11pq+d02q2+Op,q3,and finally becomes(70)w˙=F4w,z=z+∑i+j=3eijwizj+Ow,z4,z˙=G4w,z=f11wz+∑i+j=3fijwizj+Ow,z4.Taking into accountf11=b11+2a20 and the Lemma 1, the equilibrium point E7 is a cusp of codimension at least 3. Obviously, we notice that coefficients f11≠0, f30<0 and(71)5f30f21+3e30−3f11f40−e30f11=48αe+m2456α2e2−4αem2−7m22αe−m23143489072αe+m22e4≠0,f112+8f30=−8αe+m2264α2e2+16αem2−17m22αe−m226561e22αe+m22<0.By Lemma 3.1 in paper [21], there must exist a small neighbourhood of the origin O such that the system (70) is locally topologically equivalent to system.(72)x˙=y,y˙=f11xy+f30x3+f21+3e30x2y+f40−e30f11x4+4e40+f31+13e21f11+16f11f12x3y+Ox,y5.Hence the degenerate equilibrium pointE7 is a codimension 3 Bogdanov-Takens singularity (focus or center). Thus, we have the following theorem:Theorem 2.
In the cases (C4) and (C5), we have: (i) the equilibrium point E6 is a stable (unstable) node if A1<0A1>0; (ii) the equilibrium point E7 is a codimension 3 Bogdanov-Takens singularity (focus or center case).


### 2.6. Numerical Simulations and A Brief Summary

This subsection will give numerical simulations for above subsections, which is also a preliminaries of following bifurcations.


#### 2.6.1. The Case (C1)

Example 1.
For the case (C1), we set values of some parameters as follow: r1=0.6, α=0.5, a=1.5 and e=0.6. For the value λ=1/2, we have: when μ=3<μσ, the unique equilibrium point E3 is a stable multiple focus with multiplicity one since the negative first Lyapunov number σ≈−0.021655<0; when μ=μσ, the equilibrium point E3 is a weak stable focus of order 2; when μ=10>μσ, the equilibrium point E3 becomes an unstable multiple focus with multiplicity one; when μ=14>μ1, the equilibrium point E3 is a saddle point. A brief inspection in Figure 2 can reveal this phenomenon. Furthermore, it can be noticed that the condition (9) (iii) holds since a22−3a1a3<0 when μ=3.
Figures3 and 4 include curves of functions φA2λ,μ (red) and φσλ,μ (blue) in the case (C1) corresponding to λ=1/2 and λ=1/3, respectively. Furthermore, it can be noticed that when λ=1/3 with A2>0, μσ>μ1 or μσ does not exist.Figure 2
Phase diagrams aroundE3 in the case (C1) with λ=1/2: (a) A stable multiple focus with multiplicity one with μ=3; (b) An unstable multiple focus with multiplicity one with μ=10; (c) A stable weak focus of order 2 with μ=μσ. (d) A saddle point with μ=14.
(a)(b)(c)(d)Figure 3
Curves of functionsφA2λ,μ (in red) and φσλ,μ (in blue) in the case (C1) with λ=1/2.Figure 4
Curves of functionsφA2λ,μ (in red) and φσλ,μ (in blue) in the case (C1) with λ=1/3.


#### 2.6.2. Cases (C2) and (C3)

Moreover, as we takem2=αe/2λ=1/2, K1=23+22ar1/αe, m1=r1−3αe/2 and d=2−1α/8a to guarantee A2E42=0 in the case (C1), then A1E42=Δx=Δy=0 and our system coexists a degenerate equilibrium point E42:=x42,y42=3+22a,2ae2+2 and an unstable node E41:=x41,y41=1+22a,42ea. We denote this case as (C2).On the other hand, when we setm2=αe/3λ=1/3, m1=r1−4αe/3, K1=313+317ar1/8αe and d=−5+317α/72a, thus we have A1E52=A2E52=0, a nilpotent (or double-zero eigenvalue) E52≔x52,y52=13+317a/8,3ae7+17/8 and a stable node E51≔x51,y51=5+317a/4,317−1ea/2. This case is denoted as (C3).Obeying to Zhifen Zhang’s book [20] once more, for the case (C2) and the system (29), we denote it in analogy with the equations (64). Then an implicit function from the right hand side of the equation (64) is v=hu=α1−22/16ea2u3+⋯. Substituting v=hu into the equation (65), we have series(73)ψu=Ψu,hu=a2u2+⋯,a2=b20=2−2α2e32a≠0,Φuu,v+Ψvu,v∣v=hu=b1u+⋯,b1=2a20+b11=4−52α8a≠0.Hence the equilibrium pointE42 is a degenerate singular point (in the sense of [20]). The point E52 in the case (C3) is also a degenerate singular point since a2=−5+317α2e/432a≠0 and b1=α17−717/36a≠0.Example 2.
For the cases (C2) and (C3), we take values of parameters from the Example 1. Then a cusp of codimension 2 is E42≈8.742641,6.145584, and an unstable node is E41≈5.742641,5.091169, which is surrounded by an closed orbit (see Figure 5(a)); a cusp of codimension 2 is E52≈4.756747,3.754048, and a stable node is E51≈6.513494,4.216193 with characteristic direction θ=θ1≈0.234711, which appeals all trajectories in a small neighbourhood of it (see Figure 5(b)).Figure 5
(a) A cuspE42 of codimension 2 and an unstable node E41 in the case (C2). (b) A cusp E52 of codimension 2 and a stable node E51 in the case (C3).
(a)(b)


#### 2.6.3. The Cases (C4) and (C5)

Example 3.
For the cases (C4) and (C5), we take some parameters as α=0.6, e=0.9, m1=0.43, m2=0.2 and a=0.125. A codimension 3 BT singularity equilibrium point is E7≈0.470588,0.569531 with r1≈1.230417 (see Figure 6). When r1=1, a stable node is E6≈0.470588,0.405580 (see Figure 7(a)). While for r1=1.5, an unstable node is E6≈0.470588,0.761352, and the Poincare-Bendixson theorem yields that there exists a limit cycle enclosing this equilibrium (see Figure 7(b)).
Finally, at the end of this section, we will conclude stability and type of the equilibrium pointE3 as follows:(74)λ=25μ1=μσ,∄μσmultiple stable focus with multiplicity one,μ<μ1,codimension3BT singularityfocus or centerE3∈E7,μ=μ1,saddle point,μ>μ1;λ∈25,1μ1>μσ,∃μσmultiple stable focus with multiplicity one,μ∈0,μσ,weak stable focus of order2,μ=μσ,multiple unstable focus with multiplicity one,μ∈μσ,μ1,cusp of codimension2E3∈E42,μ=μ1,saddle point,μ>μ1;λ∈0,25μ1>μσ,∄μσmultiple stable focus with multiplicity one,μ∈0,μ1,cusp of codimension2E3∈E52,μ=μ1,saddle point,μ>μ1..
In whichE represents the “equivalence class” of equilibria, which have same properties in stability and type with the representative element E, symbols ∃μσ and ∄μσ is the meaning that μσ exists and μσ does not exist.Figure 6
A codimension 3 Bogdanov-Takens singularity (focus type) in the case (C5).Figure 7
Phase diagrams in the case (C4): (a) Stable node E6 with r1=1; (b) Unstable node E6 with r1=1.5.
(a)(b)


## 2.1. Preliminaries

In this subsection, we fully discuss preliminaries of the interior equilibrium point in the system (1). Here we denote it as E∗=x∗,y∗ for later use. An equilibrium point E∗ must satisfy following algebraic polynomial form equations(2)fx,y=r11−xK1a+x−αy−m1a+x=0,gx,y=αex−m2a+x−dya+x=0.In other words,x∗ and y∗ must be positive roots of following third-order polynomials (cubic equations) px=∑i=03aixi and py=∑i=03biyi, respectively, where coefficients are(3)a3=dr1,a2=dm1−r1K1+2adr1,a1=2am1−r1d+ααe−m2K1+a2dr1,a0=am1−r1d−m2αaK1,b3=K1d2,b2=−2K1dαe−m2,b1=αe−m22+ader1−m1K1+a2der1,b0=em1−r1αe−m2K1+am2r1a.We will mention an eliminant method in linear algebraic merely for completeness and shall not involve it in detail. If we sort above equations (2) as the form of fx,y=a0xy+a1x=0, gx,y=b0xy+b1x=0, the cubic equation px=0 can also be obtained again from the first Sylvester’s resultant Ryf,g=a0xb1x−a1xb0x. Similarly, if we sort the equations (2) as the form of fx,y=c0yx2+c1yx+c2y=0, gx,y=d0yx+d1y=0, the second resultant(4)Rxf,g=c0yc1yc2yd0yd1y00d0yd1y.Also deduces the cubic equationqy=0.Thanks to Scipione del Ferro, Niccolo Fontana(Tartaglia), Gerolamo Cardano and Shengjin Fan’s groundbreaking masterpieces, a third-order algebraic equation’s root(s) can be formulated by Cardano’s formula or Shengjin’s formula, and may qualitatively have one, two or three positive roots in the light of Rene Descartes’ rule of signs or Shengjin’s discriminant. Here we let these complicated expressions(5)px=a1a3−13a2a32,qx=227a2a33−13a2a3a1a3+a0a3,Δx=qx22+px33,and(6)py=b1b3−13b2b32,qy=227b2b33−13b2b3b1b3+b0b3,Δy=qy22+py33,be discriminants of above cubic equations px=0 and qy=0, respectively [5].The interior equilibrium pointE∗ does not exist when r1≤m1 or αe≤m2, thus we always assume r1>m1 and αe>m2 in the rest of this paper. The isoclines from the equations (2) and cobweb model show that an interior equilibrium point E∗ exists if condition(7)0<m2aαe−m2<K11−m1r1,holds. Noticing that the condition (7) deduces a0<0, b0<0, according to the zero theorem and(8)pm2aαe−m2=−αex2+m2a+x2e2r1a2α2dαe−m23<0,qαe−m2d=a2αe2r1>0.It is supposed that such condition (7) is well-done. In addition, the condition (7) and one of following conditions(9)ix2−a2≤m2aαe−m2,iia1,a2≥0,iiia22−3a1a3≤0,ivb22−3b1b3≤0,can lead to the uniqueness of the equilibrium point E∗, here x2=K11−m1/r1.


## 2.2. Multiple Focus with Multiplicity One

Due to the Routh-Hurwitz criterion and the Perron’s theorems, we merely need to consider the Jacobian matrix at a non-hyperbolic interior equilibrium pointE∗:(10)JE∗=Jij2×2=αx∗y∗a+x∗2−r1x∗K1−αx∗a+x∗αeay∗a+x∗2−dy∗.Furthermore, the trace, determinant and discriminant of matrixJE∗ are denoted as(11)A1:=A1E∗=trJE∗=−2r1x∗2+−αe−m1+m2+r1K1−r1ax∗+K1am2K1a+x∗,A2:=A2E∗=detJE∗ and Δ∗:=Δ∗E∗=A12−4A2, respectively. In order to obtain potential Hopf bifurcation, i.e. A1=0 and A2>0 (matrix JE∗ has a pair of pure imaginary eigenvalues), we take the threshold of parameter m1=r1−αe−m2, which can ensure that the numerator of A1 has a perfect square discriminant with respect of variable x∗. Thus the system (1) has an equilibrium point E3:=x3,y3=aλμ,λμ+1ea with parameters(12)m2=λαe,K1=μar1αe,d=−αλ−1μ+1λaλμ+12,and control variables λ∈0,1, μ>μm=1/1−λ. This special case is now denoted as (C1) for later use.Let us firstly discuss the equilibriumE3 more qualitatively, under the assumption of this case, the determinant is simplified as(13)A2E3=−λα2e2φA2λ,μλμ+12,φA2λ,μ=λ3μ2+−2μ2+2μλ2+μ−12λ−μ.The equationφA2λ,μ=0 has a quite clear positive root.(14)μ1=−2λ2+2λ+1+s2λ1−λ2>μm,s=−4λ2+4λ+1.That is to say: (i)A2E3>0 and E3 is a center or focus when μ∈μm,μ1; (ii) A2E3=0 and E3 is a potential cusp of codimension at least 2 when μ=μ1; (iii) A2E3<0 and E3 is just a saddle point when μ>μ1.This subsection will deal with above case (i). By using transformationx=u+x3,y=v+y3, the point E3 is translated to the origin O=0,0. Making a change of variables as u=−dy3X+βY, v=−αeay3/a+x32X, in which β=A2E3>0, then the system (1) becomes(15)X˙=F2X,Y=−βY+∑i+j=23aijXiYj+OX,Y4,Y˙=G2X,Y=βX+∑i+j=23bijXiYj+OX,Y4.Following [17, 18], the first Lyapunov number of the system (15) at the point E3, which is used to determine the stability of limit cycles as well, is presented as.(16)σ=3π2β3a30+b03+a12+b21−1β2a20b20−a02b02−a11a02+a20+b11b02+b20=3πα3e3λμφσλ,μ2a2λμ+14βφA2λ,μ,where(17)φσλ,μ=2+2λ4μ2+4μ−2μ2λ3+2−μλ2+μ+1λ.The transformationu=βX+dy3Y, v=αeay3/a+x32Y can also be used to deduce σ in (16).For any arbitraryλ, in view of a quadratic function φσλ,μ with respect to variable μ, it is a downward opening parabola with a symmetry axis μ=4λ2−λ+1/4λ21−λ>0 and has a positive maximum λ2+14λ+1/8λ1−λ. The equation φσλ,μ=0 straightly has a unique positive root(18)μσ=4λ2+λ2+14λ+1−λ+14λ21−λ>μm.What is more,μ1>μσ (or μσ exists) if λ>2/5; μ1<μσ if λ<2/5; μ1=μσ only if λ=2/5. Therefore, the equilibrium point E3 is respectively a stable(unstable if μσ exists) multiple focus with multiplicity one if μ<μσ>rbinμσ. When μ=μσ (if μσ exists), σ vanishes and one need to calculate the second Lyapunov coefficient l2. The generalized Hopf bifurcation of codimension 2 (Bautin bifurcation) may occur and two homocentric limit cycles may appear in some particular cases around the Bautin point E3 if l2≠0.


## 2.3. Weak Focus of Order 2

Whenμ=μσ, the first Lyapunov coefficient will be zero and we need to treat a weak (fine) focus E3 of order at least 2 in the so-called center-focus problem. This subsection starts with the system (15) or a standard form(19)x˙=−βy+φx,y,y˙=βx+ψx,y,where φx,y, ψx,y are both analytical functions and φ0,0=ψ0,0=0. Taking the polar coordinate transformation x=rcosθ, y=rsinθ, we derive a new system r˙=rRr,θ, θ.=β+Qr,θ, where Rr,θ=φcosθ+ψsinθ, Qr,θ=1/rψcosθ−φsinθ are also analytical functions. Now we rewrite the above system as(20)drdθ=rRβ+Q=∑k=2∞Rkθrk,here Rkθ are just polynomials of triangle functions cosθ and sinθ.We suppose that the special solution of the system (20) with initial condition rθ=0=c is rθ,c=∑k=1∞rkθck. Substituting it into the equation (20) and comparing all terms in ck, we have a series of coupled first-order differential equations(21)dr1dθ=0,dr2dθ=R2r12,dr3dθ=2R2r1r2+R3r13,⋯.These equations naturally yield solutionsr1=1, rk=gkθ+φkθk>1, where involved constants are(22)g2=12π∫02πR2τdτ,g3=12π∫02πR3τ+2R2τr2τdτ,and φkθ are periodic functions with a period 2π as well. For the system (19) with standard form, we obviously derive the first focal quantity g3 with g2=0, which repeats the formula (16) again, i.e. g3/σ≡constant.>0.Ifg3=0 or σ=0, from the integrals (22), we know g4=0 and the second focal quantity g5 (see the Appendix). From now on, with the system (15) at hand, the second focal quantity g5 reads(23)g5=−2097152λ8e4α42λ−16φg5λ3a4s1+3λ+1124λ2−λ+s1+141−λψg5λ7/2,in which auxiliary functions are ψg5λ=9λ3−λ2s1+3λ2+4λs1−3λ−s1−1 and(24)φg5λ=1+s1+54s1+61λ+1179s1+1533λ2+13252s1+20377λ3+79878s1+151954λ4+230064s1+604458λ5+57696s1+873264λ6+−1682544s1−2369616λ7+−6203262s1−14628750λ8+−11558132s1−36469478λ9+−6708162s1−43299230λ10+22132104s1+19510890λ11+88954156s1+206053636λ12+179705232s1+536172580λ13+268668984s1+899927400λ14+310805040s1+1188987000λ15+295112037s1+1226575797λ16+224630382s1+1060506033λ17+141939919s1+727205209λ18+69055348s1+415540245λ19+27134470s1+180678370λ20+6942224s1+61201034λ21+1283776s1+13517744λ22+19840s1+1703872λ23+−800s1−14240λ24+800λ25.Symbols1 is λ2+14λ+1. The resultant(25)Rs1φg5=10485765λ−22λ17λ−11412λ3+2λ2+2λ−1λ2+2λ+322λ2+λ+24≠0.Implies thatφg5λ>0 when λ∈2/5,1. While ψg5λ is also positive since Rs1ψg5=16λ2λ−15λ−2λ2+2λ+3≠0, then g5<0≠0. In other words, the equilibrium point E3 with σ=0 is a stable weak focus of order 2.


## 2.4. Cusp of Codimension 2

From the case (C1), undaunted by the scale of the distinction, this subsection reasonably set about finishing a puzzle when μ=μ1. Letting A2=0 in this case, we have threshold of parameters, a nilpotent E∗2=x∗2,y∗2 and an elementary equilibrium E∗1=x∗1,y∗1 as follows:(26)m1=r1−λ+1αe,d=λαμ1−λμ1−1aλμ1+12,x∗2=μ1aλ,y∗2=eaλμ1+1,x∗1=−λ2−λ+1s−λ2+λ+1aλs+1λ−1,y∗1=2aesλ+2λ2−2s−λ−2λ−1−2λ2+2λ+1+s.Firstly, we can make a transformation (I): x=X+x∗2, y=Y+y∗2 for the system (1) with above threshold of parameters:(27)X˙=F1X,Y,Y˙=G1X,Y.Secondly, we take the transformation(28)II:X=14s−2λ+1αeu+v,Y=αe2s−2λ+1s−1u8λ.Thus above system becomes Jordan’s standard form(29)u˙=F2u,v=v+a20u2+a11uv+a02v2+Ou,v3,v˙=G2u,v=b20u2+b11uv+b02v2+Ou,v3.By using the Lemma1 (or the Lemma 1 in [19]), this system (29) is equivalent to following system(30)x˙=y,y˙=d1λx2+d2λxy+Ox,y3,where discriminants manipulated by λ are(31)d1=d1λ=b20=α3e3λ−14φd1λ4a2λs+4λ−3s−5,d2=d2λ=b11+2a20=α2e2λ1−λ3s3+5−2λs2+19−12λs+23−18λas+1ψd2λ.Here temporarily auxiliary functions areφd1λ=2λs−2s2+8λ−s+1 and ψd2λ=s3+2λ−1s2+12λ−15s+14λ−17. We have added control variable s∈1,2.Lemma 1.
The following system(32)x˙=y+Ax2+Bxy+Cy2+Ox,y3,y˙=Dx2+Exy+Fy2+Ox,y3,is equivalent to the system(33)x˙=y,y˙=Dx2+E+2Axy+Ox,y3.After some nonsingular transformations in the neighbourhood of the originO.Letting functionφd1λ be zero and solving out the threshold of λ from the resultant (determinant)(34)Rsφd1=11−2λ2−1−8λ20011−2λ2−1−8λ210hλ0010hλ=4λ35λ−2,or equation φd1λ=0, it leads to a unique solution λ=2/5, where hλ=4λ2−4λ−1. And the derivative(35)dφd1λdλ=2ss2−1+8sλ+4λ1−λ>0.Implies that the functionφd1λ is a monotonic increasing function defined on interval 0,1. Combining the denominator of d1:(36)4a2λs+4λ−3s−5=4a2s+4λ−1−s+1<0.Here we know thatd1<>0 if λ><2/5. Notice that the numerator of d2 is obviously positive since the range of λ. In fact, we can assert that d2 is non zero. Indeed, for the function ψd2λ in the denominator of d2, its negativity will be found and well verified from an inequality(37)ψd2λ=s2s+λ−3s+1+λ−1s2+12s+14<s2−3s+1.That is to say,ψd2λ is always negative, namely d2 is negative here. All in all, d1d2≠0 when m2≠2/5αe, while d1d2=0 when m2=2/5αe.Next, we mainly concentrate on the stability of the equilibrium pointE∗1 whose trace, determinant and discriminant are parameterized as(38)A1E∗1=αe5−sφA1λ2s+12λs+4λ−3s−5,A2E∗1=αeαe−m2s−3φA2λ2s+12λs+4λ−3s−5,Δ∗E∗1=−m22φΔλs+14ψΔλ.In which some continuous auxiliary functions can be respectively expressed by(39)φA1λ=−2s3+2λ−7s2+−6λ+6s−16λ+19,φA2λ=s3+14λ−19s2+12λ−9s−34λ+43,φΔλ=∑j=19xjs9−j,ψΔλ=s4+4s3+8λ−6s2+28λ−34s+24λ−29,with(40)x1=1,x2=64λ−226,x3=124λ−354,x4=−1532λ+3156,x5=−2452λ+5010,x6=8248λ−11682,x7=13172λ−18786,x8=−10108λ+12208,x9=−16476λ+19889.The determinants(41)RsφA1=1−λ+723λ−38λ−192001−λ+723λ−38λ−19210hλ00010hλ00010hλ=165λ−2λ−15,RsφA2=114λ−1912λ−9−34λ+4300114λ−1912λ−9−34λ+4310hλ00010hλ00010hλ=1285λ−22λ−14,and(42)RsφΔ=1x2x3x4x5x6x7x8x9001x2x3x4x5x6x7x8x910hλ0000000010hλ0000000010hλ0000000010hλ0000000010hλ0000000010hλ0000000010hλ0000000010hλ=6553641λ2−54λ+81λ−1125λ−22,which are respectively from equations φA1λ=0, φA2λ=0 and φΔλ=0, yield the unique solution (zero point) λ=2/5 once more.For the sake of mathematical natures of the functionsφA1λ, φA2λ, φΔλ and ψΔλ defined on interval 0,1, we calculate their derivatives up to second order with respect to λ:(43)dφA1λdλ=φDφA1λs,d2φA1λdλ2=24s4+−48λ+72s3+−96λ2+168λ−72s2+96λ−1λ−1/22s3,dφA2λdλ=φDφA2λs,d2φA2λdλ2=−12s4+−336λ+264s3+48λ2−192λ+96s2−192λ3+336λ2−192λ+36s3,dφΔλdλ=φDφΔλs,d2φΔλdλ2=1s3−32s9+−5376λ+8120s8+768λ2−9696λ+11664s7+35840λ3−162400λ2+227440λ−125400s6+47616λ3−183552λ2+265536λ−153376s5+−367680λ3+1125120λ2−1146288λ+428520s4+−313856λ3+955136λ2−1035872λ+415984s3+395904λ3−956640λ2+781008λ−229448s2+161728λ−122λ−436361,dψΔλdλ=φDψΔλs,where new auxiliary continuous functions are(44)φDφA1λ=2s3+24λ−18s2+−16λ2+64λ−44s+24λ2−36λ+12,φDφA2λ=14s3+−12λ+18s2+−112λ2+208λ−110s−48λ2+60λ−18,φDφΔλ=64∑j=19yjs9−j,φDψΔλ=−16λ+16s3+−48λ+52s2+−64λ2+80λs−112λ2+192λ−68,in which(45)y1=1,y2=−32λ+14064,y3=−1792λ2+7224λ−469664,y4=−2976λ2+9984λ−670064,y5=30640λ2−78440λ+3980864,y6=39232λ2−99776λ+5325264,y7=−98976λ2+189672λ−8020064,y8=−105376λ2+202976λ−9162064,y9=40432λ2−69048λ+2441664.The following determinants(46)RsφDφA1=112λ−9−8λ2+32λ−2212λ2−18λ+600112λ−9−8λ2+32λ−2212λ2−18λ+610hλ00010hλ00010hλ=144λ−1420λ2−12λ−3,RsφDφA2=1−6λ+97−8λ2+104λ−557−24λ2+30λ−97001−6λ+97−8λ2+104λ−557−24λ2+30λ−9710hλ00010hλ00010hλ=1152495λ−25λ+1λ−14,RsφDφΔ=1y2y3y4y5y6y7y8y9001y2y3y4y5y6y7y8y910hλ0000000010hλ0000000010hλ0000000010hλ0000000010hλ0000000010hλ0000000010hλ0000000010hλ=1286560λ4−12808λ3+14305λ2−4395λ−15205λ−2λ−111,and(47)RsφDψΔ=1−48λ+52−16λ+16−64λ2+80λ−16λ+16−112λ2+192λ−68−16λ+16001−48λ+52−16λ+16−64λ2+80λ−16λ+16−112λ2+192λ−68−16λ+1610hλ00010hλ00010hλ=64λλ−15,which are respectively constructed from equations φDφA1λ=0, φDφA2λ=0, φDφΔλ=0 and φDψΔλ=0, yield that:(i)
The functionφA1λ has a (local) minimum φA13+26/10=227−1476/25 on interval 0,1 since d2/dλ2φA13+26/10=4386−1008>0, and, the detailed numerical simulation results are shown in Figure 1(a);(ii)
The functionφA2λ has a (local) minimum φA22/5=0 on interval 0,1 since d2/dλ2φA22/5=6480/49>0, and, the detailed numerical simulation results are shown in Figure 1(b);(iii)
The functionφΔλ has a local minimum φΔ2/5=0 and a local maximum φΔλ3≈534.507461 on interval 0,1 since d2/dλ2φΔ2/5=1156415616/30625>0 and d2/dλ2φΔλ3≈−13531.790941<0, where λ3 is another positive root of equation φDφΔλ=0 or the unique positive root of equation 6560λ4−12808λ3+14305λ2−4395λ−1520=0 in open interval 0,1, i.e.(48)λ3=19840M1/6N1/46−205NM2/3+−4040497M1/3+17198445275N+69931867473M1/2+4803M1/6N1/4−3N3/4≈0.764691,with(49)M=2783067866245+1917582022669769211,N=410M2/3−4040497M1/3−34396890550,and then, the detailed numerical simulation results are shown in Figure1(c); (iv) The function ψΔλ is a monotonic function on interval 0,1 since dψΔλ/dλ≠0, and, the detailed numerical simulation results are shown in Figure 1(d).Hence, maxima and minima of the functionsφA1λ, φA2λ, φΔλ and ψΔλ on unit interval 0,1 are listed as follows:(50)maxφA1λ=16,minφA1λ=227−147625<0,maxφA2λ=16,minφA2λ=0,maxφΔλ=9216,minφΔλ=0,maxψΔλ=−4,minψΔλ=−64.Combining Figures1(a)–1(d) and the determinant(51)RsψΔ=148λ−628λ−3424λ−2900148λ−628λ−3424λ−2910hλ000010hλ000010hhλλ000010hλ=256λ−18.Thus we have a rough summary:(i)
φA1λ is positive (negative) when λ<2/5λ>2/5, and it has a unique zero point 2/5;(ii)
φA2λ and φΔλ are nonnegative functions defined on interval 0,1, and they both have the unique zero point λ=2/5;(iii)
ψΔλ is a negative function defined on interval 0,1.Thus, it follows that: (i)A1E∗1 is negative (positive) when λ<2/5λ>2/5; (ii) A2E∗1 and Δ∗E∗1 are always positive when λ≠2/5. In a word, the equilibrium point E∗1 is a stable or unstable elementary node when λ≠2/5.On the whole, whenm2≠2/5αe, we have d1d2≠0, an interior equilibrium point E∗2 with A1=A2=0 from the case (C1) is a cusp of codimension 2, and an interior equilibrium point E∗1 is a stable (unstable) node when m2<2/5αem2>2/5αe. When m2=2/5αe, the unique interior equilibrium point E∗ becomes a cusp of codimension at least 3.Theorem 1.
From the case (C1), when m2≠2/5αe, an interior equilibrium point E∗2 with A1=A2=0 is a cusp of codimension 2 (Bogdanov-Takens bifurcation point) since d1d2≠0, and an interior equilibrium point E∗1 is an asymptotic stable (unstable) node when m2<2/5αem2>2/5αe. When m2=2/5αe, the unique interior equilibrium point E∗=4a,5ea is a cusp of codimension at least 3, indeed it is a codimension 3 BT singularity (focus or center case), see next subsection.Remark 1.
When we use the transformations (III):u=p+a02pq,v=q−a20p2, and (IV): p=w, q=z−c11wz, the system (29) firstly becomes(52)p˙=q+c11pq+Op,q3,q˙=d20p2+d11pq+d02q2+Op,q3,and then becomes(53)w˙=z+Ow,z3,z˙=f20w2+f11wz+f02z2+Ow,z3.
Finally, in order to obtain standard form, we construct a transformation (V):w=x1+1/2f02x12, z=f02x1y1+y1, above system becomes(54)x1.=y1+Ox1,y13,y1.=d1x12+d2x1y1+Ox1,y13.
This also ensures the Theorem1 and (31).Figure 1
(a) Figure of functionφA1λ on the interval 0,1; (b) figure of function φA2λ on the interval 0,1; (c) figure of function φΔλ on the interval 0,1; (d) figure of function ψΔλ on the interval 0,1.
(a)(b)(c)(d)


## 2.5. Bogdanov-Takens Singularity (Focus or Center) of Codimension 3

Lettingpx=Δx=0 or py=Δy=0, we obviously derive thresholds of parameters(55)K1=ar18αe+m2r1−m1αe−m2,d=8αe+m2αe−m2227αar1−m1e2,and a unique degenerate equilibrium point E6≔x6,y6=a2αe+m2/αe−m2,18r1−m1e2αa/8αe+m2αe−m2. In this case, it is quite obvious to see that A2E6=0, while the trace reduces to(56)A1E6=−16α2e2+6αem1−14αem2−6αer1+3m1m2−2m22−3m2r138αe+m2.The combinations in this case withA1E6≠0 is denoted as (C4). Furthermore, if A1E6=0 in the case (C4), thresholds of parameters r1, K1, d and a unique nilpotent equilibrium point E7 are(57)r1=2αe−m28αe+m2+3m12αe+m232αe+m2,K1=a2αe−m28αe+m2+3m12αe+m22αe−m22,d=2αe+m2αe−m218e2αa,E7≔x7,y7=a2αe+m2αe−m2,12e2αa2αe+m2.Similarly, we denote this special case as (C5). Note that the conditions (9) (iii) and (9) (iv) all hold in cases (C4) or (C5) since a22−3a1a3=b22−3b1b3=0.For the case (C4), we take transformations (I): x=X+x6, y=Y+y6 and(58)II:X=αx6ua+x6+v,Y=dy6u+3r1−m1e8αe+m2v,τ=tA1.The system (1) becomes a normal form(59)u˙=u+Φu,v,(60)v˙=Ψu,v,where Φu,v,Ψu,v=Ou,v2. The implicit function u=hv=h2v2+⋯ from right hand side of the equation (59) implies ψv=Ψhv,v=a3v3+⋯ with(61)a3=2r1−m1αe−m232αe+m2278αe+m2e2α2a2A12≠0.There is only one zero eigenvalue and a nonzero eigenvalue for the matrixJE6, by using the Theorem 7.1 in Zhifen Zhang’s book [20], it is easy to see that the equilibrium point E6 is a stable (unstable) node if A1<0A1>0.We will provide some explicitly smooth transformations to obtain a normal form with terms up to fourth order and determine the exact type of the equilibrium pointE7. Firstly, we take a transformation (I): x=X+x7, y=Y+y7 and translate the equilibrium point E7 to the origin O. Secondly, we take a linear transformation(62)II:u=−27αe−m22αe+m2Xaα8αe+m22−9αe−m22αe+m2Yaαe8αe+m22,v=−6αe−m22Xαa8αe+m2+3αe−m22αe+m2Yeaα8αe+m2,and transform matrix JE7 to its real Jordan’s canonical form, then above system becomes(63)u˙=F2u,v=v+∑i+j=24aijuivj+Ou,v5,v˙=G2u,v=∑i+j=24bijuivj+Ou,v5,where a11=b20=0 and a20a02b11b02≠0.Coming to the equation (63), we denote it in obvious notations:(64)u˙=v+Φu,v,(65)v˙=Ψu,v,where Φu,v,Ψu,v=Ou,v2. Substituting an implicit function v=hu=α2u2+α3u3+⋯ from right hand side of the first equation (64) into the second equation (65), where(66)α2=−28αe+m2αe−m281e,α3=4αe−m228αe+m22187e2.We immediately derive seriesψu=Ψu,hu=a3u3+⋯, Φuu,v+Ψvu,v∣v=hu=b1u+⋯ and a discriminant μ=b12+8a3, where(67)a3=b30−a20b11=−4αe−m228αe+m226561e2<0,b1=αe−m28αe+7m28αe+m281e2αe+m2≠0,μ=−αe−m228αe+m2264α2e2+16αem2−17m226561e22αe+m22<0.Thus the equilibrium pointE7 is a center or a focus.Thirdly, we take near identity transformations(68)III:u=p,v=q−a20p2−a02q2,IV:p=w+12d02w2,q=d02wz+z.To eliminate second order terms inu˙ and q2 term in q˙, thus above system firstly becomes(69)p˙=F3p,q=q+Op,q3,q˙=G3p,q=d11pq+d02q2+Op,q3,and finally becomes(70)w˙=F4w,z=z+∑i+j=3eijwizj+Ow,z4,z˙=G4w,z=f11wz+∑i+j=3fijwizj+Ow,z4.Taking into accountf11=b11+2a20 and the Lemma 1, the equilibrium point E7 is a cusp of codimension at least 3. Obviously, we notice that coefficients f11≠0, f30<0 and(71)5f30f21+3e30−3f11f40−e30f11=48αe+m2456α2e2−4αem2−7m22αe−m23143489072αe+m22e4≠0,f112+8f30=−8αe+m2264α2e2+16αem2−17m22αe−m226561e22αe+m22<0.By Lemma 3.1 in paper [21], there must exist a small neighbourhood of the origin O such that the system (70) is locally topologically equivalent to system.(72)x˙=y,y˙=f11xy+f30x3+f21+3e30x2y+f40−e30f11x4+4e40+f31+13e21f11+16f11f12x3y+Ox,y5.Hence the degenerate equilibrium pointE7 is a codimension 3 Bogdanov-Takens singularity (focus or center). Thus, we have the following theorem:Theorem 2.
In the cases (C4) and (C5), we have: (i) the equilibrium point E6 is a stable (unstable) node if A1<0A1>0; (ii) the equilibrium point E7 is a codimension 3 Bogdanov-Takens singularity (focus or center case).


## 2.6. Numerical Simulations and A Brief Summary

This subsection will give numerical simulations for above subsections, which is also a preliminaries of following bifurcations.


### 2.6.1. The Case (C1)

Example 1.
For the case (C1), we set values of some parameters as follow: r1=0.6, α=0.5, a=1.5 and e=0.6. For the value λ=1/2, we have: when μ=3<μσ, the unique equilibrium point E3 is a stable multiple focus with multiplicity one since the negative first Lyapunov number σ≈−0.021655<0; when μ=μσ, the equilibrium point E3 is a weak stable focus of order 2; when μ=10>μσ, the equilibrium point E3 becomes an unstable multiple focus with multiplicity one; when μ=14>μ1, the equilibrium point E3 is a saddle point. A brief inspection in Figure 2 can reveal this phenomenon. Furthermore, it can be noticed that the condition (9) (iii) holds since a22−3a1a3<0 when μ=3.
Figures3 and 4 include curves of functions φA2λ,μ (red) and φσλ,μ (blue) in the case (C1) corresponding to λ=1/2 and λ=1/3, respectively. Furthermore, it can be noticed that when λ=1/3 with A2>0, μσ>μ1 or μσ does not exist.Figure 2
Phase diagrams aroundE3 in the case (C1) with λ=1/2: (a) A stable multiple focus with multiplicity one with μ=3; (b) An unstable multiple focus with multiplicity one with μ=10; (c) A stable weak focus of order 2 with μ=μσ. (d) A saddle point with μ=14.
(a)(b)(c)(d)Figure 3
Curves of functionsφA2λ,μ (in red) and φσλ,μ (in blue) in the case (C1) with λ=1/2.Figure 4
Curves of functionsφA2λ,μ (in red) and φσλ,μ (in blue) in the case (C1) with λ=1/3.


### 2.6.2. Cases (C2) and (C3)

Moreover, as we takem2=αe/2λ=1/2, K1=23+22ar1/αe, m1=r1−3αe/2 and d=2−1α/8a to guarantee A2E42=0 in the case (C1), then A1E42=Δx=Δy=0 and our system coexists a degenerate equilibrium point E42:=x42,y42=3+22a,2ae2+2 and an unstable node E41:=x41,y41=1+22a,42ea. We denote this case as (C2).On the other hand, when we setm2=αe/3λ=1/3, m1=r1−4αe/3, K1=313+317ar1/8αe and d=−5+317α/72a, thus we have A1E52=A2E52=0, a nilpotent (or double-zero eigenvalue) E52≔x52,y52=13+317a/8,3ae7+17/8 and a stable node E51≔x51,y51=5+317a/4,317−1ea/2. This case is denoted as (C3).Obeying to Zhifen Zhang’s book [20] once more, for the case (C2) and the system (29), we denote it in analogy with the equations (64). Then an implicit function from the right hand side of the equation (64) is v=hu=α1−22/16ea2u3+⋯. Substituting v=hu into the equation (65), we have series(73)ψu=Ψu,hu=a2u2+⋯,a2=b20=2−2α2e32a≠0,Φuu,v+Ψvu,v∣v=hu=b1u+⋯,b1=2a20+b11=4−52α8a≠0.Hence the equilibrium pointE42 is a degenerate singular point (in the sense of [20]). The point E52 in the case (C3) is also a degenerate singular point since a2=−5+317α2e/432a≠0 and b1=α17−717/36a≠0.Example 2.
For the cases (C2) and (C3), we take values of parameters from the Example 1. Then a cusp of codimension 2 is E42≈8.742641,6.145584, and an unstable node is E41≈5.742641,5.091169, which is surrounded by an closed orbit (see Figure 5(a)); a cusp of codimension 2 is E52≈4.756747,3.754048, and a stable node is E51≈6.513494,4.216193 with characteristic direction θ=θ1≈0.234711, which appeals all trajectories in a small neighbourhood of it (see Figure 5(b)).Figure 5
(a) A cuspE42 of codimension 2 and an unstable node E41 in the case (C2). (b) A cusp E52 of codimension 2 and a stable node E51 in the case (C3).
(a)(b)


### 2.6.3. The Cases (C4) and (C5)

Example 3.
For the cases (C4) and (C5), we take some parameters as α=0.6, e=0.9, m1=0.43, m2=0.2 and a=0.125. A codimension 3 BT singularity equilibrium point is E7≈0.470588,0.569531 with r1≈1.230417 (see Figure 6). When r1=1, a stable node is E6≈0.470588,0.405580 (see Figure 7(a)). While for r1=1.5, an unstable node is E6≈0.470588,0.761352, and the Poincare-Bendixson theorem yields that there exists a limit cycle enclosing this equilibrium (see Figure 7(b)).
Finally, at the end of this section, we will conclude stability and type of the equilibrium pointE3 as follows:(74)λ=25μ1=μσ,∄μσmultiple stable focus with multiplicity one,μ<μ1,codimension3BT singularityfocus or centerE3∈E7,μ=μ1,saddle point,μ>μ1;λ∈25,1μ1>μσ,∃μσmultiple stable focus with multiplicity one,μ∈0,μσ,weak stable focus of order2,μ=μσ,multiple unstable focus with multiplicity one,μ∈μσ,μ1,cusp of codimension2E3∈E42,μ=μ1,saddle point,μ>μ1;λ∈0,25μ1>μσ,∄μσmultiple stable focus with multiplicity one,μ∈0,μ1,cusp of codimension2E3∈E52,μ=μ1,saddle point,μ>μ1..
In whichE represents the “equivalence class” of equilibria, which have same properties in stability and type with the representative element E, symbols ∃μσ and ∄μσ is the meaning that μσ exists and μσ does not exist.Figure 6
A codimension 3 Bogdanov-Takens singularity (focus type) in the case (C5).Figure 7
Phase diagrams in the case (C4): (a) Stable node E6 with r1=1; (b) Unstable node E6 with r1=1.5.
(a)(b)


## 2.6.1. The Case (C1)

Example 1.
For the case (C1), we set values of some parameters as follow: r1=0.6, α=0.5, a=1.5 and e=0.6. For the value λ=1/2, we have: when μ=3<μσ, the unique equilibrium point E3 is a stable multiple focus with multiplicity one since the negative first Lyapunov number σ≈−0.021655<0; when μ=μσ, the equilibrium point E3 is a weak stable focus of order 2; when μ=10>μσ, the equilibrium point E3 becomes an unstable multiple focus with multiplicity one; when μ=14>μ1, the equilibrium point E3 is a saddle point. A brief inspection in Figure 2 can reveal this phenomenon. Furthermore, it can be noticed that the condition (9) (iii) holds since a22−3a1a3<0 when μ=3.
Figures3 and 4 include curves of functions φA2λ,μ (red) and φσλ,μ (blue) in the case (C1) corresponding to λ=1/2 and λ=1/3, respectively. Furthermore, it can be noticed that when λ=1/3 with A2>0, μσ>μ1 or μσ does not exist.Figure 2
Phase diagrams aroundE3 in the case (C1) with λ=1/2: (a) A stable multiple focus with multiplicity one with μ=3; (b) An unstable multiple focus with multiplicity one with μ=10; (c) A stable weak focus of order 2 with μ=μσ. (d) A saddle point with μ=14.
(a)(b)(c)(d)Figure 3
Curves of functionsφA2λ,μ (in red) and φσλ,μ (in blue) in the case (C1) with λ=1/2.Figure 4
Curves of functionsφA2λ,μ (in red) and φσλ,μ (in blue) in the case (C1) with λ=1/3.


## 2.6.2. Cases (C2) and (C3)

Moreover, as we takem2=αe/2λ=1/2, K1=23+22ar1/αe, m1=r1−3αe/2 and d=2−1α/8a to guarantee A2E42=0 in the case (C1), then A1E42=Δx=Δy=0 and our system coexists a degenerate equilibrium point E42:=x42,y42=3+22a,2ae2+2 and an unstable node E41:=x41,y41=1+22a,42ea. We denote this case as (C2).On the other hand, when we setm2=αe/3λ=1/3, m1=r1−4αe/3, K1=313+317ar1/8αe and d=−5+317α/72a, thus we have A1E52=A2E52=0, a nilpotent (or double-zero eigenvalue) E52≔x52,y52=13+317a/8,3ae7+17/8 and a stable node E51≔x51,y51=5+317a/4,317−1ea/2. This case is denoted as (C3).Obeying to Zhifen Zhang’s book [20] once more, for the case (C2) and the system (29), we denote it in analogy with the equations (64). Then an implicit function from the right hand side of the equation (64) is v=hu=α1−22/16ea2u3+⋯. Substituting v=hu into the equation (65), we have series(73)ψu=Ψu,hu=a2u2+⋯,a2=b20=2−2α2e32a≠0,Φuu,v+Ψvu,v∣v=hu=b1u+⋯,b1=2a20+b11=4−52α8a≠0.Hence the equilibrium pointE42 is a degenerate singular point (in the sense of [20]). The point E52 in the case (C3) is also a degenerate singular point since a2=−5+317α2e/432a≠0 and b1=α17−717/36a≠0.Example 2.
For the cases (C2) and (C3), we take values of parameters from the Example 1. Then a cusp of codimension 2 is E42≈8.742641,6.145584, and an unstable node is E41≈5.742641,5.091169, which is surrounded by an closed orbit (see Figure 5(a)); a cusp of codimension 2 is E52≈4.756747,3.754048, and a stable node is E51≈6.513494,4.216193 with characteristic direction θ=θ1≈0.234711, which appeals all trajectories in a small neighbourhood of it (see Figure 5(b)).Figure 5
(a) A cuspE42 of codimension 2 and an unstable node E41 in the case (C2). (b) A cusp E52 of codimension 2 and a stable node E51 in the case (C3).
(a)(b)


## 2.6.3. The Cases (C4) and (C5)

Example 3.
For the cases (C4) and (C5), we take some parameters as α=0.6, e=0.9, m1=0.43, m2=0.2 and a=0.125. A codimension 3 BT singularity equilibrium point is E7≈0.470588,0.569531 with r1≈1.230417 (see Figure 6). When r1=1, a stable node is E6≈0.470588,0.405580 (see Figure 7(a)). While for r1=1.5, an unstable node is E6≈0.470588,0.761352, and the Poincare-Bendixson theorem yields that there exists a limit cycle enclosing this equilibrium (see Figure 7(b)).
Finally, at the end of this section, we will conclude stability and type of the equilibrium pointE3 as follows:(74)λ=25μ1=μσ,∄μσmultiple stable focus with multiplicity one,μ<μ1,codimension3BT singularityfocus or centerE3∈E7,μ=μ1,saddle point,μ>μ1;λ∈25,1μ1>μσ,∃μσmultiple stable focus with multiplicity one,μ∈0,μσ,weak stable focus of order2,μ=μσ,multiple unstable focus with multiplicity one,μ∈μσ,μ1,cusp of codimension2E3∈E42,μ=μ1,saddle point,μ>μ1;λ∈0,25μ1>μσ,∄μσmultiple stable focus with multiplicity one,μ∈0,μ1,cusp of codimension2E3∈E52,μ=μ1,saddle point,μ>μ1..
In whichE represents the “equivalence class” of equilibria, which have same properties in stability and type with the representative element E, symbols ∃μσ and ∄μσ is the meaning that μσ exists and μσ does not exist.Figure 6
A codimension 3 Bogdanov-Takens singularity (focus type) in the case (C5).Figure 7
Phase diagrams in the case (C4): (a) Stable node E6 with r1=1; (b) Unstable node E6 with r1=1.5.
(a)(b)


## 3. Local Bifurcations

In this section, more words are included here about the existence of Hopf bifurcation curve in a small neighbourhood of the equilibrium pointE3 in the case (C1), BT bifurcations of codimension 2 in the cases (C2) and (C3), and a degenerate focus type BT bifurcation of codimension 3 in the case (C5), respectively.


### 3.1. Hopf Bifurcation Curve AroundE3

For the case (C1), the system (1) undergoes a nondegenerate Hopf bifurcation around the equilibrium point E3 as σ≠0. The Hopf bifurcation is supercritical (subcritical) and limit cycles generated by the critical point are stable (unstable) if σ<0σ>0. On occasion, there may exist some parameter values such that σ=0 or the system (1) may undergo a degenerate Hopf bifurcation for some values of parameters [17].Conversely, we also study existence of Hopf bifurcation curve in this special case. Choosingm2 and d as bifurcation parameters and introducing sufficiently small parameter λ1,λ2 . For an unfolding system(75)x˙=xr11−xαeμar1+eλ+1α−r1−αya+x,y˙=yxαea+x−λαe−λ1−−λαλμ−μ+1aλμ+12+λ2y.We suppose that the above system has an equilibrium pointE∗=x∗,y∗ when λ1,λ2≠0, where x∗=x3+w with w≪1. Calculating A1=0 and A2>0, we have a solution(76)λ1=weα2aλμ+a+2waaλμ+a+wμ,λ2=−αwaaμ−wλμ+12λμ+1a+w2λ3+λμ2+2λ2+λ−1μ+λ+1λμ+1a2+2wλ3−1/2λ2+1/2λμ2+2λ2+1/2λμ+λ+1a+w21+λ−1μλ,with(77)A2=−e2α2μ2a2aλμ+a+w2aλμ+wλλ−12μ2+2λ2−2λ−1μ+λμa3+5w15+λ2−65λ+15μ2+65λ−15μa2+8w21/2+λ−1/2μa+4w3.At this point, the Hopf bifurcation curve of the system (75) is defined by(78)Hp=λ1,λ2∣λ1,λ2satisfy3.2,and the approximation of Hp is a straight line with slope(79)k=limw⟶0λ2wλ1w=−μ2λ3−2μλ2+−μ2−μ−1λ+μ−1aλμ+12e2λμ+1.In a small neighbourhood of the origin in parameter plane.Similarly, for the bifurcation parametersm1, m2 and corresponding unfolding system(80)x˙=xr11−xαeμar1+eλ+1α−r1−λ1−αya+x,y˙=yxαea+x−λαe−λ2+λαλμ−μ+1yaλμ+12.The Hopf bifurcation curve of the system (80) is defined by(81)λ1=−wαa2λ4μ3+3aa+2/3wμ2λ3+μμ2+μ+3a2−wμ−4a+w2λ2+2μ+1a2+wμ2+μ+2a−w2μ−1λ−μ−1a−2waea2λ4μ3+3aa+2/3wμ2λ3+3a2−1/3wμ−4a+1/3w2μλ2+a2+2wa−w2μ−1λ+waaμ,λ2=wαa3λ5μ4+4a+3/4wa2μ3λ4−μ2+μ−6a2+wμ−9a−3w2aμ2λ3−2μ−2a3+wμ−3/2μ+3a2+w2μ−3a−1/2w3μλ2+a+wμ2−μ+1a2−wμ+2μ−1a−w2μ−1λ+a2μweaλμ+a+wa2λ4μ3+3aa+2/3wμ2λ3+3a2−1/3wμ−4a+1/3w2μλ2+a2+2wa−w2μ−1λ+waaμ,with A2>0. Hence the slope of approximation around the origin is(82)k=limw⟶0λ2wλ1w=−λ1+λ3−λμ3+3λ2−λ+1μ2+3λ−1μλ3+λμ2+2λ2+λ−1μ+λ+1λμ+1.From the Example1 with values of parameters r1=0.6, α=0.5, a=1.5 and e=0.6, Figure 8 show the Hopf bifurcation curves with respect to (i) λ=1/2, μ=3, (ii) λ=1/2, μ=10, (iii) λ=1/3, μ=3 and (iv) λ=1/3, μ=10, respectively. Following subsections will make further efforts to illustrate these curvesFigure 8
(a) Hopf bifurcation curve withλ=1/2 and μ=3 for the case (C2); (b) Hopf bifurcation curve with λ=1/2 and μ=10 for the case (C2); (c) Hopf bifurcation curve with λ=1/3 and μ=3 for the case (C3); (d) Hopf bifurcation curve with λ=1/3 and μ=10 for the case (C3).
(a)(b)(c)(d)


### 3.2. BT Bifurcations of Codimension 2 AroundE42

In this subsection, we firstly choosem2 and d as bifurcation parameters, and then investigate BT bifurcation of codimension 2 in the case (C2) by following the techniques and steps in [22, 23] for an unfolding system(83)x˙=xr11−xK1−m1−αya+x,y˙=yαexa+x−m2+λ1−d+λ2y,where parameters λ1 and λ2 are sufficiently small. Denoting a provisional parameter vector λ=λ1,λ2 in a small neighbourhood of the origin O, we firstly select a linear transformation (I): x=X+x42, y=Y+y42 to rewrite above system as(84)X˙=F1X,Y=∑i+j=12aijλXiYj+OX,Y3,Y˙=G1X,Y=∑i+j=02bijλXiYj+OX,Y3.Secondly, we should be able to make an affine transformation (II):u=X, v=a10λX+a01λY, then the system (84) is transformed into(85)u˙=F2u,v=v+c20λu2+c11λuv+Ou,v3,u˙=G2u,v=∑i+j=02dijλuivj+Ou,v3,without demur. Thirdly, letting a transformation (III): p=u−c11λ+d02λ/2u2, q=c20λu2−d02λuv+v, we derive a new system(86)p˙=F3p,q=q+Op,q3,(87)q˙=G3p,q=∑i+j=02fijλpiqj+Op,q3,where f02λ=0. Finally, in order to eliminate all higher order terms in the equation (86), we mechanically construct a transformation (IV): w=p, z=F3p,q, then the system (86) is transformed into(88)w˙=F4w,z=z,z˙=G4w,z=∑i+j=02hijλwizj+Ow,z3,where h00λ=f00λ, h10λ=f10λ, h01λ=f01λ and h02λ=0. Denoting functions depended on parameters λ1 and λ2:(89)d1λ=f20λ,d1=d10=4−32α2e232a,d2λ=h11λ,d2=d20=αe92−148a,ξ1λ=f00λd1λ,ξ2λ=f10λ−d1λ,ηλ=−f01λ−d1λ,Φw,λ=f00λ+f10λw+f20λw2+⋯,and noticing that d1,d2<0, then we take a time transformation (V): X1=−w, X2=z, τ=−t to make it positive and use the Malgrange preparation theorem [24] to obtain a decomposition(90)−ΦX1,λ=ξ1λ+ξ2λX1+X12ΨX1,λ,Ψ0,λ=−d1λ.Applying a transformation (VI):Y1=X1, Y2=X2/ΨX1,λ, dτ=ΨX1,λdt and a parameter dependent affine transformation (VII): x=Y1+1/2ξ2λ, y=Y2, we finally transform above system into the normal form, which still use symbol t as time variable:(91)x˙=y,y˙=μ1λ+μ2λy+x2+d2−d1xy+Ox,y3,where μ1λ=ξ1λ−1/4ξ2λ2 and μ2λ=ηλ−d2/2−d1ξ2λ. Since the Jacobian(92)∂μ1,μ2∂λ1,λ2∣λ=0=−12832−415912+2250a7/2eα2≠0,We point out that the above transformation is nonsingular, and the results in [25–27] yield that the system (91) is strongly topologically equivalent to system(93)x˙=y,y˙=μ1+μ2y+x2−xy.In other words, the system (91) will become a standard form(94)u˙=v,v˙=ϵ1+ϵ2v+u2−uv+⋯.Under a transformation:x=−d1/d22u, y=−−d1/d23v, t=−d2/−d1τ, where ϵ1=d2/−d14μ1 and ϵ2=−d2/−d1μ2. Thus, the system (83) will undergo a BT bifurcation with bifurcation parameters μ1 and μ2 when parameter λ varies in as small neighbourhood of the origin. Hence, we have following theorem.Theorem 3 (Bogdanov-Takens bifurcation of codimension 2 aroundE42).
For the case (C2), in a small neighbourhood of the equilibrium point E42, the system (83) undergoes a Bogdanov-Takens bifurcation of codimension 2 when parameter λ varies in a small neighbourhood of the origin when we choose m2 and d as bifurcation parameters. At the same time, there exist values of parameters such that this system (83) has a limit cycle or a homoclinic loop surround the cusp E42.Finally, based on the above Theorem3, the local representations of bifurcation curves up to third-order approximations in a small neighbourhood of the origin with slope k=−2+2/4ae are presented as follows, including description of saddle-node (SN), Hopf (H) and homoclinic (HL) bifurcation curves [28, 29]. The slope k can be viewed as the limiting case of the slope (79) when μ⟶μ1. These bifurcation curves can divide a small neighbourhood of the origin in the parameter plane into several regions, which can exhibit dynamics of the system (83).(i)
The saddle-node bifurcation curve is formulated by(95)SN=λ∣μ1λ=0,μ2λ≠0=λ∣λ2=2−24aeλ1+Oλ12,and(96)SN+=λ∣μ1λ=0,μ2λ>0,SN−=λ∣μ1λ=0,μ2λ<0.(ii)
The Hopf bifurcation curve is formulated by(97)H=λ|μ2λ2−d22d1μ1λ=0,μ1λ<0=λ|−1681+562aαeλ1−32274+1932a2αλ2+167132+1000aα2e2λ12+12810672+2269a2eα2λ1λ2+25634852+4927a3α2λ22−643139+22132ae3α3λ13−128317532+44932a2α3e2λ12λ2−2561507062+213155a3eα3λ1λ22−10241176092+166328a4α3λ23+Oλ1,λ24=0,μ1λ<0.(iii)
The homoclinic bifurcation curve is formulated by(98)HL=λ∣μ2λ2−25d2249d1μ1λ=0,μ1λ<0=λ∣−40081+562a49αeλ1−800274+1932a249αλ2+16178972+25096a49α2e2λ12+128408592+57685a249α2eλ1λ2+256937672+132559a349α2λ22−6478787+555412a49e3α3λ13−1288075292+1142740a249α3e2λ12λ2−25639950742+5650643a349eα3λ1λ22−102433584012+4749632a449α3λ23+Oλ1,λ24=0,μ1λ<0.Example 4.
The values of parameters will be recalled from the Example2. Figure 9 gives the saddle-node, Hopf and homoclinic bifurcation curves when the value of parameter λ falls in a small neighbourhood of the origin in the parameter plane. Figures 10–13 depict the dynamics of the system (83) in a small neighbourhood of the origin in the parameter plane when we set values of λ1 and λ2.(i)
When the value ofλ lies on the saddle-node bifurcation curve SN−, there exist two interior equilibria and one is unstable node.(ii)
When the value ofλ lies on region I (the region below the saddle-node bifurcation curve), there exists a unique unstable node. The Poincre-Bendixson theorem implies that there is a large limit cycle enclosing this equilibrium point, the detailed numerical simulation results are shown in Figure 10(a) with λ1,λ2≈1×10−4,−3.255138×10−5.(iii)
When the value ofλ lies on another saddle-node bifurcation curve SN+, there exists a unique unstable node and a limit cycle.(iv)
When the value ofλ lies on region IV (the region between the saddle-node bifurcation curve SN+ and the homoclinic curve), there exist an unstable node, a saddle, a stable focus, and the homoclinic loop in case (v) is broken, in which the saddle and the focus are bifurcated from the curve SN+, a large limit cycle and enlarged phase diagrams around a stable focus can be seen in Figure 11 with λ1,λ2≈1×10−4,−8.131779×10−6.(v)
When the value ofλ lies on the homoclinic curve, there exist three interior equilibrium points and a homoclinic loop enclosing a stable focus, the detailed numerical simulation results are shown in Figure 12(a) with λ1,λ2≈1×10−4,−1.626356×10−5 and enlarged phase diagrams can be seen in Figure 12(b).(vi)
When the value ofλ crosses the homoclinic curve into region III (the region between the homoclinic curve and the Hopf curve), there exist an unstable node, a saddle and a stable focus, the detailed numerical simulation results are shown in Figure 13 with λ1,λ2≈1×10−4,−1.626653×10−5.(vii)
When the value ofλ lies on the Hopf curve, there exist an unstable node, a saddle, an non-hyperbolic equilibrium point (multiple focus or center) according to Hopf bifurcation and a large limit cycle enclosing these equilibria.(viii)
When the value ofλ crosses the Hopf curve into region II (the region between the Hopf curve and the saddle-node bifurcation curve SN−), there exist an unstable node, a saddle and an unstable focus, the detailed numerical simulation results are shown in Figure 10(b) with λ1,λ2≈1×10−4,−1.627259×10−5. At the same time, it can be noticed that the stable focus in case (vi) can become unstable.Figure 9
Saddle-node (in red), hopf (in green) and homoclinic (in blue) bifurcation curves in case (C2).Figure 10
Phase diagrams in cases (ii) and (viii).
(a)(b)Figure 11
Phase diagrams in the case (iv).
(a)(b)Figure 12
Phase diagrams in case (v).
(a)(b)Figure 13
Phase diagrams in case (vi).
(a)(b)


### 3.3. BT Bifurcation of Codimension 2 AroundE52

In this subsection, we set out to deal with BT bifurcation of codimension 2 around the equilibrium pointE52 in the case (C3). Here we choose m1 and m2 as bifurcation parameters and rewrite the system (1) as the unfolding form:(99)x˙=xr11−xK1+m1+λ1−αya+x,y˙=yαexa+x−m2+λ2−dy,where λ1 and λ2 are sufficiently small parameters and λ=λ1,λ2. Firstly, we take a mere linear transformation (I): x=X+x52,y=Y+y52, the equilibrium point E52 is translated to the origin O:(100)X˙=F1X,Y=∑i+j=02aijλXiYj+OX,Y3,Y˙=G1X,Y=∑i+j=02bijλXiYj+OX,Y3.Secondly, we take a transformation (II):u=X, v=F1X,Y, then the above system becomes(101)u˙=F2u,v=v,v˙=G2u,v=∑i+j=02dijλuivj+Ou,v3.Thirdly, we take a transformation (III):p=u+d01λ/d11λ, q=v since d110,0=−85+1917αe/36a≠0, then we have following system with f01λ=0:(102)p˙=F3p,q=q,q˙=G3p,q=∑i+j=02fijλpiqj+Op,q3.Letting (IV):w=p, z=1−f02λpq, dt=1−f02λpdτ and rewrite symbol τ as t, we derive a new system(103)w˙=F4w,z=z,z˙=G4w,z=∑i+j=02hijλwizj+Ow,z3.It can be noticed thath200,0=33−717α2e2/432a>0 and h110,0=d110,0≠0, so h20λ>0 and h11λ≠0 when λ changes in a small neighbourhood of the origin, which depends smoothly on λ. Finally, we construct a time transformation (V): m=h11λ2/h20λw, n=h11λ3/h20λ2z, dt=h11λ/h20λdτ and rewrite symbol τ as t, then the above system (103) becomes(104)m˙=F5m,n=n,n˙=G5m,n=μ1λ+μ2λm+m2+mn+Om,n3,where(105)μ1=μ1λ=5781−17αeλ1+8679−17αeλ2+Oλ1,λ22,μ2=μ2λ=2−51+1117αeλ1+6917−85αeλ2+Oλ1,λ22.Therefore, owing to the Jacobian ofμ1 and μ2(106)∂μ1,μ2∂λ1,λ2|λ=0=346885+1917α2e2≠0,or(107)−1∂μ1/∂λ2⋅∂μ1,μ2∂λ1,λ2|λ=0<0.The system (99) is a generic family unfolding at the codimension 2 cusp E52, and we have following local representations of the bifurcation curves up to second-order approximations with slope k=1+17/12 for the system (104) [30]. Thus, it can be noticed that this slope k can be viewed as the limiting case of the slope (82) when μ⟶μ1.(i)
The saddle-node bifurcation curve is formulated by(108)SN=λ|μ1=14μ22=λ|5781−17αeλ1+8679−17αeλ2+179775−63172α2e2λ12+2551003−49917α2e2λ1λ2+3063111−62317α2e2λ22+Oλ1,λ23=0.(ii)
The Hopf bifurcation curve is formulated by(109)H=λ|μ1=0,μ2<0=λ|5781−17αeλ1+8679−17αeλ2+513441−65172α2e2λ12+515723−265917α2e2λ1λ2+1224841−16717α2e2λ22+Oλ1,λ23=0.(iii)
The homoclinic bifurcation curve is formulated by(110)HL=λ∣μ1=−625μ22,μ2<0=λ|5781−17αeλ1+8679−17αeλ2+5190409−26811750α2e2λ12+51160067−704111725α2e2λ1λ2+122422543−44451725α2e2λ22+Oλ1,λ23=0.Theorem 4.
(Bogdanov-Takens bifurcation of codimension 2 aroundE52). For the case (C3), in a small neighbourhood of the equilibrium point E52, there exist values of parameters such that the system (99) undergoes an attracting Bogdanov-Takens bifurcation of codimension 2 when parameter λ varies in a small neighbourhood of the origin with bifurcation parameters m1 and m2. And this system is a generic family unfolding at the cusp E52 of codimension 2.Example 5.
Here we extract values of parameters in the Subsection3.2 and the Example 3. Figure 14 presents the saddle-node, Hopf and homoclinic bifurcation curves when the value of parameter λ falls in a small neighbourhood of the origin in the parameter plane, which can divide a small neighbourhood of the origin in the parameter plane into several regions and exhibit different dynamical behavior of the system (99).(i)
When the value ofλ crosses the λ1 axis into region I (the region between the saddle-node bifurcation curve SN2 and the homoclinic bifurcation curve), there exist an unstable focus, a saddle and a stable node, in which the focus and saddle are bifurcated from the curve SN2.(ii)
When the value ofλ lies on the homoclinic bifurcation curve, there exist a saddle, a stable node and a homoclinic loop enclosing an unstable focus.(iii)
When the value ofλ crosses the homoclinic bifurcation curve into region II (the region between the homoclinic bifurcation curve and the Hopf bifurcation curve), there exist an unstable focus, a saddle and a stable node.(iv)
When the value ofλ lies on the Hopf bifurcation curve, there exist three interior equilibria including a saddle and a stable node.(v)
When the value ofλ crosses the Hopf bifurcation curve into region III (the region between the Hopf bifurcation curve and the saddle-node bifurcation curve SN1), there exist a stable focus, a saddle and a stable node. Furthermore, ti shall be noticed that the unstable focus in case (iii) becomes stable, which can ensure potential Hopf bifurcation.(vi)
When the value ofλ lies on the saddle-node bifurcation curve SN1, there exist two interior equilibria including a stable node.(vii)
When the value ofλ crosses the saddle-node bifurcation curve into region IV (the region above the saddle-node bifurcation curve), there exists a unique stable node.(viii)
When the value ofλ lies on the saddle-node bifurcation curve SN2, the unique equilibrium point is a stable node.Figure 14
Saddle-node (in red), hopf (in green) and homoclinic (in blue) bifurcation curves in the case (C3).


### 3.4. Degenerate Focus Type BT Bifurcation of Codimension 3

We should always keep the case (C5) in mind and choose r1, K1 and d as bifurcation parameters and use the method in [14]. Thus we introduce sufficiently small variables α1, α2, α3, a parameter vector A=α1,α2,α3, the origin O=0,0,0 in parameter space and an unfolding form system(111)x˙=xr1+α11−xK1+α2−αxya+x−m1x,y˙=αexya+x−m2y−d+α3y2.Of the system (1) since the Theorem 2 (ii). As indicated in Subsection 2.5 for the system (111), by using the coordinate transformations (I), (II), (III) and (IV), actual calculations can yield a new system(112)w˙=F4w,z=∑i+j=03eijAwizj+Ow,z4,(113)z˙=G4w,z=∑i+j=03fijAwizj+Ow,z4.Here we omit complicated, but consider inconsequential expressions of smooth functionseijA and fijA for the sake of convenience, while eij0=eij,fij0=fij, for instance, e000=e100=e200=e110=e020=0 and f000=f100=f010=f200=f020=0.Secondly, in order to eliminate third-order terms in right hand side of equation (112) when A=0, we construct a transformation(114)V:w=x1+e03x1y12,z=−e12x1y12−e21x12y1−e30x13+y1.Above system becomes(115)x1.=F5x1,y1=y1+∑i+j=03gijAx1iy1j+Ox1,y14,y1.=G5x1,y1=∑i+j=03hijAx1iy1j+Ox1,y14.Similarly, coefficientsgijA, hijA are all omitted for space and can be deduced by coefficients eijA and fijA. Moreover, it shall be noticed that gij0=0i+j≤3. Finally, the transformation (VI): x2=x1, y2=F5x1,y1 can translate the above system into the following system(116)x2.=F6x2,y2=y2,y2.=G6x2,y2=∑i+j=03kijAx2iy2j+Ox2,y24,where coefficients kijA can be expressed by gijA and hijA recursively, and we also omit them here since they are much too tedious.In addition, it shall be noticed that(117)k300=−4αe−m228αe+m226561e2<0,k210=−αe−m28αe+m2264α2e2−20αem2−35m2287482αe+m22e2<0.Combining the paper [14, 31] and making a time transformation τ=−k30A/k21At in a small neighbourhood of A=0 (we still use symbol t), the new version of the system (116)(118)u˙=σAνAv,v˙=−k30AσAλ1A+λ2AνAu−νA3u3+k21Aλ3A+AAνAu+νA2u2+v2Q1u,v,A+Ou,v4,can be rewriten as(119)u˙=v,v˙=μ1A+μ2Au−u3+vμ3A+A1Au+u2+v2Q2u,v,A+Ou,v4,where coefficients are(120)λ1A=−k00Ak30A+k10Ak20A3k30A2−k20A39k30A3+k20A327k30A3,λ2A=−k10Ak30A+k20A23k30A2,λ3A=k01Ak21A−k11Ak20A3k21Ak30A+k21Ak20A29k21Ak30A2,AA=k11Ak21A+2k20A3k30A,νA=−k30Ak21A2,σA=−k30Ak21AνA,Q1u,v,A=σAk02A+k12Ak20A29k30A2+σAk03Av+νAk12Au,A1A=k21A−k30Ak30AAA,Q2u,v,A=−k21Ak30AQ1u,v,A,and a transformation is(121)μ1A=k21A3k30A−k30Aλ1A,μ2A=−k21A2k30Aλ2A,μ3A=−k21A2k30Aλ3A.With some necessary calculations, we can derive(122)A10=−8αe+7m222αe+m2<0,∂μ1,μ2,μ3∂α1,α2,α3∣A=0=−9e64α2e2−20αem2−35m22613107216α2e2+6αem1−14αem2+3m1m2−2m222αe+m212≠0.That is to say, the transformation (121) is a homeomorphism and independent in a small neighbourhood of the origin. By the results of [14, 31–33], the above system (119) is a generic three parameters family or standard family of BT singularity of codimension 3 (focus case). Hence, the system (1) in the case (C5) undergoes a degenerate focus type BT bifurcation of codimension 3 with the value of bifurcation parameters r1, K1 and d falling in a small neighbourhood of E7.Theorem 5 (Degenerate focus type Bogdanov-Takens bifurcation of codimension 3).
For the case (C5), if the variable A falls in a small neighbourhood of the origin, then the system (111) undergoes a degenerate focus type BT bifurcation of codimension 3 in a small neighbourhood of E7. In a sufficiently small neighbourhood of r1,K1,d of the bifurcation parameter space, there is a Hopf bifurcation surface, two homoclinic bifurcation surfaces, two saddle-node loop bifurcation surfaces, a multiple limit cycle bifurcation surface, and two saddle-node bifurcation surfaces for the system (111). When parameters r1,K1,d cross above surfaces, the system (111) undergoes above bifurcations, respectively.Example 6.
We take out values of parameters from the Example3, the system (111) undergoes a degenerate focus type Bogdanov-Takens bifurcation of codimension 3 in a small neighbourhood of E7 when the variable vector A=α1,α2,α3 varies in a small neighbourhood of the origin O.


## 3.1. Hopf Bifurcation Curve AroundE3

For the case (C1), the system (1) undergoes a nondegenerate Hopf bifurcation around the equilibrium point E3 as σ≠0. The Hopf bifurcation is supercritical (subcritical) and limit cycles generated by the critical point are stable (unstable) if σ<0σ>0. On occasion, there may exist some parameter values such that σ=0 or the system (1) may undergo a degenerate Hopf bifurcation for some values of parameters [17].Conversely, we also study existence of Hopf bifurcation curve in this special case. Choosingm2 and d as bifurcation parameters and introducing sufficiently small parameter λ1,λ2 . For an unfolding system(75)x˙=xr11−xαeμar1+eλ+1α−r1−αya+x,y˙=yxαea+x−λαe−λ1−−λαλμ−μ+1aλμ+12+λ2y.We suppose that the above system has an equilibrium pointE∗=x∗,y∗ when λ1,λ2≠0, where x∗=x3+w with w≪1. Calculating A1=0 and A2>0, we have a solution(76)λ1=weα2aλμ+a+2waaλμ+a+wμ,λ2=−αwaaμ−wλμ+12λμ+1a+w2λ3+λμ2+2λ2+λ−1μ+λ+1λμ+1a2+2wλ3−1/2λ2+1/2λμ2+2λ2+1/2λμ+λ+1a+w21+λ−1μλ,with(77)A2=−e2α2μ2a2aλμ+a+w2aλμ+wλλ−12μ2+2λ2−2λ−1μ+λμa3+5w15+λ2−65λ+15μ2+65λ−15μa2+8w21/2+λ−1/2μa+4w3.At this point, the Hopf bifurcation curve of the system (75) is defined by(78)Hp=λ1,λ2∣λ1,λ2satisfy3.2,and the approximation of Hp is a straight line with slope(79)k=limw⟶0λ2wλ1w=−μ2λ3−2μλ2+−μ2−μ−1λ+μ−1aλμ+12e2λμ+1.In a small neighbourhood of the origin in parameter plane.Similarly, for the bifurcation parametersm1, m2 and corresponding unfolding system(80)x˙=xr11−xαeμar1+eλ+1α−r1−λ1−αya+x,y˙=yxαea+x−λαe−λ2+λαλμ−μ+1yaλμ+12.The Hopf bifurcation curve of the system (80) is defined by(81)λ1=−wαa2λ4μ3+3aa+2/3wμ2λ3+μμ2+μ+3a2−wμ−4a+w2λ2+2μ+1a2+wμ2+μ+2a−w2μ−1λ−μ−1a−2waea2λ4μ3+3aa+2/3wμ2λ3+3a2−1/3wμ−4a+1/3w2μλ2+a2+2wa−w2μ−1λ+waaμ,λ2=wαa3λ5μ4+4a+3/4wa2μ3λ4−μ2+μ−6a2+wμ−9a−3w2aμ2λ3−2μ−2a3+wμ−3/2μ+3a2+w2μ−3a−1/2w3μλ2+a+wμ2−μ+1a2−wμ+2μ−1a−w2μ−1λ+a2μweaλμ+a+wa2λ4μ3+3aa+2/3wμ2λ3+3a2−1/3wμ−4a+1/3w2μλ2+a2+2wa−w2μ−1λ+waaμ,with A2>0. Hence the slope of approximation around the origin is(82)k=limw⟶0λ2wλ1w=−λ1+λ3−λμ3+3λ2−λ+1μ2+3λ−1μλ3+λμ2+2λ2+λ−1μ+λ+1λμ+1.From the Example1 with values of parameters r1=0.6, α=0.5, a=1.5 and e=0.6, Figure 8 show the Hopf bifurcation curves with respect to (i) λ=1/2, μ=3, (ii) λ=1/2, μ=10, (iii) λ=1/3, μ=3 and (iv) λ=1/3, μ=10, respectively. Following subsections will make further efforts to illustrate these curvesFigure 8
(a) Hopf bifurcation curve withλ=1/2 and μ=3 for the case (C2); (b) Hopf bifurcation curve with λ=1/2 and μ=10 for the case (C2); (c) Hopf bifurcation curve with λ=1/3 and μ=3 for the case (C3); (d) Hopf bifurcation curve with λ=1/3 and μ=10 for the case (C3).
(a)(b)(c)(d)


## 3.2. BT Bifurcations of Codimension 2 AroundE42

In this subsection, we firstly choosem2 and d as bifurcation parameters, and then investigate BT bifurcation of codimension 2 in the case (C2) by following the techniques and steps in [22, 23] for an unfolding system(83)x˙=xr11−xK1−m1−αya+x,y˙=yαexa+x−m2+λ1−d+λ2y,where parameters λ1 and λ2 are sufficiently small. Denoting a provisional parameter vector λ=λ1,λ2 in a small neighbourhood of the origin O, we firstly select a linear transformation (I): x=X+x42, y=Y+y42 to rewrite above system as(84)X˙=F1X,Y=∑i+j=12aijλXiYj+OX,Y3,Y˙=G1X,Y=∑i+j=02bijλXiYj+OX,Y3.Secondly, we should be able to make an affine transformation (II):u=X, v=a10λX+a01λY, then the system (84) is transformed into(85)u˙=F2u,v=v+c20λu2+c11λuv+Ou,v3,u˙=G2u,v=∑i+j=02dijλuivj+Ou,v3,without demur. Thirdly, letting a transformation (III): p=u−c11λ+d02λ/2u2, q=c20λu2−d02λuv+v, we derive a new system(86)p˙=F3p,q=q+Op,q3,(87)q˙=G3p,q=∑i+j=02fijλpiqj+Op,q3,where f02λ=0. Finally, in order to eliminate all higher order terms in the equation (86), we mechanically construct a transformation (IV): w=p, z=F3p,q, then the system (86) is transformed into(88)w˙=F4w,z=z,z˙=G4w,z=∑i+j=02hijλwizj+Ow,z3,where h00λ=f00λ, h10λ=f10λ, h01λ=f01λ and h02λ=0. Denoting functions depended on parameters λ1 and λ2:(89)d1λ=f20λ,d1=d10=4−32α2e232a,d2λ=h11λ,d2=d20=αe92−148a,ξ1λ=f00λd1λ,ξ2λ=f10λ−d1λ,ηλ=−f01λ−d1λ,Φw,λ=f00λ+f10λw+f20λw2+⋯,and noticing that d1,d2<0, then we take a time transformation (V): X1=−w, X2=z, τ=−t to make it positive and use the Malgrange preparation theorem [24] to obtain a decomposition(90)−ΦX1,λ=ξ1λ+ξ2λX1+X12ΨX1,λ,Ψ0,λ=−d1λ.Applying a transformation (VI):Y1=X1, Y2=X2/ΨX1,λ, dτ=ΨX1,λdt and a parameter dependent affine transformation (VII): x=Y1+1/2ξ2λ, y=Y2, we finally transform above system into the normal form, which still use symbol t as time variable:(91)x˙=y,y˙=μ1λ+μ2λy+x2+d2−d1xy+Ox,y3,where μ1λ=ξ1λ−1/4ξ2λ2 and μ2λ=ηλ−d2/2−d1ξ2λ. Since the Jacobian(92)∂μ1,μ2∂λ1,λ2∣λ=0=−12832−415912+2250a7/2eα2≠0,We point out that the above transformation is nonsingular, and the results in [25–27] yield that the system (91) is strongly topologically equivalent to system(93)x˙=y,y˙=μ1+μ2y+x2−xy.In other words, the system (91) will become a standard form(94)u˙=v,v˙=ϵ1+ϵ2v+u2−uv+⋯.Under a transformation:x=−d1/d22u, y=−−d1/d23v, t=−d2/−d1τ, where ϵ1=d2/−d14μ1 and ϵ2=−d2/−d1μ2. Thus, the system (83) will undergo a BT bifurcation with bifurcation parameters μ1 and μ2 when parameter λ varies in as small neighbourhood of the origin. Hence, we have following theorem.Theorem 3 (Bogdanov-Takens bifurcation of codimension 2 aroundE42).
For the case (C2), in a small neighbourhood of the equilibrium point E42, the system (83) undergoes a Bogdanov-Takens bifurcation of codimension 2 when parameter λ varies in a small neighbourhood of the origin when we choose m2 and d as bifurcation parameters. At the same time, there exist values of parameters such that this system (83) has a limit cycle or a homoclinic loop surround the cusp E42.Finally, based on the above Theorem3, the local representations of bifurcation curves up to third-order approximations in a small neighbourhood of the origin with slope k=−2+2/4ae are presented as follows, including description of saddle-node (SN), Hopf (H) and homoclinic (HL) bifurcation curves [28, 29]. The slope k can be viewed as the limiting case of the slope (79) when μ⟶μ1. These bifurcation curves can divide a small neighbourhood of the origin in the parameter plane into several regions, which can exhibit dynamics of the system (83).(i)
The saddle-node bifurcation curve is formulated by(95)SN=λ∣μ1λ=0,μ2λ≠0=λ∣λ2=2−24aeλ1+Oλ12,and(96)SN+=λ∣μ1λ=0,μ2λ>0,SN−=λ∣μ1λ=0,μ2λ<0.(ii)
The Hopf bifurcation curve is formulated by(97)H=λ|μ2λ2−d22d1μ1λ=0,μ1λ<0=λ|−1681+562aαeλ1−32274+1932a2αλ2+167132+1000aα2e2λ12+12810672+2269a2eα2λ1λ2+25634852+4927a3α2λ22−643139+22132ae3α3λ13−128317532+44932a2α3e2λ12λ2−2561507062+213155a3eα3λ1λ22−10241176092+166328a4α3λ23+Oλ1,λ24=0,μ1λ<0.(iii)
The homoclinic bifurcation curve is formulated by(98)HL=λ∣μ2λ2−25d2249d1μ1λ=0,μ1λ<0=λ∣−40081+562a49αeλ1−800274+1932a249αλ2+16178972+25096a49α2e2λ12+128408592+57685a249α2eλ1λ2+256937672+132559a349α2λ22−6478787+555412a49e3α3λ13−1288075292+1142740a249α3e2λ12λ2−25639950742+5650643a349eα3λ1λ22−102433584012+4749632a449α3λ23+Oλ1,λ24=0,μ1λ<0.Example 4.
The values of parameters will be recalled from the Example2. Figure 9 gives the saddle-node, Hopf and homoclinic bifurcation curves when the value of parameter λ falls in a small neighbourhood of the origin in the parameter plane. Figures 10–13 depict the dynamics of the system (83) in a small neighbourhood of the origin in the parameter plane when we set values of λ1 and λ2.(i)
When the value ofλ lies on the saddle-node bifurcation curve SN−, there exist two interior equilibria and one is unstable node.(ii)
When the value ofλ lies on region I (the region below the saddle-node bifurcation curve), there exists a unique unstable node. The Poincre-Bendixson theorem implies that there is a large limit cycle enclosing this equilibrium point, the detailed numerical simulation results are shown in Figure 10(a) with λ1,λ2≈1×10−4,−3.255138×10−5.(iii)
When the value ofλ lies on another saddle-node bifurcation curve SN+, there exists a unique unstable node and a limit cycle.(iv)
When the value ofλ lies on region IV (the region between the saddle-node bifurcation curve SN+ and the homoclinic curve), there exist an unstable node, a saddle, a stable focus, and the homoclinic loop in case (v) is broken, in which the saddle and the focus are bifurcated from the curve SN+, a large limit cycle and enlarged phase diagrams around a stable focus can be seen in Figure 11 with λ1,λ2≈1×10−4,−8.131779×10−6.(v)
When the value ofλ lies on the homoclinic curve, there exist three interior equilibrium points and a homoclinic loop enclosing a stable focus, the detailed numerical simulation results are shown in Figure 12(a) with λ1,λ2≈1×10−4,−1.626356×10−5 and enlarged phase diagrams can be seen in Figure 12(b).(vi)
When the value ofλ crosses the homoclinic curve into region III (the region between the homoclinic curve and the Hopf curve), there exist an unstable node, a saddle and a stable focus, the detailed numerical simulation results are shown in Figure 13 with λ1,λ2≈1×10−4,−1.626653×10−5.(vii)
When the value ofλ lies on the Hopf curve, there exist an unstable node, a saddle, an non-hyperbolic equilibrium point (multiple focus or center) according to Hopf bifurcation and a large limit cycle enclosing these equilibria.(viii)
When the value ofλ crosses the Hopf curve into region II (the region between the Hopf curve and the saddle-node bifurcation curve SN−), there exist an unstable node, a saddle and an unstable focus, the detailed numerical simulation results are shown in Figure 10(b) with λ1,λ2≈1×10−4,−1.627259×10−5. At the same time, it can be noticed that the stable focus in case (vi) can become unstable.Figure 9
Saddle-node (in red), hopf (in green) and homoclinic (in blue) bifurcation curves in case (C2).Figure 10
Phase diagrams in cases (ii) and (viii).
(a)(b)Figure 11
Phase diagrams in the case (iv).
(a)(b)Figure 12
Phase diagrams in case (v).
(a)(b)Figure 13
Phase diagrams in case (vi).
(a)(b)


## 3.3. BT Bifurcation of Codimension 2 AroundE52

In this subsection, we set out to deal with BT bifurcation of codimension 2 around the equilibrium pointE52 in the case (C3). Here we choose m1 and m2 as bifurcation parameters and rewrite the system (1) as the unfolding form:(99)x˙=xr11−xK1+m1+λ1−αya+x,y˙=yαexa+x−m2+λ2−dy,where λ1 and λ2 are sufficiently small parameters and λ=λ1,λ2. Firstly, we take a mere linear transformation (I): x=X+x52,y=Y+y52, the equilibrium point E52 is translated to the origin O:(100)X˙=F1X,Y=∑i+j=02aijλXiYj+OX,Y3,Y˙=G1X,Y=∑i+j=02bijλXiYj+OX,Y3.Secondly, we take a transformation (II):u=X, v=F1X,Y, then the above system becomes(101)u˙=F2u,v=v,v˙=G2u,v=∑i+j=02dijλuivj+Ou,v3.Thirdly, we take a transformation (III):p=u+d01λ/d11λ, q=v since d110,0=−85+1917αe/36a≠0, then we have following system with f01λ=0:(102)p˙=F3p,q=q,q˙=G3p,q=∑i+j=02fijλpiqj+Op,q3.Letting (IV):w=p, z=1−f02λpq, dt=1−f02λpdτ and rewrite symbol τ as t, we derive a new system(103)w˙=F4w,z=z,z˙=G4w,z=∑i+j=02hijλwizj+Ow,z3.It can be noticed thath200,0=33−717α2e2/432a>0 and h110,0=d110,0≠0, so h20λ>0 and h11λ≠0 when λ changes in a small neighbourhood of the origin, which depends smoothly on λ. Finally, we construct a time transformation (V): m=h11λ2/h20λw, n=h11λ3/h20λ2z, dt=h11λ/h20λdτ and rewrite symbol τ as t, then the above system (103) becomes(104)m˙=F5m,n=n,n˙=G5m,n=μ1λ+μ2λm+m2+mn+Om,n3,where(105)μ1=μ1λ=5781−17αeλ1+8679−17αeλ2+Oλ1,λ22,μ2=μ2λ=2−51+1117αeλ1+6917−85αeλ2+Oλ1,λ22.Therefore, owing to the Jacobian ofμ1 and μ2(106)∂μ1,μ2∂λ1,λ2|λ=0=346885+1917α2e2≠0,or(107)−1∂μ1/∂λ2⋅∂μ1,μ2∂λ1,λ2|λ=0<0.The system (99) is a generic family unfolding at the codimension 2 cusp E52, and we have following local representations of the bifurcation curves up to second-order approximations with slope k=1+17/12 for the system (104) [30]. Thus, it can be noticed that this slope k can be viewed as the limiting case of the slope (82) when μ⟶μ1.(i)
The saddle-node bifurcation curve is formulated by(108)SN=λ|μ1=14μ22=λ|5781−17αeλ1+8679−17αeλ2+179775−63172α2e2λ12+2551003−49917α2e2λ1λ2+3063111−62317α2e2λ22+Oλ1,λ23=0.(ii)
The Hopf bifurcation curve is formulated by(109)H=λ|μ1=0,μ2<0=λ|5781−17αeλ1+8679−17αeλ2+513441−65172α2e2λ12+515723−265917α2e2λ1λ2+1224841−16717α2e2λ22+Oλ1,λ23=0.(iii)
The homoclinic bifurcation curve is formulated by(110)HL=λ∣μ1=−625μ22,μ2<0=λ|5781−17αeλ1+8679−17αeλ2+5190409−26811750α2e2λ12+51160067−704111725α2e2λ1λ2+122422543−44451725α2e2λ22+Oλ1,λ23=0.Theorem 4.
(Bogdanov-Takens bifurcation of codimension 2 aroundE52). For the case (C3), in a small neighbourhood of the equilibrium point E52, there exist values of parameters such that the system (99) undergoes an attracting Bogdanov-Takens bifurcation of codimension 2 when parameter λ varies in a small neighbourhood of the origin with bifurcation parameters m1 and m2. And this system is a generic family unfolding at the cusp E52 of codimension 2.Example 5.
Here we extract values of parameters in the Subsection3.2 and the Example 3. Figure 14 presents the saddle-node, Hopf and homoclinic bifurcation curves when the value of parameter λ falls in a small neighbourhood of the origin in the parameter plane, which can divide a small neighbourhood of the origin in the parameter plane into several regions and exhibit different dynamical behavior of the system (99).(i)
When the value ofλ crosses the λ1 axis into region I (the region between the saddle-node bifurcation curve SN2 and the homoclinic bifurcation curve), there exist an unstable focus, a saddle and a stable node, in which the focus and saddle are bifurcated from the curve SN2.(ii)
When the value ofλ lies on the homoclinic bifurcation curve, there exist a saddle, a stable node and a homoclinic loop enclosing an unstable focus.(iii)
When the value ofλ crosses the homoclinic bifurcation curve into region II (the region between the homoclinic bifurcation curve and the Hopf bifurcation curve), there exist an unstable focus, a saddle and a stable node.(iv)
When the value ofλ lies on the Hopf bifurcation curve, there exist three interior equilibria including a saddle and a stable node.(v)
When the value ofλ crosses the Hopf bifurcation curve into region III (the region between the Hopf bifurcation curve and the saddle-node bifurcation curve SN1), there exist a stable focus, a saddle and a stable node. Furthermore, ti shall be noticed that the unstable focus in case (iii) becomes stable, which can ensure potential Hopf bifurcation.(vi)
When the value ofλ lies on the saddle-node bifurcation curve SN1, there exist two interior equilibria including a stable node.(vii)
When the value ofλ crosses the saddle-node bifurcation curve into region IV (the region above the saddle-node bifurcation curve), there exists a unique stable node.(viii)
When the value ofλ lies on the saddle-node bifurcation curve SN2, the unique equilibrium point is a stable node.Figure 14
Saddle-node (in red), hopf (in green) and homoclinic (in blue) bifurcation curves in the case (C3).


## 3.4. Degenerate Focus Type BT Bifurcation of Codimension 3

We should always keep the case (C5) in mind and choose r1, K1 and d as bifurcation parameters and use the method in [14]. Thus we introduce sufficiently small variables α1, α2, α3, a parameter vector A=α1,α2,α3, the origin O=0,0,0 in parameter space and an unfolding form system(111)x˙=xr1+α11−xK1+α2−αxya+x−m1x,y˙=αexya+x−m2y−d+α3y2.Of the system (1) since the Theorem 2 (ii). As indicated in Subsection 2.5 for the system (111), by using the coordinate transformations (I), (II), (III) and (IV), actual calculations can yield a new system(112)w˙=F4w,z=∑i+j=03eijAwizj+Ow,z4,(113)z˙=G4w,z=∑i+j=03fijAwizj+Ow,z4.Here we omit complicated, but consider inconsequential expressions of smooth functionseijA and fijA for the sake of convenience, while eij0=eij,fij0=fij, for instance, e000=e100=e200=e110=e020=0 and f000=f100=f010=f200=f020=0.Secondly, in order to eliminate third-order terms in right hand side of equation (112) when A=0, we construct a transformation(114)V:w=x1+e03x1y12,z=−e12x1y12−e21x12y1−e30x13+y1.Above system becomes(115)x1.=F5x1,y1=y1+∑i+j=03gijAx1iy1j+Ox1,y14,y1.=G5x1,y1=∑i+j=03hijAx1iy1j+Ox1,y14.Similarly, coefficientsgijA, hijA are all omitted for space and can be deduced by coefficients eijA and fijA. Moreover, it shall be noticed that gij0=0i+j≤3. Finally, the transformation (VI): x2=x1, y2=F5x1,y1 can translate the above system into the following system(116)x2.=F6x2,y2=y2,y2.=G6x2,y2=∑i+j=03kijAx2iy2j+Ox2,y24,where coefficients kijA can be expressed by gijA and hijA recursively, and we also omit them here since they are much too tedious.In addition, it shall be noticed that(117)k300=−4αe−m228αe+m226561e2<0,k210=−αe−m28αe+m2264α2e2−20αem2−35m2287482αe+m22e2<0.Combining the paper [14, 31] and making a time transformation τ=−k30A/k21At in a small neighbourhood of A=0 (we still use symbol t), the new version of the system (116)(118)u˙=σAνAv,v˙=−k30AσAλ1A+λ2AνAu−νA3u3+k21Aλ3A+AAνAu+νA2u2+v2Q1u,v,A+Ou,v4,can be rewriten as(119)u˙=v,v˙=μ1A+μ2Au−u3+vμ3A+A1Au+u2+v2Q2u,v,A+Ou,v4,where coefficients are(120)λ1A=−k00Ak30A+k10Ak20A3k30A2−k20A39k30A3+k20A327k30A3,λ2A=−k10Ak30A+k20A23k30A2,λ3A=k01Ak21A−k11Ak20A3k21Ak30A+k21Ak20A29k21Ak30A2,AA=k11Ak21A+2k20A3k30A,νA=−k30Ak21A2,σA=−k30Ak21AνA,Q1u,v,A=σAk02A+k12Ak20A29k30A2+σAk03Av+νAk12Au,A1A=k21A−k30Ak30AAA,Q2u,v,A=−k21Ak30AQ1u,v,A,and a transformation is(121)μ1A=k21A3k30A−k30Aλ1A,μ2A=−k21A2k30Aλ2A,μ3A=−k21A2k30Aλ3A.With some necessary calculations, we can derive(122)A10=−8αe+7m222αe+m2<0,∂μ1,μ2,μ3∂α1,α2,α3∣A=0=−9e64α2e2−20αem2−35m22613107216α2e2+6αem1−14αem2+3m1m2−2m222αe+m212≠0.That is to say, the transformation (121) is a homeomorphism and independent in a small neighbourhood of the origin. By the results of [14, 31–33], the above system (119) is a generic three parameters family or standard family of BT singularity of codimension 3 (focus case). Hence, the system (1) in the case (C5) undergoes a degenerate focus type BT bifurcation of codimension 3 with the value of bifurcation parameters r1, K1 and d falling in a small neighbourhood of E7.Theorem 5 (Degenerate focus type Bogdanov-Takens bifurcation of codimension 3).
For the case (C5), if the variable A falls in a small neighbourhood of the origin, then the system (111) undergoes a degenerate focus type BT bifurcation of codimension 3 in a small neighbourhood of E7. In a sufficiently small neighbourhood of r1,K1,d of the bifurcation parameter space, there is a Hopf bifurcation surface, two homoclinic bifurcation surfaces, two saddle-node loop bifurcation surfaces, a multiple limit cycle bifurcation surface, and two saddle-node bifurcation surfaces for the system (111). When parameters r1,K1,d cross above surfaces, the system (111) undergoes above bifurcations, respectively.Example 6.
We take out values of parameters from the Example3, the system (111) undergoes a degenerate focus type Bogdanov-Takens bifurcation of codimension 3 in a small neighbourhood of E7 when the variable vector A=α1,α2,α3 varies in a small neighbourhood of the origin O.


## 4. Summary

In summary, the main object of writing this paper is to consider Hopf and BT bifurcations of codimension 2 and 3 in the Bazykin’s predator-prey system. The system in our paper has complicated and rich dynamical behaviors. The particular cases (C1)–(C5) based on polynomial equations and Jacobian matrix at point E∗ are analytical obtained and investigated in detail. In the case (C1), for the equilibrium point E3, with the normal form and the Lyapunov number at hand, we analyze the stability, Hopf bifurcation with standard bifurcation theory and existence of Hopf bifurcation curve with bifurcation parameters m1 and d. In the cases (C2) and (C3), the equilibria E42 and E52 are both cusp of codimension 2, and we analytically give formulae of saddle-node, Hopf and homoclinic bifurcation curves, respectively. These bifurcation curves can divide a small neighbourhood of the origin in the parameter plane into several regions, which can exhibit dynamical behaviors of corresponding unfolding system, respectively. In the case (C5), the nilpotent equilibrium point E7 is a codimension 3 BT singularity (focus or center), which can exhibit a degenerate focus type BT bifurcation of codimension 3 in a small neighbourhood. For the system (1) or a general ODEs system, some open problems are:(i)
Whether the BT bifurcations in Subsection3.2 and Subsection 3.3 are limiting cases of the Hopf bifurcations in Subsection 3.1 when μ⟶μ1?(ii)
The existence and uniqueness of a codimensionN cusp, and the corresponding bifurcation;(iii)
Whether a codimensionN+1 cusp (bifurcation) is the threshold or limiting case of some codimension N cusps (bifurcations) or not?(iv)
The classification work of a codimensionN cusp, for instance, topologically equivalent systems, diffeomorphic systems, limit cycles or homoclinic loops;(v)
The existence of invariant quantities under some nonsingular transformations, for instance, the Lyapunov quantities, the symbolled “index”sgnd1d2=±1or0 within codimension 2 cusps. Correspondingly, a codimension 2 and 3 cusps respectively signifies sgnd1d2=±1 and sgnd1d2=0.Besides to the Theorem3, Theorem 4 and Theorem 5, in a sense, we point out that the maximum number of limit cycles in our system is an important research content, particularly the versal unfolding of a focus type BT singularity of codimension 3. These are reflected in the open problems in [34] or the Hilbert’s 16th problem. Up to now, the Poincare-Bendixson’s existence theorem, the Zhifen Zhang’s uniqueness theorem and the Wintner-Perko termination principle to determine at most two limit cycles surrounding a singular point in [35, 36] are functional tools. Aparting from these theorems, the Theorem 5.4 (ii) in [37] established a result that there exists a unique stable limit cycle in the first quadrant if required condition holds; the Theorem 4, Corollary 1 and Corollary 2 in [38] proved that the system has at least one stable limit cycle and one unstable limit cycle.Based on above open problem (iv), the paper [10] mainly gave occurrence of Bogdanov-Takens bifurcations and modified the approximate calculation of limit cycles via a perturbation procedure and canonical transformation in view of supercritical Hopf bifurcation. It further illustrate the BT singularity (focus or center) E7 of codimension 3, which also shows its uniqueness and nonexistence of codimension N≥4 cusps. Finally, it is our expectancy that the qualitative analysis of stability and bifurcations can be suitable to more predator-prey systems to reveal these phenomena, which are much more complex and have richer dynamical behaviors, even physical mechanical systems and epidemic models, etc.In the follow-up research works, we will first deepen theoretical research of bifurcation dynamics by learning from relevant results in these papers [39–41], and then further explore the dynamic behavior of patch pattern in ecosystem in the help of these papers [42–44], finally put mathematical models in specific ecological and environmental problems to study their ecological significance by means of these papers [45–47]. In a word, all these results are expected to be useful in studying dynamic behavior of the ecosystem.



---

*Source: 4844228-2022-09-22.xml*