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

Authors: Xiaoqin P. Wu
Journal: Discrete Dynamics in Nature and Society
 (2009)

---

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


---

## Body

## 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)dY(t)dt=α[I(Y(t),K(t))-S(Y(t),K(t))],dK(t)dt=I(Y(t-τ),K(t))-qK(t),
where Y is the gross product, K is the capital stock, α>0 is the adjustment coefficient in the goods market, q∈(0,1) is the depreciation rate of capital stock, I(Y,K) and S(Y,K) 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 K in the investment function I(Y,K) of the second equation of Sys. (1.1) leads to the following Kaldor-Kalecki model of business cycles:(1.2)dY(t)dt=α[I(Y(t),K(t))-S(Y(t),K(t))],dK(t)dt=I(Y(t-τ),K(t-τ))-qK(t).As in [3]; also see [10, 16, 17], using the following saving and investment functions S and I, respectively,
(1.3)S(Y,K)=γY,I(Y,K)=I(Y)-βK,
where β>0 and γ∈(0,1) are constants, we obtain the following system:(1.4)dY(t)dt=α[I(Y(t))-βK(t)-γY(t)],dK(t)dt=I(Y(t-τ))-βK(t-τ)-qK(t).
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 k and τ as bifurcation parameters (where k 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 Section2, a detailed presentation is given for the distribution of eigenvalues of the linear part of Sys. (1.4) at an equilibrium point in the (k,τ)-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(2.1)α,β>0,q,γ∈(0,1),andI(s)isanonlinearC4function,
and that (Y*,K*) is an equilibrium point of Sys. (1.4). Let I*=I(Y*), u1=Y-Y*, u2=K-K*, and i(s)=I(s+Y*)-I*. Then Sys. (1.4) can be transformed as(2.2)du1(t)dt=α[i(u1(t))-βu2(t)-γu1(t)],du2(t)dt=i(u1(t-τ))-βu2(t-τ)-qu2(t).
Let the Taylor expansion of i at 0 be(2.3)i(u)=ku+i(2)u2+i(3)u3+𝒪(|u|4),
where(2.4)k=i′(0)=I′(Y*),i(2)=12i′′(0)=12I′′(Y*),i(3)=13!i′′′(0)=13!I′′′(Y*).
The linear part of Sys. (2.2) at (0,0) is(2.5)du1(t)dt=α[(k-γ)u1(t)-βu2(t)],du2(t)dt=ku1(t-τ)-βu2(t-τ)-qu2(t),
and the corresponding characteristic equation is(2.6)Δ(λ)≡λ2+Aλ+B+(βλ+C)e-λτ=0,
where(2.7)A=q-α(k-γ),B=-αq(k-γ),C=αβγ.
For τ=0, (2.6) becomes(2.8)λ2+(A+β)λ+B+C=0.
Define(2.9)k*=βγq+γ,k**=q+βα+γ.Theorem 2.1.
Letτ=0. If k<min{k*,k**}, then all roots of (2.8) have negative real parts, and hence (Y*,K*) is asymptotically stable. If k>min{k*,k**}, then (2.8) has a positive root and a negative root, and hence, (Y*,K*) is unstable.Now assumeτ>0. Clearly Δ(0)=0 if and only if k=k*. Next we always assume that k=k*. It is easy to attain(2.10)Δ′(λ)=2λ+q-αβγq+βe-λτ-(βλ+C)τe-λτ,Δ′′(λ)=2-2βτe-λτ+βτ2λe-λτ+Cτ2e-λτ.
Define τ*=(q2+qβ-αβγ)/αβγq. Then we have that,(2.11)Δ′(0)=αβγq(τ*-τ),Δ′′(0)|τ=τ*=q4-β2q2+α2β2γ2αβγq2.
Define(2.12)f(x)=x2+βx-αβγ,g(x)=x2-β2x+α2β2γ2.
Hence if f(q)≤0, τ*≤0, and hence Δ′(0)<0, and if f(q)>0, τ*>0, and hence Δ′(0)=0 if and only if τ=τ*. Also Δ′′(0)|τ=τ*≠0 if and only if g(q2)≠0. Thus we obtain the following result.Lemma 2.2.
Suppose thatk=k*. 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 g(q2)≠0.Letωi (ω>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(2.13)ω2+αβγ=αβγcos(ωτ)+βωsin(ωτ),q2-αβγqω=αβγsin(ωτ)-βωcos(ωτ).
Adding squares of two equations yields
(2.14)ω2+g(q2)q2=0.
Then (2.14) has a nonzero solution if and only if g(q2)<0 and does not have a nonzero solution if and only if g(q2)≥0. If g(q2)<0, from (2.14), we solve ω as follows:(2.15)ω=ω0≡1q-g(q2),
and from (2.13), we solve cos(ω0τ), sin(ω0τ) as:(2.16)cos(ω0τ)=-q2ω02+αβγω02+qαγ(αβγ+ω02)qβ(α2γ2+ω02)≡a,sin(ω0τ)=q2αγω0-α2βγ2ω0+qαβγω0+qω03qβ(α2γ2+ω02)≡b.
Define(2.17)δ={arccosa,ifb≥0,2π-arccosa,ifb<0.
From (2.16), we obtain(2.18)τ=τj≡1ω0(δ+2jπ),j=0,1,2,….Clearly ifβ>2αγ, then g(x)=0 has two positive roots, and if β≤2αγ, then g(x)≥0. Now, under k=k*, we impose the following conditions:(H1)
β≤2αγ, τ*≤0,(H2)
β≤2αγ, τ*>0, τ≠τ*,(H3)
β≤2αγ, τ*>0, τ=τ*,(H4)
β>2αγ, τ*>0, τ≠τ*,g(q2)≥0,(H5)
β>2αγ, τ*>0, τ≠τ*,g(q2)<0,(H6)
β>2αγ, τ*>0, τ=τ*,g(q2)≥0,(H7)
β>2αγ, τ*>0, τ=τ*,g(q2)<0.Based on Lemma2.2, we have the following result.Lemma 2.3.
Suppose thatk=k* and 0<q<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 ±ω0i in the imaginary axis if τ=τj, j=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 ±ω0i in the imaginary axis if τ*=τj for some j.Now we use the roots off(x)=0, g(x)=0 to give a more detailed discussion for the roots of (2.6). Define(2.19)q0=12(-β+β2+4αβγ),q1=12(β2-β4-4α2β2γ2),q2=12(β2+β4-4α2β2γ2).
Clearly q0 is the positive root of f(x)=0 and q1, q2 are two positive roots of g(x)=0 if β>2αγ. Note that f(x)≤0 if 0<x≤q0,andf(x)>0 if x>q0, g(x)≥0 if 0<x≤q1, or x≥q2, then g(x)<0 if q1<x<q2. Also note that as well as if β>2αγ, q02<q1. In fact it is based on the following calculation:(2.20)q1-q02=12(β2-β4-4α2β2γ2)-14(-β+β2+4αβγ)2=β2(β2+4αβγ-β2-4α2γ2-2αγ)=2αβγ(β-β2-4α2γ2)β2+4αβγ+β2-4α2γ2+2αγ>0.
Thus for β>2αγ, we always have q0<q1<q2. Noting that q∈(0,1), we have the following result.Lemma 2.4.
Letβ>2αγ. Then the following are given.(i)
Suppose thatq0≥1. Then for 0<q<1, then (2.6) has a simple zero root and does not have roots in the imaginary axis.(ii)
Suppose thatq0<1≤q1<q2. If 0<q≤q0, then (2.6) has a simple zero root and does not have roots in the imaginary axis. And if q0<q<1, (2.6) has a double zero root and does not have roots in the imaginary axis.(iii)
Suppose thatq0<q1<1<q2. If 0<q≤q0, then (2.6) has a simple zero root and does not have roots in the imaginary axis. If q0<q≤q1, then (2.6) has a double zero root and does not have roots in the imaginary axis. And if q1<q<1, then (2.6) has a double zero root and has a pair of purely imaginary roots.(iv)
Suppose thatq2≥1. Then if 0<q≤q0, then (2.6) has a simple zero root and does not have roots in the imaginary axis. If q0<q≤q1, then (2.6) has a double zero root and does not have roots in the imaginary axis. If q1<q<q2, then (2.6) has a double zero root and has a pair of purely imaginary roots when τ*=τj for some j. And if q2≤q<1, (2.6) has a double zero root and does not have a pair of purely imaginary roots.Defineλ(τ)=σ(τ)+iω(τ) to be the root of (2.6) such that σ(τj)=0 and ω(τj)=ω0. Then we have the following result.Lemma 2.5.
Suppose thatk=k* and g(q2)<0. Then σ′(τj)>0.Proof.
Differentiating (2.6) with respect to τ yields
(2.21)(dλdτ)-1=[2λ+q-α(k-γ)]eλτ+βλβ(λ+αγ)-τλ,
and a simple calculation gives
(2.22)Re(dλdτ)-1|τ=τj=α2β2γ2+q2(-β2+q2+2ω02)β2q2(α2γ2+ω02)=q2β2-α2β2γ2-q4β2q2(α2γ2+ω02),
which gives
(2.23)SignRe(dλdτ)-1|τ=τj=Sign(-g(q2))=1,
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)P(λ,e-λτ)=p(λ)+q(λ)e-λτ,
where p, q are real polynomials such that deg(q)<deg(p) and τ≥0. As τ varies, the sum of the order of zeros of P(λ,e-λτ) on the open right half-plane can change only if a zero appears on or crosses the imaginary axis.Lemma 2.7.
Letk=k* and τ>0. Then, the following are obtained.(i)
Ifq>q0, then all roots of (2.6) except 0 and purely imaginary roots have negative real parts,(ii)
If0<q≤q0, then (2.6) has at least one positive root.Proof.
Note that, forτ=0, if q>q0 or q2+qβ>αβγ, Δ(λ)=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<q≤q0 or q2+qβ≤αβγ. For τ>0, let
(2.25)f(λ)=Δ(λ)λ=λ+A+βe-λτ+B+Ce-λτλ.
Also noting that B+C=0 when k=k*, we have that
(2.26)limλ→0+f(λ)=A+β-Cτ=1q[(q2+qβ-αβγ)-αβγτ]<0,
and limλ→∞f(λ)=∞. 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)dxdt=L(μ)xt+G(xt,μ),
where x∈C([-τ,0],ℝn), μ∈ℝp. This equation is equivalent to(3.2)dxdt=L(μ)xt+G(xt,μ),dμdt=0,
which can be written as(3.3)dXdt=ℒXt+F(Xt),
where X=(x,μ)T, F(Xt)=(G(xt),0)T, and ℒ=diag(L,0). Define X∈C:=C([-τ,0],ℝn+p) with supreme norm and Xt∈C is defined by Xt(θ)=X(t+θ), -τ≤θ≤0; ℒ:C→L(ℝn+p) is a bounded linear operator; and F:C→C is a Ck(k≥2) function with F(0)=0, DF(0)=0. Consider the following linear system:(3.4)Ẋ(t)=ℒXt.
Since ℒ is a bounded linear operator, then ℒ can be represented by a Riemann-Stieltjes integral(3.5)ℒφ=∫-τ0dη(θ)φ(θ),∀φ∈C,
by the Riesz representation theorem, where η(θ) (θ∈[-τ,0]) is an (n+p)×(n+p) matrix function of bounded variation. Let 𝒜0 be the infinitesimal generator for the solution semigroup defined by Sys. (3.4) such that(3.6)𝒜0φ=φ̇,D(𝒜0)={φ∈C1([-τ,0],ℝn+p):φ̇(0)=∫-τ0dη(θ)φ(θ)}.
Define the bilinear form between C and C*=C([0,τ],ℝ(n+p)*) (where ℝ(n+p)* is the space of all row (n+p)-vectors) by(3.7)〈ψ,φ〉=ψ(0)φ(0)-∫-τ0∫0θψ(ξ-θ)dη(θ)φ(ξ)dξ,∀ψ∈C*,∀φ∈C.
The adjoint of 𝒜0 is defined by 𝒜0* as(3.8)𝒜0*ψ=-ψ̇,D(𝒜0*)={φ∈C1([0,τ],ℝ(n+p)*):ψ̇(0)=-∫-τ0ψ(-θ)dη(θ)}.
In our setting, (3.3) has p trivial components. Assume that the characteristic equation of (3.3) has eigenvalue zero with multiplicity 2p and all other eigenvalues have negative real parts. Then ℒ has a generalized eigenspace P which is invariant under the flow (3.4). Let P* be the space adjoint with P in C*. Then C can be decomposed as C=P⊕Q where Q={φ∈C:〈ψ,φ〉=0, ∀ψ∈P*}. Choose the bases Φ and Ψ for P and P*, respectively, such that(3.9)〈Ψ,Φ〉=I,Φ̇=ΦJ,Ψ̇=-JΨ,
where J is Jordan matrix associated with the eigenvalue 0.To consider Sys. (3.3), we need to enlarge the space C to the following BC:(3.10)BC={φ:[-τ,0]→ℝn+p:φiscontinuouson[-τ,0),∃limθ→0-φ(θ)∈ℝn+p}.
The elements of BC can be expressed as ψ=φ+X0α with φ∈C,α∈ℝn+p, and(3.11)X0(θ)={0,-τ≤θ<0,I,θ=0,
where I is the n×n identity matrix. Define the projection π:BC→P by(3.12)π(φ+X0α)=Φ[(Ψ,φ)+Ψ(0)α].
Then the enlarged phase space BC can be decomposed as BC=P⊕kerπ. Let X=Φx+y with x∈ℝ2p and y∈Q1={φ∈Q:φ̇∈C}. Then (3.3) can be decomposed as(3.13)ẋ=Jx+Ψ(0)F(Φx+y),ẏ=𝒜Q1y+(I-π)X0F(Ψx+y),
where 𝒜 is an extension of the infinitesimal generator 𝒜0 from C1 to BC, defined by(3.14)𝒜0φ=φ̇+X0[Lφ-φ̇(0)]={φ̇,-1≤θ<0,∫-τ0dη(t)φ(t),θ=0,
for φ∈C1 and its adjoint by 𝒜* is defined by(3.15)𝒜*ψ={-ψ̇,0<s≤θ,∫-τ0ψ(-θ)dη(θ),s=0,
for ψ∈C1*. Let F(v)=∑j≥2(1/j!)Fj(v). Then Sys. (3.13) becomes(3.16)ẋ=Jx+∑j≥21j!fj1(x,y),ẏ=𝒜Q1y+∑j≥21j!fj2(x,y),
where(3.17)fj1(x,y)=Ψ(0)Fj(Φx+y),fj2(x,y)=(I-π)X0Fj(Φx+y).
On the center manifold, (3.16) can be approximated as(3.18)ẋ=Jx+∑j≥21j!fj1(x,0).


## 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 q>q0. Therefore (H2) is equivalent to(4.1)k=k*,q>q0,τ>0,τ≠τ*.
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 k as a bifurcation parameter near k*.SetC:=C([-τ,0],ℝ3), C*:=C([0,τ],ℝ3*). Let μ=k-k*. Then Sys. (2.5) can be rewritten as(4.2)du1dt=α[βγqu1(0)-βu2(0)+μu1(0)+i(2)u12(0)+i(0)u13(t)βγq]+𝒪(|μ||u|2+|u|4),du2dt=k*u1(-τ)-qu2(t)+μu1(-τ)-βu2(-τ)+i(2)u12(-τ)+i(3)u13(-τ)+𝒪(|μ||u|2+|u|4),dμdt=0.
The linearization of Sys. (4.2) at (0,0,0) is(4.3)du1dt=αβγqu1(0)-αβu2(0),du2dt=k*u1(-τ)-qu2(0)-βu2(-τ),dμdt=0.
Let η(θ)=𝔸δ(θ)+𝔹δ(θ+τ) where(4.4)𝔸=(αβγq-αβ00-q0000),𝔹=(000k*-β0000).
Let X=(u1,u2,μ)T and(4.5)F(Xt)=(αμu1(0)+αi(2)u12(0)+αi(3)u13(0)+𝒪(|μ||u|2+|u|4)μu1(-τ)+i(2)u12(-τ)+i(3)u13(-τ)+𝒪(|μ||u|2+|u|4)0).
Define(4.6)Lφ=∫-τ0dη(θ)φ(θ),∀φ∈C.
Then Sys. (4.2) becomes(4.7)Ẋ(t)=LXt+F(Xt).
From (3.7), the bilinear form can be expressed as(4.8)〈ψ,φ〉=ψ(0)φ(0)+∫-τ0ψ(ξ+τ)𝔹φ(ξ)dξ.
It is not hard to see that the infinitesimal generator 𝒜:C1→BC is given by(4.9)𝒜φ=φ̇+X0[Lφ-φ̇(0)]={φ̇,-τ≤θ<0,𝔸φ(0)+𝔹φ(-τ),θ=0,
for φ∈C1 and its adjoint 𝒜* by(4.10)𝒜*ψ={-ψ̇,0<s≤θ,ψ(0)𝔸+ψ(τ)𝔹,s=0,
for ψ∈C1*.Next we obtain the bases for the center spaceP and its adjoint space P*, respectively. Let 𝒜φ=0 for φ∈C1, that is,(4.11)φ̇(θ)=0for-τ≤θ<0,𝔸φ(0)+𝔹φ(-τ)=0forθ=0.
then we know that φ is a constant vector (a1,a2,a3)T∈ℝ3∖{0} such that(4.12)(𝔸+𝔹)(a1,a2,a3)T=0.
Then we have two linearly independent solutions φ1=(q,γ,0)T, φ2=(0,0,1)T which are bases for the center space P. Let Φ=(φ1,φ2).Similarly, let𝒜*ψ=0 for ψ∈C1*, that is,(4.13)-ψ̇(s)=0for0<s≤τ,ψ(0)𝔸+ψ(τ)𝔹=0fors=0,
then we know that ψ is a constant vector (b1,b2,b3)∈ℝ3*∖{0} such that
(4.14)(b1,b2,b3)(𝔸+𝔹)=0.
From this we have two linearly independent solutions ψ1=(-(q+β),αβ,0) and ψ2=(0,0,1) which are bases for the center space P*. Let Ψ=(rψ1,ψ2)T with r being determined such that 〈ψ1,φ1〉=1. In fact(4.15)r=1qαβγ(τ-τ*).
Clearly r is well defined since τ-τ*≠0. It is not hard to check that Φ̇=ΦJ, Ψ̇=-JΨ and 〈Ψ,Φ〉=I, where J=(0000).Letu=Φx+y. Then Sys. (4.2) can be decomposed as(4.16)ẋ=ΨF(Φx+y),ẏ=AQ1y+(I-π)X0F(Φx+y).
Write x=(x1,μ). Note that(4.17)Ψ(0)F(Φx)=(-rαq2(μx1+qi(2)x12+q2i(3)x13)0)+h.o.t..
Here h.o.t. represents higher-order terms. Thus, for sufficiently small μ, on the center manifold, if i(2)≠0, then Sys. (4.2) becomes(4.18)ẋ1=-rαq2μx1-rαq3i(2)x12+h.o.t.,μ̇=0.
If i(2)=0 and i(3)≠0, then Sys. (4.2) can be transformed into the following form:(4.19)ẋ1=-rαq2μx1-rαq4i(3)x13+h.o.t.,μ̇=0.
Thus we have the following results.Theorem 4.1.
Letμ be small. Then consider what follows.(i)
Suppose thatμ=0. Then if i(2)≠0, the equilibrium (Y*,K*) is unstable, and if i(2)=0 and i(3)≠0, then the equilibrium (Y*,K*) is asymptotically stable for (τ-τ*)i(3)>0 and unstable if (τ-τ*)i(3)<0.(ii)
The equilibrium(Y*,K*) is asymptotically stable if (τ-τ*)μ>0 and unstable if (τ-τ*)μ<0.(iii)
At(Y*,K*,k*), Sys. (1.4) undergoes a transcritical bifurcation if i(2)≠0 and a pitchfork bifurcation if i(2)=0 and i(3)≠0.


## 5. Double-Zero Singularity

In this section, we assume that one of the conditions (H3) and (H6) holds andg(q2)>0, or equivalently, as(5.1)k=k*,τ=τ*,q>q0,g(q2)>0.
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 k=k* and τ=τ*. We treat (k,τ) as a bifurcation parameter near (k*,τ*).By scalingt→t/τ, Sys. (2.2) can be written as(5.2)du1(t)dt=ατ(k-γ)u1(t)-αβτu2(t)+ατi(2)u12(t)+τi(3)u13(t)+𝒪(|u1|4),du2(t)dt=τku1(t-1)-(q+β)τu2(t)+τi(2)u12(t-1)+τi(3)u13(t-1)+𝒪(|u1|4).
Let C:=C([-1,0],ℝ4),C*:=C([0,1],ℝ4*). Let μ1=k-k*, μ2=τ-τ*. Then on C we have(5.3)du1dt=α[βγqτ*u1(0)-βτ*u2(0)+τ*μ1u1(0)+βγqμ2u1(0)-βμ2u2(0)+τ*i(2)u12(0)+i(2)μ2u12(0)+i(3)τ*u13(0)+i(3)μ2u13(0)βγq]+𝒪(|μ|2|u|+|μ||u|4),du2dt=βγqτ*u1(t-1)-qτ*u2(t)-βτ*u2(-1)+τ*μ1u1(-1)+βγqμ2u1(-1)-qμ2u2(0)-βμ2u2(-1)+i(2)τ*u12(-1)+i(2)μ2u12(-1)+i(3)τ*u13(-1)+i(3)μ2u13(-1)+𝒪(|μ|2|u|+|μ||u|4),dμ1dt=0,dμ2dt=0.
The linearization of Sys. (5.3) at (0,0,0,0) is(5.4)du1(t)dt=αβγqτ*u1(0)-αβτ*u2(0),du2(t)dt=k*τ*u1(-1)-qτ*u2(0)-βτ*u2(-1),dμ1dt=0,dμ2dt=0.
Let(5.5)η(θ)=𝔸δ(θ)+𝔹δ(θ+1),
where(5.6)𝔸=(αβγqτ*-αβτ*000-qτ*0000000000),𝔹=(0000k*τ*-βτ*0000000000).
Define(5.7)ℒφ=∫-10dη(θ)φ(θ),∀φ∈C.
Let C1=C1([-1,0],ℝ4). Let X=(u1,u2,μ1,μ2)T and F(Xt)=(F1,F2,0,0)T where(5.8)F1=α[τ*μ1u1(0)+βγqμ2u1(0)-βμ2u2(0)+i(2)τ*u12(0)+i(3)τ*u13(0)]+𝒪(|μ|2|u|+|μ||u|4),F2=τ*μ1u1(-1)+βγqμ2u1(-1)-qμ2u2(0)-βμ2u2(-1)+i(2)τ*u12(0)+i(3)τ*u13(-1)+𝒪(|μ|2|u|+|μ||u|4).
Then Sys. (5.3) can be transformed into(5.9)Ẋ(t)=ℒXt+F(Xt).
Let C*=C([0,1],ℝ4*). From (3.7), the bilinear inner product between C and C* can be expressed by(5.10)〈ψ,φ〉=ψ(0)φ(0)+∫-10ψ(ξ+1)𝔹φ(ξ)dξ,
for φ∈C and ψ∈C*. As in Section 4, the infinitesimal generator 𝒜:C1→BC associated with ℒ is given by(5.11)𝒜φ=φ̇+X0[ℒφ-φ̇(0)]={φ̇,-1≤θ<0,𝔸φ(0)+𝔹φ(-1),θ=0,
for φ∈C1 and its adjoint by(5.12)𝒜*ψ={-ψ̇,0<s≤1,ψ(0)𝔸+ψ(1)𝔹,s=0,
for ψ∈C1*. 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 spaceP and its adjoint space P*, respectively. Let 𝒜φ=0 for φ∈C1. This means that(5.13)φ̇(θ)=0for-1≤θ<0,𝔸φ(0)+𝔹φ(-1)=0forθ=0.
From this we obtain that φ(θ)=φ0 is a constant vector in ℝ4 satisfying(5.14)(𝔸+𝔹)φ0=0.
This equation has three linearly independent solutions: a1=(q,γ,0,0)T, a3=(0,0,1,0)T, a4=(0,0,0,1)T. Let φ10 be one of those. Suppose that 𝒜a2=φ10 for a2∈C1, namely,(5.15)ȧ2(θ)=φ10for-1≤θ<0,𝔸a2(0)+𝔹a2(-1)=φ10forθ=0.
This implies that there is a constant vector φ20 in ℝ4 such that a2(θ)=φ10θ+φ20 and(5.16)ℒ(φ10θ+φ20)=φ10.
Since(5.17)ℒ(φ10θ+φ20)=ℒ(φ10θ)+ℒ(φ20)=-Bφ10+(𝔸+𝔹)φ20,
we have that(5.18)(𝔸+𝔹)φ20=(I+𝔹)φ10.
It is easy to see that (5.18) has no solution if φ10 is either a3 or a4. For φ10=a1, setting φ20=(0,l,0,0)T in (5.18), we obtain
(5.19)l=-q2γq2+qβ-αβγ,
and hence a2(θ)=(θq,l+γθ,0,0)T. Thus we obtain bases a1, a2, a3, a4 for the center space P. Let Φ=(a1,a2,a3,a4). Then we have that Φ̇=ΦJ where J=(0100000000000000).Similarly, let𝒜*ψ2=0 for ψ2∈C1*, that is,(5.20)-ψ̇2(s)=0for0<s≤1,ψ2(0)𝔸+ψ2(-1)𝔹=0fors=0,
which means that ψ2(s)=ψ20 is a constant vector ψ20∈ℝ4*∖{0} satisfying(5.21)ψ20(𝔸+𝔹)=0.
This equation has three linearly independent solutions: b2=(m(q+β),-mαβ,0,0), b3=(0,0,1,0), b4=(0,0,0,1). Asserting that 〈b2,a2〉=1 gives(5.22)m=2(q2+qβ-αβγ)q4-q2β2+α2β2γ2.
Let ψ20 be one of b2, b3, b4. Suppose 𝒜*b1=ψ20, that is,(5.23)-ḃ1(s)=ψ20for0<s≤1,b1(0)𝔸+b1(1)𝔹=ψ20fors=0,
which implies that there is ψ10∈ℝ4* such that b1(s)=-ψ20s+ψ10 satisfying(5.24)ℒ*(-ψ20s+ψ10)=ψ20.
Since(5.25)ℒ*(-ψ20s+ψ10)=-ℒ*(ψ20s)+ℒ*(ψ10)=-ψ20𝔹+ψ10(𝔸+𝔹),
we have(5.26)ψ10(𝔸+𝔹)=ψ20(I+𝔹).
It is not hard to check that (5.26) has no solution if ψ20=b3 or b4. Letting ψ20=b2, setting ψ10=(n1,n2,0,0) in (5.26) and using 〈b1,a2〉=0, we can get n1 and n2:(5.27)n1=2(q+β)(q6-3q5β-3q4β2+q3β3+3q2α2β2γ2-3qα2β3γ2+2α3β3γ3)3(q4-q2β2+α2β2γ2)2,n2=2αβ(-q2+2qβ+αβγ)(q2+qβ-αβγ)23(q4-q2β2+α2β2γ2)2.
Hence(5.28)b1(s)=ψ10-sψ20=(-m(q+β)s+n1,mαβs+n2,0,0).
Then b1, b2, b3, b4 are bases of the center space P*. Let Ψ=(b1,b2,b3,b4)T. Then 〈Ψ,Φ〉=I, Φ̇=ΦJ and Ψ̇=-JΨ.Letu=Φx+y, namely,(5.29)u1(θ)=qx1+qθx2+y1(θ),u2(θ)=γx1+(l+γθ)x2+y2(θ).
Then Sys. (5.9) can be decomposed as(5.30)ẋ=Jx+Ψ(0)F(Φx+y),ẏ=AQ1y+(I-π)X0F(Φx+y).
Write x=(x1,x2,μ1,μ2). Then, on the center manifold, Sys. (5.30) becomes(5.31)ẋ1=(a11μ1+a21μ2)x1+(a21μ1+a22μ2)x2+αn1τ*(i(2)q2x12+i(3)q3τ*x13)+n2τ*×[i(2)q2(x1-x2)2+i(3)q3τ*(x1-x2)3]+h.o.t.,ẋ2=(b11μ1+b21μ2)x1+(b12μ1+b22μ2)x2+mα(q+β)τ*(i(2)q2x12+i(3)q3τ*x13)+mαβτ*×[i(2)q2(x1-x2)2+i(3)q3τ*(x1-x2)3]+h.o.t.,μ̇1=0,μ̇2=0,
where(5.32)a11=qτ*(αn1+n2),a21=-α((q+β)γ-qk*)n1+(-(2q+β)γ+qk*)n2,a12=-qn2τ*,a22=-(lαβn1+((q+β)(l-γ)+qk*)n2),b11=mq2ατ*,b21=mqαβγ,b12=mqαβτ*,b22=mαβ(-(q+β)γ+qk*).Next we use techniques of nonlinear transformations in [22] to transform Sys. (5.31) into normal forms. If i(2)≠0, then up to the second order, Sys. (5.31) can be written as(5.33)ẋ1=(a11μ1+a21μ2)x1+(a21μ1+a22μ2)x2+αn1τ*i(2)q2x12+n2τ*i(2)q2(x1-x2)2+h.o.t.,ẋ2=(b11μ1+b21μ2)x1+(b12μ1+b22μ2)x2+mα(q+β)τ*i(2)q2x12+mαβτ*i(2)q2(x1-x2)2+h.o.t.,μ̇1=0,μ̇2=0.
This system can be transformed into the following normal form:(5.34)ẋ1=x2+h.o.t.,ẋ2=ρ1x1+ρ2x2+a1x12+b1x1x2+h.o.t.,
where(5.35)ρ1=b11μ1+b21μ2,ρ2=(a11+b12)μ1+(a21+b22)μ2,a1=mq3ατ*i(2),b1=2q2(mαβ+αn1+n2)τ*i(2).
Since(5.36)|∂ρ∂μ|=det(∂ρ1∂μ1∂ρ1∂μ2∂ρ2∂μ1∂ρ2∂μ2)=-mq2αγ(mαβ2+αβn1+(q+β)n2)τ*=-4q3α3β2γ2τ*(q2+qβ-αβγ)(q4-q2β2+α2β2γ2)2≠0,
we have that (μ1,μ2)→(ρ1,ρ2) is regular and hence the transversality condition holds.Ifi(2)=0 and i(3)≠0, then up to the third order, Sys. (5.31) becomes(5.37)ẋ1=(a11μ1+a21μ2)x1+(a21μ1+a22μ2)x2+αn1τ*i(3)q3τ*x13+n2τ*i(3)q3τ*(x1-x2)3+h.o.t.,ẋ2=(b11μ1+b21μ2)x1+(b12μ1+b22μ2)x2+mα(q+β)τ*i(3)q3τ*x13+mαβτ*i(3)q3τ*(x1-x2)3+h.o.t.,μ̇1=0,μ̇2=0.
This system can be transformed into the following normal form:(5.38)ẋ1=x2+h.o.t.,ẋ2=ρ1x1+ρ2x2+a2x13+b2x12x2+h.o.t.,
where a2=mq4ατ*i(3), b2=3q3(mαβ+αn1+n2)τ*i(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)ẋ1=x2,ẋ2=ρ1x1+ρ2x2+a1x12+b1x1x2,
where a2 and b2 are in Section 5. Note that (a1,b1)→(-a1,-b1) under the transformation (x1,x2)→(-x1,-x2). We may assume that i(2)>0. After the change of coordinates(6.2)x1=a1b12(ξ1-ρ1b1),x2→sa12b13ξ2,t→|b1a1|τ,
we have (still using x1,x2 for simplicity) that
(6.3)ẋ1=x2,ẋ2=ν1+ν2x1+x12+sx1x2,
where ν1=-(b12/a13)(b1ρ1-a1ρ2), ν2=(b1/a12)(b1ρ1-2a1ρ1) and s=signa1b1=±1. Simple calculation shows that s=sign(δ) where(6.4)δ=q6-3q4β2-2q3β3+3q2α2β2γ2+2α3β3γ3.
Now take s=-1, namely δ<0. The complete bifurcation diagrams of Sys. (6.3) can be found in [22]. Here, we just briefly list some results. For (ν1,ν2) small enough, consider the following.(i)
Sys. (6.3) undergoes a fold bifurcation when (ν1,ν2) is on the curves
(6.5)T+={(ν1,ν2):4ν1-ν22=0,ν2>0},T-={(ν1,ν2):4ν1-ν22=0,ν2<0}.(ii)
Sys. (6.3) undergoes a Hopf bifurcation when (ν1,ν2) is on the half-line
(6.6)H={(ν1,ν2):ν1=0,ν2<0},
and the Hopf bifurcation gives rise to a stable limit cycle.(iii)
Sys. (6.3) undergoes a homoclinic loop bifurcation when (ν1,ν2) is on the curve
(6.7)P={(ν1,ν2):ν1=-625ν22,ν2<0}.Moreover, when(ν1,ν2) is in the region between the curves H and P, Sys. (6.1) has a unique stable periodic orbit.Fors=1, under the transformation t→-t, x1→-x1, 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 ν1,ν2, we obtain the following result regarding Sys. (5.3).Theorem 6.1.
Suppose thati(2)>0 and δ<0. For sufficiently small μ1,μ2, consider the following(i)
Sys. (5.3) undergoes a fold bifurcation in the half-lines
(6.8)T¯+={(μ1,μ2):μ1=-βγqτ*μ2,μ2<0},T¯-={(μ1,μ2):μ1=-βγqτ*μ2,μ2>0}.(ii)
Sys. (5.3) undergoes a Hopf bifurcation on the curve
(6.9)H¯={(μ1,μ2):μ1=-γ(2mαβ2+2αβn1+(q+2β)n2)q(mαβ+αn1+n2)τ*μ2+𝒪(μ22),μ2>0}.(iii)
Sys. (5.3) undergoes a saddle of homoclinic bifurcation on the curve
(6.10)P¯={(μ1,μ2):μ1=-γ(12mαβ2+2αβn1+(7q+12β)n2)5q(mαβ+αn1+n2)τ*μ2+𝒪(μ23/2),μ2>0}.Moreover, if (μ1,μ2) is in the region between the curves H¯ and P¯, Sys. (5.3) has a unique stable periodic orbit.Next, we consider the truncated system of (5.38):
(6.11)ẋ1=x2,ẋ2=ρ1x1+ρ2x2+a2x13+b2x12x2,
where a2,b2 are in Section 5. The bifurcation diagrams of this system are more complicated and interesting. We must consider two cases.Case 1 (i(3)>0 so that a2>0).
We can assumeb2<0 which is equivalent to δ<0 in (6.4). Then Sys. (6.11) can be transformed as
(6.12)ẋ1=x2,ẋ2=ε1x1+ε2x2+x13-x12x2,
where
(6.13)ε1=(b2a2)2ρ1,ε2=-b2a2ρ2.
The complete bifurcation diagrams of Sys. (6.12) can be found, for example, in [22–24]. For b2>0, under the transformation t→-t and x1→-x1, we can get Sys. (6.12). Here, we briefly list some results: for small ε1, ε2 as follows. (i)
When(ε1,ε2) is in the line
(6.14)T1={(ε1,ε2):ε1=0,ε2∈ℝ},
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
(6.15)H1={(ε1,ε2):ε2=0,ε1<0};(iii)
On the curve(6.16)C={(ε1,ε2):ε2=-15ε1+𝒪(|ε1|3/2),ε1<0},Sys. (6.12) undergoes a heteroclinic bifurcation. Moreover, if (ε1,ε2) is in the region between the curves H1 and C then Sys. (6.12) has a unique stable periodic orbit.Applying the above results and using the expressions ofρ1, ρ2, ε1, ε2, we obtain the following theorem regarding Sys. (5.3).Theorem 6.2.
Suppose thati(2)=0, i(3)>0 and δ<0. For sufficiently small μ1, μ2, the following are given.(i)
Sys. (5.3) undergoes a pitchfork bifurcation in the line
(6.17)T¯1={(μ1,μ2):μ1=-βγqτ*μ2,μ2∈ℝ}.(ii)
Sys. (5.3) undergoes a stable Hopf bifurcation in the half-line
(6.18)H¯1={(μ1,μ2):μ1=γn2τ*(mαβ+αn1+n2)μ2,μ2>0}.(iii)
Sys. (5.3) undergoes a branch of homoclinic bifurcation on the curve
(6.19)C¯={(μ1,μ2):μ1=γ(3mαβ2+3αβn1+(5q+3β)n2)2q(mαβ+αn1+n2)τ*μ2+𝒪(|μ2|3/2),μ2>0}.Moreover, if (μ1,μ2) is in the region between the curves H¯1 and C¯, Sys. (5.3) has a unique stable periodic orbit.Case 2 (i(3)<0 so that a2<0).
We assumeb2<0, namely δ>0. After rescaling, then Sys. (6.11) can be transformed as
(6.20)ẋ1=x2,ẋ2=ε1x1+ε2x2-x13-x12x2,
where
(6.21)ε1=(b1a1)2ρ1,ε2=b1a1ρ2.
For b2>0, under the transformation t→-t,x1→-x1, 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 ε1, ε2, as follows.(i)
When(ε1,ε2) is in the line
(6.22)T2={(ε1,ε2):ε1=0,ε2∈ℝ},
Sys. (6.20) undergoes a pitchfork bifurcation.(ii)
When(ε1,ε2) is in the half-line
(6.23)H2={(ε1,ε2):ε2=ε1,ε1>0},
Sys. (6.20) undergoes a stable Hopf bifurcation at E1,2 and the bifurcation is subcritical.(iii)
When(ε1,ε2) is on the curve
(6.24)C′={(ε1,ε2):ε2=45ε1+𝒪(|ε1|3/2),ε1>0},
Sys. (6.20) has a unique homoclinic orbit connecting E1 and E2 and two homoclinic orbits simultaneously at E0. Moreover, if (ε1,ε2) is in the region between the curves H2 and C′, Sys. (6.20) has three limit periodic orbits: a “large” one and two “small” ones.(iv)
When(ε1,ε2) is on the curve
(6.25)Cd={(ε1,ε2):ε2=cε1+𝒪(|ε1|3/2),ε1>0},
where c≈0.752, Sys. (6.20) undergoes a double limit cycle bifurcation. Moreover, if (ε1,ε2) is in the region between the curves C′ and Cd, 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 Cd.Applying the above results and using the expressions ofρ1, ρ2, ε1, ε2, we obtain the following theorem regarding Sys. (5.3).Theorem 6.3.
Suppose thati(2)=0, i(3)<0 and δ>0. For sufficiently small μ1, μ2, considered the following.(i)
Sys. (5.3) undergoes a pitchfork bifurcation in the line
(6.26)T¯2={(μ1,μ2):μ1=-βγqτ*,μ2∈ℝ}.(ii)
Sys. (5.3) undergoes a branch of stable Hopf bifurcation on the curve
(6.27)H¯2={(μ1,μ2):μ2=-γ(3mαβ2+3αβn1+(q+3β)n2)2q(mαβ+αn1+n2)τ*μ2,μ2<0}.(iii)
Sys. (5.3) has two small homoclinic orbits simultaneously at (Y*,K*) and a large homoclinic orbit on the curve
(6.28)C¯′={(μ1,μ2):μ2=-γ(12mαβ2+12αβn1+(5q+12β)n2)7q(mαβ+αn1+n2)τ*μ2+𝒪(|μ2|3/2),μ2<0}.
Moreover, if (μ1,μ2) is in the region between the curves H¯2 and C¯′, 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
(6.29)C¯d={(μ1,μ2):μ2=-γ(3cmαβ2+3cαβn1+(q+3cβ)n2)(-1+3c)q(mαβ+αn1+n2)τ*μ2+𝒪(|μ2|3/2),μ2<0},
where the constant c≈0.752. Moreover, if (μ1,μ2) is in the region between the curves C¯′ and C¯d, 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 C¯d.


## 7. Numerical Simulations

In this section, we give some examples to verify the theoretical results obtained in Section6. For simplicity, we assume that (0,0) is one of the equilibrium points.Example 7.1.
This example demonstrates the result of Theorem6.1. Let α=1, β=0.8, γ=0.5625, q=0.9. Then k*=1.0625, τ*=2.66667. Take
(7.1)I(s)=tanh(ks)+0.1s2.
Then (0,0) is the trivial equilibrium point and i(2)=0.1>0 and δ=-0.800442<0. Simple calculation shows that
(7.2)H¯={(μ1,μ2):μ1=-0.079918μ2,μ2>0},P¯={(μ1,μ2):μ1=-0.0368852μ2,μ2>0}.
Take μ1=-0.00058, μ2=0.01. Then k=k*+μ1=1.06192,τ=τ*+μ2=2.67777, and (μ1,μ2) is in the region between H¯ and P¯. 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 H¯ and P¯.Example 7.2.
This example supports the result of Theorem6.2. Let α=1, β=0.8, γ=0.5625, q=0.9. Then k*=1.0625, τ*=2.66667, and δ=-0.800442<0. Take
(7.3)I(s)=ks+0.001s3.
Then i(2)=0, i(3)=0.001<0, and hence, the condition of Theorem 6.2 is satisfied. After using the algorithm in Section 6, we have
(7.4)H¯1={(μ1,μ2):μ1=-0.295082μ2,μ2>0},C¯={(μ1,μ2):μ1=-0.456455μ2+𝒪(μ23/2),μ2>0}.
Take μ1=-0.000035, μ2=0.0001, and hence (μ1,μ2) is in the region between H¯1 and C¯. 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 H¯ and P¯.Example 7.3.
This example verifies the result of Theorem6.3. Let α=1, β=0.5, γ=0.5625, q=0.8. Then k*=0.914063, τ*=3.37222 and δ=0.0233>0. Take
(7.5)I(s)=ks-0.001s3,
and hence, I′′(0)=0, I′′′(0)<0. Thus the condition of Theorem 6.3 holds. Simple calculation shows
(7.6)H¯2={(μ1,μ2):μ1=-0.724353μ2,μ2<0},C¯′={(μ1,μ2):μ1=-0.99011μ2+𝒪(|μ2|3/2),μ2<0},C¯d={(μ1,μ2):μ1=-1.09167μ2+𝒪(|μ2|3/2),μ2<0}.
If we take μ1=0.0000861, μ2=-0.0001, then (μ1,μ2) is in the region between H¯2 and C¯′, and hence there are two small limit cycles and a large limit cycle (Figure 3). If we take μ1=0.0001, μ2=-0.0001, then (μ1,μ2) is in the region between C¯′ and C¯d, and hence, there are two large limit cycles (Figure 4).Three limit (two small and one large) cycles are generated when(μ1,μ2) is in the region between H¯2 and C¯′.
(a)(b)(c)Figure 4
Two large limit cycles are generated when(μ1,μ2) is in the region between C¯′ and C¯d.



---

*Source: 923809-2010-02-03.xml*