Bogdanov–Takens bifurcation of a Holling IV prey–predator model with constant-effort harvesting

Abstract

A prey–predator model with constant-effort harvesting on the prey and predators is investigated in this paper. First, we discuss the number and type of the equilibria by analyzing the equations of equilibria and the distribution of eigenvalues. Second, with the rescaled harvesting efforts as bifurcation parameters, a subcritical Hopf bifurcation is exhibited near the multiple focus and a Bogdanov–Takens bifurcation is also displayed near the BT singularity by analyzing the versal unfolding of the model. With the variation of bifurcation parameters, the system shows multi-stable structure, and the attractive domains for different attractors are constituted by the stable and unstable manifolds of saddles and the limit cycles bifurcated from Hopf and Bogdanov–Takens bifurcations. Finally, a cusp point and two generalized Hopf points are found on the saddle-node bifurcation curve and the Hopf bifurcation curves, respectively. Several phase diagrams for parameters near one of the generalized Hopf points are exhibited through the generalized Hopf bifurcation.

1 Introduction

The prey–predator model based on Lotka–Volterra model is one of the most popular models in mathematical ecology and has been widely applied in understanding population dynamics of the species, which is characterized by the complicated interaction among the species and the interaction between the species and their surroundings. Functional response is the rate of prey consumption by the predators, which can characterize the interaction between the prey and the predators. The initial functional responses are presented by Holling in [11, 12] and have been classified into three types called Holling I, II and III. Generally speaking, Holling type I response is seen mostly at filter-feeding predators, such as mollusks, algae, and cells. Holling types II and III responses are applicable to the invertebrates and complicated vertebrates, respectively. The common feature of these functions is all monotonic increasing and bounded, which means that the more the prey is in the environment, the better for the predators [6]. However, many experiments and observations indicate that monotonicity does not always hold and a non-monotonic response may exist when the nutrient concentration reaches a high level [3, 9, 31]. It is also explained that when the number of the prey is large enough, their group defense ability and camouflage ability are strengthened to decrease the predation ability of the predators. For example, large swarms of insects make individual identification difficult for their predators. A lone musk ox can be successfully attacked by wolves, while large groups of musk oxen are less likely to be attacked successfully [28]. Filamentous algae are often qualified as inedible by herbivorous zooplankton [8]. To model such an inhibitory effect, a Holling IV response $\frac{mx}{a+bx+x^{2}}$ and a simplified Holling IV response $\frac{mx}{a+x^{2}}$ have been proposed by Andrews [1] and Sokoland [25], respectively. Since then the prey–predator system with Holling IV functional response has been widely studied by many researchers [4, 15, 20, 30, 32].

For commercial and economical purposes, the exploitation of biological resources and harvesting of populations are commonly practiced in fishery, forestry and wildlife management [22]. It is known that harvesting, as a scientific management of renewable resources, has a strong effect on the dynamic evolution of the biological models, which can determine the survival and sustainable development of the species. In order to reasonably govern renewable resources, the study of the exploitation of biological resources and harvesting of populations is becoming a hot topic in the field of biological economics, which is related to the optimal management of renewable resources [5].

Recently, it has become very popular among many researchers in using bifurcation theory to analyze the harvested prey–predator model to reach an optimal management in the field of mathematical biology [10, 13, 15–18, 21, 29, 33, 34], which can predict the species’ evolutionary direction, including the constant state, periodic state and chaotic state, as the critical biological parameters vary.

From the present work, we find that most of the harvested prey–predator models which have been considered by many authors are confined to three aspects: (i) constant-yield harvesting, (ii) harvesting subjected to only one species or (iii) the numerical response is proportional to the functional response, where the numerical response is the change in predator density as a function of change in prey density [26]. It is the reproduction rate of a consumer as a function of food density. That is, when the predators consume their prey, the numerical response is dynamic relation between the increasing number of the predators and the consumed prey. From the convention, we know that harvesting frequently varies with season, market demand, the species quantity, harvesting cost and so on; the prey and the predators can be both valuable from the business point; the numerical response can be uncorrelated with the functional response. Based on above analysis and the virtue of the Holling IV response, in this paper we extend previous work from three aspects and establish a prey–predator model with Leslie–Gower and Holling IV schemes with constant-effort harvesting, which is given by

$$ \begin{aligned}& \dot{x}=r_{1}x \biggl(1- \frac{x}{K} \biggr)-\frac{mx}{b+x^{2}}y-c_{1}x, \\ &\dot{y}=r_{2}y \biggl(1-\frac{y}{sx} \biggr)-c_{2}y. \end{aligned} $$

(1.1)

The logistic model $r_{1}x(1-\frac{x}{K})$ is used to describe the growth of the prey in the absence of predators in the limited environmental carrying capacity. The second term $\frac{mx}{b+x^{2}}y$ describes the change in the density of the prey attacked per unit time as the prey density changes. The last term $c_{1}x$ is the constant-effort harvesting of the prey. $r_{2}y (1-\frac{y}{sx} )$ is the predator growth term, which is described by the Leslie–Gower function. $\frac{y}{sx}$ measures the loss in the predator population due to rarity (per capita of $\frac{y}{x}$) of its favorite term [19]. The term $c_{2}y$ is the constant-effort harvesting of the predators. The meanings of the parameters are given as follows:

(i) x and y denote the densities of the prey and the predators.

(ii) K is the environmental carrying capacity for the prey.

(iii) m denotes the maximal predation rate.

(iv) b is the so-called half-saturation constant.

(v) s is a measure of food quality that the prey provides for conversion into the predator birth.

(vi) $r_{1}$ and $r_{2}$ are intrinsic growth rates of the prey and predators, respectively.

(vii) $c_{1}$ and $c_{2}$ measure the effort being spent by a harvesting agency (harvesting efforts).

The model (1.1) without harvesting (i.e., $c_{1}=c_{2}=0$) has been discussed in [7, 14, 20]. Li and Xiao [20] proved that the model could simultaneously undergo a Bogdanov–Takens bifurcation and a subcritical Hopf bifurcation in the small neighborhoods of two different equilibria, respectively. It was shown that for different parameters the model could have a stable limit cycle enclosing two equilibria, or an unstable limit cycle enclosing a hyperbolic equilibria, or two limit cycles enclosing a hyperbolic equilibrium. In Ref. [14], for the same model, Huang et al. investigated the degenerate Bogdanov–Takens bifurcation of codimension 3 around the degenerate positive equilibrium. Some other abundant dynamical phenomena could emerge, such as the coexistence of three hyperbolic positive equilibria, a homoclinic loop enclosing an unstable limit cycle, or a stable limit cycle enclosing three unstable hyperbolic positive equilibria. From the above-mentioned work, Hopf cyclicity and the global dynamics of the same model were investigated by Dai and Zhao [7] for the cases of one non-degenerate positive equilibrium and three distinct positive equilibria, respectively. Some explicit conditions for the globally stable equilibrium were established by applying Dulac’s criterion and constructing the Lyapunov function. It was shown that the Hopf bifurcation could occur at each weak focus and more complicated and new dynamics were observed. When the system had a unique positive equilibrium, there existed parameter values such that the system had two limit cycles around it. When the system had three positive equilibria, one limit cycle could bifurcate from each of the two positive anti-saddles simultaneously.

In order to discuss if the added harvesting terms can affect the bifurcation structure of the prey–predator model, we also do a series of bifurcation analysis for the model with liner harvesting. We find that, besides some original dynamical behaviors, some new phenomena and bifurcations appear in the harvested model, such as a globally stable equilibrium, a cusp bifurcation and a generalized Hopf bifurcation. From the cusp point, two saddle-node bifurcation curves are bifurcated, through which three equilibria collide to form an equilibrium or one equilibrium splits into three equilibria. From the generalized Hopf bifurcation point, a fold bifurcation curve of limit cycles is bifurcated, on which two limit cycles with different stability collide to form a semi-stable cycle. These dynamical behaviors have not been found or discussed in the original model in Ref. [7, 14, 20].

Before entering the topic, we carry out the following scaling transformations:

$$\begin{aligned} &\bar{t}=r_{1}t,\qquad \bar{x}=\frac{x}{K}, \qquad\bar{y}= \frac{my}{r_{1}K^{2}},\qquad a= \frac{b}{K^{2}}, \qquad\delta = \frac{r_{2}}{r_{1}}, \qquad \beta = \frac{r_{2}K}{sm}, \\ & h_{1}= \frac{c_{1}}{r_{1}},\qquad h_{2}= \frac{c_{2}}{r_{2}}. \end{aligned}$$

Dropping the bars system (1.1) becomes

$$ \begin{aligned} &\dot{x}=x(1-x)-\frac{x}{a+x^{2}}y-h_{1}x, \\ &\dot{y}=y \biggl(\delta -\beta \frac{y}{x} \biggr)-h_{2}y, \end{aligned} $$

(1.2)

where $a,\delta, \beta, h_{1}$ and $h_{2}$ are positive parameters. $h_{1}$ and $h_{2}$ are called rescaled harvesting efforts. From the perspective of biology, system (1.2) is defined on the set $D= \{ (x,y) \in R^{2} {|} x>0, y\geq 0 \}$.

By analyzing (1.2), we can prove that solution trajectories starting from the initial value $(x_{0},0)$ remain within the positive x-axis for all the time. Through integrating two equations in (1.2) and discussing the initial value $x(0)>1-h_{1}$ or $x(0)<1-h_{1}$, we can get the following lemma (the detailed proof is similar to Appendix A in [24]).

Lemma 1.1

