Simple-Zero and Double-Zero Singularities of a Kaldor-Kalecki Model of Business Cycles with Delay

Abstract

We study the Kaldor-Kalecki model of business cycles with delay in both the gross product and the capital stock. Simple-zero and double-zero singularities are investigated when bifurcation parameters change near certain critical values. By performing center manifold reduction, the normal forms on the center manifold are derived to obtain the bifurcation diagrams of the model such as Hopf, homoclinic and double limit cycle bifurcations. Some examples are given to confirm the theoretical results.

1. Introduction

In the last decade, the study of delayed differential equations that arose in business cycles has received much attention. The first model of business cycles can be traced back to Kaldor [1] who used a system of ordinary differential equations to study business cycles in 1940 by proposing nonlinear investment and saving functions so that the system may have cyclic behaviors or limit cycles, which are important from the point of view of economics. Kalecki [2] introduced the idea that there is a time delay for investment before a business decision. Krawiec and Szydłowski [3–5] incorporated the idea of Kalecki into the model of Kaldor by proposing the following Kaldor-Kalecki model of business cycles:

𝑑𝑌(𝑡)𝑑𝑡=𝛼[𝐼(𝑌(𝑡),𝐾(𝑡))−𝑆(𝑌(𝑡),𝐾(𝑡))],𝑑𝐾(𝑡)𝑑𝑡=𝐼(𝑌(𝑡−𝜏),𝐾(𝑡))−𝑞𝐾(𝑡), (1.1) where 𝑌 is the gross product, 𝐾 is the capital stock, 𝛼>0 is the adjustment coefficient in the goods market, 𝑞∈(0,1) is the depreciation rate of capital stock, 𝐼(𝑌,𝐾) and 𝑆(𝑌,𝐾) are investment and saving functions, and 𝜏≥0 is a time lag representing delay for the investment due to the past investment decision. This model has been studied extensively by many authors; see [6–11]. Several authors also discussed similar models [12–14] and established the existence of limit cycles.

Considering that past investment decisions [6] also influence the change in the capital stock, Kaddar and Talibi Alaoui [15] extended the model (1.1) by imposing delays in both the gross product and capital stock. Thus adding the same delay to the capital stock 𝐾 in the investment function 𝐼(𝑌,𝐾) of the second equation of Sys. (1.1) leads to the following Kaldor-Kalecki model of business cycles:

𝑑𝑌(𝑡)𝑑𝑡=𝛼[𝐼(𝑌(𝑡),𝐾(𝑡))−𝑆(𝑌(𝑡),𝐾(𝑡))],𝑑𝐾(𝑡)𝑑𝑡=𝐼(𝑌(𝑡−𝜏),𝐾(𝑡−𝜏))−𝑞𝐾(𝑡). (1.2)

As in [3]; also see [10, 16, 17], using the following saving and investment functions 𝑆 and 𝐼, respectively, 𝑆(𝑌,𝐾)=𝛾𝑌,𝐼(𝑌,𝐾)=𝐼(𝑌)−𝛽𝐾, (1.3) where 𝛽>0 and 𝛾∈(0,1) are constants, we obtain the following system:

𝑑𝑌(𝑡)𝑑𝑡=𝛼[𝐼(𝑌(𝑡))−𝛽𝐾(𝑡)−𝛾𝑌(𝑡)],𝑑𝐾(𝑡)𝑑𝑡=𝐼(𝑌(𝑡−𝜏))−𝛽𝐾(𝑡−𝜏)−𝑞𝐾(𝑡). (1.4) Kaddar and Talibi Alaoui [15] studied the characteristic equation of the linear part of Sys. (1.4) at an equilibrium point and used the delay 𝜏 as a bifurcation parameter to show that the Hopf bifurcation may occur under some conditions as 𝜏 passes some critical values. However, they did not obtain the stability of the bifurcating limit cycles and the direction of the Hopf bifurcation. Wang and Wu [18] further studied Sys. (1.4) and gave a more detailed discussion of the distribution of the eigenvalues of the characteristic equation which has a pair of purely imaginary roots. They derived the normal forms on the center manifold for sys. (1.4) to give the direction of the Hopf bifurcation and the stability of the bifurcating limit cycles for some critical values of 𝜏.

However, under certain conditions, the characteristic equation of the linear part of Sys. (1.4) may have a simple-zero root, a double-zero root, or a simple zero root and a pair of purely imaginary roots. In this paper, simple-zero (fold) and double-zero (Bogdanov-Takens) singularities for Sys. (1.4) and their corresponding dynamical behaviors are investigated by using 𝑘 and 𝜏 as bifurcation parameters (where 𝑘 is defined in Section 2). The discussion of zero-Hopf singularity will be addressed in a coming paper.

The rest of this manuscript is organized as follows. In Section 2, a detailed presentation is given for the distribution of eigenvalues of the linear part of Sys. (1.4) at an equilibrium point in the (𝑘,𝜏)-parameter space. In Section 3, the theory of center manifold reduction for general delayed differential equations (DDEs) is briefly introduced. In Sections 4 and 5, center manifold reduction is performed for Sys. (1.4); and hence, the normal forms for simple-zero and double-zero singularities are obtained on the center manifold, respectively. In Section 6, the normal forms for the double-zero singularity are used to predict the bifurcation diagrams such as Hopf, homoclinic, and double limit cycle bifurcations for the original Sys. of (1.4). Finally in Section 7, some numerical simulations are presented to confirm the theoretical results.

2. Distribution of Eigenvalues

Throughout the rest of this paper, we assume that

𝛼,𝛽>0,𝑞,𝛾∈(0,1),and𝐼(𝑠)isanonlinear𝐶⁴function, (2.1) and that (𝑌^∗,𝐾^∗) is an equilibrium point of Sys. (1.4). Let 𝐼^∗=𝐼(𝑌^∗), 𝑢₁=𝑌−𝑌^∗, 𝑢₂=𝐾−𝐾^∗, and 𝑖(𝑠)=𝐼(𝑠+𝑌^∗)−𝐼^∗. Then Sys. (1.4) can be transformed as

𝑑𝑢₁(𝑡)𝑑𝑡=𝛼[𝑖(𝑢₁(𝑡))−𝛽𝑢₂(𝑡)−𝛾𝑢₁(𝑡)],𝑑𝑢₂(𝑡)𝑑𝑡=𝑖(𝑢₁(𝑡−𝜏))−𝛽𝑢₂(𝑡−𝜏)−𝑞𝑢₂(𝑡). (2.2) Let the Taylor expansion of 𝑖 at 0 be

𝑖(𝑢)=𝑘𝑢+𝑖⁽²⁾𝑢²+𝑖⁽³⁾𝑢³+𝒪(|𝑢|⁴), (2.3) where

𝑘=𝑖^′(0)=𝐼^′(𝑌^∗),𝑖⁽²⁾=12𝑖^′′1(0)=2𝐼^′′(𝑌^∗),𝑖⁽³⁾=1𝑖3!^′′′1(0)=𝐼3!^′′′(𝑌^∗). (2.4) The linear part of Sys. (2.2) at (0,0) is

𝑑𝑢₁(𝑡)𝑑𝑡=𝛼[(𝑘−𝛾)𝑢₁(𝑡)−𝛽𝑢₂(𝑡)],𝑑𝑢₂(𝑡)𝑑𝑡=𝑘𝑢₁(𝑡−𝜏)−𝛽𝑢₂(𝑡−𝜏)−𝑞𝑢₂(𝑡), (2.5) and the corresponding characteristic equation is

Δ(𝜆)≡𝜆²+𝐴𝜆+𝐵+(𝛽𝜆+𝐶)𝑒^−𝜆𝜏=0, (2.6) where

𝐴=𝑞−𝛼(𝑘−𝛾),𝐵=−𝛼𝑞(𝑘−𝛾),𝐶=𝛼𝛽𝛾. (2.7) For 𝜏=0, (2.6) becomes

𝜆²+(𝐴+𝛽)𝜆+𝐵+𝐶=0. (2.8) Define

𝑘^∗=𝛽𝛾𝑞+𝛾,𝑘^∗∗=𝑞+𝛽𝛼+𝛾. (2.9)

Theorem 2.1. Let 𝜏=0. If 𝑘<min{𝑘^∗,𝑘^∗∗}, then all roots of (2.8) have negative real parts, and hence (𝑌^∗,𝐾^∗) is asymptotically stable. If 𝑘>min{𝑘^∗,𝑘^∗∗}, then (2.8) has a positive root and a negative root, and hence, (𝑌^∗,𝐾^∗) is unstable.

Now assume 𝜏>0. Clearly Δ(0)=0 if and only if 𝑘=𝑘^∗. Next we always assume that 𝑘=𝑘^∗. It is easy to attain

Δ^′(𝜆)=2𝜆+𝑞−𝛼𝛽𝛾𝑞+𝛽𝑒^−𝜆𝜏−(𝛽𝜆+𝐶)𝜏𝑒^−𝜆𝜏,Δ^′′(𝜆)=2−2𝛽𝜏𝑒^−𝜆𝜏+𝛽𝜏²𝜆𝑒^−𝜆𝜏+𝐶𝜏²𝑒^−𝜆𝜏. (2.10) Define 𝜏^∗=(𝑞²+𝑞𝛽−𝛼𝛽𝛾)/𝛼𝛽𝛾𝑞. Then we have that,

Δ^′(0)=𝛼𝛽𝛾𝑞(𝜏^∗−𝜏),Δ^′′(0)|_𝜏=𝜏∗=𝑞⁴−𝛽²𝑞²+𝛼²𝛽²𝛾²𝛼𝛽𝛾𝑞². (2.11) Define