When $h_{1}<1$ and $h_{2}<\delta$, each positive solution $(x(t),y(t))$ of system (1.2) is bounded. Furthermore, there exists a $T\geq 0$ such that $0< x(t)<1-h_{1}$ and $0\leq y<\frac{\delta -h_{2}}{\beta }(1-h_{1})$ for $t\geq T$.

The rest of this paper is organized as follows. In Sect. 2, we analyze the number and property of the equilibria by using the root formula of the cubic equation. In Sect. 3, we discuss the subcritical Hopf bifurcation near the multiple focus and the Bogdanov–Takens bifurcation near a BT singularity. The characteristics of the cusp point and the generalized Hopf points are analyzed simply. Numerical simulations support our results of theoretical analysis. Finally, we end this paper with a conclusion in Sect. 4.

2 Types of equilibria

In this section, for system (1.2) we will study the number and type of equilibria in the region D. It is clear that if $h_{1}<1$ system (1.2) has a boundary equilibrium $E_{0}=(1-h_{1},0)$, which is a sink for $\delta < h_{2}$ and a saddle for $\delta >h_{2}$. The corresponding phase portraits are displayed in Fig. 1.

Figure 1

The phase portraits for system (1.2). (a) $E_{0}$ is a sink if $a=0.2, \beta =0.8,\delta =0.3, h_{1}=0.4$ and $h_{2}=0.6$; (b) $E_{0}$ is a saddle if $a=0.2, \beta =0.8,\delta =0.9, h_{1}=0.4$ and $h_{2}=0.6$

Obviously, a positive equilibrium $E^{*}(x^{*},y^{*})$ of system (1.2) should satisfy the following equations:

$$ \begin{aligned} &1-h_{1}-x^{*}- \frac{y^{*}}{a+{x^{*}}^{2}}=0, \\ &y^{*}-\frac{\delta -h_{2}}{\beta } x^{*}=0. \end{aligned} $$

(2.1)

From Eqs. (2.1), we find that, if positive equilibria exist, then $h_{2}<\delta$ and $x^{*}$ is a root of the equation

$$ F(x)={x}^{3}+ (h_{1}-1) {x}^{2}+ \biggl( a+{\frac{\delta -h_{2}}{\beta }} \biggr) x+a(h_{1}-1)=0. $$

(2.2)

To determine the type of $E^{*}(x^{*},y^{*})$, the Jacobian matrix J evaluated at $E^{*}(x^{*},y^{*})$ is given by

$$J\left(E^{\ast }\right)=\left(\begin{array}{cc}1-h_1-2x^{\ast }+\frac{(\delta -h_2)({x^{\ast }}^3-a{x}^{\ast })}{\beta {(a+{x^{\ast }}^2)}^2}& \frac{-x^{\ast }}{a+{x^{\ast }}^2}\\ \frac{{(\delta -h_2)}^2}{\beta }& -(\delta -h_2)\end{array}\right).$$

Thus

$$ \begin{aligned} &\operatorname{Det} \bigl(J \bigl(E^{*} \bigr) \bigr)=(h_{2}-\delta ) \bigl(1-h_{1}-2x^{*} \bigr) + \frac{2a(\delta -h_{2})^{2} x^{*}}{\beta (a+{x^{*}}^{2} )^{2}}=(\delta -h_{2})\frac{x^{*}}{a+{x^{*}}^{2}}F' \bigl(x^{*} \bigr), \\ &\operatorname{Tr} \bigl(J \bigl(E^{*} \bigr) \bigr)=1-h_{1}+h_{2}- \delta -2x^{*}+ \frac{(\delta -h_{2}) ({x^{*}}^{3}-ax^{*} )}{\beta (a+{x^{*}}^{2} )^{2}} \\ &\phantom{\operatorname{Tr} (J (E^{*} ) )}=\frac{x^{*}}{a+{x^{*}}^{2}} \biggl( \frac{\delta -h_{2}}{\beta }-F' \bigl(x^{*} \bigr) \biggr)-(\delta -h_{2}). \end{aligned} $$

According to the sign of $\operatorname{Det}(J(E^{*}))$, the equilibrium $E^{*}(x^{*},y^{*})$ can be divided into three types: an anti-saddle equilibrium, a hyperbolic saddle or a degenerated equilibrium if $\operatorname{Det}(J(E^{*}))>0, \operatorname{Det}(J(E^{*}))<0$ or $\operatorname{Det}(J(E^{*}))=0$, respectively. Obviously, the type of $E^{*}(x^{*},y^{*})$ is determined by the sign of $F'(x^{*})$.

Using the root formula of the third-order equation and estimating the sign of $F'(x^{*})$, we can get the following results.

Lemma 2.1

Suppose $h_{1} < 1$, $h_{2} < \delta$, $A=(h_{1}-1)^{2}-3(a+\frac{\delta -h_{2}}{\beta })$ and $\Delta = -4A^{3} + (h_{1}-1)^{2} [3A+27a - (h_{1} - 1)^{2} ]^{2}$.

• (a) If $\Delta >0$, then system (1.2) has a unique positive equilibrium $E^{*}(x^{*},y^{*})$, which is an anti-saddle equilibrium whatever the sign of A is.

• (b) If $\Delta =0, A=0$, then system (1.2) has a unique positive equilibrium $E^{*}(x^{*},y^{*})= (\frac{1-h_{1}}{3}, \frac{(\delta -h_{2})(1-h_{1})}{3\beta } )$, which is degenerated.

• (c) If $\Delta =0, A>0$ and $\operatorname{max} \{ \frac{(h_{1}-1)^{2}-4A}{27},0 \} < a< \frac{(h_{1}-1)^{2}-A}{27}$, then the system (1.2) has two positive equilibria: an anti-saddle equilibrium

  $$ E_{1}^{*} \bigl(x_{1}^{*},y_{1}^{*} \bigr)= \biggl( \frac{(1-h_{1}) [4A+27a-(h_{1}-1)^{2} ]}{3A}, \frac{(\delta -h_{2})(1-h_{1}) [4A+27a-(h_{1}-1)^{2} ]}{3\beta A} \biggr) $$

  and a degenerated equilibrium

  $$ E^{*} \bigl(x^{*},y^{*} \bigr)= \biggl( \frac{(1-h_{1}) [(h_{1}-1)^{2}-A-27a ]}{6A}, \frac{(\delta -h_{2})(1-h_{1}) [(h_{1}-1)^{2}-A-27a ]}{6\beta A} \biggr). $$

• (d) If $\Delta <0$, then system (1.2) has three positive equilibria $E_{1}^{*}(x_{1}^{*},y_{1}^{*}),E_{2}^{*}(x_{2}^{*},y_{2}^{*})$ and $E_{3}^{*}(x_{3}^{*},y_{3}^{*})$. $E_{1}^{*}$ and $E_{3}^{*}$ are anti-saddle equilibria and $E_{2}^{*}$ is a hyperbolic saddle.

The detailed proof of Lemma 2.1 is given in the appendix. The phase portraits for the four cases are exhibited in Fig. 2.

Figure 2

(a) An anti-saddle $E^{*}=(0.116,0.078)$ with $a=0.1, \beta =0.3, \delta =0.6, h_{1}=0.2, h_{2}=0.4, \Delta =23.75$ and $A=-1.66$; (b) a degenerated equilibrium $E^{*}=(0.299,0.072)$ with $a=0.03$, $\beta =1.25$, $\delta =0.8$, $h_{1}=0.1$, $h_{2}=0.5$, $\Delta =0$ and $A=0$; (c) a degenerated equilibrium $E^{*}=(0.21,0.047)$ and an anti-saddle $E_{1}^{*}=(0.48,0.11)$ with $a=0.024$, $\beta =0.9$, $\delta =0.5$, $h_{1}=0.1$, $h_{2}=0.3$, $\Delta =0$ and $A=0.073$; (d) two anti-saddles $E_{1}^{*}=(0.06,0.012)$ and $E_{3}^{*}=(0.599,0.11)$, a hyperbolic saddle $E_{2}^{*}=(0.24, 0.04)$ with $a=0.01, \beta =0.54, \delta =0.6, h_{1}=0.1, h_{2}=0.5$ and $\Delta =-0.031$

Next, we discuss the case (c) of Lemma 2.1 and look for some parameter values such that system (1.2) has a degenerated equilibrium $E^{*}(x^{*},y^{*})$ with $\operatorname{Det}(J(E^{*}))=0$ and $\operatorname{Tr}(J(E^{*}))=0$ and a non-hyperbolic equilibrium $E_{1}^{*}(x_{1}^{*},y_{1}^{*})$ with $\operatorname{Det}(J(E_{1}^{*}))>0$ and $\operatorname{Tr}(J(E_{1}^{*}))=0$.

From $\operatorname{Det}(J(E^{*}))=0$ and $\operatorname{Tr}(J(E^{*}))=0$, we can get

$$ \begin{aligned} &x^{*}=1-h_{1}- \delta +h_{2}, \qquad y^{*}=\frac{\delta -h_{2}}{\beta }x^{*}, \\ & a= \frac{(1-h_{1}-\delta +h_{2})^{2}[2(\delta -h_{2})+h_{1}-1]}{1-h_{1}},\qquad\beta =\frac{1-h_{1}}{2(\delta -h_{2})(1-h_{1}-\delta +h_{2})}, \end{aligned} $$

(2.3)

where $(\delta -h_{2})<1-h_{1}<2(\delta -h_{2})$.

Thus when $a,\beta$ satisfy (2.3) and $3(\delta -h_{2})\neq 2(1-h_{1})$, system (1.2) has two positive equilibria

$$\begin{aligned} &E_{1}^{*}= \biggl(2(\delta -h_{2})+h_{1}-1, \frac{2(\delta -h_{2})^{2}(2\delta -2h_{2}+h_{1}-1)(1-h_{1}-\delta +h_{2})}{1-h_{1}} \biggr) \quad\text{{and}} \\ & E^{*}= \biggl(1-h_{1}- \delta +h_{2}, \frac{2(\delta -h_{2})^{2}(1-h_{1}-\delta +h_{2})^{2}}{1-h_{1}} \biggr). \end{aligned}$$

Furthermore, when $E_{1}^{*}$ satisfies $\operatorname{Det}(J(E_{1}^{*}))>0$ and $\operatorname{Tr}(J(E_{1}^{*}))=0$, the following conclusions hold.

Theorem 2.2

If $0< h_{1}<1$ and $(a,\beta,\delta )= (\frac{41\sqrt{17}-169}{2}(1-h_{1})^{2}, \frac{4+\sqrt{17}}{4(1-h_{1})}, \frac{\sqrt{17}-3}{2}(1-h_{1})+h_{2} )$, then (1.2) has two positive equilibria

$$\begin{aligned} &E_{1}^{*}= \bigl( (\sqrt{17}-4 ) (1-h_{1}), (114\sqrt{17}-470 ) (1-h_{1})^{3} \bigr) \quad\textit{{and}} \\ & E^{*}= \biggl(\frac{5-\sqrt{17}}{2}(1-h_{1}), (264-64 \sqrt{17} ) (1-h_{1})^{3} \biggr). \end{aligned}$$

Moreover, (i) $E_{1}^{*}$ is an unstable multiple focus with multiplicity one. (ii) $E^{*}$ is a codimension 2 BT singularity. The phase diagram is shown in Fig. 3.

Figure 3

There exists an unstable multiple focus $E_{1}^{*}=(0.1108,0.0248)$ with multiplicity one and a BT singularity $E^{*}=(0.3946,0.0884)$ of codimension 2 for $h_{1}=0.1$ and $h_{2}=0.3$

Proof

When $a,\beta$ and δ satisfy the conditions in Theorem 2.2, system (1.2) becomes

$$ \begin{aligned} &\dot{x}=x(1-x)-\frac{2xy}{(41\sqrt{17}-169)(1-h_{1})^{2}+2x^{2}}-h_{1}x, \\ &\dot{y}=y \biggl( \frac{\sqrt{17}-3}{2}(1-h_{1})- \frac{4+\sqrt{17}}{4(1-h_{1})} \frac{y}{x} \biggr). \end{aligned} $$

(2.4)

Obviously, (2.4) has two equilibria $E_{1}^{*}$ and $E^{*}$ as stated above.

(1) To prove (i), let $u=x- (\sqrt{17}-4 )(1-h_{1})$, $v=y- (114\sqrt{17}-470 )(1-h_{1})^{3}$ and Taylor expand (2.4), thus we get

$$ \begin{aligned} \dot{u}={}&\frac{\sqrt{17}-3}{2}(1-h_{1})u- \frac{13+3\sqrt{17}}{8(1-h_{1})}v+3u^{2} - \frac{45+11\sqrt{17}}{16(1-h_{1})^{2}}uv \\ &{}- \frac{85+19\sqrt{17}}{8(1-h_{1})}u^{3} + \frac{235+57\sqrt{17}}{4(1-h_{1})^{3}}u^{2}v+O \bigl( \vert u,v \vert ^{4} \bigr), \\ \dot{v}={}& (50\sqrt{17}-206 ) (1-h_{1})^{3}u+ \frac{3-\sqrt{17}}{2}(1-h_{1})v \\ &{}+ (6\sqrt{17}-26 ) (1-h_{1})^{2}u^{2}+ ( \sqrt{17}+5 )uv \\ &{}-\frac{33+8\sqrt{17}}{4(1-h_{1})^{2}}v^{2}+ (2\sqrt{17}+2 ) (1-h_{1})u^{3} \\ &{}- \frac{37+9\sqrt{17}}{1-h_{1}}u^{2}v + \frac{268+65\sqrt{17}}{4(1-h_{1})^{3}}uv^{2}+O \bigl( \vert u,v \vert ^{4} \bigr). \end{aligned} $$

(2.5)

Suppose

$$ x=\frac{\sqrt{3526\sqrt{17}-14\text{,}538}}{2}(1-h_{1})^{2} u,\qquad y= \frac{45-11\sqrt{17}}{2}(1-h_{1})^{2} u+v, $$

then (2.5) can be changed into

$$ \dot{x}=-\omega _{0} y+f(x,y), \qquad \dot{y}= \omega _{0} x+g(x,y), $$

(2.6)

where $\omega _{0}=\frac{\sqrt{22\sqrt{17}-90}}{2}(1-h_{1})$ and

$$\begin{aligned} f(x,y)={}&\frac{\sqrt{7269+1763\sqrt{17}}}{8(1-h_{1})^{2}}x^{2} - \frac{45+11\sqrt{17}}{16(1-h_{1})^{2}}xy+ \frac{11\text{,}351+2753\sqrt{17}}{1024(1-h_{1})^{5}}x^{3} \\ &{}+ \frac{\sqrt{401\text{,}460\text{,}573+97\text{,}368\text{,}491\sqrt{17}}}{32(1-h_{1})^{5}}x^{2}y+O \bigl( \vert x,y \vert ^{4} \bigr), \\ g(x,y)={}&{-}\frac{235+57\sqrt{17}}{16(1-h_{1})^{2}}x^{2} + \frac{\sqrt{133\text{,}205+32\text{,}307\sqrt{17}}}{16(1-h_{1})^{2}}xy- \frac{33+8\sqrt{17}}{4(1-h_{1})^{2}}y^{2} \\ &{}+ \frac{\sqrt{151\text{,}456\text{,}733\text{,}853+36\text{,}733\text{,}653\text{,}611\sqrt{17}}}{1024(1-h_{1})^{5}}x^{3} \\ &{}-\frac{86\text{,}243+20\text{,}917\sqrt{17}}{64(1-h_{1})^{5}}x^{2}y+ \frac{\sqrt{2\text{,}088\text{,}374\text{,}221+50\text{,}650\text{,}514\sqrt{17}}}{64(1-h_{1})^{5}}xy^{2} \\ &{}+O \bigl( \vert x,y \vert ^{4} \bigr). \end{aligned}$$

The first Lyapunov number [23] can be expressed as

$$\begin{aligned} \operatorname{Re}c_{1}={}&\frac{1}{16} \biggl\{ (f_{xxx}+f_{xyy}+g_{xxy}+g_{yyy} ) \\ &{}+ \frac{1}{\omega _{0}} \bigl[f_{xy} (f_{xx} +f_{yy} ) - g_{xy} (g_{xx} + g_{yy} ) - f_{xx}g_{xx} + f_{yy}g_{yy} \bigr] \biggr\} {\bigg|}_{x=y=0} \\ ={}&\frac{1\text{,}428\text{,}479+346\text{,}457\sqrt{17}}{8192(1-h_{1})^{5}}>0, \end{aligned}$$

thus $E_{1}^{*}$ is an unstable multiple focus with multiplicity one.

(2) To prove (ii), let $x_{1}=x-\frac{5-\sqrt{17}}{2}(1-h_{1})$, $x_{2}=y- (264-64\sqrt{17} )(1-h_{1})^{3}$ and expand (2.4) in a power series at the origin, thus (2.4) becomes

$$ \begin{aligned} \dot{x_{1}}={}& \frac{\sqrt{17}-3}{2}(1-h_{1})x_{1}- \frac{\sqrt{17}+4}{4(1-h_{1})}x_{2}- \frac{\sqrt{17}+9}{8}x_{1}^{2} \\ &{}+ \frac{7\sqrt{17}+29}{16(1-h_{1})^{2}}x_{1}x_{2}+O \bigl( \vert x_{1},x_{2} \vert ^{3} \bigr), \\ \dot{x_{2}}={}& (50\sqrt{17}-206 ) (1-h_{1})^{3}x_{1}- \frac{(\sqrt{17}-3)(1-h_{1})}{2}x_{2} \\ &{}+ (45-11\sqrt{17} ) (1-h_{1})^{2}x_{1}^{2}+ \frac{\sqrt{17}+1}{2}x_{1}x_{2} \\ &{}-\frac{9\sqrt{17}+37}{16(1-h_{1})^{2}}x_{2}^{2}+O \bigl( \vert x_{1},x_{2} \vert ^{3} \bigr). \end{aligned} $$

(2.7)

Making the affine transformation

$$ y_{1}=\frac{25\sqrt{17}+103}{32(1-h_{1})^{3}}x_{2},\qquad y_{2}=x_{1}- \frac{7\sqrt{17}+29}{16(1-h_{1})^{2}}x_{2}, $$

we obtain

$$ \begin{aligned}& \dot{y_{1}}=y_{2}- \frac{\sqrt{17}+5}{4(1-h_{1})}y_{2}^{2}+O \bigl( \vert y_{1},y_{2} \vert ^{3} \bigr), \\ &\dot{y_{2}}=\frac{19-5\sqrt{17}}{8}(1-h_{1})^{2}y_{1}^{2}- \frac{\sqrt{17}+1}{4}(1-h_{1})y_{1}y_{2} + \frac{\sqrt{17}-7}{8}y_{2}^{2}+O \bigl( \vert y_{1},y_{2} \vert ^{3} \bigr). \end{aligned} $$

(2.8)

By a $C^{\infty }$ change of variables