𝑓(𝑥)=𝑥²+𝛽𝑥−𝛼𝛽𝛾,𝑔(𝑥)=𝑥²−𝛽²𝑥+𝛼²𝛽²𝛾². (2.12) Hence if 𝑓(𝑞)≤0, 𝜏^∗≤0, and hence Δ^′(0)<0, and if 𝑓(𝑞)>0, 𝜏^∗>0, and hence Δ^′(0)=0 if and only if 𝜏=𝜏^∗. Also Δ^′′(0)|_𝜏=𝜏∗≠0 if and only if 𝑔(𝑞²)≠0. Thus we obtain the following result.

Lemma 2.2. Suppose that 𝑘=𝑘^∗. Then the following are considered. (i)If 𝜏^∗≤0, then (2.6) has a simple root 0 for all 𝜏>0. (ii)Let 𝜏^∗>0. Then the following are given. (a)Equation (2.6) has a simple root 0 if and only if 𝜏≠𝜏^∗, (b)Equation (2.6) has a double root 0 if and only if 𝜏=𝜏^∗ and 𝑔(𝑞²)≠0.

Let 𝜔𝑖 (𝜔>0) be a purely imaginary root of (2.6). After plugging it into (2.6) and separating the real and imaginary parts, we have that

𝜔²𝑞+𝛼𝛽𝛾=𝛼𝛽𝛾cos(𝜔𝜏)+𝛽𝜔sin(𝜔𝜏),₂−𝛼𝛽𝛾𝑞𝜔=𝛼𝛽𝛾sin(𝜔𝜏)−𝛽𝜔cos(𝜔𝜏). (2.13) Adding squares of two equations yields 𝜔²+𝑔(𝑞²)𝑞²=0. (2.14) Then (2.14) has a nonzero solution if and only if 𝑔(𝑞²)<0 and does not have a nonzero solution if and only if 𝑔(𝑞²)≥0. If 𝑔(𝑞²)<0, from (2.14), we solve 𝜔 as follows:

𝜔=𝜔₀≡1𝑞−𝑔(𝑞²), (2.15) and from (2.13), we solve cos(𝜔₀𝜏), sin(𝜔₀𝜏) as:

cos(𝜔₀𝜏)=−𝑞²𝜔²₀+𝛼𝛽𝛾𝜔²₀+𝑞𝛼𝛾(𝛼𝛽𝛾+𝜔²₀)𝑞𝛽(𝛼²𝛾²+𝜔²₀)≡𝑎,sin(𝜔₀𝑞𝜏)=²𝛼𝛾𝜔₀−𝛼²𝛽𝛾²𝜔₀+𝑞𝛼𝛽𝛾𝜔₀+𝑞𝜔³₀𝑞𝛽(𝛼²𝛾²+𝜔²₀)≡𝑏. (2.16) Define

𝛿=arccos𝑎,if𝑏≥0,2𝜋−arccos𝑎,if𝑏<0. (2.17) From (2.16), we obtain

𝜏=𝜏_𝑗≡1𝜔₀(𝛿+2𝑗𝜋),𝑗=0,1,2,…. (2.18)

Clearly if 𝛽>2𝛼𝛾, then 𝑔(𝑥)=0 has two positive roots, and if 𝛽≤2𝛼𝛾, then 𝑔(𝑥)≥0. Now, under 𝑘=𝑘^∗, we impose the following conditions:

(H1)𝛽≤2𝛼𝛾, 𝜏^∗≤0, (H2)𝛽≤2𝛼𝛾, 𝜏^∗>0, 𝜏≠𝜏^∗, (H3)𝛽≤2𝛼𝛾, 𝜏^∗>0, 𝜏=𝜏^∗, (H4)𝛽>2𝛼𝛾, 𝜏^∗>0, 𝜏≠𝜏^∗,𝑔(𝑞²)≥0,(H5)𝛽>2𝛼𝛾, 𝜏^∗>0, 𝜏≠𝜏^∗,𝑔(𝑞²)<0,(H6)𝛽>2𝛼𝛾, 𝜏^∗>0, 𝜏=𝜏^∗,𝑔(𝑞²)≥0,(H7)𝛽>2𝛼𝛾, 𝜏^∗>0, 𝜏=𝜏^∗,𝑔(𝑞²)<0.

Based on Lemma 2.2, we have the following result.

Lemma 2.3. Suppose that 𝑘=𝑘^∗ and 0<𝑞<1. Then the following are obtained.(i)Under one of the conditions (H1), (H2), and (H4), (2.6) has a simple zero root and does not have other roots in the imaginary axis.(ii)Under the condition (H5), (2.6) has a simple zero root and a pair of purely imaginary roots ±𝜔₀𝑖 in the imaginary axis if 𝜏=𝜏_𝑗, 𝑗=0,1,2,….(iii)Under one of the conditions (H3) and (H6), then (2.6) has a double root 0 and does not have other roots in the imaginary axis. (iv)Under the condition (H7), (2.6) has a double zero root and a pair of purely imaginary roots ±𝜔₀𝑖 in the imaginary axis if 𝜏^∗=𝜏_𝑗 for some 𝑗.

Now we use the roots of 𝑓(𝑥)=0, 𝑔(𝑥)=0 to give a more detailed discussion for the roots of (2.6). Define

𝑞₀=12(−𝛽+𝛽²𝑞+4𝛼𝛽𝛾),₁=12(𝛽²−𝛽⁴−4𝛼²𝛽²𝛾²𝑞),₂=12(𝛽²+𝛽⁴−4𝛼²𝛽²𝛾²). (2.19) Clearly 𝑞₀ is the positive root of 𝑓(𝑥)=0 and 𝑞₁, 𝑞₂ are two positive roots of 𝑔(𝑥)=0 if 𝛽>2𝛼𝛾. Note that 𝑓(𝑥)≤0 if 0<𝑥≤𝑞₀,and𝑓(𝑥)>0 if 𝑥>𝑞₀, 𝑔(𝑥)≥0 if 0<𝑥≤𝑞₁, or 𝑥≥𝑞₂, then 𝑔(𝑥)<0 if 𝑞₁<𝑥<𝑞₂. Also note that as well as if 𝛽>2𝛼𝛾, 𝑞²₀<𝑞₁. In fact it is based on the following calculation:

𝑞₁−𝑞²₀=12(𝛽²−𝛽⁴−4𝛼²𝛽²𝛾²1)−4(−𝛽+𝛽²+4𝛼𝛽𝛾)²=𝛽2(𝛽²+4𝛼𝛽𝛾−𝛽²−4𝛼²𝛾²=−2𝛼𝛾)2𝛼𝛽𝛾(𝛽−𝛽²−4𝛼²𝛾²)𝛽²+4𝛼𝛽𝛾+𝛽²−4𝛼²𝛾²+2𝛼𝛾>0. (2.20) Thus for 𝛽>2𝛼𝛾, we always have 𝑞₀<√𝑞₁<√𝑞₂. Noting that 𝑞∈(0,1), we have the following result.

Lemma 2.4. Let 𝛽>2𝛼𝛾. Then the following are given.(i)Suppose that 𝑞₀≥1. Then for 0<𝑞<1, then (2.6) has a simple zero root and does not have roots in the imaginary axis. (ii)Suppose that 𝑞₀<1≤√𝑞₁<√𝑞₂. If 0<𝑞≤𝑞₀, then (2.6) has a simple zero root and does not have roots in the imaginary axis. And if 𝑞₀<𝑞<1, (2.6) has a double zero root and does not have roots in the imaginary axis. (iii)Suppose that 𝑞₀<√𝑞₁<1<√𝑞₂. If 0<𝑞≤𝑞₀, then (2.6) has a simple zero root and does not have roots in the imaginary axis. If 𝑞₀<𝑞≤√𝑞₁, then (2.6) has a double zero root and does not have roots in the imaginary axis. And if √𝑞₁<𝑞<1, then (2.6) has a double zero root and has a pair of purely imaginary roots. (iv)Suppose that √𝑞₂≥1. Then if 0<𝑞≤𝑞₀, then (2.6) has a simple zero root and does not have roots in the imaginary axis. If 𝑞₀<𝑞≤√𝑞₁, then (2.6) has a double zero root and does not have roots in the imaginary axis. If √𝑞₁<𝑞<√𝑞₂, then (2.6) has a double zero root and has a pair of purely imaginary roots when 𝜏^∗=𝜏_𝑗 for some 𝑗. And if √𝑞₂≤𝑞<1, (2.6) has a double zero root and does not have a pair of purely imaginary roots.

Define 𝜆(𝜏)=𝜎(𝜏)+𝑖𝜔(𝜏) to be the root of (2.6) such that 𝜎(𝜏_𝑗)=0 and 𝜔(𝜏_𝑗)=𝜔₀. Then we have the following result.

Lemma 2.5. Suppose that 𝑘=𝑘^∗ and 𝑔(𝑞²)<0. Then 𝜎^′(𝜏_𝑗)>0.

Proof. Differentiating (2.6) with respect to 𝜏 yields (𝑑𝜆)𝑑𝜏⁻¹=[2𝜆+𝑞−𝛼(𝑘−𝛾)]𝑒^𝜆𝜏+𝛽−𝜏𝜆𝛽(𝜆+𝛼𝛾)𝜆, (2.21) and a simple calculation gives Re(𝑑𝜆)𝑑𝜏⁻¹|_{𝜏=𝜏𝑗}=𝛼²𝛽²𝛾²+𝑞²(−𝛽²+𝑞²+2𝜔²₀)𝛽²𝑞²(𝛼²𝛾²+𝜔²₀)=𝑞²𝛽²−𝛼²𝛽²𝛾²−𝑞⁴𝛽²𝑞²(𝛼²𝛾²+𝜔²₀), (2.22) which gives SignRe(𝑑𝜆)𝑑𝜏⁻¹|_{𝜏=𝜏𝑗}=Sign(−𝑔(𝑞²))=1, (2.23) thus completing the proof.