$$ z_{1}=y_{1}-\frac{\sqrt{17}-7}{16}y_{1}^{2}+ \frac{\sqrt{17}+5}{4(1-h_{1})}y_{1}y_{2},\qquad z_{2}=y_{2}- \frac{\sqrt{17}-7}{8}y_{1}y_{2}, $$

system (2.8) becomes

$$ \begin{aligned}& \dot{z_{1}}=z_{2}+O \bigl( \vert z_{1},z_{2} \vert ^{3} \bigr), \\ &\dot{z_{2}}=\frac{19-5\sqrt{17}}{8}(1-h_{1})^{2}z_{1}^{2}- \frac{\sqrt{17}+1}{4}(1-h_{1})z_{1}z_{2} +O \bigl( \vert z_{1},z_{2} \vert ^{3} \bigr). \end{aligned} $$

(2.9)

Suppose ${\rho _{1}}=\frac{19-5\sqrt{17}}{8}(1-h_{1})^{2}$ and ${\rho _{2}}=-\frac{\sqrt{17}+1}{4}(1-h_{1})$, then ${\rho _{1}\rho _{2}}=\frac{33-7\sqrt{17}}{16}(1-h_{1})^{3}\neq 0$ for $0< h_{1}<1$. Hence $E^{*}$ is a BT singularity of codimension 2. □

3 Bifurcations

In this section, we study the subcritical Hopf bifurcation in a neighborhood of $E_{1}^{*}$ and the Bogdanov–Takens bifurcation in a neighborhood of $E^{*}$. Now we carry out bifurcation analysis for system (1.2) by choosing $h_{1}$ and $h_{2}$ as bifurcation parameters. When the conditions of Theorem 2.2 are satisfied, the unfolding system of (1.2) is given by

$$ \begin{aligned}& \dot{x}=x(1-x)-\frac{2xy}{(41\sqrt{17}-169)(1-h_{1})^{2}+2x^{2}}-(h_{1}+ \mu _{1})x, \\ &\dot{y}=y \biggl( \frac{\sqrt{17}-3}{2}(1-h_{1})+h_{2}- \frac{4+\sqrt{17}}{4(1-h_{1})} \frac{y}{x} \biggr)-(h_{2}+\mu _{2})y, \end{aligned} $$

(3.1)

where $\mu _{1}$ and $\mu _{2}$ are small parameters.

Through analysis, we have the following conclusions.

Theorem 3.1

When $\mu _{1}$ and $\mu _{2}$ vary near the origin, system (3.1) undergoes a subcritical Hopf bifurcation in a small neighborhood of $E_{1}^{*}$ and a Bogdanov–Takens bifurcation in a small neighborhood of $E^{*}$. Hence, when $h_{1}$ and $h_{2}$ vary, system (1.2) may have an unstable closed cycle around $E_{1}^{*}$, and an unstable closed cycle or an unstable homoclinic loop around $E^{*}$.

Proof

If $\mu _{1}=\mu _{2}=0$, then (3.1) has two equilibria $E_{1}^{*}$ and $E^{*}$, whose types are discussed in Theorem 2.2.

(1) First we verify that system (3.1) experiences a subcritical Hopf bifurcation around $E_{1}^{*}$ when $(\mu _{1},\mu _{2})$ varies in a small neighborhood of the origin. If $(\mu _{1},\mu _{2})\neq (0,0)$, (3.1) has an equilibrium with the following form:

$$ \begin{aligned} E_{1}&=(x_{1},y_{1})\\ &= \biggl( (\sqrt{17}-4 ) (1-h_{1})+\omega, \frac{7\sqrt{17}-29}{2} \bigl[ (3+ \sqrt{17} )\mu _{2}-4(1-h_{1}) \bigr](1-h_{1})x_{1} \biggr). \end{aligned}$$

(3.2)

Here ω is an infinitely small quantity for $\mu _{1}$ and $\mu _{2}$ small enough.

Substituting (3.2) into (3.1), one can find that $x_{1}$ should satisfy the following equations:

$$\begin{aligned} \begin{aligned} & x^{3} + (h_{1} + \mu _{1} - 1)x^{2}+ \biggl[\frac{1}{16}(13\sqrt{17}-53) (8h_{1}+9 \mu _{2}+\mu _{2}\sqrt{17}-8) (h_{1} - 1) \biggr]x \\ &\qquad{}+ \frac{1}{2}(41\sqrt{17}-169) (h_{1} + \mu _{1} - 1) (h_{1} - 1)^{2} = 0, \end{aligned} \\ &\operatorname{Tr} \bigl(J(E_{1}) \bigr) \\ &\quad = \frac{5-\sqrt{17}}{2}(1-h_{1})-2x_{1} \\ &\qquad{}+ \frac{(44+12\sqrt{17})(1-h_{1})^{2}[(169+41\sqrt{17})x_{1}^{3}-8(1-h_{1})^{2}x_{1}]}{[(169+41\sqrt{17})x_{1}^{2}+8(1-h_{1})^{2}]^{2}} -\mu _{1} \\ &\qquad{}+ \biggl(1- \frac{(84+20\sqrt{17})(1-h_{1})[(169+41\sqrt{17})x_{1}^{3}-8(1-h_{1})^{2} x_{1}]}{[(169+41\sqrt{17})x_{1}^{2}+8(1-h_{1})^{2}]^{2}} \biggr)\mu _{2}, \\ &\operatorname{Det} \bigl(J(E_{1}) \bigr) \\ &\quad = \biggl( 1 - h_{1} - \mu _{1} - 2x_{1} \\ &\qquad{}- \frac{(11+3\sqrt{17})(1 - h_{1})[(3+\sqrt{17})\mu _{2}-4(1-h_{1})][(169+ 41\sqrt{17}) x_{1}^{3}-8(1 - h_{1})^{2}x_{1}]}{[(169 +41\sqrt{17})x_{1}^{2}+8(1 - h_{1})^{2}]^{2}} \biggr) \\ &\qquad{}\times \biggl(\mu _{2} + \frac{3 - \sqrt{17}}{2}(1 - h_{1}) \biggr) \\ &\qquad{}+ \frac{(9+\sqrt{17})(1-h_{1})[(3+\sqrt{17})\mu _{2}-4(1-h_{1})]^{2}x_{1}}{4[(169+ 41\sqrt{17})x_{1}^{2}+8(1-h_{1})^{2}]}. \end{aligned}$$

(3.3)

Letting $\operatorname{Tr}(J(E_{1}))=0,\operatorname{Det}(J(E_{1}))>0$ and according to (3.3), we can get

$$ \begin{aligned} \mu _{1}={}& \bigl(-2\omega ^{5}+2(7\sqrt{17}-27) (1-h_{1})\omega ^{4} -2(349 \sqrt{17}-1439) (1-h_{1})^{2}\omega ^{3}\\ &{}+8(1137 \sqrt{17}-4688) (1-h_{1})^{3}\omega ^{2} -(37 \text{,}001\sqrt{17}-152\text{,}559) (1-h_{1})^{4}\omega \bigr) \\ &{}/\bigl(2\omega ^{4}-8(\sqrt{17}-4) (1-h_{1})\omega ^{3} +2(185\sqrt{17}-763) (1-h_{1})^{2}\omega ^{2}\\ &{}-4(983\sqrt{17}-4053) (1-h_{1})^{3}\omega +(5873\sqrt{17}-24\text{,}215) (1-h_{1})^{4}\bigr), \\ \mu _{2}={}& \bigl(2\omega ^{5}+(11\sqrt{17}-43) (1-h_{1})\omega ^{4}-2(25\sqrt{17}-103) (1-h_{1})^{2} \omega ^{3} \\ &{}+8(252\sqrt{17}-1039) (1-h_{1})^{3} \omega ^{2} +7(2575\sqrt{17}-10\text{,}617) (1-h_{1})^{4} \omega \bigr) \\ &{}/\bigl(2\omega ^{4}-8(\sqrt{17}-4) (1-h_{1})\omega ^{3}+2(185\sqrt{17}-763) (1-h_{1})^{2}\omega ^{2} \\ &{}+4(983\sqrt{17}-4053) (1-h_{1})^{3}\omega +(5873\sqrt{17}-24\text{,}215) (1-h_{1})^{4}\bigr). \end{aligned} $$

(3.4)

Substituting (3.2) and (3.4) into $\operatorname{Det}(J(E_{1}))$, we obtain $\operatorname{Det}(J(E_{1}))>0$ if and only if

$$ \frac{463\sqrt{17}+1783}{2}\omega ^{2}+ 2 (3 \sqrt{17}-22 ) (1-h_{1})\omega + \frac{11\sqrt{17}-45}{2}(1-h_{1})^{2}>0. $$

(3.5)

The discriminant of inequality (3.5) is always negative for $h_{1}<1$, so $\operatorname{Det}(J(E_{1}))>0$ for all small ω. The Hopf bifurcation curve of system (3.1) at $E_{1}^{*}$ is defined by

$$ H_{1}= \bigl\{ (\mu _{1},\mu _{2}) | \mu _{1} \text{{ and }} \mu _{2} \text{{ are given by }} \text{(3.4)} \text{{ and }} \omega { \text{{ is sufficiently small}}} \bigr\} . $$

Due to ${\lim_{\omega \rightarrow 0}} \frac{\mu _{2}}{\mu _{1}}= \frac{35+7\sqrt{17}}{41+13\sqrt{17}}$, thus the approximate expression of $H_{1}$ is given by $\mu _{2}=\frac{35+7\sqrt{17}}{41+13\sqrt{17}}\mu _{1}$ in a small neighborhood of the origin.