Next we discuss the distribution of other roots of (2.6). We need the following lemma due to Ruan and Wei [19].

Lemma 2.6. Consider the exponential polynomial 𝑃(𝜆,𝑒^−𝜆𝜏)=𝑝(𝜆)+𝑞(𝜆)𝑒^−𝜆𝜏, (2.24) where 𝑝, 𝑞 are real polynomials such that deg(𝑞)<deg(𝑝) and 𝜏≥0. As 𝜏 varies, the sum of the order of zeros of 𝑃(𝜆,𝑒^−𝜆𝜏) on the open right half-plane can change only if a zero appears on or crosses the imaginary axis.

Lemma 2.7. Let 𝑘=𝑘^∗ and 𝜏>0. Then, the following are obtained.(i)If 𝑞>𝑞₀, then all roots of (2.6) except 0 and purely imaginary roots have negative real parts,(ii)If 0<𝑞≤𝑞₀, then (2.6) has at least one positive root.

Proof. Note that, for 𝜏=0, if 𝑞>𝑞₀ or 𝑞²+𝑞𝛽>𝛼𝛽𝛾, Δ(𝜆)=0 has a zero root and a negative root. Using Lemmas 2.2 and 2.6, we obtain claim (i). For 𝜏=0, Δ(𝜆)=0 has a zero root and a positive root if 0<𝑞≤𝑞₀ or 𝑞²+𝑞𝛽≤𝛼𝛽𝛾. For 𝜏>0, let 𝑓(𝜆)=Δ(𝜆)𝜆=𝜆+𝐴+𝛽𝑒^−𝜆𝜏+𝐵+𝐶𝑒^−𝜆𝜏𝜆. (2.25) Also noting that 𝐵+𝐶=0 when 𝑘=𝑘^∗, we have that lim_𝜆→0+1𝑓(𝜆)=𝐴+𝛽−𝐶𝜏=𝑞[(𝑞²+𝑞𝛽−𝛼𝛽𝛾)−𝛼𝛽𝛾𝜏]<0, (2.26) and lim_𝜆→∞𝑓(𝜆)=∞. This proves the second part of the lemma and completes the proof of the lemma.

3. Center Manifold Reduction

In this section, we briefly summarize the theory of center manifold reduction for general DDEs. The material is mainly taken from [20, 21]. Consider the following DDE:

𝑑𝑥𝑑𝑡=𝐿(𝜇)𝑥_𝑡+𝐺(𝑥_𝑡,𝜇), (3.1) where 𝑥∈𝐶([−𝜏,0],ℝ^𝑛), 𝜇∈ℝ^𝑝. This equation is equivalent to

𝑑𝑥𝑑𝑡=𝐿(𝜇)𝑥_𝑡+𝐺(𝑥_𝑡,𝜇),𝑑𝜇𝑑𝑡=0, (3.2) which can be written as

𝑑𝑋𝑑𝑡=ℒ𝑋_𝑡+𝐹(𝑋_𝑡), (3.3) where 𝑋=(𝑥,𝜇)^𝑇, 𝐹(𝑋_𝑡)=(𝐺(𝑥_𝑡),0)^𝑇, and ℒ=diag(𝐿,0). Define 𝑋∈𝐶∶=𝐶([−𝜏,0],ℝ^𝑛+p) with supreme norm and 𝑋_𝑡∈𝐶 is defined by 𝑋_𝑡(𝜃)=𝑋(𝑡+𝜃), −𝜏≤𝜃≤0; ℒ∶𝐶→𝐿(ℝ^𝑛+𝑝) is a bounded linear operator; and 𝐹∶𝐶→𝐶 is a 𝐶^𝑘(𝑘≥2) function with 𝐹(0)=0, 𝐷𝐹(0)=0. Consider the following linear system:

𝑋̇(𝑡)=ℒ𝑋_𝑡. (3.4) Since ℒ is a bounded linear operator, then ℒ can be represented by a Riemann-Stieltjes integral

ℒ𝜑=⁰_−𝜏𝑑𝜂(𝜃)𝜑(𝜃),∀𝜑∈𝐶, (3.5) by the Riesz representation theorem, where 𝜂(𝜃) (𝜃∈[−𝜏,0]) is an (𝑛+𝑝)×(𝑛+𝑝) matrix function of bounded variation. Let 𝒜₀ be the infinitesimal generator for the solution semigroup defined by Sys. (3.4) such that

𝒜₀𝜑=𝜑̇,𝐷(𝒜₀)={𝜑∈𝐶¹([−𝜏,0],ℝ^𝑛+𝑝)∶𝜑̇(0)=⁰_−𝜏𝑑𝜂(𝜃)𝜑(𝜃)}. (3.6) Define the bilinear form between 𝐶 and 𝐶^∗=𝐶([0,𝜏],ℝ^{(𝑛+𝑝)∗}) (where ℝ^{(𝑛+𝑝)∗} is the space of all row (𝑛+𝑝)-vectors) by

⟨𝜓,𝜑⟩=𝜓(0)𝜑(0)−⁰_−𝜏^𝜃₀𝜓(𝜉−𝜃)𝑑𝜂(𝜃)𝜑(𝜉)𝑑𝜉,∀𝜓∈𝐶^∗,∀𝜑∈𝐶. (3.7) The adjoint of 𝒜₀ is defined by 𝒜^∗₀ as

𝒜^∗₀𝜓=−𝜓̇,𝐷(𝒜^∗₀)={𝜑∈𝐶¹([0,𝜏],ℝ^{(𝑛+𝑝)∗})∶𝜓̇(0)=−⁰_−𝜏𝜓(−𝜃)𝑑𝜂(𝜃)}. (3.8) In our setting, (3.3) has 𝑝 trivial components. Assume that the characteristic equation of (3.3) has eigenvalue zero with multiplicity 2𝑝 and all other eigenvalues have negative real parts. Then ℒ has a generalized eigenspace 𝑃 which is invariant under the flow (3.4). Let 𝑃^∗ be the space adjoint with 𝑃 in 𝐶^∗. Then 𝐶 can be decomposed as 𝐶=𝑃⊕𝑄 where 𝑄={𝜑∈𝐶∶⟨𝜓,𝜑⟩=0, ∀𝜓∈𝑃^∗}. Choose the bases Φ and Ψ for 𝑃 and 𝑃^∗, respectively, such that

⟨Ψ,Φ⟩=𝐼,Φ̇=Φ𝐽,Ψ̇=−𝐽Ψ, (3.9) where 𝐽 is Jordan matrix associated with the eigenvalue 0.

To consider Sys. (3.3), we need to enlarge the space 𝐶 to the following 𝐵𝐶:

𝐵𝐶={𝜑∶[−𝜏,0]→ℝ^𝑛+𝑝∶𝜑iscontinuouson[−𝜏,0),∃lim_𝜃→0−𝜑(𝜃)∈ℝ^𝑛+𝑝}. (3.10) The elements of 𝐵𝐶 can be expressed as 𝜓=𝜑+𝑋₀𝛼 with 𝜑∈𝐶,𝛼∈ℝ^𝑛+𝑝, and

𝑋₀(𝜃)=0,−𝜏≤𝜃<0,𝐼,𝜃=0, (3.11) where 𝐼 is the 𝑛×𝑛 identity matrix. Define the projection 𝜋∶𝐵𝐶→𝑃 by

𝜋(𝜑+𝑋₀𝛼)=Φ[(Ψ,𝜑)+Ψ(0)𝛼]. (3.12) Then the enlarged phase space 𝐵𝐶 can be decomposed as 𝐵𝐶=𝑃⊕ker𝜋. Let 𝑋=Φ𝑥+𝑦 with 𝑥∈ℝ^2𝑝 and 𝑦∈𝑄¹={𝜑∈𝑄∶𝜑̇∈𝐶}. Then (3.3) can be decomposed as

̇̇𝑥=𝐽𝑥+Ψ(0)𝐹(Φ𝑥+𝑦),𝑦=𝒜_𝑄¹𝑦+(𝐼−𝜋)𝑋₀𝐹(Ψ𝑥+𝑦), (3.13) where 𝒜 is an extension of the infinitesimal generator 𝒜₀ from 𝐶¹ to 𝐵𝐶, defined by

𝒜₀𝜑=𝜑̇+𝑋₀̇̇[𝐿𝜑−𝜑(0)]=𝜑,−1≤𝜃<0,_0−𝜏𝑑𝜂(𝑡)𝜑(𝑡),𝜃=0, (3.14) for 𝜑∈𝐶¹ and its adjoint by 𝒜^∗ is defined by

𝒜^∗−̇𝜓=𝜓,0<𝑠≤𝜃,_0−𝜏𝜓(−𝜃)𝑑𝜂(𝜃),𝑠=0, (3.15) for 𝜓∈𝐶^1∗. Let 𝐹(𝑣)=∑_𝑗≥2(1/𝑗!)𝐹_𝑗(𝑣). Then Sys. (3.13) becomes

𝑥̇=𝐽𝑥+_𝑗≥21𝑓𝑗!^1𝑗̇(𝑥,𝑦),𝑦=𝒜_𝑄¹𝑦+_𝑗≥21𝑓𝑗!^2𝑗(𝑥,𝑦), (3.16) where

𝑓¹_𝑗(𝑥,𝑦)=Ψ(0)𝐹_𝑗(Φ𝑥+𝑦),𝑓²_𝑗(𝑥,𝑦)=(𝐼−𝜋)𝑋₀𝐹_𝑗(Φ𝑥+𝑦). (3.17) On the center manifold, (3.16) can be approximated as

𝑥̇=𝐽𝑥+_𝑗≥21𝑓𝑗!^1𝑗(𝑥,0). (3.18)

4. Simple-Zero Singularity