(2) Now we prove that system (3.1) undergoes a Bogdanov–Takens bifurcation at the BT singularity $E^{*}$. First translate $E^{*}$ to the origin by letting $x_{1}=x-\frac{5-\sqrt{17}}{2}(1-h_{1})$ and $x_{2}=y- (264-64\sqrt{17} )(1-h_{1})^{3}$ and Taylor expand system (3.1), then we have

$$ \begin{aligned} \dot{x_{1}}={}& \frac{\sqrt{17}-5}{2}(1-h_{1})\mu _{1} - \frac{\sqrt{17}-3}{8} ( \sqrt{17}\mu _{1}+4h_{1}+3\mu _{1}-4 )x_{1} - \frac{4+\sqrt{17}}{4(1-h_{1})}x_{2} \\ &{}- \frac{9+\sqrt{17}}{8}x_{1}^{2} + \frac{29+7\sqrt{17}}{16(1-h_{1})^{2}}x_{1}x_{2} + O \bigl( \vert x_{1},x_{2} \vert ^{3} \bigr), \\ \dot{x_{2}}={}&8\mu _{2} (8\sqrt{17}-33 ) (1-h_{1})^{3}+ (50\sqrt{17}-206 ) (1-h_{1})^{3}x_{1} \\ &{}- \frac{\sqrt{17}-3}{8} (\sqrt{17}\mu _{2}-4h_{1}+3\mu _{2}+4 )x_{2} \\ &{}+ (45-11\sqrt{17} ) (1-h_{1})^{2}x_{1}^{2}+ \frac{1+\sqrt{17}}{2}x_{1}x_{2}-\frac{37+9\sqrt{17}}{16(1-h_{1})^{2}}x_{2}^{2} +O \bigl( \vert x_{1},x_{2} \vert ^{3} \bigr). \end{aligned} $$

(3.6)

With an affine transformation $y_{1}=x_{1}$, $y_{2}=\frac{\sqrt{17}-3}{2}(1-h_{1})x_{1}- \frac{4+\sqrt{17}}{4(1-h_{1})}x_{2}$, system (3.6) becomes

$$ \begin{aligned} \dot{y_{1}}={}& \frac{\sqrt{17}-5}{2}(1-h_{1})\mu _{1}-\mu _{1}y_{1}+y_{2}- \frac{1+\sqrt{17}}{8}y_{1}^{2}- \frac{3+\sqrt{17}}{4(1-h_{1})}y_{1}y_{2}+ O \bigl( \vert y_{1},y_{2} \vert ^{3} \bigr), \\ \dot{y_{2}}= {}&(8-2\sqrt{17} ) (\mu _{1}- \mu _{2}) (1-h_{1})^{2}-\frac{(\sqrt{17}-3)(\mu _{1}-\mu _{2})}{2}(1-h_{1})y_{1}\\ &{}- \mu _{2}y_{2} +\frac{\sqrt{17}-7}{8}(1-h_{1})y_{1}^{2} \\ &{}-y_{1}y_{2}+\frac{5+\sqrt{17}}{4(1-h_{1})}y_{2}^{2} +O \bigl( \vert y_{1},y_{2} \vert ^{3} \bigr). \end{aligned} $$

(3.7)

Under a $C^{\infty }$ transformation of coordinates

$$ z_{1}=y_{1},\qquad z_{2}= \frac{\sqrt{17}-5}{2}(1-h_{1})\mu _{1}-\mu _{1}y_{1}+y_{2}- \frac{1+\sqrt{17}}{8}y_{1}^{2}- \frac{3+\sqrt{17}}{4(1-h_{1})}y_{1}y_{2}, $$

system (3.7) can be transformed into

$$ \begin{aligned}& \dot{z_{1}}=z_{2}+O \bigl( \vert z_{1},z_{2} \vert ^{3} \bigr), \\ &\dot{z_{2}}= \beta _{0} + \beta _{1}z_{1} + \frac{\sqrt{17}-5}{8} \bigl[(5+\sqrt{17})\mu _{2} - 2\mu _{1} \bigr] z_{2} + \beta _{2}z_{1}^{2}\\ &\phantom{\dot{z_{2}}=}{} + \frac{(5+\sqrt{17}) [(5+3\sqrt{17})\mu _{1}-8(1-h_{1}) ]}{32(1-h_{1})}z_{1}z_{2} +\frac{1}{2(1-h_{1})}z_{2}^{2}+O \bigl( \vert z_{1},z_{2} \vert ^{3} \bigr), \end{aligned} $$

(3.8)

where

$$ \begin{aligned} &\beta _{0}=\frac{\sqrt{17}-4}{2}(1-h_{1}) (\mu _{1}-\mu _{2}) \bigl[ (3+\sqrt{17} )\mu _{1}-4(1-h_{1}) \bigr], \\ &\beta _{1}=\frac{\sqrt{17}-3 }{16} \bigl[ (19+5\sqrt{17} )\mu _{1}^{2}- (6+2\sqrt{17} )\mu _{1}\mu _{2} -8(1-h_{1})\mu _{1}\\ &\phantom{\beta _{1}=}{}+ (10-2\sqrt{17} ) (1-h_{1})\mu _{2} \bigr], \\ &\beta _{2}=\frac{1}{128(1-h_{1})} (\sqrt{17}-7 ) \bigl[- (161+39\sqrt{17} )\mu _{1}^{2}- (24+8\sqrt{17} ) (1-h_{1})\mu _{1}\\ &\phantom{\beta _{2}=}{}+ (40+8\sqrt{17} ) (1-h_{1})\mu _{2}+16(1-h_{1})^{2} \bigr]. \end{aligned} $$

By another $C^{\infty }$ change $X_{1}=z_{1}-\frac{1}{4(1-h_{1})}z_{1}^{2}$, $X_{2}=z_{2}-\frac{1}{2(1-h_{1})}z_{1}z_{2}$, system (3.8) becomes

$$ \begin{aligned} \dot{X_{1}}={}&X_{2}+O \bigl( \vert X_{1},X_{2} \vert ^{3} \bigr), \\ \dot{X_{2}}={}&\beta _{0} - \frac{2\beta _{1}h_{1}+\beta _{0}-2\beta _{1}}{2(1-h_{1})}X_{1}+ \frac{\sqrt{17}-5}{8} \bigl[ (5+\sqrt{17} )\mu _{2} - 2\mu _{1} \bigr]X_{2} \\ &{}+ \psi X_{1}^{2} + \frac{(5+\sqrt{17}) [(5+3\sqrt{17})\mu _{1}-8(1-h_{1}) ]}{32(1-h_{1})}X_{1}X_{2} \\ &{}+O \bigl( \vert X_{1},X_{2} \vert ^{3} \bigr), \end{aligned} $$

(3.9)

where $\psi = \frac{8\beta _{2}h_{1}^{2}+2(\beta _{1}-8\beta _{2})h_{1}- \beta _{0}-2\beta _{1}+8\beta _{2}}{8(1-h_{1})^{2}}$.

Consider the $C^{\infty }$ change $Y_{1}=X_{1}, Y_{2}=X_{2}+O(|X_{1},X_{2}|^{3})$ and system (3.9) can be changed into

$$ \begin{aligned} \dot{Y_{1}}={}&Y_{2}, \\ \dot{Y_{2}}={}& \beta _{0} - \frac{2\beta _{1}h_{1} + \beta _{0}-2\beta _{1}}{2(1-h_{1})}Y_{1} + \frac{\sqrt{17}-5}{8} \bigl[ (5+\sqrt{17} )\mu _{2}-2\mu _{1} \bigr]Y_{2}\\ &{}+\psi Y_{1}^{2} + \frac{(5+\sqrt{17}) [(5+3\sqrt{17})\mu _{1}-8(1-h_{1}) ]}{32(1-h_{1})}Y_{1}Y_{2} + O \bigl( \vert Y_{1},Y_{2} \vert ^{3} \bigr). \end{aligned} $$

(3.10)

Substituting the expressions of $\beta _{0}, \beta _{1}$ and $\beta _{2}$ into ψ, we have $\psi <0$ for any small $\mu _{1}$ and $\mu _{2}$. Under the following change:

$$ Z_{1}=Y_{1}, \qquad Z_{2}= \frac{Y_{2}}{\sqrt{-\psi }},\qquad \tau =\sqrt{-\psi }t, $$

system (3.10) becomes

$$ \begin{aligned}& \dot{Z_{1}}=Z_{2}, \\ &\dot{Z_{2}}=\frac{\beta _{0}}{-\psi }+ \frac{2\beta _{1}h_{1}+\beta _{0}-2\beta _{1}}{2(1-h_{1})\psi }Z_{1}+ \frac{\sqrt{17}-5}{8\sqrt{-\psi }} \bigl[ (5+\sqrt{17} )\mu _{2}-2\mu _{1} \bigr]Z_{2}\\ &\phantom{\dot{Z_{2}}=}{}- Z_{1}^{2} + \frac{(5+\sqrt{17}) [(5+3\sqrt{17})\mu _{1}-8(1-h_{1}) ]}{32(1-h_{1})\sqrt{-\psi }}Z_{1}Z_{2} \\ &\phantom{\dot{Z_{2}}=}{}+O \bigl( \vert Z_{1},Z_{2} \vert ^{3} \bigr). \end{aligned} $$

(3.11)

The transformation $U=Z_{1}- \frac{2\beta _{1}h_{1}+\beta _{0}-2\beta _{1}}{4(1-h_{1})\psi }, V=Z_{2}$ brings (3.11) into the form

$$ \begin{aligned} &\dot{U}=V, \\ &\dot{V}= \frac{(2\beta _{1}h_{1}+\beta _{0}-2\beta _{1})^{2}}{16(1-h_{1})^{2}\psi ^{2}}- \frac{\beta _{0}}{\psi } +\phi V-U^{2}\\ &\phantom{\dot{V}=}{}+ \frac{(5+\sqrt{17}) [(5+3\sqrt{17})\mu _{1}-8(1-h_{1}) ]}{32(1-h_{1})\sqrt{-\psi }}UV+O \bigl( \vert U,V \vert ^{3} \bigr), \end{aligned} $$

(3.12)

where $\phi = \frac{(5+\sqrt{17})(2\beta _{1}h_{1}+\beta _{0}-2\beta _{1}) [(5+3\sqrt{17})\mu _{1}-8(1-h_{1}) ]}{128(1-h_{1})^{2}\psi \sqrt{-\psi }}+ \frac{(\sqrt{17}-5) [(5+\sqrt{17})\mu _{2}-2\mu _{1} ]}{8\sqrt{-\psi }}$.

The derivatives of (3.11) and (3.12) are about time τ. For $\mu _{1}$ small enough we make the following transformation:

$$ \begin{aligned} &x = - \biggl(\frac{(5+\sqrt{17}) [(5+3\sqrt{17})\mu _{1}-8(1-h_{1}) ]}{32(1-h_{1})\sqrt{-\psi }} \biggr)^{2}U, \\ &y = \biggl(\frac{(5+\sqrt{17}) [(5+3\sqrt{17})\mu _{1}-8(1-h_{1}) ]}{32(1-h_{1})\sqrt{-\psi }} \biggr)^{3}V, \\ & t = - \frac{32(1-h_{1})\sqrt{-\psi }}{(5+\sqrt{17}) [(5+3\sqrt{17})\mu _{1}-8(1-h_{1}) ]} \tau, \end{aligned} $$

then system (3.12) becomes

$$ \begin{aligned}& \dot{x}=y, \\ &\dot{y}=\eta _{1}(\mu _{1},\mu _{2})+\eta _{2}(\mu _{1},\mu _{2})y+x^{2}+xy+O \bigl( \vert x,y \vert ^{3} \bigr), \end{aligned} $$

(3.13)

where the derivatives are about time t and the expressions of $\eta _{1}(\mu _{1},\mu _{2})$ and $\eta _{2}(\mu _{1},\mu _{2})$ are given by

$$ \begin{aligned}& \eta _{1}(\mu _{1},\mu _{2}) \\ &\quad = \frac{(53+13\sqrt{17})^{2} ( (5+3\sqrt{17} )\mu _{1}-8(1-h_{1}) )^{4}}{134\text{,}217\text{,}728(1-h_{1})^{4}Q(\mu _{1},\mu _{2})^{4}} \\ &\qquad{}\times \biggl[ 8 (13\sqrt{17}-53 ) (1-h_{1}) ( \mu _{1}-\mu _{2}) (3\mu _{1}+\sqrt{17}\mu _{1} +4h_{1}-4 )Q(\mu _{1},\mu _{2}) \\ &\qquad{}+ \frac{33\sqrt{17}-137}{4} (3\sqrt{17}\mu _{1}+11\mu _{1}-\sqrt{17}\mu _{2}-3\mu _{2}+4h_{1}-4 )^{2}\mu _{1}^{2} \biggr], \\ &\eta _{2}(\mu _{1},\mu _{2}) \\ &\quad= \frac{(5\mu _{1}+3\sqrt{17}\mu _{1}+8h_{1}-8)(5+\sqrt{17})}{1024(1-h_{1})^{2}Q(\mu _{1}, \mu _{2})^{2}}\\ &\qquad{}\times \biggl[ 4 (\sqrt{17}-5 ) (1-h_{1}) (-2\mu _{1}\sqrt{17}\mu _{2}+5\mu _{2} )Q(\mu _{1},\mu _{2}) \\ &\qquad{} -\frac{1}{2} (3\sqrt{17}\mu _{1}+11\mu _{1}-\sqrt{17}\mu _{2}-3\mu _{2}+4h_{1}-4 ) (3\sqrt{17}\mu _{1}+5\mu _{1}+8h_{1}-8 ) \mu _{1} \biggr], \\ &Q(\mu _{1},\mu _{2}) \\ &\quad =\frac{1}{16(1-h_{1})} \bigl[ (46+14 \sqrt{17} )\mu _{1}^{2}- (\sqrt{17}-9 )\mu _{1}\mu _{2}+2 (5\sqrt{17}-9 ) (1-h_{1})\mu _{1}\\ &\qquad{}-2 (5 \sqrt{17}-7 ) (1-h_{1}) \mu _{2} \\ &\qquad{} + 2 (\sqrt{17}-7 ) (1-h_{1})^{2} \bigr]. \end{aligned} $$

We can check

$$ \biggl\vert \frac{\partial (\eta _{1},\eta _{2})}{\partial (\mu _{1},\mu _{2})} \biggr\vert _{\mu _{1}=\mu _{2}=0}= \frac{4299\sqrt{17}+17\text{,}725}{64(1-h_{1})^{2}}\neq 0, $$

for $0< h_{1}<1$, which indicates that above parameter transformation is regular at $(0,0)$.

From the results in [2] and [27], we know that system (3.1) undergoes a Bogdanov–Takens bifurcation when $(\mu _{1},\mu _{2})$ varies in a small neighborhood of the origin, and the local expressions of the bifurcation curves in a small neighborhood of the origin are given by:

• (a) the saddle-node bifurcation curve

  $$\begin{aligned} SN={}& \biggl\{ (\mu _{1},\mu _{2})| \eta _{1}=0,\text{ i.e.}, \mu _{1}-\mu _{2}+ \frac{(19\sqrt{17}-43)\mu _{1}^{2}}{32(1-h_{1})}\\ &{} - \frac{(53\sqrt{17}+35)\mu _{1}\mu _{2}}{16(1-h_{1})} + \frac{3(7\sqrt{17}+9)\mu _{2}^{2}}{8(1-h_{1})}+ O \bigl( \vert \mu _{1},\mu _{2} \vert ^{3} \bigr)=0 \biggr\} , \end{aligned}$$

  which includes $SN^{+}=\{(\mu _{1},\mu _{2})| \eta _{1}=0,\eta _{2}>0\}$ and $SN^{-}=\{(\mu _{1},\mu _{2})| \eta _{1}=0,\eta _{2}<0\}$;

• (b) the Hopf bifurcation curve

  $$\begin{aligned} H={}& \bigl\{ (\mu _{1},\mu _{2})| \eta _{2}=\sqrt{-\eta _{1}}, \eta _{1}< 0 \bigr\} \\ = {}&\biggl\{ (\mu _{1},\mu _{2}) | \mu _{1}-\mu _{2} + \frac{5(7\sqrt{17}-23)\mu _{1}^{2}}{16(1-h_{1})} - \frac{(73\sqrt{17}-33)\mu _{1}\mu _{2}}{16(1-h_{1})} \\ &{}+ \frac{(23\sqrt{17}+21)\mu _{2}^{2}}{8(1-h_{1})}+ O \bigl( \vert \mu _{1},\mu _{2} \vert ^{3} \bigr)=0, \eta _{1}< 0 \biggr\} ; \end{aligned}$$

• (c) the homoclinic bifurcation curve

  $$\begin{aligned} HL={}& \biggl\{ (\mu _{1},\mu _{2})| \eta _{2}=\frac{5}{7}\sqrt{-\eta _{1}}, \eta _{1}< 0 \biggr\} \\ ={}& \biggl\{ (\mu _{1},\mu _{2})| \mu _{1} - \mu _{2} + \frac{(1487\sqrt{17}-5119)\mu _{1}^{2}}{400(1-h_{1})} - \frac{(2305\sqrt{17}-2457) \mu _{1}\mu _{2} }{400(1-h_{1})} \\ &{}+ \frac{(623\sqrt{17}+381) \mu _{2}^{2}}{200(1-h_{1})} + O \bigl( \vert \mu _{1},\mu _{2} \vert ^{3} \bigr)=0,\eta _{1}< 0 \biggr\} ; \end{aligned}$$

• (d) the Hopf bifurcation curve of system (3.1) at $E_{1}^{*}$ is

  $$H_{1}= \biggl\{ (\mu _{1},\mu _{2}) | \mu _{2}= \frac{35+7\sqrt{17}}{41+13\sqrt{17}}\mu _{1} \biggr\} . $$

The sketches of these bifurcation curves are displayed in Fig. 4.

Figure 4

The bifurcation diagram of system (3.1) for $h_{1}=0.1$ and $h_{2}=0.3$

Due to $\eta _{2}<0$ around the origin, H and HL are nonexistent in the third quadrant of Fig. 4, so only $SN^{-}$ and $H_{1}$ are considered in this part. Obviously, the neighborhood of the origin in the $(\mu _{1},\mu _{2})$-plane is divided into six parts by $H_{1}, SN^{+}, H, HL$ and $SN^{-}$.

When $(\mu _{1},\mu _{2})$ lies in region I, system (3.1) has an unstable focus $E_{1}^{*}$. By Lemma 1.1 and the Poincaré–Bendixson theorem, it can be proved that there exists a stable limit cycle. The corresponding phase portrait is depicted in Fig. 5. From the biological point of view, we know that the prey and the predators can coexist and will periodically fluctuate along this cycle eventually.

Figure 5

For $(\mu _{1}, \mu _{2})=(0.01,0.007)$ in region I, there exists an unstable focus $E_{1}^{*}$ surrounded by a stable limit cycle