In this section, we assume that the condition (H2) holds. From the definition of 𝜏^∗, we know that 𝜏^∗>0 if and only if 𝑞>𝑞₀. Therefore (H2) is equivalent to

𝑘=𝑘^∗,𝑞>𝑞₀,𝜏>0,𝜏≠𝜏^∗. (4.1) From (ii) of Lemma 2.4 and (ii) of Lemma 2.7, we know that, at (0,0), the characteristic equation of the linear part of Sys. (2.5) has a simple zero root and the rest of roots have negative parts. We treat 𝑘 as a bifurcation parameter near 𝑘^∗.

Set 𝐶∶=𝐶([−𝜏,0],ℝ³), 𝐶^∗∶=𝐶([0,𝜏],ℝ^3∗). Let 𝜇=𝑘−𝑘^∗. Then Sys. (2.5) can be rewritten as

𝑑𝑢₁𝑑𝑡=𝛼[𝛽𝛾𝑞𝑢₁(0)−𝛽𝑢₂(0)+𝜇𝑢₁(0)+𝑖⁽²⁾𝑢²₁(0)+𝑖⁽⁰⁾𝑢³₁(𝑡)]+𝒪(|𝜇||𝑢|²+|𝑢|⁴),𝑑𝑢₂𝑑𝑡=𝑘^∗𝑢₁(−𝜏)−𝑞𝑢₂(𝑡)+𝜇𝑢₁(−𝜏)−𝛽𝑢₂(−𝜏)+𝑖⁽²⁾𝑢²₁(−𝜏)+𝑖⁽³⁾𝑢³₁(−𝜏)+𝒪(|𝜇||𝑢|²+|𝑢|⁴),𝑑𝜇𝑑𝑡=0. (4.2) The linearization of Sys. (4.2) at (0,0,0) is

𝑑𝑢₁=𝑑𝑡𝛼𝛽𝛾𝑞𝑢₁(0)−𝛼𝛽𝑢₂(0),𝑑𝑢₂𝑑𝑡=𝑘^∗𝑢₁(−𝜏)−𝑞𝑢₂(0)−𝛽𝑢₂(−𝜏),𝑑𝜇𝑑𝑡=0. (4.3) Let 𝜂(𝜃)=𝔸𝛿(𝜃)+𝔹𝛿(𝜃+𝜏) where

𝔸=𝛼𝛽𝛾𝑞𝑘−𝛼𝛽00−𝑞0000,𝔹=000_∗−𝛽0000. (4.4) Let 𝑋=(𝑢₁,𝑢₂,𝜇)^𝑇 and

𝐹(𝑋_𝑡)=𝛼𝜇𝑢₁(0)+𝛼𝑖⁽²⁾𝑢²₁(0)+𝛼𝑖⁽³⁾𝑢³₁(0)+𝒪(|𝜇||𝑢|²+|𝑢|⁴)𝜇𝑢₁(−𝜏)+𝑖⁽²⁾𝑢²₁(−𝜏)+𝑖⁽³⁾𝑢³₁(−𝜏)+𝒪(|𝜇||𝑢|²+|𝑢|⁴)0. (4.5) Define

𝐿𝜑=⁰_−𝜏𝑑𝜂(𝜃)𝜑(𝜃),∀𝜑∈𝐶. (4.6) Then Sys. (4.2) becomes

𝑋̇(𝑡)=𝐿𝑋_𝑡+𝐹(𝑋_𝑡). (4.7) From (3.7), the bilinear form can be expressed as

⟨𝜓,𝜑⟩=𝜓(0)𝜑(0)+⁰_−𝜏𝜓(𝜉+𝜏)𝔹𝜑(𝜉)𝑑𝜉. (4.8) It is not hard to see that the infinitesimal generator 𝒜∶𝐶¹→𝐵𝐶 is given by

𝒜𝜑=𝜑̇+𝑋₀̇̇[𝐿𝜑−𝜑(0)]=𝜑,−𝜏≤𝜃<0,𝔸𝜑(0)+𝔹𝜑(−𝜏),𝜃=0, (4.9) for 𝜑∈𝐶¹ and its adjoint 𝒜^∗ by

𝒜^∗−̇𝜓=𝜓,0<𝑠≤𝜃,𝜓(0)𝔸+𝜓(𝜏)𝔹,𝑠=0, (4.10) for 𝜓∈𝐶^1∗.

Next we obtain the bases for the center space 𝑃 and its adjoint space 𝑃^∗, respectively. Let 𝒜𝜑=0 for 𝜑∈𝐶¹, that is,

𝜑̇(𝜃)=0for−𝜏≤𝜃<0,𝔸𝜑(0)+𝔹𝜑(−𝜏)=0for𝜃=0. (4.11) then we know that 𝜑 is a constant vector (𝑎₁,𝑎₂,𝑎₃)^𝑇∈ℝ³⧵{0} such that

(𝔸+𝔹)(𝑎₁,𝑎₂,𝑎₃)^𝑇=0. (4.12) Then we have two linearly independent solutions 𝜑₁=(𝑞,𝛾,0)^𝑇, 𝜑₂=(0,0,1)^𝑇 which are bases for the center space 𝑃. Let Φ=(𝜑₁,𝜑₂).

Similarly, let 𝒜^∗𝜓=0 for 𝜓∈𝐶^1∗, that is,

−𝜓̇(𝑠)=0for0<𝑠≤𝜏,𝜓(0)𝔸+𝜓(𝜏)𝔹=0for𝑠=0, (4.13) then we know that 𝜓 is a constant vector (𝑏₁,𝑏₂,𝑏₃)∈ℝ^3∗⧵{0} such that (𝑏₁,𝑏₂,𝑏₃)(𝔸+𝔹)=0. (4.14) From this we have two linearly independent solutions 𝜓₁=(−(𝑞+𝛽),𝛼𝛽,0) and 𝜓₂=(0,0,1) which are bases for the center space 𝑃^∗. Let Ψ=(𝑟𝜓₁,𝜓₂)^𝑇 with 𝑟 being determined such that ⟨𝜓₁,𝜑₁⟩=1. In fact

1𝑟=𝑞𝛼𝛽𝛾(𝜏−𝜏^∗). (4.15) Clearly 𝑟 is well defined since 𝜏−𝜏^∗≠0. It is not hard to check that Φ̇=Φ𝐽, Ψ̇=−𝐽Ψ and ⟨Ψ,Φ⟩=𝐼, where 𝐽=(⁰⁰₀₀).

Let 𝑢=Φ𝑥+𝑦. Then Sys. (4.2) can be decomposed as

̇̇𝑥=Ψ𝐹(Φ𝑥+𝑦),𝑦=𝐴_𝑄¹𝑦+(𝐼−𝜋)𝑋₀𝐹(Φ𝑥+𝑦). (4.16) Write 𝑥=(𝑥₁,𝜇). Note that

Ψ(0)𝐹(Φ𝑥)=−𝑟𝛼𝑞²(𝜇𝑥₁+𝑞𝑖⁽²⁾𝑥²₁+𝑞²𝑖⁽³⁾𝑥³₁)0+h.o.t.. (4.17) Here h.o.t. represents higher-order terms. Thus, for sufficiently small 𝜇, on the center manifold, if 𝑖⁽²⁾≠0, then Sys. (4.2) becomes

𝑥̇₁=−𝑟𝛼𝑞²𝜇𝑥₁−𝑟𝛼𝑞³𝑖⁽²⁾𝑥²₁̇+h.o.t.,𝜇=0. (4.18) If 𝑖⁽²⁾=0 and 𝑖⁽³⁾≠0, then Sys. (4.2) can be transformed into the following form:

𝑥̇₁=−𝑟𝛼𝑞²𝜇𝑥₁−𝑟𝛼𝑞⁴𝑖⁽³⁾𝑥³₁̇+h.o.t.,𝜇=0. (4.19) Thus we have the following results.

Theorem 4.1. Let 𝜇 be small. Then consider what follows.(i)Suppose that 𝜇=0. Then if 𝑖⁽²⁾≠0, the equilibrium (𝑌^∗,𝐾^∗) is unstable, and if 𝑖⁽²⁾=0 and 𝑖⁽³⁾≠0, then the equilibrium (𝑌^∗,𝐾^∗) is asymptotically stable for (𝜏−𝜏^∗)𝑖⁽³⁾>0 and unstable if (𝜏−𝜏^∗)𝑖⁽³⁾<0.(ii)The equilibrium (𝑌^∗,𝐾^∗) is asymptotically stable if (𝜏−𝜏^∗)𝜇>0 and unstable if (𝜏−𝜏^∗)𝜇<0.(iii)At (𝑌^∗,𝐾^∗,𝑘^∗), Sys. (1.4) undergoes a transcritical bifurcation if 𝑖⁽²⁾≠0 and a pitchfork bifurcation if 𝑖⁽²⁾=0 and 𝑖⁽³⁾≠0.

5. Double-Zero Singularity

In this section, we assume that one of the conditions (H3) and (H6) holds and 𝑔(𝑞²)>0, or equivalently, as

𝑘=𝑘^∗,𝜏=𝜏^∗,𝑞>𝑞₀,𝑔(𝑞²)>0. (5.1) From Section 2, we can see that, at (0,0), the characteristic equation of Sys. (2.5) has a double root 0 and all other roots have negative real parts if 𝑘=𝑘^∗ and 𝜏=𝜏^∗. We treat (𝑘,𝜏) as a bifurcation parameter near (𝑘^∗,𝜏^∗).

By scaling 𝑡→𝑡/𝜏, Sys. (2.2) can be written as

𝑑𝑢₁(𝑡)𝑑𝑡=𝛼𝜏(𝑘−𝛾)𝑢₁(𝑡)−𝛼𝛽𝜏𝑢₂(𝑡)+𝛼𝜏𝑖⁽²⁾𝑢²₁(𝑡)+𝜏𝑖⁽³⁾𝑢³₁(𝑡)+𝒪(|𝑢₁|⁴),𝑑𝑢₂(𝑡)𝑑𝑡=𝜏𝑘𝑢₁(𝑡−1)−(𝑞+𝛽)𝜏𝑢₂(𝑡)+𝜏𝑖⁽²⁾𝑢²₁(𝑡−1)+𝜏𝑖⁽³⁾𝑢³₁(𝑡−1)+𝒪(|𝑢₁|⁴). (5.2) Let 𝐶∶=𝐶([−1,0],ℝ⁴),𝐶^∗∶=𝐶([0,1],ℝ^4∗). Let 𝜇₁=𝑘−𝑘^∗, 𝜇₂=𝜏−𝜏^∗. Then on 𝐶 we have

𝑑𝑢₁𝑑𝑡=𝛼[𝛽𝛾𝑞𝜏^∗𝑢₁(0)−𝛽𝜏^∗𝑢₂(0)+𝜏^∗𝜇₁𝑢₁+(0)𝛽𝛾𝑞𝜇₂𝑢₁(0)−𝛽𝜇₂𝑢₂(0)+𝜏^∗𝑖⁽²⁾𝑢²₁(0)+𝑖⁽²⁾𝜇₂𝑢²₁(0)+𝑖⁽³⁾𝜏^∗𝑢³₁(0)+𝑖⁽³⁾𝜇₂𝑢³₁(0)]+𝒪(|𝜇|²|𝑢|+|𝜇||𝑢|⁴),𝑑𝑢₂=𝑑𝑡𝛽𝛾𝑞𝜏^∗𝑢₁(𝑡−1)−𝑞𝜏^∗𝑢₂(𝑡)−𝛽𝜏^∗𝑢₂(−1)+𝜏^∗𝜇₁𝑢₁(−1)+𝛽𝛾𝑞𝜇₂𝑢₁(−1)−𝑞𝜇₂𝑢₂(0)−𝛽𝜇₂𝑢₂(−1)+𝑖⁽²⁾𝜏^∗𝑢²₁(−1)+𝑖⁽²⁾𝜇₂𝑢²₁(−1)+𝑖⁽³⁾𝜏^∗𝑢³₁(−1)+𝑖⁽³⁾𝜇₂𝑢³₁(−1)+𝒪(|𝜇|²|𝑢|+|𝜇||𝑢|⁴),𝑑𝜇₁𝑑𝑡=0,𝑑𝜇₂𝑑𝑡=0. (5.3) The linearization of Sys. (5.3) at (0,0,0,0) is

𝑑𝑢₁(𝑡)=𝑑𝑡𝛼𝛽𝛾𝑞𝜏^∗𝑢₁(0)−𝛼𝛽𝜏^∗𝑢₂(0),𝑑𝑢₂(𝑡)𝑑𝑡=𝑘^∗𝜏^∗𝑢₁(−1)−𝑞𝜏^∗𝑢₂(0)−𝛽𝜏^∗𝑢₂(−1),𝑑𝜇₁𝑑𝑡=0,𝑑𝜇₂𝑑𝑡=0. (5.4) Let

𝜂(𝜃)=𝔸𝛿(𝜃)+𝔹𝛿(𝜃+1), (5.5) where

𝔸=𝛼𝛽𝛾𝑞𝜏^∗−𝛼𝛽𝜏^∗000−𝑞𝜏^∗00𝑘00000000,𝔹=0000_∗𝜏^∗−𝛽𝜏^∗0000000000. (5.6) Define

ℒ𝜑=⁰₋₁𝑑𝜂(𝜃)𝜑(𝜃),∀𝜑∈𝐶. (5.7) Let 𝐶¹=𝐶¹([−1,0],ℝ⁴). Let 𝑋=(𝑢₁,𝑢₂,𝜇₁,𝜇₂)^𝑇 and 𝐹(𝑋_𝑡)=(𝐹¹,𝐹²,0,0)^𝑇 where

𝐹¹=𝛼[𝜏^∗𝜇₁𝑢₁(0)+𝛽𝛾𝑞𝜇₂𝑢₁(0)−𝛽𝜇₂𝑢₂(0)+𝑖⁽²⁾𝜏^∗𝑢²₁(0)+𝑖⁽³⁾𝜏^∗𝑢³₁(0)]+𝒪(|𝜇|²|𝑢|+|𝜇||𝑢|⁴𝐹),₂=𝜏^∗𝜇₁𝑢₁(−1)+𝛽𝛾𝑞𝜇₂𝑢₁(−1)−𝑞𝜇₂𝑢₂(0)−𝛽𝜇₂𝑢₂(−1)+𝑖⁽²⁾𝜏^∗𝑢²₁(0)+𝑖⁽³⁾𝜏^∗𝑢³₁(−1)+𝒪(|𝜇|²|𝑢|+|𝜇||𝑢|⁴). (5.8) Then Sys. (5.3) can be transformed into

𝑋̇(𝑡)=ℒ𝑋_𝑡+𝐹(𝑋_𝑡). (5.9) Let 𝐶^∗=𝐶([0,1],ℝ^4∗). From (3.7), the bilinear inner product between 𝐶 and 𝐶^∗ can be expressed by

⟨𝜓,𝜑⟩=𝜓(0)𝜑(0)+⁰₋₁𝜓(𝜉+1)𝔹𝜑(𝜉)𝑑𝜉, (5.10) for 𝜑∈𝐶 and 𝜓∈𝐶^∗. As in Section 4, the infinitesimal generator 𝒜∶𝐶¹→𝐵𝐶 associated with ℒ is given by

𝒜𝜑=𝜑̇+𝑋₀̇̇[ℒ𝜑−𝜑(0)]=𝜑,−1≤𝜃<0,𝔸𝜑(0)+𝔹𝜑(−1),𝜃=0, (5.11) for 𝜑∈𝐶¹ and its adjoint by

𝒜^∗−̇𝜓=𝜓,0<𝑠≤1,𝜓(0)𝔸+𝜓(1)𝔹,𝑠=0, (5.12) for 𝜓∈𝐶^1∗. From Section 2, we know that 0 is an eigenvalue of 𝒜 and 𝒜^∗ with multiplicity 4. Now we compute eigenvectors of 𝒜 and 𝒜^∗ associated with 0, respectively.

Next we obtain the bases for the center space 𝑃 and its adjoint space 𝑃^∗, respectively. Let 𝒜𝜑=0 for 𝜑∈𝐶¹. This means that

𝜑̇(𝜃)=0for−1≤𝜃<0,𝔸𝜑(0)+𝔹𝜑(−1)=0for𝜃=0. (5.13) From this we obtain that 𝜑(𝜃)=𝜑⁰ is a constant vector in ℝ⁴ satisfying

(𝔸+𝔹)𝜑⁰=0. (5.14) This equation has three linearly independent solutions: 𝑎₁=(𝑞,𝛾,0,0)^𝑇, 𝑎₃=(0,0,1,0)^𝑇, 𝑎₄=(0,0,0,1)^𝑇. Let 𝜑⁰₁ be one of those. Suppose that 𝒜𝑎₂=𝜑⁰₁ for 𝑎₂∈𝐶¹, namely,

𝑎̇₂(𝜃)=𝜑⁰₁for−1≤𝜃<0,𝔸𝑎₂(0)+𝔹𝑎₂(−1)=𝜑⁰₁for𝜃=0. (5.15) This implies that there is a constant vector 𝜑⁰₂ in ℝ⁴ such that 𝑎₂(𝜃)=𝜑⁰₁𝜃+𝜑⁰₂ and

ℒ(𝜑⁰₁𝜃+𝜑⁰₂)=𝜑⁰₁. (5.16) Since

ℒ(𝜑⁰₁𝜃+𝜑⁰₂)=ℒ(𝜑⁰₁𝜃)+ℒ(𝜑⁰₂)=−𝐵𝜑⁰₁+(𝔸+𝔹)𝜑⁰₂, (5.17) we have that

(𝔸+𝔹)𝜑⁰₂=(𝐼+𝔹)𝜑⁰₁. (5.18) It is easy to see that (5.18) has no solution if 𝜑⁰₁ is either 𝑎₃ or 𝑎₄. For 𝜑⁰₁=𝑎₁, setting 𝜑⁰₂=(0,𝑙,0,0)^𝑇 in (5.18), we obtain 𝑞𝑙=−²𝛾𝑞²,+𝑞𝛽−𝛼𝛽𝛾 (5.19) and hence 𝑎₂(𝜃)=(𝜃𝑞,𝑙+𝛾𝜃,0,0)^𝑇. Thus we obtain bases 𝑎₁, 𝑎₂, 𝑎₃, 𝑎₄ for the center space 𝑃. Let Φ=(𝑎₁,𝑎₂,𝑎₃,𝑎₄). Then we have that Φ̇=Φ𝐽 where 𝐽=(^01000000_00000000).

Similarly, let 𝒜^∗𝜓₂=0 for 𝜓₂∈𝐶^1∗, that is,

−𝜓̇₂(𝑠)=0for0<𝑠≤1,𝜓₂(0)𝔸+𝜓₂(−1)𝔹=0for𝑠=0, (5.20) which means that 𝜓₂(𝑠)=𝜓⁰₂ is a constant vector 𝜓⁰₂∈ℝ^4∗⧵{0} satisfying

𝜓⁰₂(𝔸+𝔹)=0. (5.21) This equation has three linearly independent solutions: 𝑏₂=(𝑚(𝑞+𝛽),−𝑚𝛼𝛽,0,0), 𝑏₃=(0,0,1,0), 𝑏₄=(0,0,0,1). Asserting that ⟨𝑏₂,𝑎₂⟩=1 gives