When $(\mu _{1},\mu _{2})$ lies on the curve $H_{1}$, the stable limit cycle still exists and the equilibrium $E_{1}^{*}$ turns into an unstable multiple focus with multiplicity one. When $(\mu _{1},\mu _{2})$ crosses $H_{1}$ and enters into the region VI, $E_{1}^{*}$ becomes a stable focus and an unstable limit cycle bifurcates, which indicates a subcritical Hopf bifurcation happens at $E_{1}^{*}$ when $(\mu _{1},\mu _{2})$ goes through $H_{1}$. Figure 6 shows that there is a bistable state (a stable equilibrium $E_{1}^{*}$ and a big stable limit cycle) and the unstable limit cycle bifurcated from Hopf bifurcation acts as the separatrices of attractive domains for different attractors. Eventually, trajectories can either tend to an equilibrium or a limit cycle depending on initial values of the species.

Figure 6

When $(\mu _{1}, \mu _{2})=(0.01,0.001)$ lies in region VI, a stable focus $E_{1}^{*}$ is surrounded by an unstable limit cycle, which is circled by a stable cycle

When $(\mu _{1},\mu _{2})$ lies on the curve $SN^{-}$, a new saddle-node $E^{*}$ emerges and the other equilibrium $E_{1}^{*}$ is still a stable focus. When $(\mu _{1},\mu _{2})$ passes through $SN^{-}$ and enters into region V, there are two new equilibria $E_{2}^{*}$ and $E_{3}^{*}$ bifurcated from the saddle-node bifurcation. $E_{2}^{*}$ is a hyperbolic saddle and $E_{3}^{*}$ is a stable focus. The type of $E_{1}^{*}$ is not changed and there is still an unstable limit cycle surrounding $E_{1}^{*}$. The Poincaré–Bendixson theorem implies that there is a big stable limit cycle surrounding the three equilibria and the small limit cycle, which is displayed in Fig. 7. It is shown that there is a tristable state in system (3.1). The stable and unstable manifolds of the saddle $E_{2}^{*}$ and the unstable limit cycle act as separatrices of different attractive domains. With different initial values, the species will eventually tend to a constant state or oscillate along a big limit cycle.

Figure 7

When $(\mu _{1},\mu _{2})=(-0.01,-0.008)$ lies in region V, there are three equilibria $E_{1}^{*}, E_{2}^{*}$ and $E_{3}^{*}$. $E_{1}^{*}$ and $E_{3}^{*}$ are stable foci and $E_{2}^{*}$ is a hyperbolic saddle

When $(\mu _{1},\mu _{2})$ lies on the curve $H_{1}$, system (3.1) has three equilibria $E_{1}^{*}, E_{2}^{*}$ and $E_{3}^{*}$, where $E_{1}^{*}$ becomes an unstable multiple focus with multiplicity one and the types of $E_{2}^{*}$ and $E_{3}^{*}$ are unchanged. When $(\mu _{1},\mu _{2})$ crosses $H_{1}$ into region IV, the unstable limit cycle enclosing $E_{1}^{*}$ disappears and $E_{1}^{*}$ turns into an unstable focus, which shows that system (3.1) undergoes a subcritical Hopf bifurcation when $(\mu _{1},\mu _{2})$ passes through $H_{1}$. The corresponding phase portrait is displayed in Fig. 8, in which $E_{3}^{*}$ is a unique globally asymptotical attractor.

Figure 8

When $(\mu _{1},\mu _{2})=(0.01,0.015)$ lies in region IV, system (3.1) has three equilibria $E_{1}^{*}, E_{2}^{*}$ and $E_{3}^{*}$. $E_{2}^{*}$ is a hyperbolic saddle. $E_{1}^{*}$ is an unstable focus and $E_{3}^{*}$ is a globally stable focus

When $(\mu _{1},\mu _{2})$ lies on the curve HL, there appears an unstable homoclinic loop to the saddle $E_{2}^{*}$ which is displayed in Fig. 9(a), and a large stable limit cycle encloses all the equilibria and the homoclinic loop shown in Fig. 9(b). Moreover, when $(\mu _{1},\mu _{2})$ deviates from the curve HL, the homoclinic loop breaks, which means that system (3.1) undergoes a homoclinic bifurcation when $(\mu _{1},\mu _{2})$ passes through the curve HL.

Figure 9

(a) A local amplified diagram of the homoclinic loop for $(\mu _{1},\mu _{2})=(0.01,0.009962)$ on the curve HL. (b) A big stable limit cycle encloses all the equilibria

When $(\mu _{1},\mu _{2})$ lies in region III, there occurs an unstable limit cycle surrounding the stable focus $E_{3}^{*}$, which means that system (3.1) undergoes a subcritical Hopf bifurcation when $(\mu _{1},\mu _{2})$ passes through the curve H. For $(\mu _{1},\mu _{2})=(0.01,0.009952)$, there are three equilibria, a small limit cycle and a large stable limit cycle enclosing all the equilibria and the small cycle, which is shown in Figs. 10(a) and 10(b). In this case there is also a bistable state (a stable equilibrium $E_{3}^{*}$ and a big stable limit cycle).

Figure 10

(a) A local amplification diagram of an unstable limit cycle for $(\mu _{1},\mu _{2})=(0.01.0.009952)$ in region III. (b) Two unstable equilibria $E_{1}^{*}$ and $E_{2}^{*}$, and a stable focus $E_{3}^{*}$ are surrounded by a stable limit cycle

When $(\mu _{1},\mu _{2})$ goes through H and enters into region II (the region between H and $SN^{+}$), the unstable limit cycle disappears and $E_{3}^{*}$ becomes an unstable focus while the types of $E_{1}^{*}$ and $E_{2}^{*}$ are not changed. When $(\mu _{1},\mu _{2})$ lies on the curve $SN^{+}$, two equilibria $E_{2}^{*}$ and $E_{3}^{*}$ overlap and become a saddle-node $E^{*}$, which means that system (3.1) undergoes a saddle-node bifurcation when $(\mu _{1},\mu _{2})$ passes through the curve $SN^{+}$. Whatever $(\mu _{1},\mu _{2})$ lies in region II or on $SN^{+}$, there both exists a stable limit cycle enclosing all the equilibria by the Poincaré–Bendixson theorem. The phase portrait of system (3.1) for $(\mu _{1},\mu _{2})$ on the curve $SN^{+}$ is exhibited in Fig. 11. From the strategy of the optimal management of renewable resources, we find that the saddle-node bifurcation curve acts as the feasible upper bound for the rescaled harvesting efforts which guarantees the coexistence and sustainable development of the species.

Figure 11

The phase portrait of system (3.1) for $(\mu _{1},\mu _{2})=(0.01,0.0099393)$ on the curve $SN^{+}$

From Fig. 12(a), we find that when $(\mu _{1},\mu _{2})$ varies along the saddle-node bifurcation curve $SN^{+}$ as shown in Fig. 4, a cusp point CP can be encountered and a cusp bifurcation will happen when $(\mu _{1},\mu _{2})$ varies near the CP point. The coordinates and parameter values for the CP point and the normal form coefficient c for the cusp bifurcation are given as follows:

$$\begin{aligned} &\mathrm{Label} = CP, (x,y,\mu _{1}, \mu _{2}) = (0.239806, 0.036775, 0.180583, 0.159373)\quad \text{and}\\ & c=7.915010. \end{aligned}$$

When $(\mu _{1},\mu _{2})$ varies along two Hopf bifurcation curves H and $H_{1}$ as indicated in Fig. 4, respectively, two generalized Hopf points $GH^{+}$ and $GH^{-}$ can be encountered, at which both of the first Lyapunov coefficients are zero. The associated coordinates and parameter values for the two points and the second Lyapunov coefficients $l^{\pm }_{2}$ of the generalized Hopf bifurcations are given as follows:

$$\begin{aligned} &\mathrm{Label} = GH^{+}, (x,y,\mu _{1},\mu _{2}) = ( 0.260349, 0.028104, 0.316434, 0.261824) \quad\text{and}\\ &l^{+}_{2}=-550.7068. \\ &\mathrm{Label} = GH^{-}, (x,y,\mu _{1},\mu _{2}) = (0.126604, 0.015720, 0.326771, 0.225224) \quad\text{and}\\ &l^{-}_{2}=-2300.709. \end{aligned}$$

Figure 12

(a) The point BT is corresponding to the origin of Fig. 4. When $(\mu _{1},\mu _{2})$ varies along curves $SN^{+}$, H and $H_{1}$, respectively, the points CP, $GH^{+}$ and $GH^{-}$ can be encountered. (b) The bifurcation portrait for $h_{1}=0.1$ and $h_{2}=0.3$, $SN_{c}$ indicates the saddle-node bifurcation curve originating from the point CP. $H^{+}$($H^{-}$) and $T^{+}$($T^{-}$) express the Hopf bifurcation curve and the fold bifurcation curve of the cycles originating from $GH^{+}$($GH^{-}$), respectively

From Fig. 12(b), we find that the saddle-node bifurcation curve $SN_{c}$ originating from the CP point and the saddle-node bifurcation curve $SN^{+}$ originating from the BT point are nearly coincident or tangent at the point CP. Through the curve $SN_{c}$ three equilibria collide to form an equilibrium or one equilibrium splits into three equilibria. Similarly, the Hopf bifurcation curve $H\ (H_{1})$ originating from the BT point and the Hopf bifurcation curve $H^{+}$ $(H^{-})$ originating from the $GH^{+}$ ($GH^{-}$) point are almost coincident and tangent at the $GH^{+}$ ($GH^{-}$) point. In addition, near the Hopf bifurcation curves $H^{+}$ and $H^{-}$, two fold bifurcation curves of the cycles $T^{+}$ and $T^{-}$ also appear, on which two limit cycles with different stability collide to form a semi-stable cycle. Several typical phase portraits for $(\mu _{1},\mu _{2})$ near the generalized Hopf point $GH^{+}$ are displayed in Fig. 13.

Figure 13

(a) For $(\mu _{1},\mu _{2})=(0.28,0.238999)$ between the curves $H^{+}$ and $T^{+}$, there are two limit cycles. The outer one is stable and the inner one is unstable. (b) For $(\mu _{1},\mu _{2})=(0.2813,0.237)$ below the curve $H^{+}$, there is only a stable cycle. (c) For $(\mu _{1},\mu _{2})=(0.2819,0.23969)$ above the curve $T^{+}$, there is no limit cycle. (d) For $(\mu _{1},\mu _{2})=(0.1593,0.1445)$ between the wedge of $SN_{c}$, three equilibria $E_{1}, E_{2}$ and $E_{3}$ are a spiral source, a nodal source and a saddle, respectively

 □

4 Conclusions

The prey–predator model in [20] added with constant-effort harvesting has been investigated in this paper. The Bogdanov–Takens bifurcation, the cusp bifurcation and the generalized Hopf bifurcation are discussed. By analysis, it is shown that this model can generate many novel dynamic behaviors compared with the model with no harvesting. For example, a globally asymptotically stable equilibrium can appear, while the stable equilibria in the original model are all locally asymptotically stable. Through the generalized Hopf bifurcation, the second Lyaponov coefficient can determine the relative position between the stable limit cycle and the unstable one. A cusp point CP has also been detected, from which two saddle-node bifurcation curves of equilibria emanate, through which three equilibria collide to form an equilibrium or one equilibrium splits into three equilibria. These bifurcations have not been discussed in the original model. In addition, compared with the original model, the harvested prey–predator model can exhibit different bifurcation structure in the parameter plane. These phenomena show that the added harvesting terms play an important role in directing the evolutional directions of the species. From the Hopf and Bogdanov–Takens bifurcations, we find that two small limit cycles bifurcating from two different equilibria are both unstable, while the big limit cycle is always stable, which shows that the ultimate numbers of the species circulate periodically along the big cycle rather than the small cycle.

Appendix

Denoting coefficients of Eq. (2.2) as $b= (h_{1}-1), c= ( a+{\frac{\delta -h_{2}}{\beta }} ), d=a(h_{1}-1)$ and letting $A_{1}=b^{2}-3c, B=bc-9d, C=c^{2}-3bd$, the discriminant $\Delta _{1}$ of (2.2) is given by $\Delta _{1}=B^{2}-4A_{1}C$, from which we can get $A_{1}=A, \Delta _{1}=\frac{\Delta }{9}$. Suppose that $x_{1}, x_{2}$ and $x_{3}$ are three roots of Eq. (2.2).

$F'(x)=3x^{2}+2(h_{1}-1)x+a+{\frac{\delta -h_{2}}{\beta }}$ and the discriminant $\Delta _{2}$ of $F'(x)$ is given by $\Delta _{2}=4(h_{1}-1)^{2}- 12 ( a+{\frac{\delta -h_{2}}{\beta }} )=4A$. Now we prove Lemma 2.1.

• (a) If $\Delta >0$, then $\Delta _{1}>0$. (2.2) has one real root $x_{1}$ and a pair of conjugate complex roots $x_{2}, x_{3}$ according to the root formula of the cubic equation. Owing to $x_{1}x_{2}x_{3}=a(1-h_{1})>0$, $x_{1}$ is the positive root of (2.2). Letting $(x^{*},y^{*})= (x_{1},\frac{\delta -h_{2}}{\beta }x_{1} )$, now we will discuss three cases.

  • (1) If $A= \frac{\Delta _{2}}{4}<0$, then $F'(x^{*})>0$.

  • (2) If $A= \frac{\Delta _{2}}{4}=0$, then $x^{*}=\frac{1}{3} (1-h_{1}+\sqrt[3]{(1-h_{1})(27a-(1-h_{1})^{2})} )$ and $F'(x^{*})=\frac{1}{3} \sqrt[3]{(1-h_{1})^{2}(27a-(1-h_{1})^{2})^{2}}=\frac{1}{3}\sqrt[3]{\Delta }>0$.

  • (3) If $A= \frac{\Delta _{2}}{4}>0$, then $x^{*}= \frac{1-h_{1}-\sqrt[3]{Y_{1}}-\sqrt[3]{Y_{2}}}{3}$ and $F'(x^{*})= \frac{1}{3} (\sqrt[3]{Y_{1}^{2}}+\sqrt[3]{Y_{2}^{2}}+2\sqrt[3]{Y_{1}}\sqrt[3]{Y_{2}}-A ) = \frac{1}{3} (\sqrt[3]{Y_{1}^{2}}+\sqrt[3]{Y_{2}^{2}}+A )>0$, where $Y_{1,2}=(h_{1}-1)A+\frac{3(-B\pm \sqrt{\Delta _{1}})}{2}$.

  For the above three cases, $F'(x^{*})>0$ always holds. So $E^{*}(x^{*},y^{*})$ is an anti-saddle of system (1.2).

• (b) If $\Delta =0, A=0$, then $x_{1}=x_{2}=x_{3}=x^{*}=\frac{1-h_{1}}{3}$ is a triple root of Eq. (2.2). Now system (1.2) has a unique equilibrium $E^{*}(x^{*},y^{*})= (\frac{1-h_{1}}{3}, \frac{(1-h_{1})(\delta -h_{2})}{3\beta } )$, which is degenerated according to $F'(x^{*})=0$.

• (c) If $\Delta =0, A>0$, then $x_{1}=1-h_{1}+\frac{B}{A}$ and $x_{2}=x_{3}=-\frac{B}{2A}$ are two positive roots of Eq. (2.2) for $\operatorname{max} \{ \frac{(h_{1}-1)^{2}-4A}{27},0 \} < a<\frac{(h_{1}-1)^{2}-A}{27}$. So system (1.2) has two equilibria:

  $$\begin{aligned} &E_{1}^{*} \bigl(x_{1}^{*},y_{1}^{*} \bigr)\\ &\quad= \biggl(x_{1}, \frac{\delta -h_{2}}{\beta }x_{1} \biggr)\\ &\quad= \biggl(\frac{(1-h_{1}) [4A+27a-(h_{1}-1)^{2} ]}{3A},\frac{(\delta -h_{2})(1-h_{1}) [4A+27a-(h_{1}-1)^{2} ]}{3\beta A} \biggr), \\ &E^{*} \bigl(x^{*},y^{*} \bigr)\\ &\quad= \biggl(x_{2},\frac{\delta -h_{2}}{\beta }x_{2} \biggr)\\ &\quad= \biggl( \frac{(1-h_{1}) [(h_{1}-1)^{2}-A-27a ]}{6A},\frac{(\delta -h_{2})(1-h_{1}) [(h_{1}-1)^{2}-A-27a ]}{6\beta A} \biggr). \end{aligned}$$

  Because of $F'(x_{1}^{*})=\frac{1}{3A^{2}} [(1-h_{1})(3A+27a-(1-h_{1})^{2}) ]^{2}=A>0$ and $F'(x^{*})=0$, $E_{1}^{*}(x_{1}^{*},y_{1}^{*})$ is an anti-saddle equilibrium and $E^{*}(x^{*},y^{*})$ is a degenerated equilibrium.

• (d) If $\Delta <0$, then $A>0$ and (2.2) has three different roots:

  $$ x_{1}=\frac{1-h_{1}-2\sqrt{A}\operatorname{cos}\frac{\theta }{3}}{3}, \qquad x_{2,3}=\frac{1-h_{1}+\sqrt{A} (\operatorname{cos}\frac{\theta }{3}\mp \sqrt{3}\operatorname{sin}\frac{\theta }{3} )}{3}, $$

  where $\theta =\operatorname{arccos}T, -1< T=\frac{2A(h_{1}-1)-3B}{2\sqrt{A^{3}}}<1$. By simple analysis, we get

  $$ 0< 1-h_{1}+\sqrt{A}< x_{2}< 1-h_{1}+2\sqrt{A},\qquad 0< 1-h_{1}-\sqrt{A}< x_{3}< 1-h_{1}+\sqrt{A}, $$

  so $x_{1}>0$. Letting $E_{1}^{*}(x_{1}^{*},y_{1}^{*})= (x_{1},\frac{\delta -h_{2}}{\beta }x_{1} ), E_{2}^{*}(x_{2}^{*},y_{2}^{*})= (x_{2},\frac{\delta -h_{2}}{\beta }x_{2} ), E_{3}^{*}(x_{3}^{*},y_{3}^{*})= (x_{3},\frac{\delta -h_{2}}{\beta }x_{3} )$ and due to

  $$\begin{aligned} &0< F' \bigl(x_{1}^{*} \bigr)= \frac{4A\operatorname{cos}^{2}\frac{\theta }{3}-A}{3} < A,\qquad 0< F' \bigl(x_{3}^{*} \bigr)= \frac{A}{9} \biggl(12\operatorname{cos}^{2} \frac{\pi -\theta }{3} - 3 \biggr) < A,\\ & \frac{-A}{3} < F' \bigl(x_{2}^{*} \bigr)=\frac{A}{9} \biggl(12 \operatorname{cos}^{2}\frac{2\pi -\theta }{3} - 3 \biggr) < 0, \end{aligned}$$

  hence $E_{1}^{*},E_{3}^{*}$ are two anti-saddle equilibria and $E_{2}^{*}$ is a hyperbolic saddle.