𝑚=2(𝑞²+𝑞𝛽−𝛼𝛽𝛾)𝑞⁴−𝑞²𝛽²+𝛼²𝛽²𝛾². (5.22) Let 𝜓⁰₂ be one of 𝑏₂, 𝑏₃, 𝑏₄. Suppose 𝒜^∗𝑏₁=𝜓⁰₂, that is,

−𝑏̇₁(𝑠)=𝜓⁰₂for0<𝑠≤1,𝑏₁(0)𝔸+𝑏₁(1)𝔹=𝜓⁰₂for𝑠=0, (5.23) which implies that there is 𝜓⁰₁∈ℝ^4∗ such that 𝑏₁(𝑠)=−𝜓⁰₂𝑠+𝜓⁰₁ satisfying

ℒ^∗(−𝜓⁰₂𝑠+𝜓⁰₁)=𝜓⁰₂. (5.24) Since

ℒ^∗(−𝜓⁰₂𝑠+𝜓⁰₁)=−ℒ^∗(𝜓⁰₂𝑠)+ℒ^∗(𝜓⁰₁)=−𝜓⁰₂𝔹+𝜓⁰₁(𝔸+𝔹), (5.25) we have

𝜓⁰₁(𝔸+𝔹)=𝜓⁰₂(𝐼+𝔹). (5.26) It is not hard to check that (5.26) has no solution if 𝜓⁰₂=𝑏₃ or 𝑏₄. Letting 𝜓⁰₂=𝑏₂, setting 𝜓⁰₁=(𝑛₁,𝑛₂,0,0) in (5.26) and using ⟨𝑏₁,𝑎₂⟩=0, we can get 𝑛₁ and 𝑛₂:

𝑛₁=2(𝑞+𝛽)(𝑞⁶−3𝑞⁵𝛽−3𝑞⁴𝛽²+𝑞³𝛽³+3𝑞²𝛼²𝛽²𝛾²−3𝑞𝛼²𝛽³𝛾²+2𝛼³𝛽³𝛾³)3(𝑞⁴−𝑞²𝛽²+𝛼²𝛽²𝛾²)²,𝑛₂=2𝛼𝛽(−𝑞²+2𝑞𝛽+𝛼𝛽𝛾)(𝑞²+𝑞𝛽−𝛼𝛽𝛾)²3(𝑞⁴−𝑞²𝛽²+𝛼²𝛽²𝛾²)². (5.27) Hence

𝑏₁(𝑠)=𝜓⁰₁−𝑠𝜓⁰₂=(−𝑚(𝑞+𝛽)𝑠+𝑛₁,𝑚𝛼𝛽𝑠+𝑛₂,0,0). (5.28) Then 𝑏₁, 𝑏₂, 𝑏₃, 𝑏₄ are bases of the center space 𝑃^∗. Let Ψ=(𝑏₁,𝑏₂,𝑏₃,𝑏₄)^𝑇. Then ⟨Ψ,Φ⟩=𝐼, Φ̇=Φ𝐽 and Ψ̇=−𝐽Ψ.

Let 𝑢=Φ𝑥+𝑦, namely,

𝑢₁(𝜃)=𝑞𝑥₁+𝑞𝜃𝑥₂+𝑦₁𝑢(𝜃),₂(𝜃)=𝛾𝑥₁+(𝑙+𝛾𝜃)𝑥₂+𝑦₂(𝜃). (5.29) Then Sys. (5.9) can be decomposed as

̇̇𝑥=𝐽𝑥+Ψ(0)𝐹(Φ𝑥+𝑦),𝑦=𝐴_𝑄¹𝑦+(𝐼−𝜋)𝑋₀𝐹(Φ𝑥+𝑦). (5.30) Write 𝑥=(𝑥₁,𝑥₂,𝜇₁,𝜇₂). Then, on the center manifold, Sys. (5.30) becomes

𝑥̇₁=(𝑎₁₁𝜇₁+𝑎₂₁𝜇₂)𝑥₁+(𝑎₂₁𝜇₁+𝑎₂₂𝜇₂)𝑥₂+𝛼𝑛₁𝜏^∗(𝑖⁽²⁾𝑞²𝑥²₁+𝑖⁽³⁾𝑞³𝜏^∗𝑥³₁)+𝑛₂𝜏^∗×[𝑖⁽²⁾𝑞²(𝑥₁−𝑥₂)²+𝑖⁽³⁾𝑞³𝜏^∗(𝑥₁−𝑥₂)³𝑥̇]+h.o.t.,₂=(𝑏₁₁𝜇₁+𝑏₂₁𝜇₂)𝑥₁+(𝑏₁₂𝜇₁+𝑏₂₂𝜇₂)𝑥₂+𝑚𝛼(𝑞+𝛽)𝜏^∗(𝑖⁽²⁾𝑞²𝑥²₁+𝑖⁽³⁾𝑞³𝜏^∗𝑥³₁)+𝑚𝛼𝛽𝜏^∗×[𝑖⁽²⁾𝑞²(𝑥₁−𝑥₂)²+𝑖⁽³⁾𝑞³𝜏^∗(𝑥₁−𝑥₂)³𝜇̇]+h.o.t.,₁=0,𝜇̇₂=0, (5.31) where

𝑎₁₁=𝑞𝜏^∗(𝛼𝑛₁+𝑛₂),𝑎₂₁=−𝛼((𝑞+𝛽)𝛾−𝑞𝑘^∗)𝑛₁+(−(2𝑞+𝛽)𝛾+𝑞𝑘^∗)𝑛₂,𝑎₁₂=−𝑞𝑛₂𝜏^∗,𝑎₂₂=−(𝑙𝛼𝛽𝑛₁+((𝑞+𝛽)(𝑙−𝛾)+𝑞𝑘^∗)𝑛₂𝑏),₁₁=𝑚𝑞²𝛼𝜏^∗,𝑏₂₁=𝑚𝑞𝛼𝛽𝛾,𝑏₁₂=𝑚𝑞𝛼𝛽𝜏^∗,𝑏₂₂=𝑚𝛼𝛽(−(𝑞+𝛽)𝛾+𝑞𝑘^∗). (5.32)

Next we use techniques of nonlinear transformations in [22] to transform Sys. (5.31) into normal forms. If 𝑖⁽²⁾≠0, then up to the second order, Sys. (5.31) can be written as

𝑥̇₁=(𝑎₁₁𝜇₁+𝑎₂₁𝜇₂)𝑥₁+(𝑎₂₁𝜇₁+𝑎₂₂𝜇₂)𝑥₂+𝛼𝑛₁𝜏^∗𝑖⁽²⁾𝑞²𝑥²₁+𝑛₂𝜏^∗𝑖⁽²⁾𝑞²(𝑥₁−𝑥₂)²𝑥̇+h.o.t.,₂=(𝑏₁₁𝜇₁+𝑏₂₁𝜇₂)𝑥₁+(𝑏₁₂𝜇₁+𝑏₂₂𝜇₂)𝑥₂+𝑚𝛼(𝑞+𝛽)𝜏^∗𝑖⁽²⁾𝑞²𝑥²₁+𝑚𝛼𝛽𝜏^∗𝑖⁽²⁾𝑞²(𝑥₁−𝑥₂)²𝜇̇+h.o.t.,₁=0,𝜇̇₂=0. (5.33) This system can be transformed into the following normal form:

𝑥̇₁=𝑥₂𝑥̇+h.o.t.,₂=𝜌₁𝑥₁+𝜌₂𝑥₂+𝑎₁𝑥²₁+𝑏₁𝑥₁𝑥₂+h.o.t., (5.34) where

𝜌₁=𝑏₁₁𝜇₁+𝑏₂₁𝜇₂,𝜌₂=(𝑎₁₁+𝑏₁₂)𝜇₁+(𝑎₂₁+𝑏₂₂)𝜇₂,𝑎₁=𝑚𝑞³𝛼𝜏^∗𝑖⁽²⁾,𝑏₁=2𝑞²(𝑚𝛼𝛽+𝛼𝑛₁+𝑛₂)𝜏^∗𝑖⁽²⁾. (5.35) Since

|𝜕𝜌𝜕𝜇|=det𝜕𝜌₁𝜕𝜇₁𝜕𝜌₁𝜕𝜇₂𝜕𝜌₂𝜕𝜇₁𝜕𝜌₂𝜕𝜇₂=−𝑚𝑞²𝛼𝛾(𝑚𝛼𝛽²+𝛼𝛽𝑛₁+(𝑞+𝛽)𝑛₂)𝜏^∗=−4𝑞³𝛼³𝛽²𝛾²𝜏^∗(𝑞²+𝑞𝛽−𝛼𝛽𝛾)(𝑞⁴−𝑞²𝛽²+𝛼²𝛽²𝛾²)²≠0, (5.36) we have that (𝜇₁,𝜇₂)→(𝜌₁,𝜌₂) is regular and hence the transversality condition holds.

If 𝑖⁽²⁾=0 and 𝑖⁽³⁾≠0, then up to the third order, Sys. (5.31) becomes

𝑥̇₁=(𝑎₁₁𝜇₁+𝑎₂₁𝜇₂)𝑥₁+(𝑎₂₁𝜇₁+𝑎₂₂𝜇₂)𝑥₂+𝛼𝑛₁𝜏^∗𝑖⁽³⁾𝑞³𝜏^∗𝑥³₁+𝑛₂𝜏^∗𝑖⁽³⁾𝑞³𝜏^∗(𝑥₁−𝑥₂)³𝑥̇+h.o.t.,₂=(𝑏₁₁𝜇₁+𝑏₂₁𝜇₂)𝑥₁+(𝑏₁₂𝜇₁+𝑏₂₂𝜇₂)𝑥₂+𝑚𝛼(𝑞+𝛽)𝜏^∗𝑖⁽³⁾𝑞³𝜏^∗𝑥³₁+𝑚𝛼𝛽𝜏^∗𝑖⁽³⁾𝑞³𝜏^∗(𝑥₁−𝑥₂)³𝜇̇+h.o.t.,₁=0,𝜇̇₂=0. (5.37) This system can be transformed into the following normal form:

𝑥̇₁=𝑥₂𝑥̇+h.o.t.,₂=𝜌₁𝑥₁+𝜌₂𝑥₂+𝑎₂𝑥³₁+𝑏₂𝑥²₁𝑥₂+h.o.t., (5.38) where 𝑎₂=𝑚𝑞⁴𝛼𝜏^∗𝑖⁽³⁾, 𝑏₂=3𝑞³(𝑚𝛼𝛽+𝛼𝑛₁+𝑛₂)𝜏^∗𝑖⁽³⁾.

6. Bifurcation Diagrams

In this section, we will use the truncated systems (5.34) and (5.38) to obtain bifurcation diagrams of Sys. (5.3).

First, we consider the truncated system of (5.34): 𝑥̇₁=𝑥₂,𝑥̇₂=𝜌₁𝑥₁+𝜌₂𝑥₂+𝑎₁𝑥²₁+𝑏₁𝑥₁𝑥₂, (6.1) where 𝑎₂ and 𝑏₂ are in Section 5. Note that (𝑎₁,𝑏₁)→(−𝑎₁,−𝑏₁) under the transformation (𝑥₁,𝑥₂)→(−𝑥₁,−𝑥₂). We may assume that 𝑖⁽²⁾>0. After the change of coordinates

𝑥₁=𝑎₁𝑏²₁(𝜉₁−𝜌₁𝑏₁),𝑥₂𝑎→𝑠²¹𝑏³₁𝜉₂𝑏,𝑡→|¹𝑎₁|𝜏, (6.2) we have (still using 𝑥₁,𝑥₂ for simplicity) that 𝑥̇₁=𝑥₂,𝑥̇₂=𝜈₁+𝜈₂𝑥₁+𝑥²₁+𝑠𝑥₁𝑥₂, (6.3) where 𝜈₁=−(𝑏²₁/𝑎³₁)(𝑏₁𝜌₁−𝑎₁𝜌₂), 𝜈₂=(𝑏₁/𝑎²₁)(𝑏₁𝜌₁−2𝑎₁𝜌₁) and 𝑠=sign𝑎₁𝑏₁=±1. Simple calculation shows that 𝑠=sign(𝛿) where

𝛿=𝑞⁶−3𝑞⁴𝛽²−2𝑞³𝛽³+3𝑞²𝛼²𝛽²𝛾²+2𝛼³𝛽³𝛾³. (6.4) Now take 𝑠=−1, namely 𝛿<0. The complete bifurcation diagrams of Sys. (6.3) can be found in [22]. Here, we just briefly list some results. For (𝜈₁,𝜈₂) small enough, consider the following.

(i)Sys. (6.3) undergoes a fold bifurcation when (𝜈₁,𝜈₂) is on the curves 𝑇⁺={(𝜈₁,𝜈₂)∶4𝜈₁−𝜈²₂=0,𝜈₂>0},𝑇⁻={(𝜈₁,𝜈₂)∶4𝜈₁−𝜈²₂=0,𝜈₂<0}. (6.5)(ii)Sys. (6.3) undergoes a Hopf bifurcation when (𝜈₁,𝜈₂) is on the half-line 𝐻={(𝜈₁,𝜈₂)∶𝜈₁=0,𝜈₂<0}, (6.6) and the Hopf bifurcation gives rise to a stable limit cycle.(iii)Sys. (6.3) undergoes a homoclinic loop bifurcation when (𝜈₁,𝜈₂) is on the curve 𝑃={(𝜈₁,𝜈₂)∶𝜈₁6=−𝜈25²²,𝜈₂<0}. (6.7)

Moreover, when (𝜈₁,𝜈₂) is in the region between the curves 𝐻 and 𝑃, Sys. (6.1) has a unique stable periodic orbit.

For 𝑠=1, under the transformation 𝑡→−𝑡, 𝑥₁→−𝑥₁, we can get Sys. (6.12) whose parametric portrait remains as it was but the cycle becomes unstable. Applying the above results and using the expressions of 𝜈₁,𝜈₂, we obtain the following result regarding Sys. (5.3).

Theorem 6.1. Suppose that 𝑖⁽²⁾>0 and 𝛿<0. For sufficiently small 𝜇₁,𝜇₂, consider the following(i)Sys. (5.3) undergoes a fold bifurcation in the half-lines 𝑇₊={(𝜇₁,𝜇₂)∶𝜇₁=−𝛽𝛾𝑞𝜏^∗𝜇₂,𝜇₂<0},𝑇₋={(𝜇₁,𝜇₂)∶𝜇₁=−𝛽𝛾𝑞𝜏^∗𝜇₂,𝜇₂>0}. (6.8)(ii)Sys. (5.3) undergoes a Hopf bifurcation on the curve 𝐻={(𝜇₁,𝜇₂)∶𝜇₁=−𝛾(2𝑚𝛼𝛽²+2𝛼𝛽𝑛₁+(𝑞+2𝛽)𝑛₂)𝑞(𝑚𝛼𝛽+𝛼𝑛₁+𝑛₂)𝜏^∗𝜇₂+𝒪(𝜇²₂),𝜇₂>0}. (6.9)(iii)Sys. (5.3) undergoes a saddle of homoclinic bifurcation on the curve 𝑃={(𝜇₁,𝜇₂)∶𝜇₁=−𝛾(12𝑚𝛼𝛽²+2𝛼𝛽𝑛₁+(7𝑞+12𝛽)𝑛₂)5𝑞(𝑚𝛼𝛽+𝛼𝑛₁+𝑛₂)𝜏^∗𝜇₂+𝒪(𝜇₂^3/2),𝜇₂>0}. (6.10)Moreover, if (𝜇₁,𝜇₂) is in the region between the curves 𝐻 and 𝑃, Sys. (5.3) has a unique stable periodic orbit.

Next, we consider the truncated system of (5.38): 𝑥̇₁=𝑥₂,𝑥̇₂=𝜌₁𝑥₁+𝜌₂𝑥₂+𝑎₂𝑥³₁+𝑏₂𝑥²₁𝑥₂, (6.11) where 𝑎₂,𝑏₂ are in Section 5. The bifurcation diagrams of this system are more complicated and interesting. We must consider two cases.

Case 1 (𝑖⁽³⁾>0 so that a₂>0). We can assume 𝑏₂<0 which is equivalent to 𝛿<0 in (6.4). Then Sys. (6.11) can be transformed as 𝑥̇₁=𝑥₂,𝑥̇₂=𝜀₁𝑥₁+𝜀₂𝑥₂+𝑥³₁−𝑥²₁𝑥₂, (6.12) where 𝜀₁𝑏=(²𝑎₂)²𝜌₁,𝜀₂𝑏=−²𝑎₂𝜌₂. (6.13) The complete bifurcation diagrams of Sys. (6.12) can be found, for example, in [22–24]. For 𝑏₂>0, under the transformation 𝑡→−𝑡 and 𝑥₁→−𝑥₁, we can get Sys. (6.12). Here, we briefly list some results: for small 𝜀₁, 𝜀₂ as follows. (i)When (𝜀₁,𝜀₂) is in the line 𝑇₁={(𝜀₁,𝜀₂)∶𝜀₁=0,𝜀₂∈ℝ}, (6.14) Sys. (6.12) undergoes a pitchfork bifurcation.(ii)Sys. (6.12) undergoes a stable Hopf bifurcation for the trivial equilibrium point on the half-line 𝐻₁={(𝜀₁,𝜀₂)∶𝜀₂=0,𝜀₁<0}; (6.15)(iii)On the curve 𝐶={(𝜀₁,𝜀₂)∶𝜀₂1=−5𝜀₁+𝒪(|𝜀₁|^3/2),𝜀₁<0}, (6.16)Sys. (6.12) undergoes a heteroclinic bifurcation. Moreover, if (𝜀₁,𝜀₂) is in the region between the curves 𝐻₁ and 𝐶 then Sys. (6.12) has a unique stable periodic orbit.

Applying the above results and using the expressions of 𝜌₁, 𝜌₂, 𝜀₁, 𝜀₂, we obtain the following theorem regarding Sys. (5.3).

Theorem 6.2. Suppose that 𝑖⁽²⁾=0, 𝑖⁽³⁾>0 and 𝛿<0. For sufficiently small 𝜇₁, 𝜇₂, the following are given.(i)Sys. (5.3) undergoes a pitchfork bifurcation in the line 𝑇₁={(𝜇₁,𝜇₂)∶𝜇₁=−𝛽𝛾𝑞𝜏^∗𝜇₂,𝜇₂∈ℝ}. (6.17)(ii)Sys. (5.3) undergoes a stable Hopf bifurcation in the half-line 𝐻₁={(𝜇₁,𝜇₂)∶𝜇₁=𝛾𝑛₂𝜏^∗(𝑚𝛼𝛽+𝛼𝑛₁+𝑛₂)𝜇₂,𝜇₂>0}. (6.18)(iii)Sys. (5.3) undergoes a branch of homoclinic bifurcation on the curve 𝐶={(𝜇₁,𝜇₂)∶𝜇₁=𝛾(3𝑚𝛼𝛽²+3𝛼𝛽𝑛₁+(5𝑞+3𝛽)𝑛₂)2𝑞(𝑚𝛼𝛽+𝛼𝑛₁+𝑛₂)𝜏^∗𝜇₂+𝒪(|𝜇₂|^3/2),𝜇₂>0}. (6.19)Moreover, if (𝜇₁,𝜇₂) is in the region between the curves 𝐻₁ and 𝐶, Sys. (5.3) has a unique stable periodic orbit.

Case 2 (𝑖⁽³⁾<0 so that 𝑎₂<0). We assume 𝑏₂<0, namely 𝛿>0. After rescaling, then Sys. (6.11) can be transformed as 𝑥̇₁=𝑥₂,𝑥̇₂=𝜀₁𝑥₁+𝜀₂𝑥₂−𝑥³₁−𝑥²₁𝑥₂, (6.20) where 𝜀₁𝑏=(¹𝑎₁)²𝜌₁,𝜀₂=𝑏₁𝑎₁𝜌₂. (6.21) For 𝑏₂>0, under the transformation 𝑡→−𝑡,𝑥₁→−𝑥₁, we can get Sys. (6.20). The complete bifurcation diagrams of this system can be found, for example, in [22–24]. Here, we briefly list some results: for small 𝜀₁, 𝜀₂, as follows.(i)When (𝜀₁,𝜀₂) is in the line 𝑇₂={(𝜀₁,𝜀₂)∶𝜀₁=0,𝜀₂∈ℝ}, (6.22) Sys. (6.20) undergoes a pitchfork bifurcation.(ii)When (𝜀₁,𝜀₂) is in the half-line 𝐻₂={(𝜀₁,𝜀₂)∶𝜀₂=𝜀₁,𝜀₁>0}, (6.23) Sys. (6.20) undergoes a stable Hopf bifurcation at 𝐸_1,2 and the bifurcation is subcritical.(iii)When (𝜀₁,𝜀₂) is on the curve 𝐶^′={(𝜀₁,𝜀₂)∶𝜀₂=45𝜀₁+𝒪(|𝜀₁|^3/2),𝜀₁>0}, (6.24) Sys. (6.20) has a unique homoclinic orbit connecting 𝐸₁ and 𝐸₂ and two homoclinic orbits simultaneously at 𝐸₀. Moreover, if (𝜀₁,𝜀₂) is in the region between the curves 𝐻₂ and 𝐶^′, Sys. (6.20) has three limit periodic orbits: a “large” one and two “small” ones.(iv)When (𝜀₁,𝜀₂) is on the curve 𝐶_𝑑={(𝜀₁,𝜀₂)∶𝜀₂=𝑐𝜀₁+𝒪(|𝜀₁|^3/2),𝜀₁>0}, (6.25) where 𝑐≈0.752, Sys. (6.20) undergoes a double limit cycle bifurcation. Moreover, if (𝜀₁,𝜀₂) is in the region between the curves 𝐶^′ and 𝐶_𝑑, then Sys. (6.20) has two large limit cycles: the outer one which is stable and the inner one which is unstable, and these two cycles collide on 𝐶_𝑑.

Applying the above results and using the expressions of 𝜌₁, 𝜌₂, 𝜀₁, 𝜀₂, we obtain the following theorem regarding Sys. (5.3).

Theorem 6.3. Suppose that 𝑖⁽²⁾=0, 𝑖⁽³⁾<0 and 𝛿>0. For sufficiently small 𝜇₁, 𝜇₂, considered the following.(i)Sys. (5.3) undergoes a pitchfork bifurcation in the line 𝑇₂={(𝜇₁,𝜇₂)∶𝜇₁=−𝛽𝛾𝑞𝜏^∗,𝜇₂∈ℝ}. (6.26)(ii)Sys. (5.3) undergoes a branch of stable Hopf bifurcation on the curve 𝐻₂={(𝜇₁,𝜇₂)∶𝜇₂=−𝛾(3𝑚𝛼𝛽²+3𝛼𝛽𝑛₁+(𝑞+3𝛽)𝑛₂)2𝑞(𝑚𝛼𝛽+𝛼𝑛₁+𝑛₂)𝜏^∗𝜇₂,𝜇₂<0}. (6.27)(iii)Sys. (5.3) has two small homoclinic orbits simultaneously at (Y^∗,𝐾^∗) and a large homoclinic orbit on the curve 𝐶^′={(𝜇₁,𝜇₂)∶𝜇₂=−𝛾(12𝑚𝛼𝛽²+12𝛼𝛽𝑛₁+(5𝑞+12𝛽)𝑛₂)7𝑞(𝑚𝛼𝛽+𝛼𝑛₁+𝑛₂)𝜏^∗𝜇₂+𝒪(|𝜇₂|^3/2),𝜇₂<0}. (6.28) Moreover, if (𝜇₁,𝜇₂) is in the region between the curves 𝐻₂ and 𝐶^′, then Sys. (5.3) has three limit periodic orbits: a “large” one and two “small” ones.(iv)Sys. (5.3) undergoes a branch of a double limit cycle bifurcation on the curve 𝐶_𝑑={(𝜇₁,𝜇₂)∶𝜇₂=−𝛾(3𝑐𝑚𝛼𝛽²+3𝑐𝛼𝛽𝑛₁+(𝑞+3𝑐𝛽)𝑛₂)(−1+3𝑐)𝑞(𝑚𝛼𝛽+𝛼𝑛₁+𝑛₂)𝜏^∗𝜇₂+𝒪(|𝜇₂|^3/2),𝜇₂<0}, (6.29) where the constant 𝑐≈0.752. Moreover, if (𝜇₁,𝜇₂) is in the region between the curves 𝐶^′ and 𝐶_𝑑, Sys. (5.3) has two large different limit cycles. The outer one is stable, the inner one is unstable, and these two cycles collide on 𝐶_𝑑.

7. Numerical Simulations

In this section, we give some examples to verify the theoretical results obtained in Section 6. For simplicity, we assume that (0,0) is one of the equilibrium points.

Example 7.1. This example demonstrates the result of Theorem 6.1. Let 𝛼=1, 𝛽=0.8, 𝛾=0.5625, 𝑞=0.9. Then 𝑘^∗=1.0625, 𝜏^∗=2.66667. Take 𝐼(𝑠)=tanh(𝑘𝑠)+0.1𝑠². (7.1) Then (0,0) is the trivial equilibrium point and 𝑖⁽²⁾=0.1>0 and 𝛿=−0.800442<0. Simple calculation shows that 𝐻={(𝜇₁,𝜇₂)∶𝜇₁=−0.079918𝜇₂,𝜇₂>0},𝑃={(𝜇₁,𝜇₂)∶𝜇₁=−0.0368852𝜇₂,𝜇₂>0}. (7.2) Take 𝜇₁=−0.00058, 𝜇₂=0.01. Then 𝑘=𝑘^∗+𝜇₁=1.06192,𝜏=𝜏^∗+𝜇₂=2.67777, and (𝜇₁,𝜇₂) is in the region between 𝐻 and 𝑃. According to Theorem 6.1, Sys. (5.3) has a unique stable periodic orbit (see Figure 1).

Figure 1 

A stable limit cycle is generated when (𝜇1,𝜇2) is in the region between 𝐻 and 𝑃.

Example 7.2. This example supports the result of Theorem 6.2. Let 𝛼=1, 𝛽=0.8, 𝛾=0.5625, 𝑞=0.9. Then 𝑘^∗=1.0625, 𝜏^∗=2.66667, and 𝛿=−0.800442<0. Take 𝐼(𝑠)=𝑘𝑠+0.001𝑠³. (7.3) Then 𝑖⁽²⁾=0, 𝑖⁽³⁾=0.001<0, and hence, the condition of Theorem 6.2 is satisfied. After using the algorithm in Section 6, we have 𝐻₁={(𝜇₁,𝜇₂)∶𝜇₁=−0.295082𝜇₂,𝜇₂>0},𝐶={(𝜇₁,𝜇₂)∶𝜇₁=−0.456455𝜇₂+𝒪(𝜇₂^3/2),𝜇₂>0}. (7.4) Take 𝜇₁=−0.000035, 𝜇₂=0.0001, and hence (𝜇₁,𝜇₂) is in the region between 𝐻₁ and 𝐶. Figure 2 shows that there is a limit cycle which is stable according to Theorem 6.2.

Figure 2 

A stable periodic orbit is generated when (𝜇1,𝜇2) is located in the region between 𝐻 and 𝑃.

Example 7.3. This example verifies the result of Theorem 6.3. Let 𝛼=1, 𝛽=0.5, 𝛾=0.5625, 𝑞=0.8. Then 𝑘^∗=0.914063, 𝜏^∗=3.37222 and 𝛿=0.0233>0. Take 𝐼(𝑠)=𝑘𝑠−0.001𝑠³, (7.5) and hence, 𝐼^′′(0)=0, 𝐼^′′′(0)<0. Thus the condition of Theorem 6.3 holds. Simple calculation shows 𝐻₂={(𝜇₁,𝜇₂)∶𝜇₁=−0.724353𝜇₂,𝜇₂<0},𝐶^′={(𝜇₁,𝜇₂)∶𝜇₁=−0.99011𝜇₂+𝒪(|𝜇₂|^3/2),𝜇₂<0},𝐶_𝑑={(𝜇₁,𝜇₂)∶𝜇₁=−1.09167𝜇₂+𝒪(|𝜇₂|^3/2),𝜇₂<0}. (7.6) If we take 𝜇₁=0.0000861, 𝜇₂=−0.0001, then (𝜇₁,𝜇₂) is in the region between 𝐻₂ and 𝐶^′, and hence there are two small limit cycles and a large limit cycle (Figure 3). If we take 𝜇₁=0.0001, 𝜇₂=−0.0001, then (𝜇₁,𝜇₂) is in the region between 𝐶^′ and 𝐶_𝑑, and hence, there are two large limit cycles (Figure 4).

(a)

(b)

(c)

Figure 3 

Three limit (two small and one large) cycles are generated when (𝜇1,𝜇2) is in the region between 𝐻2 and 𝐶′.

Figure 4 

Two large limit cycles are generated when (𝜇1,𝜇2) is in the region between 𝐶′ and 𝐶𝑑.

Acknowledgment

The author would like to thank the referees for their valuable suggestions.