DYNAMICS AND DOUBLE HOPF BIFURCATIONS OF THE ROSE ...
December 9, 2009
11:8
02508
International Journal of Bifurcation and Chaos, Vol. 19, No. 11 (2009) 3733–3751 c World Scientific Publishing Company
DYNAMICS AND DOUBLE HOPF BIFURCATIONS OF THE ROSE–HINDMARSH MODEL WITH TIME DELAY* SUQI MA Department of Mathematics, Chinese Agricultural University, Beijing 100083, P. R. China ZHAOSHENG FENG† Department of Mathematics, University of Texas-Pan American, Edinburg, TX 78539, USA [email protected] QISHAI LU School of Science, Beijing University of Aeronautics and Astronautics, Beijing 100083, P. R. China Received January 11, 2008; Revised April 22, 2009 In this paper, we are concerned with the Rose–Hindmarsh model with time delay. By applying the generalized Sturm criterion, a number of imaginary roots of the characteristic equation are classified. The absolutely stable regions for any value of time delay are detected. By the continuous software DDE-Biftool, both the Hopf bifurcation curves and double Hopf bifurcation points are illustrated in parametric spaces. The normal form and universal unfolding at double Hopf bifurcation points are considered by the center manifold method. Some examples also indicate that the corresponding unique attractor near each double Hopf point is asymptotically stable. Keywords: Hopf bifurcation; Rose–Hindmarsh model; normal form; center manifold method; autonomous system; elliptic function.
1. Introduction
and the Rose–Hindmarsh model is described as
The electric activity of neurons is of great concern. A typical example is given by the Rose–Hindmarsh model of action potential which has continuously attracted considerable attention in the past decades [Hindmarsh et al., 1982, 1984]. In [Hindmarsh et al., 1984], in order to simulate the bursting behavior of neurons, an adaption current z was introduced to the two-variable model of the excitable membrane,
x = y − ax3 + bx(t)2 − cz + Iapp , y = c − dx2 − y, z = r(S(x − χ) − z), where the variable x represents membrane potential, y represents a recovery variable, z represents an adaption current, and Iapp represents an applied current. Here r is the parameter of control and
∗ Main contents have been presented at the Seventh AIMS International Conference on Dynamical Systems, Differential Equations and Applications, University of Texas, Arlington, 18–21 May 2008. This work is supported by UTPA Faculty Research Council Grant 119100 and NSF of China under 10771212. † Author for correspondence
3733
December 9, 2009
3734
11:8
02508
S. Ma et al.
all other parameters a, b, c, d, S and χ are real constants. Dynamical behaviors of the Rose–Hindmarsh model have been intensively studied. In [Holden et al., 1992a, 1992b] both periodic and aperiodic “bursting” oscillations were considered which describe the behavior of a period of rapid spiking followed by a quiescent period. Numerical simulations of evolution from simple periodic bursting to complex periodic bursting were considered in [Holden et al., 1992a]. Transition from periodicity to aperiodicity via an alternation between the perioddoubling bifurcation sequence and the saddle-node bifurcation was presented in [Holden et al., 1992b]. It was shown [Pomeau et al., 1980] and [Grebogi et al., 1986, 1987] that transition from intermittent chaos to complex oscillatory behavior was caused by a crisis of the chaotic attractor, and after the crisis periodicity was maintained by period-adding bifurcations. Mathematical models with time delays have been broadly applied in many scientific fields such as neuronscience [Shi et al., 2007], biology [Zhao et al., 2001] and [Xu et al., 1999], and mechanics [Ma et al., 2007]. We see that in many cases controlled systems with delay feedback may have double Hopf bifurcations on the margin of “amplitude death origin” [Ma et al., 2007, 2008]. Such kind of phenomena has also been observed in [Campell et al., 1995b]. Normal form and universal unfolding in a vector field at nonresonant double Hopf bifurcation points for particular retarded functional differential equations (RFDEs) were presented in [Buono, 2003; Campell, 1995a; Elphick, 1987; Guckenheimer, 1983] with some restrictions on the flows on the center manifold for certain singularities, based on the well-known center manifold theory for RFDEs (see [Ma et al., 2007, 2008] and references therein). Recently, it was found that delay factors are inherent in many biological systems due to finite propagation speed of signals and finite processing time in synapses. The Rose–Hindmarsh model with time delay takes the form x = y − ax3 + bx(t − τ )2 − cz + Iapp , y = c − dx2 − y, z = r(S(x − χ) − z).
few references available so far. So, in this paper we will focus our attention on the dynamical behavior of system (1) by changing values of time delay τ , parameters S and Iapp , and fixing values for other parameters as a = 1.0,
When the time delay τ is equal to zero, Rose and Hindmarsh observed chaos in system (1) for parameter Iapp ∈ [3.28, 3.405]. When the time delay τ is not equal to zero, it seems that there are very
c = 1.0, d = 5.0, r = 0.001.
The rest of this paper is organized as follows. In Sec. 2, an analytic method is used to derive the condition under which Hopf bifurcation occurs for system (1). We apply the formula of the Sturm sequence (see [Yang et al., 1996a, 1996b; Wang et al., 1999, 2000]) to consider the stability when imaginary roots iω of the corresponding polynomial equation has high multiplicity. Hopf bifurcation curves are plotted in the (Iapp − τ ) plane by the numerical software DDE-Biftool. Double Hopf bifurcation points are obtained at the intersection points of Hopf bifurcation curves. In Sec. 3, the normal form and universal unfolding at double Hopf bifurcation points are studied using the theory of center manifold theorem. Numerical simulations of some examples are illustrated to indicate that numerical simulations agree well with theoretical analysis. Section 4 is a brief conclusion.
2. Stability and Hopf Bifurcation Analysis Suppose that (x∗ , y ∗ , z ∗ ) is the equilibrium solution of system (1), so it satisfies z ∗ = S(x∗ − χ), y ∗ = c − dx∗2 , c − dx∗2 − ax∗3 + bx∗2 − cS(x∗ − χ) + Iapp = 0.
(2)
To solve system (2), we set x∗ = ((3va + b − d)/3a) and have v 3 + 3pv + 2q = 0,
(3)
with p=− q=
(1)
b = 3.0, χ = −1.6,
cs (b − d)2 + , 2 9a 3a
(9acS(b − d) − 2(b − d)3 ) c + Iapp + cSχ . − 54a3 2a
After analyzing Eq. (3), we have four cases: (a) If p3 + q 2 > 0, then Eq. (3) has one real solution ( p3 + q 2 − q)2/3 − p ; v= ( p3 + q 2 − q)1/3
December 9, 2009
11:8
02508
Dynamics and Double Hopf Bifurcations of the Rose–Hindmarsh Model
(b) If p3 + q 2 = 0 and p = q = 0, then Eq. (3) has one real solution v = 0; (c) If p3 + q 2 = 0 and p, q = 0, then Eq. (3) has two real solutions v1 = −2q 1/3
and
v2 = q 1/3 ;
(d) If p3 + q 2 < 0, then Eq. (3) has three real solutions φ v1 = −2r cos , 3 φ+π , v2 = 2r cos 3 π−φ , v3 = 2r cos 3 3 respectively, where ∗cos φ = q/r and r = sign(q) · |p|. So x can be obtained explicitly and the associated equilibrium solution (x∗ , y ∗ , z ∗ ) can also be found accordingly.
3735
To analyze the stability of equilibrium solutions, we set x = x − x∗ , y = y − y ∗ , z = z − z ∗ , then system (1) becomes x = −3ax∗2x + 2bx∗x(t − τ ) + y − cz −ax 3 − 3ax∗x 2 + bx 2 (t − τ ), y = −2dx∗x − y − dx 2 , z = rSx − rz .
(4)
The characteristic equation of system (4) is −3ax∗2 − λ + 2bx∗ e−λτ 1 −c −1 − λ 0 = 0. −2dx∗ rS 0 −r − λ (5) Setting λ = iω in the characteristic equation (5) and separating the real part from the imaginary part, we have
sin(ωτ ) · 2bx∗ ω(1 + r) + cos(ωτ ) · 2bx∗ (r − ω 2 ) = 2dx∗ r + rSc + 3ax∗2 r − 3ax∗2 ω 2 − ω 2 − ω 2 r, (6) sin(ωτ ) · (−2b)x∗ (r − ω 2 ) + cos(ωτ ) · 2bx∗ ω(1 + r) = 2dx∗ ω + 3ax∗2 ω + ωr − ω 3 + rScω + 3ax∗2 ωr. System (6) is equivalent to sin(ωτ ) =
ω(2dx∗ ω 2 + rSc − ω 2 + 2dx∗ r 2 − r 2 − ω 2 r 2 + rScω 2 − ω 4 ) , 2bx∗ (ω 2 + r 2 + ω 4 + ω 2 r 2 )
cos(ωτ ) =
2dx∗ r 2 + r 2 Sc + r 2 Scω 2 + 3ax∗2 r 2 ω 2 + 3r 2 ax∗2 + 2dx∗ ω 2 + 3ax∗2 ω 4 + 3ax∗2 ω 2 , 2bx∗ (ω 2 + r 2 + ω 4 + ω 2 r 2 )
with G(ω 2 ) := ω 6 − a1 ω 4 − a2 ω 2 − a3 = 0,
(7)
where a1 = (2rSc + 4dx∗ − 9a2 x∗4 − r 2 + 4b2 x∗2 − 1), a2 = (−9a2 x∗4 + 4b2 x∗2 − 4d2 x∗2 + 4b2 x∗2 r 2 + 2rSc − 6r 2 Scax∗2 − 12dx∗3 a − r 2 S 2 c2 + 4dx∗ r 2 − 9a2 x∗4 r 2 − 4dx∗ rSc − r 2 ), a3 = −r 2 (Sc + 2dx∗ + 3ax∗2 − 2x∗ b) × (Sc + 2dx∗ + 3ax∗2 + 2x∗ b). The key step of the above discussion in testing the stability of system (1) is to determine the number of real roots of the polynomial as given on the left side of Eq. (7). Given a real polynomial with constant coefficients, the classical Sturm criterion gives a full answer to this problem. A complete discrimination system for polynomials
recently developed in [Yang et al., 1996a, 1996b], which is called the generalized Sturm criterion, offers a powerful tool to determine the number of a polynomial with unknown parameters. Wang and Hu [1999] presented some basic properties and provide an interpretation on the generalized Sturm criterion. In the following, we just introduce the generalized Sturm criterion, but omit the proof and basic concepts of discrimination sequence and the modified sign table (see [Yang et al. 1996a, 1996b; Wang et al., 1999, 2000]) for details). Proposition 1. Suppose that f (x) is a polynomial of order n. Let D1 (f ), D2 (f ), . . . , Dn (f ) be the corresponding discrimination sequence, and σi = Dqi (i = 1, 2, . . . , k) be the ith nonzero term of the discrimination sequence, where σ0 = 1. Let q0 ≡ 0 and ri = qi+1 − qi − 1 (i = 0, 1, 2, . . . , k − 1). Assume that the number of variation of signs in the modified sign table is s. If Dl (f ) = 0 and Dm (f ) = 0, m > l,
December 9, 2009
11:8
02508
S. Ma et al.
3736
then the following properties hold: (a) The number of distinct pairs of conjugate complex roots of f (x) is s. (b) The number of distinct real roots of f (x) is l − 2s, which satisfies k−1
l − 2s =
ri /2
(−1)
sgn
i=0, ri is even
σi+1 σi
(d) Except for some positive factors, D1 (f ), D2 (f ), . . . , Dl (f ) is the discrimination sequence of polynomial f /g.c.d(f, f ), which has no repeated roots. According to Proposition 1, we find that the discrimination sequence of G(ω 2 ) is: S := [1, a1 , a1 d1 , −d1 d2 , d2 d3 , d23 a3 ]
.
(c) α is a root of f (x) with multiplicity p if and only if it is a root of ∆n−l (f ) = 0 with multiplicity p − 1.
with d1 = (3a2 + a21 ), d2 = (a2 a21 − 3a3 a1 + 4a22 ), d3 = (4a31 a3 − a22 a21 + 18a1 a3 a2 − 4a32 + 27a23 ).
6
6
s
5
s I
4
4
IV
I (b)
I (a)
I (b)
3
3 III
III
II 2
2
1
1 II 0 −2
5
II
IV 0
−1.5
−1
−0.5
0
0.5 *
1
x
−1.6
−1.4
−1.2
(a)
−1
−0.8
−0.6
x*
(b)
1.5
6 5.5
s
s 5
I
4.5
1
4 3.5
II (b)
II (b)
II (a)
3
0.5
2.5 2
II (a)
1.5
0
−1.3
−1.2
−1.1
−1
−0.9
−0.8
−0.7
x (c)
*
−0.6
1
0.05
0.1
0.15
0.2
0.25
0.3
0.35
*
0.4
x (d)
Fig. 1. Classification diagram of the number of real roots of G(ω 2 ) in the (x∗ − S) plane. G(ω 2 ) is defined by Eq. (7). (a) The whole classification diagram, (b)–(d) the amplification cases.
December 9, 2009
11:8
02508
Dynamics and Double Hopf Bifurcations of the Rose–Hindmarsh Model Table 1.
Imaginary roots of the characteristic equation (5).
S
x∗
Region
Imaginary Roots ±iω
Modified Sign Table
3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 0.5 0.5 0.5
−1.8 −1.4 −1.2 −0.8 −0.5 0.0 0.1 0.2 −1.1 −0.8 −0.6
III I(b) I(a) I(b) III IV II(b) II(a) II(a) II(b) IV
±1.655996i ±2.75883i ±0.002059i, ±0.990345i, ±2.683806i ±0.000988i ±0.0009077i none none ±0.9251473i, ±1.89095i ±1.41178i, ±2.321846i none none
[1, −1, −1, −1, −1, 1] [1, 1, 1, −1, −1, 1] [1, 1, 1, 1, 1, 1] [1, 1, −1, 1, 1, 1] [1, −1, 1, 1, 1, 1] [1, −1, −1, −1, 1, −1] [1, 1, −1, 1, 1, −1] [1, 1, 1, 1, 1, −1] [1, 1, 1, 1, 1, −1] [1, 1, −1, 1, 1, −1] [1, −1, 1, 1, 1, −1]
We classify the number of real roots of G(ω 2 ) by the following four cases: Case I . Suppose a1 > 0, a3 > 0, then (a) if d1 > 0, d2 < 0, d3 < 0, G(ω 2 ) has three real roots; (b) otherwise, G(ω 2 ) has one real root. Case II . Suppose a1 > 0, a3 < 0, then (a) if d1 > 0, d2 < 0 or d1 > 0, d2 > 0, d3 < 0 or d1 < 0, d2 < 0, d3 > 0, G(ω 2 ) has two real roots; (b) if d2 > 0, d3 > 0 or d1 < 0, d3 < 0, G(ω 2 ) has no real roots. Case III . If a1 < 0, a3 > 0, G(ω 2 ) has one real root. Case IV . If a1 < 0, a3 < 0, G(ω 2 ) has two real roots provided that d1 > 0, d2 > 0, d3 < 0; otherwise, it has no real roots.
Set a = 1.0, b = 3.0, c = 1.0 and d = 5.0. The classification diagram of the number of real roots of G(ω 2 ) is plotted in the (x∗ − S) plane as shown in Fig. 1(a) and amplifying diagrams as shown in Figs. 1(b)–1(d). For example, transit (x∗ , S) between different regimes from left to right, imaginary roots of the characteristic equation (5) are listed in Table 1. The linearized system of system (4) is stable for any value of τ if there exists no imaginary root. Therefore, from Table 1 we obtain Theorem 1. The linearized system of (4) is abso-
lutely stable for any value of time delay τ if and only if parameter S and the value of x∗ lie in regions II(b) and IV . By applying the DDE-Biftool, we plot Hopf bifurcation curves in the (x∗ − τ ) plane while taking S = 1. As shown in Fig. 2, two pairs of imaginary
15
τ
14
3737
6 DH
τ
13 12
5 DH
4
11
DH
3
10 9
2
DH
8
6 0.1
1
DH
7 0.15
Fig. 2.
0.2
x*
0.25
0 1.7
1.75
1.8
Hopf bifurcation curves of the linearized system of (4) when S = 1.0.
1.85
x*
1.9
December 9, 2009
3738
11:8
02508
S. Ma et al. Table 2.
Values of double Hopf bifurcation points.
x∗
τ
Imaginary Roots
x∗
τ
Imaginary Roots
0.1768 1.7523
6.8693 1.7741
ω1 = 0.9107; ω2 = 1.6810 ω1 = 3.5918; ω2 = 6.8435
1.8207 1.8383
2.9934 4.2114
ω1 = 4.1953; ω2 = 6.1936 ω1 = 4.4642; ω2 = 5.9052
roots are acquired simultaneously at the intersection point of two Hopf bifurcation curves. Table 2 lists the values of imaginary roots and the related Double Hopf bifurcation points.
+ 2bˆ x∗ [x(t − τˆ − εˆ τε ) − x(t − τˆ)] x∗ x + 2εˆ xε bx(t − τˆ) − aεx3 + o(ε), − 6εˆ xε aˆ √ x∗ dx − y − 2εˆ xε dx − εdx2 , y = −2ˆ z = rSx − rz. (8)
3. Normal Form for Double Hopf Bifurcations
√ Set x∗ = √ x ˆ∗ + εˆ xε , τ = τˆ + εˆ τε and x = εx, y = √ εy, z = εz, then system (4) is equivalent to x∗2 x + 2bˆ x∗ x(t − τˆ) + y − cz x = −3aˆ √ √ x∗ x2 + εbx(t − τˆ)2 − 3 εaˆ
Choose the phase space C = C([−ˆ τ , 0), R3 ) as the Banach space of continuous functions from [−ˆ τ , 0) to R3 with the supremum norm. For φ ∈ C, we define 0 [dη(θ)]φ(θ), (9) L(0)φ = −ˆ τ
where η : [−ˆ τ , 0] → R3 ×R3 is a real-valued function of bounded variation in the interval [−ˆ τ , 0] with x∗ δ(θ + τˆ) δ(θ) −cδ(θ) −3aˆ x∗2 δ(θ) + 2bˆ −δ(θ) 0 dθ. −2ˆ x∗ dδ(θ) dη(θ) = rSδ(θ) 0 −rδ(θ)
Then, we let L(ˆ xε , τˆε )φ =
0
−ˆ τ
and
where
−6aˆ x∗ x ˆε δ(θ) + 2bˆ xε δ(θ + τˆ) −2ˆ xε dδ(θ) dη1 (θ) = 0
and
[εdη1 (θ, x ˆε ) + dη2 (θ, ε, τˆε )]φ(θ), 0 0 0 0 dθ, 0 0
2bˆ x∗ (δ(θ + τˆ + εˆ τε ) − δ(θ + τˆ)) 0 0 0 0 0 dθ. dη2 (θ) = 0 0 0 For φ ∈ C, the linear operator defined by system (9) generates an infinitesimal generator of the semi-flow of bounded linear operators with dφ , θ ∈ [−ˆ τ , 0), (10) A(0)φ = dθ L(0)φ, θ = 0. Hence, we define
0, A(ˆ xε , τˆε )φ = L(ˆ xε , τˆε )φ,
θ ∈ [−τ, 0), θ = 0,
Q(φ) = where
0, F (φ),
θ ∈ [−τ, 0), θ = 0,
√ √ x∗ φ21 (0) + εbφ21 (−ˆ τ ) − εaφ31 (0) −3 εaˆ √ − εdφ21 (0) F (φ) = . 0
System (8) can be re-expressed as an operator differential equation xε , τˆε )ut + εQ(φ), u (t) = A(0)ut + A(ˆ
(11)
τ ≤ where u = (x, y, z)T and ut = u(t + θ) for −ˆ θ < 0. For ψ ∈ C ∗ = C([0, τˆ], R2 ), the adjoint operators A∗ (0), A∗ (σ1 , σ2 ) of A(0) and A(σ1 , σ2 ), respectively, are defined as dψ − 0 < s ≤ τˆ, ds , A∗ (0)ψ(s) = 0 dη T (s)ψ(−s), s = 0, −ˆ τ
December 9, 2009
11:8
02508
Dynamics and Double Hopf Bifurcations of the Rose–Hindmarsh Model
and A∗ (σ1 , σ2 )ψ(s) 0, = 0 dη T (s, σ1 , σ2 )ψ(−s),
0 < s ≤ τˆ, s = 0.
−ˆ τ
For ϕ ∈ C and ψ ∈ C ∗ , the bilinear form is given as 0 θ ψ T (ξ − θ)dη(θ)ϕ(ξ)dξ. ψ, ϕ = ψ T (0)ϕ(0) − −ˆ τ
0
From the discussion shown in the preceding section, we know that the characteristic equation (5) has two pairs of pure imaginary characteristic values Λ = {iω1 , iω2 } at the double Hopf bifurcation point (ˆ x∗ , τˆ), and other characteristic values with
1
2x∗ d − iω θ 1 , q1 (θ) = 1 + iω1 e rS r + iω1 and
r − iω1 − c 2 c(1 − iω1 ) (r − iω1 ) r − iω1 eiω1 s , p1 (s) = l1 + im1 c(−1 + iω1 ) 1 r − iω2 − c 2 iω s c(1 − iω2 ) (r − iω2 ) r − iω e 2 , 2 p2 (s) = l2 + im2 c(−1 + iω2 ) 1
where li and mi (i = 1, 2) are given as follows: l1 = (−3τ ax∗2 − 2τ − 2rτ − 1)ω14
negative real parts. Hence, C can be decomposed as C = PΛ ⊕ QΛ , where PΛ is the characteristic space corresponding to eigenvalues Λ, and QΛ is the complementary subspace of PΛ . System (11) is an ordinary differential system on the Banach space C of functions from [−ˆ τ , 0] to R3 , which is bounded and continuous on (−ˆ τ , 0) with a possible jump discontinuity at 0 [Buono et al., 2003]. We now choose Φ(θ) = (q1 (θ),q 1 (θ), q2 (θ),q 2 (θ)), (12) Ψ(θ) = (p1 (s), p 1 (s), p2 (s), p 2 (s)), as the bases of PΛ and PΛ∗ , where p1 , q1 = 1, p1 , q 1 = 0, p1 , q2 = 0, p1 , q 2 = 0, p2 , q1 = 0, p2 ,q 1 = 0, p2 , q2 = 1, p2 , q 2 = 0. Then, we may choose 1 2x∗ d − iω θ 2 , q2 (θ) = 1 + iω2 e rS r + iω2 + (4rτ x∗ d + 6r 2 τ ax∗2 + 2r 2 τ Sc + r 2 τ + 2r 2 + 2r 2 τ x∗ d − 4x∗ dr + 6rτ ax∗2 + rτ Sc + 2r − 2rSc)ω1 , l2 = (−3τ ax∗2 − 2τ − 2rτ − 1)ω24 + (2r 2 τ + r 2 τ Sc + r 2 + 3τ ax∗2 + 12rτ ax∗2 + 2rτ Sc − rSc + 4r + 4rτ x∗ d + 2rτ + 3r 2 τ ax∗2 + 2τ x∗ d − 2dx∗ + 1)ω22 − r 2 − 2r 2 τ x∗ d − 3r 2 τ ax∗2 + rSc − r 2 τ Sc + 2x∗ dr 2 , m2 = ω25 τ + (−rτ Sc − 6rτ ax∗2 − r 2 τ − 4rτ − 2r − τ − 2 − 6τ ax∗2 − 2τ x∗ d)ω23
+ (2r 2 τ + r 2 τ Sc + r 2 + 3τ ax∗2 + 12rτ ax∗2
+ (4rτ x∗ d + 6r 2 τ ax∗2 + 2r 2 τ Sc
+ 2rτ Sc − rSc + 4r + 4rτ x∗ d + 2rτ
+ r 2 τ + 2r 2 + 2r 2 τ x∗ d − 4x∗ dr
+ 3r 2 τ ax∗2 + 2τ x∗ d − 2dx∗ + 1)ω12
+ 6rτ ax∗2 + rτ Sc + 2r − 2rSc)ω2 .
− r 2 − 2r 2 τ x∗ d − 3r 2 τ ax∗2 + rSc − r 2 τ Sc + 2x∗ dr 2 , m1 = ω15 τ + (−rτ Sc − 6rτ ax∗2 − r 2 τ − 4rτ − 2r − τ − 2 − 6τ ax∗2 − 2τ x∗ d)ω13
3739
Thus, explicit expressions for Φ and Ψ can be obtained by substituting pi and qi (i = 1, 2) into system (12). Let z = (z1 ,z 1 , z2 ,z 2 )T and vt ∈ QΛ ∩ τ , 0), R), and denote ut = Φz + vt . System C ([−ˆ
December 9, 2009
3740
11:8
02508
S. Ma et al.
(11) is equivalent to z = Bz + Ψ T (0)(A(xε , τε )(Φz + vt ) + εF (Φz + vt )), vt = A(0)vt + (I − Π)X0 (A(pε , τε ) × (Φz + vt ) + εF (Φz + vt )), where
iω1 0 B= 0 0
0 −iω1 0 0
0 0 iω2 0
The matrix B generates the torus group T2 = S1 × S1 whose action on C2 is given by
Then, the T2 -equivalent normal form, which is truncated to the quadratic order, becomes
0 0 . 0 −iω2
ˆε z1 + εb12 τˆε z1 , z1 = iω1 z1 + εb11 x z2 = iω2 z2 + εb21 x ˆε z2 + εb22 τˆε z2 ,
b11 = −
2(−bω1 e−iω1 τˆ − ibe−iω1 τˆ + id + 3aω1 x ˆ∗ + 3iaˆ x∗ )(r − iω1 )2 (i + ω1 ) , l1 + im1
b12 = −
−2ibˆ x∗ e−iω1 τˆ ω1 (r − iω1 )2 (i + ω1 )2 , l1 + im1
b21 = −
2(−bω2 e−iω2 τˆ − ibe−iω2 τˆ + id + 3aω2 x ˆ∗ + 3iaˆ x∗ )(r − iω2 )2 (i + ω2 ) , l2 + im2
b22 = −
−2ibˆ x∗ e−iω2 τˆ ω2 (r − iω2 )2 (i + ω2 )2 . l2 + im2
xε + Re(b12 )ˆ τε )r1 , r1 = (Re(b11 )ˆ xε + Re(b22 )ˆ τε )r2 . r2 = (Re(b21 )ˆ
(15)
If we remove the dependence on the unfolding parameters x ˆε and τˆε in system (13), we obtain =
iω1 z1 + p T1 (0)F (Φz
=
iω1 z1 + p T1 (0)fˆ(z1 ,z 1 , z2 , z 2 ), iω2 z2 + p T2 (0)F (Φz + vt )
=
+ vt ) (16)
= iω2 z2 + p T2 (0)fˆ(z1 ,z 1 , z2 , z 2 ), = Avt + (I − π)X0 (F (Φz + vt )),
where fˆ(z1 ,z 1 , z2 , z 2 ) = F (2 Re{z1 q1 (θ)} + 2 Re{z2 q2 (θ)} + vt ). We now expand vt = w(z, z z1 ,z 1 , z2 ,z 2 ) into the series wi,j,k,l z1i z j1 z2kz l2 , vt = w(z1 ,z 1 , z2 ,z 2 ) = i+j+k+l≥2
(14)
where
In polar coordinates z1 = r1 eiρ1 and z2 = r2 eiρ2 , the associated amplitude equations generated from system (14) are
z1 z2 vt
(θ1 , θ2 )(z1 , z2 ) = (eiθ1 z1 , eiθ2 z2 ).
(13)
where the coefficient vectors {wi,j,k,l } only depend 1 2 3 , wi,j,k,l , wi,j,k,l )T . Simion θ and wi,j,k,l = (wi,j,k,l larly, we are able to re-express the first two equations in system (16) as z1 = iω1 z1 + g1 (z1 ,z 1 , z2 ,z 2 ), (17) z2 = iω2 z2 + g2 (z1 ,z 1 , z2 ,z 2 ), with the expansions g1 (z1 , z 1 , z2 ,z 2 ) = g2 (z1 , z 1 , z2 ,z 2 ) =
i+j+k+l≥2
i+j+k+l≥2
1 gi,j,k,l z1i z j1 z2kz l2 , 2 gi,j,k,l z1i z j1 z2kz l2 ,
where all coefficients {gi,j,k,l } are real constants. From system (16), we can derive Avt − 2 Re{p T1 (0)fˆ(z1 , z 1 , z2 ,z 2 )q1 (θ)} − 2 Re{p T2 (0)fˆ(z1 ,z 1 , z2 ,z 2 )q2 (θ)}, (−τ ≤ θ < 0), dvt = dt Avt − 2 Re{p T1 (0)fˆ(z1 , z 1 , z2 ,z 2 )q1 (0)} − 2 Re{p T2 (0)fˆ(z1 ,z 1 , z2 ,z 2 )q2 (0)} + fˆ(z1 , z 1 , z2 ,z 2 ), (θ = 0). (18)
December 9, 2009
11:8
02508
Dynamics and Double Hopf Bifurcations of the Rose–Hindmarsh Model
The first equation of (18) can be rewritten as dw (z1 ,z 1 , z2 , z 2 ) dt = Aw(z1 ,z 1 , z2 , z 2 ) + H(z1 ,z 1 , z2 , z 2 ),
−τ ≤ θ < 0,
(19)
where H(z1 ,z 1 , z2 ,z 2 ) = −2 Re{p T1 (0)fˆ(z1 , z 1 , z2 ,z 2 )q1 (θ)} − 2 Re{p T2 (0)fˆ(z1 ,z 1 , z2 ,z 2 )q2 (θ)} Hi,j,k,lz1i z j1 z2kz l2 , = i+j+k+l≥2
k g2000
k g1100
k g0200
k g1010
k g1001
k g0110
k g0101
k g0020
k g0011
k g0002
3741
1 2 with the coefficients Hi,j,k,l = (Hi,j,k,l (θ), Hi,j,k,l (θ), 3 T Hi,j,k,l(θ)) . In order to determine the behavior of bifurcating solution, it is necessary to compute vt = w(z1 ,z 1 , z2 , z 2 ) through system (18). Then {gij } and {Hij } should be found first. Let us define that
φ(0) = z1 q1 (0) + z 1q 1 (0) + z2 q2 (0) + z 2q 2 (0) + w(z1 , z 1 , z2 , z 2 )|θ=0 , φ(−τ ) = z1 q1 (−τ ) + z 1q 1 (−τ ) + z2 q2 (−τ ) + z 2q 2 (−τ ) + w(z1 ,z 1 , z2 ,z 2 )|θ=−τ . From systems (16) and (17), by direct calculations for k = 1, 2, we have
√
ε(r + iωk )2 (−i + ωk ) (−bωk e−2iω1 τ + ibe−2iω1 τ − id + 3ax∗ ωk − 3iax∗ ), lk − imk √ 2 ε(r + iωk )2 (−i + ωk ) =− (−bωk − id + ib + 3ax∗ ωk − 3iax∗ ), lk − imk √ i ε(r + iωk )2 (−i + ωk ) = (−ibωk e2iω1 τ − be2iω1 τ + d + 3iax∗ ωk + 3ax∗ ), lk − imk √ 2i ε(r + iωk )2 (−i + ωk ) = (−ibωk e−iτ (ω1 +ω2 ) − be−iτ (ω1 +ω2 ) + d + 3iaωk x∗ + 3ax∗ ), lk − imk √ 2i ε(r + iωk )2 (−i + ωk ) =− (−ibω1 e−iτ (ω1 −ω2 ) − be−iτ (ω1 −ω2 ) + d + 3iaω1 x∗ + 3ax∗ ), lk − imk √ 2i ε(r + iωk )2 (−i + ωk ) = (−ibωk eiτ (ω1 −ω2 ) − beiτ (ω1 −ω2 ) + d + 3iaωk x∗ + 3ax∗ ), lk − imk √ 2 ε(r + iωk )2 (−i + ωk ) =− (−bωk eiτ (ω1 +ω2 ) + ibeiτ (ω1 +ω2 ) − id + 3aωk x∗ − 3iax∗ ), lk − imk √ − ε(r + iωk )2 (−i + ωk ) = (−bωk e−2iω2 τ + ibe−2iω2 τ − id + 3ax∗ ωk − 3iax∗ ), lk − imk √ 2 ε(r + iωk )2 (−i + ωk ) =− (−bωk + ib − id + 3aωk x∗ − 3iax∗ ), lk − imk √ i ε(r + iωk )2 (−i + ωk ) = (−ibωk e2iω2 τ − be2iω2 τ + d + 3iaωk x∗ + 3ax∗ ). lk − imk =−
Thus, for −τ ≤ θ < 0, we have H(z1 , z 1 , z2 ,z 2 ) = −2 Re{p T1 (0)fˆ(z1 ,z 1 , z2 , z 2 )q1 (θ)} − 2 Re{p T2 (0)fˆ(z1 ,z 1 , z2 , z 2 )q2 (θ)}, = −g1 (z1 ,z 1 , z2 , z 2 )q1 (θ) − g 1 (z1 , z 1 , z2 , z 2 )q 1 (θ) − g2 (z1 , z 1 , z2 ,z 2 )q2 (θ) − g 2 (z1 , z 1 , z2 ,z 2 )q 2 (θ).
December 9, 2009
3742
11:8
02508
S. Ma et al.
Comparing this formula with the expansion of H(z1 ,z 1 , z2 , z 2 ), we have 1 2 q1 (θ) − g 10200q 1 (θ) − g2000 q2 (θ) − g 20200q 2 (θ), H2000 (θ) = −g2000 1 2 q1 (θ) − g 11100q 1 (θ) − g1100 q2 (θ) − g 21100q 2 (θ), H1100 (θ) = −g1100 1 2 q1 (θ) − g 12000q 1 (θ) − g0200 q2 (θ) − g 22000q 2 (θ), H0200 (θ) = −g0200 1 2 q1 (θ) − g 10101q 1 (θ) − g1010 q2 (θ) − g 20101q 2 (θ), H1010 (θ) = −g1010 1 2 q1 (θ) − g 10110q 1 (θ) − g1001 q2 (θ) − g 20110q 2 (θ), H1001 (θ) = −g1001 1 2 q1 (θ) − g 11001q 1 (θ) − g0110 q2 (θ) − g 21001q 2 (θ), H0110 (θ) = −g0110 1 2 q1 (θ) − g 11010q 1 (θ) − g0101 q2 (θ) − g 21010q 2 (θ), H0101 (θ) = −g0101 1 2 q1 (θ) − g 10002q 1 (θ) − g0020 q2 (θ) − g 20002q 2 (θ), H0020 (θ) = −g0020 1 2 q1 (θ) − g 10011q 1 (θ) − g0011 q2 (θ) − g 20011q 2 (θ), H0011 (θ) = −g0011 1 2 q1 (θ) − g 10020q 1 (θ) − g0002 q2 (θ) − g 20020q 2 (θ). H0002 (θ) = −g0002
It can be derived from Eq. (19) that (A − 2iω1 )w2000 = −H2000 (θ), Aw1100 = −H1100 (θ), (A + 2iω1 )w0200 = −H0200 (θ), (A − iω1 − iω2 )w1010 = −H1010 (θ), (A − iω1 + iω2 )w1001 = −H1001 (θ), (A + iω1 − iω2 )w0110 = −H0110 (θ), (A + iω1 + iω2 )w0101 = −H0101 (θ), (A − 2iω2 )w0020 = −H0020 (θ), Aw0011 = −H0011 (θ), (A + 2iω2 )w0002 = −H0002 (θ). Taking account of formula (10) of operator A, we have dw2000 = 2iω1 w2000 (θ) − H2000 (θ), dθ dw1100 = −H1100 (θ), dθ dw1010 = i(ω1 + ω2 )w1010 (θ) − H1010 (θ), dθ dw1001 = i(ω1 − ω2 )w1001 (θ) − H1001 (θ), dθ dw0110 = i(−ω1 + ω2 )w0110 (θ) − H0110 (θ), dθ dw0020 = 2iω2 w0020 (θ) − H0020 (θ), dθ dw0011 = −H0011 (θ). dθ
(20)
December 9, 2009
11:8
02508
Dynamics and Double Hopf Bifurcations of the Rose–Hindmarsh Model
3743
Similarly, from the second equation of system (18), it can be derived that 0 [dη(θ)]w2000 −ˆ τ 0 [dη(θ)]w1100 −ˆ τ 0 [dη(θ)]w1010 −ˆ τ 0 [dη(θ)]w1001 −ˆ τ 0 [dη(θ)]w0110 −ˆ τ 0 [dη(θ)]w0020 −ˆ τ 0 [dη(θ)]w0011 −ˆτ
be−2iω1 τ − 3ax∗ √ −d = 2iω1 w2000 (0) − H2000 (0) − ε , 0 2b − 6ax∗ √ = −H1100 (0) − ε −2d , 0 −6ax∗ + 2be−i(ω1 +ω2 )τ √ = i(ω1 + ω2 )w1010 (0) − H1010 (0) − ε −2d , 0 −6ax∗ + 2be−i(ω1 −ω2 )τ √ = i(ω1 − ω2 )w1001 (0) − H1001 (0) − ε −2d , 0 −6ax∗ + 2bei(ω1 −ω2 )τ √ = i(−ω1 + ω2 )w0110 (0) − H0110 (0) − ε −2d , 0 −2iω τ 2 − 3ax∗ be √ −d = 2iω2 w0020 (0) − H0020 (0) − ε , 0 −6ax∗ + 2b √ −2d = −H0011 (0) − ε . 0
Integrating system (20) under conditions (21), we obtain 1 (θ) = w2000
1 (6iω12g 20200 e−iω2 θ − 3iω1 ω2g 20200 e−iω2 θ 3ω1 (2ω1 + ω2 )(2ω1 − ω2 ) 1 1 + 12iω12 g2000 eiω1 θ − 3iω22 g2000 eiω1 θ + 4iω12g 10200 e−iω1 θ − iω22g 10200 e−iω1 θ 2 2 + 6iω12 g2000 eiω2 θ + 3iω1 ω2 g2000 eiω2 θ + 12E1 ω13 e2iω1 θ − 3E1 ω1 ω22 e2iω1 θ ),
1 (θ) = w1100
1 1 2 (−iω2 g1100 eiω1 θ + ig 11100 ω2 e−iω1 θ − ig1100 ω1 eiω2 θ + ig 21100 ω1 e−iω2 θ ) + E2 , ω1 ω2
1 (θ) = w1010
1 (2iω12 ω2g 20101 e−iω2 θ + iω1 ω22g 20101 e−iω2 θ ω1 ω2 (ω1 + 2ω2 )(2ω1 + ω2 ) 1 1 1 + 2iω13 eiω1 θ g1010 + 5iω12 ω2 g1010 eiω1 θ + 2iω1 ω22 g1010 eiω1 θ + iω2 ω12g 10101 e−iω1 θ 2 2 2 + 2iω22 ω1g 10101 e−iω1 θ + 2iω2 ω12 g1010 eiω2 θ + 5iω22 ω1 g1010 eiω2 θ + 2iω23 g1010 eiω2 θ
+ 2E3 ω2 ω13 ei(ω1 +ω2 )θ + 5E3 ω22 ω12 ei(ω1 +ω2 )θ + 2E3 ω23 ω1 ei(ω1 +ω2 )θ ),
(21)
December 9, 2009
3744
11:8
02508
S. Ma et al. 1 w1001 (θ) =
i (2ω2 ω12g 20110 e−iω2 θ − 5ω22 ω1g 20110 e−iω2 θ ω1 ω2 (2ω1 − ω2 )(ω1 − 2ω2 ) 1 1 1 + 2ω23g 20110 e−iω2 θ − 2ω13 g1001 eiω1 θ + 5ω12 ω2 g1001 eiω1 θ − 2ω1 ω22 g1001 eiω1 θ 2 2 + ω12 ω2g 10110 e−iω1 θ − 2ω1 ω22g 10110 e−iω1 θ + 2ω12 ω2 g1001 eiω2 θ − ω1 ω22 g1001 eiω2 θ
− 2iE4 ω13 ω2 ei(ω1 −ω2 )θ + 5iE4 ω12 ω22 ei(ω1 −ω2 )θ − 2iE4 ω1 ω23 ei(ω1 −ω2 )θ ), 1 w0110 (θ) =
1 (−2iω2 ω12g 21001 e−iω2 θ + iω22 ω1g 21001 e−iω2 θ ω1 ω2 (ω1 − 2ω2 )(2ω1 − ω2 ) 1 1 − iω2 ω12 g0110 eiω1 θ + 2iω22 ω1 g0110 eiω1 θ + 2iω13g 11001 e−iω1 θ − 5iω12 ω2g 11001 e−iω1 θ 2 2 2 + 2iω1 ω22g 11001 e−iω1 θ − 2iω2 ω12 g0110 eiω2 θ + 5iω22 ω1 g0110 eiω2 θ − 2iω23 g0110 eiω2 θ
+ 2E5 ω2 ω13 e−i(ω1 −ω2 )θ − 5E5 ω22 ω12 e−i(ω1 −ω2 )θ + 2E5 ω23 ω1 e−i(ω1 −ω2 )θ ), 1 (θ) = w0020
i (−4ω22g 20002 e−iω2 θ + ω12g 20002 e−iω2 θ 3ω2 (−2ω2 + ω1 )(2ω2 + ω1 ) 1 1 − 6ω22 g0020 eiω1 θ − 3ω2 ω1 g0020 eiω1 θ − 6ω22g 10002 e−iω1 θ + 3ω2 ω1g 10002 e−iω1 θ 2 2 − 12ω22 g0020 eiω2 θ + 3ω12 g0020 eiω2 θ + 12iE6 ω23 e2iω2 θ − 3iE6 ω2 ω12 e2iω2 θ ),
1 (θ) = w0011
1 1 2 (−iω2 g0011 eiω1 θ + ig 10011 ω2 e−iω1 θ − ig0011 ω1 eiω2 θ + ig 20011 ω1 e−iω2 θ ) + E7 , ω1 ω2 (22)
where Ei (i = 1, . . . , 7) are determined by the following equations, crS 2x∗ d 1 1 − (0) + 2bx∗ w2000 (−τ ) w2000 −2iω1 − 3ax∗2 − r + 2iω1 1 + 2iω1 =−
2 3 √ √ d cH2000 (0) √ (0) H2000 1 + ε + − H2000 (0) − εbe−2iω1 τ + 3 εax∗ , 1 + 2iω1 1 + 2iω1 r + 2iω1
1 1 (0) + 2bx∗ w1100 (−τ ) (−3ax∗2 − 2x∗ d − cS)w1100 3 √ √ √ cH1100 (0) 1 2 − H1100 (0) + 2 εd + (0) − 2 εb + 6 εax∗ , = −H1100 r crS 2x∗ d ∗2 1 1 − (0) + 2bx∗ w1010 (−τ ) w1010 −iω1 − iω2 − 3ax − r + iω1 + iω2 1 + iω1 + iω2
√ 2 3 √ √ 2d ε cH1010 (0) (0) H1010 1 + + − H1010 (0) − 2 εbe−i(ω1 +ω2 )τ + 6 εax∗ , =− 1 + iω1 + iω2 1 + iω1 + iω2 r + iω1 + iω2 −iω1 + iω2 − 3ax∗2 −
crS 2x∗ d − r + iω1 − iω2 1 + iω1 − iω2
1 1 (0) + 2bx∗ w1001 (−τ ) w1001
√ 2 3 √ √ 2d ε cH1001 (0) (0) H1001 1 + + − H1001 (0) − 2 εbe−i(ω1 −ω2 )τ + 6 εax∗ , =− 1 + iω1 − iω2 1 + iω1 − iω2 r + iω1 − iω2
December 9, 2009
11:8
02508
Dynamics and Double Hopf Bifurcations of the Rose–Hindmarsh Model
3745
crS 2x∗ d 1 1 − − (0) + 2bx∗ w0110 (−τ ) iω1 − iω2 − w0110 r − iω1 + iω2 1 − iω1 + iω2 √ 2 3 √ √ 2d ε cH0110 (0) (0) H0110 1 + + − H0110 (0) − 2 εbei(ω1 −ω2 )τ + 6 εax∗ , =− 1 − iω1 + iω2 1 − iω1 + iω2 r − iω1 + iω2 2x∗ d crS ∗2 1 1 − − 2iω2 w0020 (0) + 2bx∗ w0020 (−τ ) −3ax − 1 + 2iω2 r + 2iω2 3ax∗2
=−
2 3 √ √ d cH0020 (0) √ (0) H0020 1 + ε + − H0020 (0) + 3 εax∗ − εbe−2iω2 τ , 1 + 2iω2 1 + 2iω2 r + 2iω2
1 1 (0) + 2bx∗ w0011 (−τ ) (−3ax∗2 − 2x∗ d − cS)w0011 3 √ √ √ cH0011 (0) 1 2 − H0011 (0) + 2 εd + (0) + 6 εax∗ − 2 εb. = −H0011 r
From (22), we deduce 1 = g2100
√ √ (r + iω1 )2 (−i + ω1 ) √ 1 1 1 (2 εbω1 e−iω1 τ w1100 (−τ ) − 2i εbe−iω1 τ w1100 (−τ ) + 2 εbeiω1 τ ω1 w2000 (−τ ) l1 − im1 √ √ √ √ √ 1 1 1 1 1 (−τ ) + 2i εdw2000 (0) + 2i εdw1100 (0) − 6 εax∗ ω1 w2000 (0) + 6i εax∗ w2000 (0) − 2i εbeiω1 τ bw2000 √ √ 1 1 (0) + 6i εax∗ w1100 (0) − 3εaω1 + 3iεa), − 6 εax∗ ω1 w1100 √ √ 2(r + iω1 )2 (−i + ω1 ) √ (− εbω1 eiω2 τ w1010 (−τ ) − εbω1 e−iω1 τ w0011 (−τ ) − εbω1 e−iω2 τ w1001 (−τ ) l1 − im1 √ √ √ √ − i εdw1001 (0) − 3ia εx∗ w1010 (0) + 3εaω1 + i εbe−iω1 τ w0011 (−τ ) + i εbeiω2 τ w1010 (−τ ) √ √ √ √ − 3ia εx∗ w1001 (0) − i εdw1010 (0) + 3 εaω1 x∗ w0011 (0) + 3 εaω1 x∗ w1010 (0) − 3iεa √ √ √ √ + 3 εaω1 x∗ w1001 (0) − 3i εax∗ w0011 (0) + i εbe−iω2 τ w1001 (−τ ) − i εdw0011 (0)),
1 =− g1011
√ √ 2(r + iω2 )2 (−i + ω2 ) √ (− εbω2 eiω1 τ w1010 (−τ ) − 3i εax∗ w1100 (0) − εbω2 e−iω1 τ w0110 (−τ ) l2 − im2 √ √ √ √ −iω + i εbe 2 τ w1100 (−τ ) − εbω2 e−iω2 τ w1100 (−τ ) + i εbe−iω1 τ w0110 (−τ ) + i εbeiω1 τ w1010 (−τ ) √ √ √ √ √ − i εdw0110 (0) − 3i εax∗ w0110 (0) + 3 εaω2 x∗ w1010 (0) + 3 εaω2 x∗ w0110 (0) + 3 εaω2 x∗ w1100 (0) √ √ √ − 3i εax∗ w1010 (0) − i εdw1010 (0) − i εdw1100 (0) + 3εaω2 − 3iεa),
2 =− g1110
√ √ √ (r + iω2 )2 (−i + ω2 ) (−2 εbω2 e−iω2 τ w0011 (−τ ) + 2i εbe−iω2 τ w0011 (−τ ) − 2 εbω2 eiω2 τ w0020 (−τ ) l2 − im2 √ √ √ √ √ iω2 τ w0020 (−τ ) − 2i εdw0011 (0) − 2i εdw0020 (0) + 6 εaω2 x∗ w0011 (0) − 6i εax∗ w0011 (0) + 2i εbe √ √ + 6 εaω2 x∗ w0020 (0) − 6i εax∗ w0020 (0) + 3εaω2 − 3iεa).
2 =− g0021
Using the standard normalization technique as described in [Ma et al., 2008a, 2008b], we set z1 = η1 + P˜1 (η1 ,η 1 , η2 ,η 2 ), η1 = iω1 η1 + B11 η12η 1 + B12 η1 η2η 2 ,
z2 = η2 + P˜2 (η1 ,η 1 , η2 ,η 2 ), η2 = iω2 η2 + B21 η1η 1 η2 + B22 η22η 2 ,
(23)
December 9, 2009
3746
11:8
02508
S. Ma et al.
where
P˜1 (η1 ,η 1 , η2 ,η 2 ) =
j+k+l+m≥2
P˜2 (η1 ,η 1 , η2 ,η 2 ) =
j+k+l+m≥2
1 η1jη k1 η2l η m P˜jklm 2 , 2 η1jη k1 η2l η m P˜jklm 2 .
Substituting (23) into system (17), then comparing the corresponding coefficients of both sides to determine the appropriate functions P˜1 and P˜2 , so as to round off the nonresonant terms, we obtain 1 + B11 = g2100
−
i 1 i 1 i 1 2 1 g1100 g2000 + (g1010 g1100 − g1001 g 21100 ) − g1 g 2 ω1 ω2 2ω1 + ω2 0101 0200
i i 1 2 2i 1 2 1 2 g0110 g2000 − |g1100 | − |g | , 2ω1 − ω2 ω1 3ω1 0200
1 + B12 = g1011
i 1 i 2 1 1 1 1 1 (g1010 g0011 − g1001 g 20011 ) + (2g2000 g0011 − g1100 g 10011 − g0011 g 20110 ω2 ω1
1 2 −g0011 g1010 )− 2 + B21 = g1110
i 1 i 2 2 2 2 1 2 (g g2 − g 11100 g0110 ) + (2g0020 g1100 − g0011 g 21100 − g1010 g1100 ω1 1100 1010 ω2
2 −g 11001 g1100 )+ 2 + B22 = g0021
−
2i 2i i i g1 g 2 − g1 g2 − |g1 |2 − |g1 |2 , ω1 + 2ω2 0002 0101 ω1 − 2ω2 0020 1001 2ω1 − ω2 0110 2ω1 + ω2 0101
2i 2i i i 1 2 2 2 g0110 g2000 − g 10101 g0200 − |g1001 |2 − |g2 |2 , 2ω1 − ω2 2ω1 + ω2 2ω2 − ω1 ω1 + 2ω2 0101
i 1 i 2 i 2 2 2 (g0011 g1010 − g 10011 g0110 ) + g0011 g0020 − g1 g2 ω1 ω2 2ω2 − ω1 0020 1001
i i 2 2 2i 2 2 2 g 10002 g0101 − |g0011 | − |g | . 2ω2 + ω1 ω2 3ω2 0002
Let η1 = r1 eiθ1 ,
η2 = r2 eiθ2 .
In polar coordinates, the amplitude equation generated from system (17) is r1 = r1 [Re(B11 )r12 + Re(B12 )r22 ], r2 = r2 [Re(B21 )r12 + Re(B22 )r22 ]. Combining the results of (15) and (24), the formal normal form arising from system (8) is ˆε η1 + εb12 τˆε η1 η = iω1 η1 + εb11 x 1 + B11 η12η 1 + B12 η1 η2η 2 , ˆε η2 + εb22 τˆε η2 η2 = iω2 η2 + εb21 x + B21 η1η 1 η2 + B22 η22η 2 . Using polar coordinates η1 = r1 eiθ1 , η2 = r2 eiθ2 , system (25) can be transformed into r1 = r1 [µ1 + Re(B11 )r12 + Re(B12 )r22 ], r = r2 [µ2 + Re(B21 )r 2 + Re(B22 )r 2 ], 2 1 2 θ1 = ω1 + ν1 + Im(B11 )r12 + Im(B12 )r22 , θ2 = ω2 + ν2 + Im(B21 )r12 + Im(B22 )r22 , xε + Re(bi2 )ˆ τε , νi = Im(bi1 )ˆ xε + Im(bi2 )ˆ τε for i = 1, 2. where µi = Re(bi1 )ˆ
(24)
(25)
(26)
December 9, 2009
11:8
02508
Dynamics and Double Hopf Bifurcations of the Rose–Hindmarsh Model
3747
In the following, we give two examples to illustrate that under certain conditions, no Hopf bifurcations occur for the limit cycle of system (4). Therefore, the unique attractor for system (4) is asymptotically stable. Example 1. Choose S = 1.0, x∗ = 0.1768, τ = 6.8693, then the characteristic equation has two pairs of
imaginary roots with ω1 = 0.9107, ω2 = 1.6810, and l1 = −0.572128, m1 = −12.544266, l2 = −88.8017, m2 = −13.87852, b11 = 0.02409516 + 0.6241099i, b21 = 0.8705159 + 0.743506i, √ 1 = ε(0.01632 − 0.3217559i), g2000 √ 1 = ε(0.01830 − 0.28222i), g0200 √ 1 = ε(0.39171 − 0.014319i), g1001 √ 1 = ε(0.42643 + 0.00472i), g0101 √ 1 = ε(0.033536 − 0.6039247i), g0011 √ 2 = ε(0.314199 − 0.036238i), g2000 √ 2 = ε(0.277067 − 0.049426i), g0200 √ 2 = ε(−0.08694 + 0.019076i), g1001 √ 2 = ε(−0.117646 + 0.04379i), g0101 √ 2 = ε(0.5916266 − 0.086679i), g0011
b12 = −0.013262 − 0.1159697i, b22 = 0.0212975 − 0.2134128i, √ 1 g1100 = ε(0.033536 − 0.6039247i), √ 1 g1010 = ε(0.36367 − 1.24877i), √ 1 g0110 = ε(0.33103 − 1.2263575i), √ 1 g0020 = ε(0.52571 − 0.651452i), √ 1 g0002 = ε(0.558073 − 0.005054i), √ 2 g1100 = ε(0.591627 − 0.086679i), √ 2 g1010 = ε(1.059699 + 0.461955i), √ 2 g0110 = ε(1.05146 + 0.423406i), √ 2 g0020 = ε(0.424199 + 0.55770i), √ 2 g0002 = ε(−0.182930 + 0.342064i).
Following the procedure introduced as above, we have √ √ E2 = − ε3.584712844, E1 = ε(1.006098900 − 2.909499977i), √ √ E4 = ε(−10.41936727 − 2.174529644i), E3 = ε(8.145198479 − 2.422005382i), √ √ E6 = ε(1.696799807 + 1.776121238i), E5 = ε(−10.41936727 + 2.174529644i), √ E7 = − ε3.267946002, and
√ 1 (0) = ε(0.804969 − 0.595855i), w2000 √ 1 (0) = − ε4.79444, w1100 √ 1 (0) = ε(8.3924166 − .947833i), w1010 √ 1 (0) = ε(−18.6899 + 1.140244i), w1001 √ 1 (0) = ε(−18.6899 − 1.14024i), w0110 √ 1 (0) = ε(1.697437 + 2.33727i), w0020 √ 1 (0) = − ε4.69736, w0011
√ 1 w2000 (−τ ) = ε(−0.638051 − 1.97351i), √ 1 w1100 (−τ ) = − ε4.249727, √ 1 w1010 (−τ ) = ε(5.641275 + 6.350949i), √ 1 w1001 (−τ ) = ε(−14.9901 + 10.1277i), √ 1 w0110 (−τ ) = ε(−14.9901 − 10.1277i), √ 1 w0020 (−τ ) = ε(−2.47207 + .83983i), √ 1 w0011 (−τ ) = − ε4.04632.
Hence, we can deduce that 1 = (−3.00532 − 4.099777i)ε, g2100 2 g1110 = (−10.76744 − 6.98252i)ε,
1 g1011 = (−0.06603 − 15.8757i)ε, 2 g0021 = (−3.94492 − 0.2161815i)ε,
which lead to B11 = (−5.5847 − 5.156339i)ε, B21 = (−4.544096 − 5.67959i)ε,
B12 = (−1.259567 − 28.3088i)ε, B22 = (−2.687968 + 0.3059236i)ε.
Furthermore, we find b11 = 0.018523 − 0.58490648i, b21 = −0.10710 + 0.0015i,
b12 = −0.01957 + 0.11507i, b22 = −0.20014 − 0.077108i.
December 9, 2009
02508
S. Ma et al.
3748
x(t)
11:8
3
0
2
y(t)−10 −20
1
−30
0
−40
−1
−50 −60
−2
−70
−3
−80
−4 −5
−90 320
340
360
380
−100
400
320
340
360
380
400
420
t
440
t
5.5
z(t) 5
5.5
z(t)
4.5
5
4.5
4
4 3.5 3.5 0
3
−50 2.5 0
100
Fig. 3.
200
300
400
500
t
y(t)
600
−100
−6
−4
−2
0
2
x(t)
The time series solution of system (4) and phase portraits at τ = 6.9 and x∗ = 0.1768.
We now discuss the restriction on phase portraits near the double Hopf bifurcation point (x∗ , τ ) and rewrite system (26) as r1 = r1 (µ1 + ar12 + br22 ), r2 = r2 (µ2 + cr12 + dr22 ), xε + Re(bi2 )ˆ τε , and where µi = Re(bi1 )ˆ Re(B11 ) = −1, |B11 |
b=
Re(B12 ) = −0.4993579, |B22 |
Re(B21 ) = −0.033449, |B11 |
d=
Re(B22 ) = −1. |B22 |
a= c=
We observe that no Hopf bifurcation occurs for the limit cycle in this case. Therefore, the unique attractor for system (4) is asymptotically stable, as shown in Fig. 3.
December 9, 2009
11:8
02508
Dynamics and Double Hopf Bifurcations of the Rose–Hindmarsh Model
3749
Example 2. Choose S = 1.0, x∗ = 1.8207, τ = 2.9934, then the characteristic equation has two pairs of
imaginary roots with ω1 = 4.1953, ω2 = 6.1936, and l1 = −10203.38178, m1 = −4904.083215, l2 = −51514.46280, m2 = −1015.974613, b11 = −0.156559 + 0.06985i, b21 = −0.12339 + 0.14264i, √ 1 g2000 = ε(−0.0210699 + 0.08338i), √ 1 = ε(−0.02331 + 0.0816597i), g0200 √ 1 = √ε(−0.090015 + 0.139878i), g1001 1 g0101 = ε(−0.09268 + 0.13892i), √ 1 = ε(−0.04437 + .165024i), g0011 √ 2 g2000 = ε(−0.06336 + 0.0484i), √ 2 = ε(−0.06432 + 0.04570i), g0200 √ 2 = ε(−0.15256 + 0.04753i), g1001 √ 2 g0101 = ε(−.154296 + 0.045248i), √ 2 = ε(−.127664 + 0.094103i), g0011
b12 = −0.016267 − 1.325248i, b22 = 0.0176 − 1.98267i, √ 1 g1100 = ε(−0.04437 + 0.165024i), √ 1 g1010 = ε(−0.0066 + 0.2049i), √ 1 g0110 = √ε(−0.00825 + 0.202589i), 1 g0020 = ε(0.00797 + 0.125942i), √ 1 g0002 = ε(−0.07195 + 0.0646i), √ 2 g1100 = ε(−.12766 + 0.0941i), √ 2 g1010 = ε(−0.1176 + 0.14883i), √ 2 g0110 = ε(−0.1177 + 0.145965i), √ 2 g0020 = ε(−0.06219 + 0.10060i), √ 2 g0002 = ε(−0.0962 + 0.004395i).
Following the procedure introduced as above, we have √ √ E2 = ε0.62986142, E1 = ε(0.1604619 + 0.398838i), √ √ E4 = ε(−0.91337 + 0.52039i), E3 = ε(0.6055876 + 0.98356i), √ √ E6 = ε(0.28760 + 0.73738i), E5 = ε(−0.91337 − 0.5203899i), √ E7 = − ε0.864976, and 1 (0) = w2000
√
ε(0.12817 + 0.35871i), √ 1 (0) = − ε0.73892, w1100 √ 1 w1010 (0) = ε(0.549281 + .938794i), √ 1 (0) = ε(−0.75798 + 0.52173i), w1001 √ 1 (0) = ε(−0.75798 − 0.52173i), w0110 √ 1 (0) = ε(0.26012 + 0.71879i), w0020 √ 1 (0) = − ε0.75592, w0011
1 w2000 (−τ ) =
√
ε(0.13001 + 0.354975i), √ 1 w1100 (−τ ) = − ε0.724738, √ 1 w1010 (−τ ) = ε(0.21885 + 1.06856i), √ 1 w1001 (−τ ) = ε(−0.57858 + 0.76051i), √ 1 w0110 (−τ ) = ε(−0.57858 − 0.76051i), √ 1 w0020 (−τ ) = ε(−0.2191 + 0.74435i), √ 1 w0011 (−τ ) = − ε0.7701.
Consequently, we have 1 = (0.017687 − 0.102537i)ε, g2100 2 = (−0.11286 − 0.13893i)ε, g1110
1 g1011 = (−0.1316795 − 0.285369i)ε, 2 g0021 = (−0.128979 − 0.130854i)ε,
which lead to B11 = (−0.00728 − 0.103909i)ε,
B12 = (−0.15409 − 0.318484i)ε,
B21 = (−0.08765 − 0.17090i)ε,
B22 = (−0.11857 − 0.1477i)ε.
Moreover, we obtain b11 = −0.156617 − 0.06702i, b21 = −0.15513 − 0.04030i,
b12 = −0.03789 + 1.3248i, b22 = −1.13648 + 1.6247i.
December 9, 2009
11:8
02508
S. Ma et al.
3750
10
x(t)
1
y(t)
9 8
0.5
7 6
0
5 −0.5
4
−1
3 2
−1.5 −2 700
1 705
710
715
0 700
720
705
710
715
720
t
t
5.5
z(t) 5 4.5 2.8
4
2.6
z(t)
3.5
2.4
1
3
2.2
2.5
2 10
0
2 0
200
Fig. 4.
400
600
800
t
8
−1
6
y(t)
1000
4
2
0
−2
x(t)
−2
The time series solution of system (4) and phase portraits at τ = 3.1 and x∗ = 1.8207.
Now, we discuss the restriction on phase portraits near the double Hopf bifurcation point (x∗ , τ ) and rewrite system (26) as r1 = r1 (µ1 + ar12 + br22 ), r2 = r2 (µ2 + cr12 + dr22 ), xε + Re(bi2 )ˆ τε , and where µi = Re(bi1 )ˆ Re(B11 ) = −1, |B11 |
b=
Re(B12 ) = −0.813467, |B22 |
Re(B21 ) = −0.8414566, |B11 |
d=
Re(B22 ) = −1. |B22 |
a= c=
Clearly, one can see that no Hopf bifurcation for the limit cycle occurs in this case. Therefore, the unique attractor for system (4) is asymptotically stable, as illustrated in Fig. 4.
4. Conclusion The dynamical behavior of the Rose–Hindmarsh model with time delay was intensively investigated. By applying the generalized Sturm criterion, which is a powerful tool for real polynomials with constant
December 9, 2009
11:8
02508
Dynamics and Double Hopf Bifurcations of the Rose–Hindmarsh Model
coefficients, classifications of the imaginary roots of the characteristic equation were presented. It was proved that the Rose–Hindmarsh system was absolutely stable for any time delay if corresponding parameters lie in the regime of II(b) and IV. Otherwise, Hopf bifurcation curves were drawn in the (x∗ , τ ) plane by using the software DDE-Biftool. The double Hopf bifurcation points were acquired at the intersection points of two Hopf bifurcation curves. The normal form and universal unfolding at double Hopf bifurcation points were studied by the theory of center manifold theory. Two examples were demonstrated and computational results of normal form at double Hopf bifurcation points reveal that the associated unique attractor was asymptotically stable.
References Buono, P. L. & B´elair, J. [2003] “Restrictions and unfolding of double Hopf bifurcation in functional differential equations,” J. Diff. Eqs. 189, 234–266. Campell, S. A. & B´elair, J. [1995a] “Analytical and symbolically-assisted investigation of Hopf bifurcations in delay-differential equations,” Canad. Appl. Math. Quart. 3, 137–154. Campell, S. A., B´elair, J., Ohira, T. & Milton, J. [1995b] “Limit cycles, tori, and complex dynamics in a second-order differential equations with delayed negative feedback,” J. Dyn. Diff. Eqs. 7, 213–236. Elphick, C., Tiraopegui, E., Brachet, M. E., Coullet, P. & Iooss, G. [1987] “A simple global characterization for normal forms of singular vector fields,” Phys. D 29, 95–127. Grebogi, C., Ott, E., Romeiras, F. & Yorke, J. A. [1986] “Critical exponent of chaotic transients in nonlinear dynamical systems,” Phys. Rev. Lett. 57, 1284–1287. Grebogi, C., Ott, E., Romeiras, F. & Yorke, J. A. [1987] “Critical exponents for crisis-induced intermittency,” Phys. Rev. A 36, 5365–5380. Guckenheimer, J. M. & Holmes, P. J. [1983] Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields (Springer, NY). Hindmarsh, J. L. & Rose, R. M. [1982] “A model of the nerve impulse using two first-order differential equations,” Nature 296, 162–164.
3751
Hindmarsh, J. L. & Rose, R. M. [1984] “A model of neuronal bursting using three coupled first-order differential equations,” Proc. Roy. Soc. Lond. B 221, 87–102. Holden, A. V. & Fan, Y. S. [1992a] “From simple to simple bursting oscillatory behavior via chaos in the Rose–Hindmarsh model for neuronal activity,” Chaos Solit. Fract. 2, 221–236. Holden, A. V. & Fan, Y. S. [1992b] “From simple to complex oscillatory behavior via intermittent chaos in the Rose–Hindmarsh model for neuronal activity,” Chaos Solit. Fract. 2, 346–369. Ma, S. Q., Lu, Q. S. & Hogan, S. J. [2007] “Double Hopf bifurcation for Stuart-Landau system with nonlinear delay feedback and delay dependant parameters,” Adv. Compl. Syst. 10, 1–27. Ma, S. Q., Feng, Z. & Lu, Q. S. [2008a] “A two-patch ecological system with nonlinear transfer rate and noise effect,” Dyn. PDEs 5, 281–298. Ma, S. Q., Lu, Q. S. & Feng, Z. [2008b] “Double Hopf bifurcation for van der Pol-Duffing oscillator with parametric delay feedback control,” J. Math. Anal. Appl. 338, 993–1007. Pomeau, Y. & Manneville, P. “Intermittent transition to turbulence in dissipative dynamic system,” Commun. Math. Phys. 74, 189–197. Shi, X., Lu, Q. S. & Chen, G. [2007] “Anti-phase synchronization of inhibitorily coupled neurons,” Int. J. Bifurcation and Chaos 17, 4355–4364. Wang, Z. & Hu, H. Y. [1999] “Delay independent stability of retarded dynamic systems of multiple degrees of freedom,” J. Sound Vibr. 226, 57–81. Wang, Z. & Hu, H. Y. [2000] “Stability switches of timedelayed dynamic systems with unknown parameters,” J. Sound Vibr. 233, 215–233. Xu, J. & Lu, Q. S. [1999] “Hopf bifurcation of time-delay Li´enard equations,” Int. J. Bifurcation and Chaos 9, 939–951. Yang, L., Hou, X. R. & Zeng, Z. B. [1996a] “A complete discrimination system for polynomial,” Sci. in China (Series E) 26, 424–441. Yang, L., Hou, X. R. & Zeng, Z. B. [1996b] Nonlinear Algebraic Equation Systems and Automated Theory (in Chinese) (Shanghai Scientific and Technological Education Press, Shanghai). Zhao, X. Q. & Zou, X. [2001] “Threshold dynamics in a delayed SIS epidemic model,” J. Math. Anal. Appl. 257, 282–291.
11:8
02508
International Journal of Bifurcation and Chaos, Vol. 19, No. 11 (2009) 3733–3751 c World Scientific Publishing Company
DYNAMICS AND DOUBLE HOPF BIFURCATIONS OF THE ROSE–HINDMARSH MODEL WITH TIME DELAY* SUQI MA Department of Mathematics, Chinese Agricultural University, Beijing 100083, P. R. China ZHAOSHENG FENG† Department of Mathematics, University of Texas-Pan American, Edinburg, TX 78539, USA [email protected] QISHAI LU School of Science, Beijing University of Aeronautics and Astronautics, Beijing 100083, P. R. China Received January 11, 2008; Revised April 22, 2009 In this paper, we are concerned with the Rose–Hindmarsh model with time delay. By applying the generalized Sturm criterion, a number of imaginary roots of the characteristic equation are classified. The absolutely stable regions for any value of time delay are detected. By the continuous software DDE-Biftool, both the Hopf bifurcation curves and double Hopf bifurcation points are illustrated in parametric spaces. The normal form and universal unfolding at double Hopf bifurcation points are considered by the center manifold method. Some examples also indicate that the corresponding unique attractor near each double Hopf point is asymptotically stable. Keywords: Hopf bifurcation; Rose–Hindmarsh model; normal form; center manifold method; autonomous system; elliptic function.
1. Introduction
and the Rose–Hindmarsh model is described as
The electric activity of neurons is of great concern. A typical example is given by the Rose–Hindmarsh model of action potential which has continuously attracted considerable attention in the past decades [Hindmarsh et al., 1982, 1984]. In [Hindmarsh et al., 1984], in order to simulate the bursting behavior of neurons, an adaption current z was introduced to the two-variable model of the excitable membrane,
x = y − ax3 + bx(t)2 − cz + Iapp , y = c − dx2 − y, z = r(S(x − χ) − z), where the variable x represents membrane potential, y represents a recovery variable, z represents an adaption current, and Iapp represents an applied current. Here r is the parameter of control and
∗ Main contents have been presented at the Seventh AIMS International Conference on Dynamical Systems, Differential Equations and Applications, University of Texas, Arlington, 18–21 May 2008. This work is supported by UTPA Faculty Research Council Grant 119100 and NSF of China under 10771212. † Author for correspondence
3733
December 9, 2009
3734
11:8
02508
S. Ma et al.
all other parameters a, b, c, d, S and χ are real constants. Dynamical behaviors of the Rose–Hindmarsh model have been intensively studied. In [Holden et al., 1992a, 1992b] both periodic and aperiodic “bursting” oscillations were considered which describe the behavior of a period of rapid spiking followed by a quiescent period. Numerical simulations of evolution from simple periodic bursting to complex periodic bursting were considered in [Holden et al., 1992a]. Transition from periodicity to aperiodicity via an alternation between the perioddoubling bifurcation sequence and the saddle-node bifurcation was presented in [Holden et al., 1992b]. It was shown [Pomeau et al., 1980] and [Grebogi et al., 1986, 1987] that transition from intermittent chaos to complex oscillatory behavior was caused by a crisis of the chaotic attractor, and after the crisis periodicity was maintained by period-adding bifurcations. Mathematical models with time delays have been broadly applied in many scientific fields such as neuronscience [Shi et al., 2007], biology [Zhao et al., 2001] and [Xu et al., 1999], and mechanics [Ma et al., 2007]. We see that in many cases controlled systems with delay feedback may have double Hopf bifurcations on the margin of “amplitude death origin” [Ma et al., 2007, 2008]. Such kind of phenomena has also been observed in [Campell et al., 1995b]. Normal form and universal unfolding in a vector field at nonresonant double Hopf bifurcation points for particular retarded functional differential equations (RFDEs) were presented in [Buono, 2003; Campell, 1995a; Elphick, 1987; Guckenheimer, 1983] with some restrictions on the flows on the center manifold for certain singularities, based on the well-known center manifold theory for RFDEs (see [Ma et al., 2007, 2008] and references therein). Recently, it was found that delay factors are inherent in many biological systems due to finite propagation speed of signals and finite processing time in synapses. The Rose–Hindmarsh model with time delay takes the form x = y − ax3 + bx(t − τ )2 − cz + Iapp , y = c − dx2 − y, z = r(S(x − χ) − z).
few references available so far. So, in this paper we will focus our attention on the dynamical behavior of system (1) by changing values of time delay τ , parameters S and Iapp , and fixing values for other parameters as a = 1.0,
When the time delay τ is equal to zero, Rose and Hindmarsh observed chaos in system (1) for parameter Iapp ∈ [3.28, 3.405]. When the time delay τ is not equal to zero, it seems that there are very
c = 1.0, d = 5.0, r = 0.001.
The rest of this paper is organized as follows. In Sec. 2, an analytic method is used to derive the condition under which Hopf bifurcation occurs for system (1). We apply the formula of the Sturm sequence (see [Yang et al., 1996a, 1996b; Wang et al., 1999, 2000]) to consider the stability when imaginary roots iω of the corresponding polynomial equation has high multiplicity. Hopf bifurcation curves are plotted in the (Iapp − τ ) plane by the numerical software DDE-Biftool. Double Hopf bifurcation points are obtained at the intersection points of Hopf bifurcation curves. In Sec. 3, the normal form and universal unfolding at double Hopf bifurcation points are studied using the theory of center manifold theorem. Numerical simulations of some examples are illustrated to indicate that numerical simulations agree well with theoretical analysis. Section 4 is a brief conclusion.
2. Stability and Hopf Bifurcation Analysis Suppose that (x∗ , y ∗ , z ∗ ) is the equilibrium solution of system (1), so it satisfies z ∗ = S(x∗ − χ), y ∗ = c − dx∗2 , c − dx∗2 − ax∗3 + bx∗2 − cS(x∗ − χ) + Iapp = 0.
(2)
To solve system (2), we set x∗ = ((3va + b − d)/3a) and have v 3 + 3pv + 2q = 0,
(3)
with p=− q=
(1)
b = 3.0, χ = −1.6,
cs (b − d)2 + , 2 9a 3a
(9acS(b − d) − 2(b − d)3 ) c + Iapp + cSχ . − 54a3 2a
After analyzing Eq. (3), we have four cases: (a) If p3 + q 2 > 0, then Eq. (3) has one real solution ( p3 + q 2 − q)2/3 − p ; v= ( p3 + q 2 − q)1/3
December 9, 2009
11:8
02508
Dynamics and Double Hopf Bifurcations of the Rose–Hindmarsh Model
(b) If p3 + q 2 = 0 and p = q = 0, then Eq. (3) has one real solution v = 0; (c) If p3 + q 2 = 0 and p, q = 0, then Eq. (3) has two real solutions v1 = −2q 1/3
and
v2 = q 1/3 ;
(d) If p3 + q 2 < 0, then Eq. (3) has three real solutions φ v1 = −2r cos , 3 φ+π , v2 = 2r cos 3 π−φ , v3 = 2r cos 3 3 respectively, where ∗cos φ = q/r and r = sign(q) · |p|. So x can be obtained explicitly and the associated equilibrium solution (x∗ , y ∗ , z ∗ ) can also be found accordingly.
3735
To analyze the stability of equilibrium solutions, we set x = x − x∗ , y = y − y ∗ , z = z − z ∗ , then system (1) becomes x = −3ax∗2x + 2bx∗x(t − τ ) + y − cz −ax 3 − 3ax∗x 2 + bx 2 (t − τ ), y = −2dx∗x − y − dx 2 , z = rSx − rz .
(4)
The characteristic equation of system (4) is −3ax∗2 − λ + 2bx∗ e−λτ 1 −c −1 − λ 0 = 0. −2dx∗ rS 0 −r − λ (5) Setting λ = iω in the characteristic equation (5) and separating the real part from the imaginary part, we have
sin(ωτ ) · 2bx∗ ω(1 + r) + cos(ωτ ) · 2bx∗ (r − ω 2 ) = 2dx∗ r + rSc + 3ax∗2 r − 3ax∗2 ω 2 − ω 2 − ω 2 r, (6) sin(ωτ ) · (−2b)x∗ (r − ω 2 ) + cos(ωτ ) · 2bx∗ ω(1 + r) = 2dx∗ ω + 3ax∗2 ω + ωr − ω 3 + rScω + 3ax∗2 ωr. System (6) is equivalent to sin(ωτ ) =
ω(2dx∗ ω 2 + rSc − ω 2 + 2dx∗ r 2 − r 2 − ω 2 r 2 + rScω 2 − ω 4 ) , 2bx∗ (ω 2 + r 2 + ω 4 + ω 2 r 2 )
cos(ωτ ) =
2dx∗ r 2 + r 2 Sc + r 2 Scω 2 + 3ax∗2 r 2 ω 2 + 3r 2 ax∗2 + 2dx∗ ω 2 + 3ax∗2 ω 4 + 3ax∗2 ω 2 , 2bx∗ (ω 2 + r 2 + ω 4 + ω 2 r 2 )
with G(ω 2 ) := ω 6 − a1 ω 4 − a2 ω 2 − a3 = 0,
(7)
where a1 = (2rSc + 4dx∗ − 9a2 x∗4 − r 2 + 4b2 x∗2 − 1), a2 = (−9a2 x∗4 + 4b2 x∗2 − 4d2 x∗2 + 4b2 x∗2 r 2 + 2rSc − 6r 2 Scax∗2 − 12dx∗3 a − r 2 S 2 c2 + 4dx∗ r 2 − 9a2 x∗4 r 2 − 4dx∗ rSc − r 2 ), a3 = −r 2 (Sc + 2dx∗ + 3ax∗2 − 2x∗ b) × (Sc + 2dx∗ + 3ax∗2 + 2x∗ b). The key step of the above discussion in testing the stability of system (1) is to determine the number of real roots of the polynomial as given on the left side of Eq. (7). Given a real polynomial with constant coefficients, the classical Sturm criterion gives a full answer to this problem. A complete discrimination system for polynomials
recently developed in [Yang et al., 1996a, 1996b], which is called the generalized Sturm criterion, offers a powerful tool to determine the number of a polynomial with unknown parameters. Wang and Hu [1999] presented some basic properties and provide an interpretation on the generalized Sturm criterion. In the following, we just introduce the generalized Sturm criterion, but omit the proof and basic concepts of discrimination sequence and the modified sign table (see [Yang et al. 1996a, 1996b; Wang et al., 1999, 2000]) for details). Proposition 1. Suppose that f (x) is a polynomial of order n. Let D1 (f ), D2 (f ), . . . , Dn (f ) be the corresponding discrimination sequence, and σi = Dqi (i = 1, 2, . . . , k) be the ith nonzero term of the discrimination sequence, where σ0 = 1. Let q0 ≡ 0 and ri = qi+1 − qi − 1 (i = 0, 1, 2, . . . , k − 1). Assume that the number of variation of signs in the modified sign table is s. If Dl (f ) = 0 and Dm (f ) = 0, m > l,
December 9, 2009
11:8
02508
S. Ma et al.
3736
then the following properties hold: (a) The number of distinct pairs of conjugate complex roots of f (x) is s. (b) The number of distinct real roots of f (x) is l − 2s, which satisfies k−1
l − 2s =
ri /2
(−1)
sgn
i=0, ri is even
σi+1 σi
(d) Except for some positive factors, D1 (f ), D2 (f ), . . . , Dl (f ) is the discrimination sequence of polynomial f /g.c.d(f, f ), which has no repeated roots. According to Proposition 1, we find that the discrimination sequence of G(ω 2 ) is: S := [1, a1 , a1 d1 , −d1 d2 , d2 d3 , d23 a3 ]
.
(c) α is a root of f (x) with multiplicity p if and only if it is a root of ∆n−l (f ) = 0 with multiplicity p − 1.
with d1 = (3a2 + a21 ), d2 = (a2 a21 − 3a3 a1 + 4a22 ), d3 = (4a31 a3 − a22 a21 + 18a1 a3 a2 − 4a32 + 27a23 ).
6
6
s
5
s I
4
4
IV
I (b)
I (a)
I (b)
3
3 III
III
II 2
2
1
1 II 0 −2
5
II
IV 0
−1.5
−1
−0.5
0
0.5 *
1
x
−1.6
−1.4
−1.2
(a)
−1
−0.8
−0.6
x*
(b)
1.5
6 5.5
s
s 5
I
4.5
1
4 3.5
II (b)
II (b)
II (a)
3
0.5
2.5 2
II (a)
1.5
0
−1.3
−1.2
−1.1
−1
−0.9
−0.8
−0.7
x (c)
*
−0.6
1
0.05
0.1
0.15
0.2
0.25
0.3
0.35
*
0.4
x (d)
Fig. 1. Classification diagram of the number of real roots of G(ω 2 ) in the (x∗ − S) plane. G(ω 2 ) is defined by Eq. (7). (a) The whole classification diagram, (b)–(d) the amplification cases.
December 9, 2009
11:8
02508
Dynamics and Double Hopf Bifurcations of the Rose–Hindmarsh Model Table 1.
Imaginary roots of the characteristic equation (5).
S
x∗
Region
Imaginary Roots ±iω
Modified Sign Table
3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 0.5 0.5 0.5
−1.8 −1.4 −1.2 −0.8 −0.5 0.0 0.1 0.2 −1.1 −0.8 −0.6
III I(b) I(a) I(b) III IV II(b) II(a) II(a) II(b) IV
±1.655996i ±2.75883i ±0.002059i, ±0.990345i, ±2.683806i ±0.000988i ±0.0009077i none none ±0.9251473i, ±1.89095i ±1.41178i, ±2.321846i none none
[1, −1, −1, −1, −1, 1] [1, 1, 1, −1, −1, 1] [1, 1, 1, 1, 1, 1] [1, 1, −1, 1, 1, 1] [1, −1, 1, 1, 1, 1] [1, −1, −1, −1, 1, −1] [1, 1, −1, 1, 1, −1] [1, 1, 1, 1, 1, −1] [1, 1, 1, 1, 1, −1] [1, 1, −1, 1, 1, −1] [1, −1, 1, 1, 1, −1]
We classify the number of real roots of G(ω 2 ) by the following four cases: Case I . Suppose a1 > 0, a3 > 0, then (a) if d1 > 0, d2 < 0, d3 < 0, G(ω 2 ) has three real roots; (b) otherwise, G(ω 2 ) has one real root. Case II . Suppose a1 > 0, a3 < 0, then (a) if d1 > 0, d2 < 0 or d1 > 0, d2 > 0, d3 < 0 or d1 < 0, d2 < 0, d3 > 0, G(ω 2 ) has two real roots; (b) if d2 > 0, d3 > 0 or d1 < 0, d3 < 0, G(ω 2 ) has no real roots. Case III . If a1 < 0, a3 > 0, G(ω 2 ) has one real root. Case IV . If a1 < 0, a3 < 0, G(ω 2 ) has two real roots provided that d1 > 0, d2 > 0, d3 < 0; otherwise, it has no real roots.
Set a = 1.0, b = 3.0, c = 1.0 and d = 5.0. The classification diagram of the number of real roots of G(ω 2 ) is plotted in the (x∗ − S) plane as shown in Fig. 1(a) and amplifying diagrams as shown in Figs. 1(b)–1(d). For example, transit (x∗ , S) between different regimes from left to right, imaginary roots of the characteristic equation (5) are listed in Table 1. The linearized system of system (4) is stable for any value of τ if there exists no imaginary root. Therefore, from Table 1 we obtain Theorem 1. The linearized system of (4) is abso-
lutely stable for any value of time delay τ if and only if parameter S and the value of x∗ lie in regions II(b) and IV . By applying the DDE-Biftool, we plot Hopf bifurcation curves in the (x∗ − τ ) plane while taking S = 1. As shown in Fig. 2, two pairs of imaginary
15
τ
14
3737
6 DH
τ
13 12
5 DH
4
11
DH
3
10 9
2
DH
8
6 0.1
1
DH
7 0.15
Fig. 2.
0.2
x*
0.25
0 1.7
1.75
1.8
Hopf bifurcation curves of the linearized system of (4) when S = 1.0.
1.85
x*
1.9
December 9, 2009
3738
11:8
02508
S. Ma et al. Table 2.
Values of double Hopf bifurcation points.
x∗
τ
Imaginary Roots
x∗
τ
Imaginary Roots
0.1768 1.7523
6.8693 1.7741
ω1 = 0.9107; ω2 = 1.6810 ω1 = 3.5918; ω2 = 6.8435
1.8207 1.8383
2.9934 4.2114
ω1 = 4.1953; ω2 = 6.1936 ω1 = 4.4642; ω2 = 5.9052
roots are acquired simultaneously at the intersection point of two Hopf bifurcation curves. Table 2 lists the values of imaginary roots and the related Double Hopf bifurcation points.
+ 2bˆ x∗ [x(t − τˆ − εˆ τε ) − x(t − τˆ)] x∗ x + 2εˆ xε bx(t − τˆ) − aεx3 + o(ε), − 6εˆ xε aˆ √ x∗ dx − y − 2εˆ xε dx − εdx2 , y = −2ˆ z = rSx − rz. (8)
3. Normal Form for Double Hopf Bifurcations
√ Set x∗ = √ x ˆ∗ + εˆ xε , τ = τˆ + εˆ τε and x = εx, y = √ εy, z = εz, then system (4) is equivalent to x∗2 x + 2bˆ x∗ x(t − τˆ) + y − cz x = −3aˆ √ √ x∗ x2 + εbx(t − τˆ)2 − 3 εaˆ
Choose the phase space C = C([−ˆ τ , 0), R3 ) as the Banach space of continuous functions from [−ˆ τ , 0) to R3 with the supremum norm. For φ ∈ C, we define 0 [dη(θ)]φ(θ), (9) L(0)φ = −ˆ τ
where η : [−ˆ τ , 0] → R3 ×R3 is a real-valued function of bounded variation in the interval [−ˆ τ , 0] with x∗ δ(θ + τˆ) δ(θ) −cδ(θ) −3aˆ x∗2 δ(θ) + 2bˆ −δ(θ) 0 dθ. −2ˆ x∗ dδ(θ) dη(θ) = rSδ(θ) 0 −rδ(θ)
Then, we let L(ˆ xε , τˆε )φ =
0
−ˆ τ
and
where
−6aˆ x∗ x ˆε δ(θ) + 2bˆ xε δ(θ + τˆ) −2ˆ xε dδ(θ) dη1 (θ) = 0
and
[εdη1 (θ, x ˆε ) + dη2 (θ, ε, τˆε )]φ(θ), 0 0 0 0 dθ, 0 0
2bˆ x∗ (δ(θ + τˆ + εˆ τε ) − δ(θ + τˆ)) 0 0 0 0 0 dθ. dη2 (θ) = 0 0 0 For φ ∈ C, the linear operator defined by system (9) generates an infinitesimal generator of the semi-flow of bounded linear operators with dφ , θ ∈ [−ˆ τ , 0), (10) A(0)φ = dθ L(0)φ, θ = 0. Hence, we define
0, A(ˆ xε , τˆε )φ = L(ˆ xε , τˆε )φ,
θ ∈ [−τ, 0), θ = 0,
Q(φ) = where
0, F (φ),
θ ∈ [−τ, 0), θ = 0,
√ √ x∗ φ21 (0) + εbφ21 (−ˆ τ ) − εaφ31 (0) −3 εaˆ √ − εdφ21 (0) F (φ) = . 0
System (8) can be re-expressed as an operator differential equation xε , τˆε )ut + εQ(φ), u (t) = A(0)ut + A(ˆ
(11)
τ ≤ where u = (x, y, z)T and ut = u(t + θ) for −ˆ θ < 0. For ψ ∈ C ∗ = C([0, τˆ], R2 ), the adjoint operators A∗ (0), A∗ (σ1 , σ2 ) of A(0) and A(σ1 , σ2 ), respectively, are defined as dψ − 0 < s ≤ τˆ, ds , A∗ (0)ψ(s) = 0 dη T (s)ψ(−s), s = 0, −ˆ τ
December 9, 2009
11:8
02508
Dynamics and Double Hopf Bifurcations of the Rose–Hindmarsh Model
and A∗ (σ1 , σ2 )ψ(s) 0, = 0 dη T (s, σ1 , σ2 )ψ(−s),
0 < s ≤ τˆ, s = 0.
−ˆ τ
For ϕ ∈ C and ψ ∈ C ∗ , the bilinear form is given as 0 θ ψ T (ξ − θ)dη(θ)ϕ(ξ)dξ. ψ, ϕ = ψ T (0)ϕ(0) − −ˆ τ
0
From the discussion shown in the preceding section, we know that the characteristic equation (5) has two pairs of pure imaginary characteristic values Λ = {iω1 , iω2 } at the double Hopf bifurcation point (ˆ x∗ , τˆ), and other characteristic values with
1
2x∗ d − iω θ 1 , q1 (θ) = 1 + iω1 e rS r + iω1 and
r − iω1 − c 2 c(1 − iω1 ) (r − iω1 ) r − iω1 eiω1 s , p1 (s) = l1 + im1 c(−1 + iω1 ) 1 r − iω2 − c 2 iω s c(1 − iω2 ) (r − iω2 ) r − iω e 2 , 2 p2 (s) = l2 + im2 c(−1 + iω2 ) 1
where li and mi (i = 1, 2) are given as follows: l1 = (−3τ ax∗2 − 2τ − 2rτ − 1)ω14
negative real parts. Hence, C can be decomposed as C = PΛ ⊕ QΛ , where PΛ is the characteristic space corresponding to eigenvalues Λ, and QΛ is the complementary subspace of PΛ . System (11) is an ordinary differential system on the Banach space C of functions from [−ˆ τ , 0] to R3 , which is bounded and continuous on (−ˆ τ , 0) with a possible jump discontinuity at 0 [Buono et al., 2003]. We now choose Φ(θ) = (q1 (θ),q 1 (θ), q2 (θ),q 2 (θ)), (12) Ψ(θ) = (p1 (s), p 1 (s), p2 (s), p 2 (s)), as the bases of PΛ and PΛ∗ , where p1 , q1 = 1, p1 , q 1 = 0, p1 , q2 = 0, p1 , q 2 = 0, p2 , q1 = 0, p2 ,q 1 = 0, p2 , q2 = 1, p2 , q 2 = 0. Then, we may choose 1 2x∗ d − iω θ 2 , q2 (θ) = 1 + iω2 e rS r + iω2 + (4rτ x∗ d + 6r 2 τ ax∗2 + 2r 2 τ Sc + r 2 τ + 2r 2 + 2r 2 τ x∗ d − 4x∗ dr + 6rτ ax∗2 + rτ Sc + 2r − 2rSc)ω1 , l2 = (−3τ ax∗2 − 2τ − 2rτ − 1)ω24 + (2r 2 τ + r 2 τ Sc + r 2 + 3τ ax∗2 + 12rτ ax∗2 + 2rτ Sc − rSc + 4r + 4rτ x∗ d + 2rτ + 3r 2 τ ax∗2 + 2τ x∗ d − 2dx∗ + 1)ω22 − r 2 − 2r 2 τ x∗ d − 3r 2 τ ax∗2 + rSc − r 2 τ Sc + 2x∗ dr 2 , m2 = ω25 τ + (−rτ Sc − 6rτ ax∗2 − r 2 τ − 4rτ − 2r − τ − 2 − 6τ ax∗2 − 2τ x∗ d)ω23
+ (2r 2 τ + r 2 τ Sc + r 2 + 3τ ax∗2 + 12rτ ax∗2
+ (4rτ x∗ d + 6r 2 τ ax∗2 + 2r 2 τ Sc
+ 2rτ Sc − rSc + 4r + 4rτ x∗ d + 2rτ
+ r 2 τ + 2r 2 + 2r 2 τ x∗ d − 4x∗ dr
+ 3r 2 τ ax∗2 + 2τ x∗ d − 2dx∗ + 1)ω12
+ 6rτ ax∗2 + rτ Sc + 2r − 2rSc)ω2 .
− r 2 − 2r 2 τ x∗ d − 3r 2 τ ax∗2 + rSc − r 2 τ Sc + 2x∗ dr 2 , m1 = ω15 τ + (−rτ Sc − 6rτ ax∗2 − r 2 τ − 4rτ − 2r − τ − 2 − 6τ ax∗2 − 2τ x∗ d)ω13
3739
Thus, explicit expressions for Φ and Ψ can be obtained by substituting pi and qi (i = 1, 2) into system (12). Let z = (z1 ,z 1 , z2 ,z 2 )T and vt ∈ QΛ ∩ τ , 0), R), and denote ut = Φz + vt . System C ([−ˆ
December 9, 2009
3740
11:8
02508
S. Ma et al.
(11) is equivalent to z = Bz + Ψ T (0)(A(xε , τε )(Φz + vt ) + εF (Φz + vt )), vt = A(0)vt + (I − Π)X0 (A(pε , τε ) × (Φz + vt ) + εF (Φz + vt )), where
iω1 0 B= 0 0
0 −iω1 0 0
0 0 iω2 0
The matrix B generates the torus group T2 = S1 × S1 whose action on C2 is given by
Then, the T2 -equivalent normal form, which is truncated to the quadratic order, becomes
0 0 . 0 −iω2
ˆε z1 + εb12 τˆε z1 , z1 = iω1 z1 + εb11 x z2 = iω2 z2 + εb21 x ˆε z2 + εb22 τˆε z2 ,
b11 = −
2(−bω1 e−iω1 τˆ − ibe−iω1 τˆ + id + 3aω1 x ˆ∗ + 3iaˆ x∗ )(r − iω1 )2 (i + ω1 ) , l1 + im1
b12 = −
−2ibˆ x∗ e−iω1 τˆ ω1 (r − iω1 )2 (i + ω1 )2 , l1 + im1
b21 = −
2(−bω2 e−iω2 τˆ − ibe−iω2 τˆ + id + 3aω2 x ˆ∗ + 3iaˆ x∗ )(r − iω2 )2 (i + ω2 ) , l2 + im2
b22 = −
−2ibˆ x∗ e−iω2 τˆ ω2 (r − iω2 )2 (i + ω2 )2 . l2 + im2
xε + Re(b12 )ˆ τε )r1 , r1 = (Re(b11 )ˆ xε + Re(b22 )ˆ τε )r2 . r2 = (Re(b21 )ˆ
(15)
If we remove the dependence on the unfolding parameters x ˆε and τˆε in system (13), we obtain =
iω1 z1 + p T1 (0)F (Φz
=
iω1 z1 + p T1 (0)fˆ(z1 ,z 1 , z2 , z 2 ), iω2 z2 + p T2 (0)F (Φz + vt )
=
+ vt ) (16)
= iω2 z2 + p T2 (0)fˆ(z1 ,z 1 , z2 , z 2 ), = Avt + (I − π)X0 (F (Φz + vt )),
where fˆ(z1 ,z 1 , z2 , z 2 ) = F (2 Re{z1 q1 (θ)} + 2 Re{z2 q2 (θ)} + vt ). We now expand vt = w(z, z z1 ,z 1 , z2 ,z 2 ) into the series wi,j,k,l z1i z j1 z2kz l2 , vt = w(z1 ,z 1 , z2 ,z 2 ) = i+j+k+l≥2
(14)
where
In polar coordinates z1 = r1 eiρ1 and z2 = r2 eiρ2 , the associated amplitude equations generated from system (14) are
z1 z2 vt
(θ1 , θ2 )(z1 , z2 ) = (eiθ1 z1 , eiθ2 z2 ).
(13)
where the coefficient vectors {wi,j,k,l } only depend 1 2 3 , wi,j,k,l , wi,j,k,l )T . Simion θ and wi,j,k,l = (wi,j,k,l larly, we are able to re-express the first two equations in system (16) as z1 = iω1 z1 + g1 (z1 ,z 1 , z2 ,z 2 ), (17) z2 = iω2 z2 + g2 (z1 ,z 1 , z2 ,z 2 ), with the expansions g1 (z1 , z 1 , z2 ,z 2 ) = g2 (z1 , z 1 , z2 ,z 2 ) =
i+j+k+l≥2
i+j+k+l≥2
1 gi,j,k,l z1i z j1 z2kz l2 , 2 gi,j,k,l z1i z j1 z2kz l2 ,
where all coefficients {gi,j,k,l } are real constants. From system (16), we can derive Avt − 2 Re{p T1 (0)fˆ(z1 , z 1 , z2 ,z 2 )q1 (θ)} − 2 Re{p T2 (0)fˆ(z1 ,z 1 , z2 ,z 2 )q2 (θ)}, (−τ ≤ θ < 0), dvt = dt Avt − 2 Re{p T1 (0)fˆ(z1 , z 1 , z2 ,z 2 )q1 (0)} − 2 Re{p T2 (0)fˆ(z1 ,z 1 , z2 ,z 2 )q2 (0)} + fˆ(z1 , z 1 , z2 ,z 2 ), (θ = 0). (18)
December 9, 2009
11:8
02508
Dynamics and Double Hopf Bifurcations of the Rose–Hindmarsh Model
The first equation of (18) can be rewritten as dw (z1 ,z 1 , z2 , z 2 ) dt = Aw(z1 ,z 1 , z2 , z 2 ) + H(z1 ,z 1 , z2 , z 2 ),
−τ ≤ θ < 0,
(19)
where H(z1 ,z 1 , z2 ,z 2 ) = −2 Re{p T1 (0)fˆ(z1 , z 1 , z2 ,z 2 )q1 (θ)} − 2 Re{p T2 (0)fˆ(z1 ,z 1 , z2 ,z 2 )q2 (θ)} Hi,j,k,lz1i z j1 z2kz l2 , = i+j+k+l≥2
k g2000
k g1100
k g0200
k g1010
k g1001
k g0110
k g0101
k g0020
k g0011
k g0002
3741
1 2 with the coefficients Hi,j,k,l = (Hi,j,k,l (θ), Hi,j,k,l (θ), 3 T Hi,j,k,l(θ)) . In order to determine the behavior of bifurcating solution, it is necessary to compute vt = w(z1 ,z 1 , z2 , z 2 ) through system (18). Then {gij } and {Hij } should be found first. Let us define that
φ(0) = z1 q1 (0) + z 1q 1 (0) + z2 q2 (0) + z 2q 2 (0) + w(z1 , z 1 , z2 , z 2 )|θ=0 , φ(−τ ) = z1 q1 (−τ ) + z 1q 1 (−τ ) + z2 q2 (−τ ) + z 2q 2 (−τ ) + w(z1 ,z 1 , z2 ,z 2 )|θ=−τ . From systems (16) and (17), by direct calculations for k = 1, 2, we have
√
ε(r + iωk )2 (−i + ωk ) (−bωk e−2iω1 τ + ibe−2iω1 τ − id + 3ax∗ ωk − 3iax∗ ), lk − imk √ 2 ε(r + iωk )2 (−i + ωk ) =− (−bωk − id + ib + 3ax∗ ωk − 3iax∗ ), lk − imk √ i ε(r + iωk )2 (−i + ωk ) = (−ibωk e2iω1 τ − be2iω1 τ + d + 3iax∗ ωk + 3ax∗ ), lk − imk √ 2i ε(r + iωk )2 (−i + ωk ) = (−ibωk e−iτ (ω1 +ω2 ) − be−iτ (ω1 +ω2 ) + d + 3iaωk x∗ + 3ax∗ ), lk − imk √ 2i ε(r + iωk )2 (−i + ωk ) =− (−ibω1 e−iτ (ω1 −ω2 ) − be−iτ (ω1 −ω2 ) + d + 3iaω1 x∗ + 3ax∗ ), lk − imk √ 2i ε(r + iωk )2 (−i + ωk ) = (−ibωk eiτ (ω1 −ω2 ) − beiτ (ω1 −ω2 ) + d + 3iaωk x∗ + 3ax∗ ), lk − imk √ 2 ε(r + iωk )2 (−i + ωk ) =− (−bωk eiτ (ω1 +ω2 ) + ibeiτ (ω1 +ω2 ) − id + 3aωk x∗ − 3iax∗ ), lk − imk √ − ε(r + iωk )2 (−i + ωk ) = (−bωk e−2iω2 τ + ibe−2iω2 τ − id + 3ax∗ ωk − 3iax∗ ), lk − imk √ 2 ε(r + iωk )2 (−i + ωk ) =− (−bωk + ib − id + 3aωk x∗ − 3iax∗ ), lk − imk √ i ε(r + iωk )2 (−i + ωk ) = (−ibωk e2iω2 τ − be2iω2 τ + d + 3iaωk x∗ + 3ax∗ ). lk − imk =−
Thus, for −τ ≤ θ < 0, we have H(z1 , z 1 , z2 ,z 2 ) = −2 Re{p T1 (0)fˆ(z1 ,z 1 , z2 , z 2 )q1 (θ)} − 2 Re{p T2 (0)fˆ(z1 ,z 1 , z2 , z 2 )q2 (θ)}, = −g1 (z1 ,z 1 , z2 , z 2 )q1 (θ) − g 1 (z1 , z 1 , z2 , z 2 )q 1 (θ) − g2 (z1 , z 1 , z2 ,z 2 )q2 (θ) − g 2 (z1 , z 1 , z2 ,z 2 )q 2 (θ).
December 9, 2009
3742
11:8
02508
S. Ma et al.
Comparing this formula with the expansion of H(z1 ,z 1 , z2 , z 2 ), we have 1 2 q1 (θ) − g 10200q 1 (θ) − g2000 q2 (θ) − g 20200q 2 (θ), H2000 (θ) = −g2000 1 2 q1 (θ) − g 11100q 1 (θ) − g1100 q2 (θ) − g 21100q 2 (θ), H1100 (θ) = −g1100 1 2 q1 (θ) − g 12000q 1 (θ) − g0200 q2 (θ) − g 22000q 2 (θ), H0200 (θ) = −g0200 1 2 q1 (θ) − g 10101q 1 (θ) − g1010 q2 (θ) − g 20101q 2 (θ), H1010 (θ) = −g1010 1 2 q1 (θ) − g 10110q 1 (θ) − g1001 q2 (θ) − g 20110q 2 (θ), H1001 (θ) = −g1001 1 2 q1 (θ) − g 11001q 1 (θ) − g0110 q2 (θ) − g 21001q 2 (θ), H0110 (θ) = −g0110 1 2 q1 (θ) − g 11010q 1 (θ) − g0101 q2 (θ) − g 21010q 2 (θ), H0101 (θ) = −g0101 1 2 q1 (θ) − g 10002q 1 (θ) − g0020 q2 (θ) − g 20002q 2 (θ), H0020 (θ) = −g0020 1 2 q1 (θ) − g 10011q 1 (θ) − g0011 q2 (θ) − g 20011q 2 (θ), H0011 (θ) = −g0011 1 2 q1 (θ) − g 10020q 1 (θ) − g0002 q2 (θ) − g 20020q 2 (θ). H0002 (θ) = −g0002
It can be derived from Eq. (19) that (A − 2iω1 )w2000 = −H2000 (θ), Aw1100 = −H1100 (θ), (A + 2iω1 )w0200 = −H0200 (θ), (A − iω1 − iω2 )w1010 = −H1010 (θ), (A − iω1 + iω2 )w1001 = −H1001 (θ), (A + iω1 − iω2 )w0110 = −H0110 (θ), (A + iω1 + iω2 )w0101 = −H0101 (θ), (A − 2iω2 )w0020 = −H0020 (θ), Aw0011 = −H0011 (θ), (A + 2iω2 )w0002 = −H0002 (θ). Taking account of formula (10) of operator A, we have dw2000 = 2iω1 w2000 (θ) − H2000 (θ), dθ dw1100 = −H1100 (θ), dθ dw1010 = i(ω1 + ω2 )w1010 (θ) − H1010 (θ), dθ dw1001 = i(ω1 − ω2 )w1001 (θ) − H1001 (θ), dθ dw0110 = i(−ω1 + ω2 )w0110 (θ) − H0110 (θ), dθ dw0020 = 2iω2 w0020 (θ) − H0020 (θ), dθ dw0011 = −H0011 (θ). dθ
(20)
December 9, 2009
11:8
02508
Dynamics and Double Hopf Bifurcations of the Rose–Hindmarsh Model
3743
Similarly, from the second equation of system (18), it can be derived that 0 [dη(θ)]w2000 −ˆ τ 0 [dη(θ)]w1100 −ˆ τ 0 [dη(θ)]w1010 −ˆ τ 0 [dη(θ)]w1001 −ˆ τ 0 [dη(θ)]w0110 −ˆ τ 0 [dη(θ)]w0020 −ˆ τ 0 [dη(θ)]w0011 −ˆτ
be−2iω1 τ − 3ax∗ √ −d = 2iω1 w2000 (0) − H2000 (0) − ε , 0 2b − 6ax∗ √ = −H1100 (0) − ε −2d , 0 −6ax∗ + 2be−i(ω1 +ω2 )τ √ = i(ω1 + ω2 )w1010 (0) − H1010 (0) − ε −2d , 0 −6ax∗ + 2be−i(ω1 −ω2 )τ √ = i(ω1 − ω2 )w1001 (0) − H1001 (0) − ε −2d , 0 −6ax∗ + 2bei(ω1 −ω2 )τ √ = i(−ω1 + ω2 )w0110 (0) − H0110 (0) − ε −2d , 0 −2iω τ 2 − 3ax∗ be √ −d = 2iω2 w0020 (0) − H0020 (0) − ε , 0 −6ax∗ + 2b √ −2d = −H0011 (0) − ε . 0
Integrating system (20) under conditions (21), we obtain 1 (θ) = w2000
1 (6iω12g 20200 e−iω2 θ − 3iω1 ω2g 20200 e−iω2 θ 3ω1 (2ω1 + ω2 )(2ω1 − ω2 ) 1 1 + 12iω12 g2000 eiω1 θ − 3iω22 g2000 eiω1 θ + 4iω12g 10200 e−iω1 θ − iω22g 10200 e−iω1 θ 2 2 + 6iω12 g2000 eiω2 θ + 3iω1 ω2 g2000 eiω2 θ + 12E1 ω13 e2iω1 θ − 3E1 ω1 ω22 e2iω1 θ ),
1 (θ) = w1100
1 1 2 (−iω2 g1100 eiω1 θ + ig 11100 ω2 e−iω1 θ − ig1100 ω1 eiω2 θ + ig 21100 ω1 e−iω2 θ ) + E2 , ω1 ω2
1 (θ) = w1010
1 (2iω12 ω2g 20101 e−iω2 θ + iω1 ω22g 20101 e−iω2 θ ω1 ω2 (ω1 + 2ω2 )(2ω1 + ω2 ) 1 1 1 + 2iω13 eiω1 θ g1010 + 5iω12 ω2 g1010 eiω1 θ + 2iω1 ω22 g1010 eiω1 θ + iω2 ω12g 10101 e−iω1 θ 2 2 2 + 2iω22 ω1g 10101 e−iω1 θ + 2iω2 ω12 g1010 eiω2 θ + 5iω22 ω1 g1010 eiω2 θ + 2iω23 g1010 eiω2 θ
+ 2E3 ω2 ω13 ei(ω1 +ω2 )θ + 5E3 ω22 ω12 ei(ω1 +ω2 )θ + 2E3 ω23 ω1 ei(ω1 +ω2 )θ ),
(21)
December 9, 2009
3744
11:8
02508
S. Ma et al. 1 w1001 (θ) =
i (2ω2 ω12g 20110 e−iω2 θ − 5ω22 ω1g 20110 e−iω2 θ ω1 ω2 (2ω1 − ω2 )(ω1 − 2ω2 ) 1 1 1 + 2ω23g 20110 e−iω2 θ − 2ω13 g1001 eiω1 θ + 5ω12 ω2 g1001 eiω1 θ − 2ω1 ω22 g1001 eiω1 θ 2 2 + ω12 ω2g 10110 e−iω1 θ − 2ω1 ω22g 10110 e−iω1 θ + 2ω12 ω2 g1001 eiω2 θ − ω1 ω22 g1001 eiω2 θ
− 2iE4 ω13 ω2 ei(ω1 −ω2 )θ + 5iE4 ω12 ω22 ei(ω1 −ω2 )θ − 2iE4 ω1 ω23 ei(ω1 −ω2 )θ ), 1 w0110 (θ) =
1 (−2iω2 ω12g 21001 e−iω2 θ + iω22 ω1g 21001 e−iω2 θ ω1 ω2 (ω1 − 2ω2 )(2ω1 − ω2 ) 1 1 − iω2 ω12 g0110 eiω1 θ + 2iω22 ω1 g0110 eiω1 θ + 2iω13g 11001 e−iω1 θ − 5iω12 ω2g 11001 e−iω1 θ 2 2 2 + 2iω1 ω22g 11001 e−iω1 θ − 2iω2 ω12 g0110 eiω2 θ + 5iω22 ω1 g0110 eiω2 θ − 2iω23 g0110 eiω2 θ
+ 2E5 ω2 ω13 e−i(ω1 −ω2 )θ − 5E5 ω22 ω12 e−i(ω1 −ω2 )θ + 2E5 ω23 ω1 e−i(ω1 −ω2 )θ ), 1 (θ) = w0020
i (−4ω22g 20002 e−iω2 θ + ω12g 20002 e−iω2 θ 3ω2 (−2ω2 + ω1 )(2ω2 + ω1 ) 1 1 − 6ω22 g0020 eiω1 θ − 3ω2 ω1 g0020 eiω1 θ − 6ω22g 10002 e−iω1 θ + 3ω2 ω1g 10002 e−iω1 θ 2 2 − 12ω22 g0020 eiω2 θ + 3ω12 g0020 eiω2 θ + 12iE6 ω23 e2iω2 θ − 3iE6 ω2 ω12 e2iω2 θ ),
1 (θ) = w0011
1 1 2 (−iω2 g0011 eiω1 θ + ig 10011 ω2 e−iω1 θ − ig0011 ω1 eiω2 θ + ig 20011 ω1 e−iω2 θ ) + E7 , ω1 ω2 (22)
where Ei (i = 1, . . . , 7) are determined by the following equations, crS 2x∗ d 1 1 − (0) + 2bx∗ w2000 (−τ ) w2000 −2iω1 − 3ax∗2 − r + 2iω1 1 + 2iω1 =−
2 3 √ √ d cH2000 (0) √ (0) H2000 1 + ε + − H2000 (0) − εbe−2iω1 τ + 3 εax∗ , 1 + 2iω1 1 + 2iω1 r + 2iω1
1 1 (0) + 2bx∗ w1100 (−τ ) (−3ax∗2 − 2x∗ d − cS)w1100 3 √ √ √ cH1100 (0) 1 2 − H1100 (0) + 2 εd + (0) − 2 εb + 6 εax∗ , = −H1100 r crS 2x∗ d ∗2 1 1 − (0) + 2bx∗ w1010 (−τ ) w1010 −iω1 − iω2 − 3ax − r + iω1 + iω2 1 + iω1 + iω2
√ 2 3 √ √ 2d ε cH1010 (0) (0) H1010 1 + + − H1010 (0) − 2 εbe−i(ω1 +ω2 )τ + 6 εax∗ , =− 1 + iω1 + iω2 1 + iω1 + iω2 r + iω1 + iω2 −iω1 + iω2 − 3ax∗2 −
crS 2x∗ d − r + iω1 − iω2 1 + iω1 − iω2
1 1 (0) + 2bx∗ w1001 (−τ ) w1001
√ 2 3 √ √ 2d ε cH1001 (0) (0) H1001 1 + + − H1001 (0) − 2 εbe−i(ω1 −ω2 )τ + 6 εax∗ , =− 1 + iω1 − iω2 1 + iω1 − iω2 r + iω1 − iω2
December 9, 2009
11:8
02508
Dynamics and Double Hopf Bifurcations of the Rose–Hindmarsh Model
3745
crS 2x∗ d 1 1 − − (0) + 2bx∗ w0110 (−τ ) iω1 − iω2 − w0110 r − iω1 + iω2 1 − iω1 + iω2 √ 2 3 √ √ 2d ε cH0110 (0) (0) H0110 1 + + − H0110 (0) − 2 εbei(ω1 −ω2 )τ + 6 εax∗ , =− 1 − iω1 + iω2 1 − iω1 + iω2 r − iω1 + iω2 2x∗ d crS ∗2 1 1 − − 2iω2 w0020 (0) + 2bx∗ w0020 (−τ ) −3ax − 1 + 2iω2 r + 2iω2 3ax∗2
=−
2 3 √ √ d cH0020 (0) √ (0) H0020 1 + ε + − H0020 (0) + 3 εax∗ − εbe−2iω2 τ , 1 + 2iω2 1 + 2iω2 r + 2iω2
1 1 (0) + 2bx∗ w0011 (−τ ) (−3ax∗2 − 2x∗ d − cS)w0011 3 √ √ √ cH0011 (0) 1 2 − H0011 (0) + 2 εd + (0) + 6 εax∗ − 2 εb. = −H0011 r
From (22), we deduce 1 = g2100
√ √ (r + iω1 )2 (−i + ω1 ) √ 1 1 1 (2 εbω1 e−iω1 τ w1100 (−τ ) − 2i εbe−iω1 τ w1100 (−τ ) + 2 εbeiω1 τ ω1 w2000 (−τ ) l1 − im1 √ √ √ √ √ 1 1 1 1 1 (−τ ) + 2i εdw2000 (0) + 2i εdw1100 (0) − 6 εax∗ ω1 w2000 (0) + 6i εax∗ w2000 (0) − 2i εbeiω1 τ bw2000 √ √ 1 1 (0) + 6i εax∗ w1100 (0) − 3εaω1 + 3iεa), − 6 εax∗ ω1 w1100 √ √ 2(r + iω1 )2 (−i + ω1 ) √ (− εbω1 eiω2 τ w1010 (−τ ) − εbω1 e−iω1 τ w0011 (−τ ) − εbω1 e−iω2 τ w1001 (−τ ) l1 − im1 √ √ √ √ − i εdw1001 (0) − 3ia εx∗ w1010 (0) + 3εaω1 + i εbe−iω1 τ w0011 (−τ ) + i εbeiω2 τ w1010 (−τ ) √ √ √ √ − 3ia εx∗ w1001 (0) − i εdw1010 (0) + 3 εaω1 x∗ w0011 (0) + 3 εaω1 x∗ w1010 (0) − 3iεa √ √ √ √ + 3 εaω1 x∗ w1001 (0) − 3i εax∗ w0011 (0) + i εbe−iω2 τ w1001 (−τ ) − i εdw0011 (0)),
1 =− g1011
√ √ 2(r + iω2 )2 (−i + ω2 ) √ (− εbω2 eiω1 τ w1010 (−τ ) − 3i εax∗ w1100 (0) − εbω2 e−iω1 τ w0110 (−τ ) l2 − im2 √ √ √ √ −iω + i εbe 2 τ w1100 (−τ ) − εbω2 e−iω2 τ w1100 (−τ ) + i εbe−iω1 τ w0110 (−τ ) + i εbeiω1 τ w1010 (−τ ) √ √ √ √ √ − i εdw0110 (0) − 3i εax∗ w0110 (0) + 3 εaω2 x∗ w1010 (0) + 3 εaω2 x∗ w0110 (0) + 3 εaω2 x∗ w1100 (0) √ √ √ − 3i εax∗ w1010 (0) − i εdw1010 (0) − i εdw1100 (0) + 3εaω2 − 3iεa),
2 =− g1110
√ √ √ (r + iω2 )2 (−i + ω2 ) (−2 εbω2 e−iω2 τ w0011 (−τ ) + 2i εbe−iω2 τ w0011 (−τ ) − 2 εbω2 eiω2 τ w0020 (−τ ) l2 − im2 √ √ √ √ √ iω2 τ w0020 (−τ ) − 2i εdw0011 (0) − 2i εdw0020 (0) + 6 εaω2 x∗ w0011 (0) − 6i εax∗ w0011 (0) + 2i εbe √ √ + 6 εaω2 x∗ w0020 (0) − 6i εax∗ w0020 (0) + 3εaω2 − 3iεa).
2 =− g0021
Using the standard normalization technique as described in [Ma et al., 2008a, 2008b], we set z1 = η1 + P˜1 (η1 ,η 1 , η2 ,η 2 ), η1 = iω1 η1 + B11 η12η 1 + B12 η1 η2η 2 ,
z2 = η2 + P˜2 (η1 ,η 1 , η2 ,η 2 ), η2 = iω2 η2 + B21 η1η 1 η2 + B22 η22η 2 ,
(23)
December 9, 2009
3746
11:8
02508
S. Ma et al.
where
P˜1 (η1 ,η 1 , η2 ,η 2 ) =
j+k+l+m≥2
P˜2 (η1 ,η 1 , η2 ,η 2 ) =
j+k+l+m≥2
1 η1jη k1 η2l η m P˜jklm 2 , 2 η1jη k1 η2l η m P˜jklm 2 .
Substituting (23) into system (17), then comparing the corresponding coefficients of both sides to determine the appropriate functions P˜1 and P˜2 , so as to round off the nonresonant terms, we obtain 1 + B11 = g2100
−
i 1 i 1 i 1 2 1 g1100 g2000 + (g1010 g1100 − g1001 g 21100 ) − g1 g 2 ω1 ω2 2ω1 + ω2 0101 0200
i i 1 2 2i 1 2 1 2 g0110 g2000 − |g1100 | − |g | , 2ω1 − ω2 ω1 3ω1 0200
1 + B12 = g1011
i 1 i 2 1 1 1 1 1 (g1010 g0011 − g1001 g 20011 ) + (2g2000 g0011 − g1100 g 10011 − g0011 g 20110 ω2 ω1
1 2 −g0011 g1010 )− 2 + B21 = g1110
i 1 i 2 2 2 2 1 2 (g g2 − g 11100 g0110 ) + (2g0020 g1100 − g0011 g 21100 − g1010 g1100 ω1 1100 1010 ω2
2 −g 11001 g1100 )+ 2 + B22 = g0021
−
2i 2i i i g1 g 2 − g1 g2 − |g1 |2 − |g1 |2 , ω1 + 2ω2 0002 0101 ω1 − 2ω2 0020 1001 2ω1 − ω2 0110 2ω1 + ω2 0101
2i 2i i i 1 2 2 2 g0110 g2000 − g 10101 g0200 − |g1001 |2 − |g2 |2 , 2ω1 − ω2 2ω1 + ω2 2ω2 − ω1 ω1 + 2ω2 0101
i 1 i 2 i 2 2 2 (g0011 g1010 − g 10011 g0110 ) + g0011 g0020 − g1 g2 ω1 ω2 2ω2 − ω1 0020 1001
i i 2 2 2i 2 2 2 g 10002 g0101 − |g0011 | − |g | . 2ω2 + ω1 ω2 3ω2 0002
Let η1 = r1 eiθ1 ,
η2 = r2 eiθ2 .
In polar coordinates, the amplitude equation generated from system (17) is r1 = r1 [Re(B11 )r12 + Re(B12 )r22 ], r2 = r2 [Re(B21 )r12 + Re(B22 )r22 ]. Combining the results of (15) and (24), the formal normal form arising from system (8) is ˆε η1 + εb12 τˆε η1 η = iω1 η1 + εb11 x 1 + B11 η12η 1 + B12 η1 η2η 2 , ˆε η2 + εb22 τˆε η2 η2 = iω2 η2 + εb21 x + B21 η1η 1 η2 + B22 η22η 2 . Using polar coordinates η1 = r1 eiθ1 , η2 = r2 eiθ2 , system (25) can be transformed into r1 = r1 [µ1 + Re(B11 )r12 + Re(B12 )r22 ], r = r2 [µ2 + Re(B21 )r 2 + Re(B22 )r 2 ], 2 1 2 θ1 = ω1 + ν1 + Im(B11 )r12 + Im(B12 )r22 , θ2 = ω2 + ν2 + Im(B21 )r12 + Im(B22 )r22 , xε + Re(bi2 )ˆ τε , νi = Im(bi1 )ˆ xε + Im(bi2 )ˆ τε for i = 1, 2. where µi = Re(bi1 )ˆ
(24)
(25)
(26)
December 9, 2009
11:8
02508
Dynamics and Double Hopf Bifurcations of the Rose–Hindmarsh Model
3747
In the following, we give two examples to illustrate that under certain conditions, no Hopf bifurcations occur for the limit cycle of system (4). Therefore, the unique attractor for system (4) is asymptotically stable. Example 1. Choose S = 1.0, x∗ = 0.1768, τ = 6.8693, then the characteristic equation has two pairs of
imaginary roots with ω1 = 0.9107, ω2 = 1.6810, and l1 = −0.572128, m1 = −12.544266, l2 = −88.8017, m2 = −13.87852, b11 = 0.02409516 + 0.6241099i, b21 = 0.8705159 + 0.743506i, √ 1 = ε(0.01632 − 0.3217559i), g2000 √ 1 = ε(0.01830 − 0.28222i), g0200 √ 1 = ε(0.39171 − 0.014319i), g1001 √ 1 = ε(0.42643 + 0.00472i), g0101 √ 1 = ε(0.033536 − 0.6039247i), g0011 √ 2 = ε(0.314199 − 0.036238i), g2000 √ 2 = ε(0.277067 − 0.049426i), g0200 √ 2 = ε(−0.08694 + 0.019076i), g1001 √ 2 = ε(−0.117646 + 0.04379i), g0101 √ 2 = ε(0.5916266 − 0.086679i), g0011
b12 = −0.013262 − 0.1159697i, b22 = 0.0212975 − 0.2134128i, √ 1 g1100 = ε(0.033536 − 0.6039247i), √ 1 g1010 = ε(0.36367 − 1.24877i), √ 1 g0110 = ε(0.33103 − 1.2263575i), √ 1 g0020 = ε(0.52571 − 0.651452i), √ 1 g0002 = ε(0.558073 − 0.005054i), √ 2 g1100 = ε(0.591627 − 0.086679i), √ 2 g1010 = ε(1.059699 + 0.461955i), √ 2 g0110 = ε(1.05146 + 0.423406i), √ 2 g0020 = ε(0.424199 + 0.55770i), √ 2 g0002 = ε(−0.182930 + 0.342064i).
Following the procedure introduced as above, we have √ √ E2 = − ε3.584712844, E1 = ε(1.006098900 − 2.909499977i), √ √ E4 = ε(−10.41936727 − 2.174529644i), E3 = ε(8.145198479 − 2.422005382i), √ √ E6 = ε(1.696799807 + 1.776121238i), E5 = ε(−10.41936727 + 2.174529644i), √ E7 = − ε3.267946002, and
√ 1 (0) = ε(0.804969 − 0.595855i), w2000 √ 1 (0) = − ε4.79444, w1100 √ 1 (0) = ε(8.3924166 − .947833i), w1010 √ 1 (0) = ε(−18.6899 + 1.140244i), w1001 √ 1 (0) = ε(−18.6899 − 1.14024i), w0110 √ 1 (0) = ε(1.697437 + 2.33727i), w0020 √ 1 (0) = − ε4.69736, w0011
√ 1 w2000 (−τ ) = ε(−0.638051 − 1.97351i), √ 1 w1100 (−τ ) = − ε4.249727, √ 1 w1010 (−τ ) = ε(5.641275 + 6.350949i), √ 1 w1001 (−τ ) = ε(−14.9901 + 10.1277i), √ 1 w0110 (−τ ) = ε(−14.9901 − 10.1277i), √ 1 w0020 (−τ ) = ε(−2.47207 + .83983i), √ 1 w0011 (−τ ) = − ε4.04632.
Hence, we can deduce that 1 = (−3.00532 − 4.099777i)ε, g2100 2 g1110 = (−10.76744 − 6.98252i)ε,
1 g1011 = (−0.06603 − 15.8757i)ε, 2 g0021 = (−3.94492 − 0.2161815i)ε,
which lead to B11 = (−5.5847 − 5.156339i)ε, B21 = (−4.544096 − 5.67959i)ε,
B12 = (−1.259567 − 28.3088i)ε, B22 = (−2.687968 + 0.3059236i)ε.
Furthermore, we find b11 = 0.018523 − 0.58490648i, b21 = −0.10710 + 0.0015i,
b12 = −0.01957 + 0.11507i, b22 = −0.20014 − 0.077108i.
December 9, 2009
02508
S. Ma et al.
3748
x(t)
11:8
3
0
2
y(t)−10 −20
1
−30
0
−40
−1
−50 −60
−2
−70
−3
−80
−4 −5
−90 320
340
360
380
−100
400
320
340
360
380
400
420
t
440
t
5.5
z(t) 5
5.5
z(t)
4.5
5
4.5
4
4 3.5 3.5 0
3
−50 2.5 0
100
Fig. 3.
200
300
400
500
t
y(t)
600
−100
−6
−4
−2
0
2
x(t)
The time series solution of system (4) and phase portraits at τ = 6.9 and x∗ = 0.1768.
We now discuss the restriction on phase portraits near the double Hopf bifurcation point (x∗ , τ ) and rewrite system (26) as r1 = r1 (µ1 + ar12 + br22 ), r2 = r2 (µ2 + cr12 + dr22 ), xε + Re(bi2 )ˆ τε , and where µi = Re(bi1 )ˆ Re(B11 ) = −1, |B11 |
b=
Re(B12 ) = −0.4993579, |B22 |
Re(B21 ) = −0.033449, |B11 |
d=
Re(B22 ) = −1. |B22 |
a= c=
We observe that no Hopf bifurcation occurs for the limit cycle in this case. Therefore, the unique attractor for system (4) is asymptotically stable, as shown in Fig. 3.
December 9, 2009
11:8
02508
Dynamics and Double Hopf Bifurcations of the Rose–Hindmarsh Model
3749
Example 2. Choose S = 1.0, x∗ = 1.8207, τ = 2.9934, then the characteristic equation has two pairs of
imaginary roots with ω1 = 4.1953, ω2 = 6.1936, and l1 = −10203.38178, m1 = −4904.083215, l2 = −51514.46280, m2 = −1015.974613, b11 = −0.156559 + 0.06985i, b21 = −0.12339 + 0.14264i, √ 1 g2000 = ε(−0.0210699 + 0.08338i), √ 1 = ε(−0.02331 + 0.0816597i), g0200 √ 1 = √ε(−0.090015 + 0.139878i), g1001 1 g0101 = ε(−0.09268 + 0.13892i), √ 1 = ε(−0.04437 + .165024i), g0011 √ 2 g2000 = ε(−0.06336 + 0.0484i), √ 2 = ε(−0.06432 + 0.04570i), g0200 √ 2 = ε(−0.15256 + 0.04753i), g1001 √ 2 g0101 = ε(−.154296 + 0.045248i), √ 2 = ε(−.127664 + 0.094103i), g0011
b12 = −0.016267 − 1.325248i, b22 = 0.0176 − 1.98267i, √ 1 g1100 = ε(−0.04437 + 0.165024i), √ 1 g1010 = ε(−0.0066 + 0.2049i), √ 1 g0110 = √ε(−0.00825 + 0.202589i), 1 g0020 = ε(0.00797 + 0.125942i), √ 1 g0002 = ε(−0.07195 + 0.0646i), √ 2 g1100 = ε(−.12766 + 0.0941i), √ 2 g1010 = ε(−0.1176 + 0.14883i), √ 2 g0110 = ε(−0.1177 + 0.145965i), √ 2 g0020 = ε(−0.06219 + 0.10060i), √ 2 g0002 = ε(−0.0962 + 0.004395i).
Following the procedure introduced as above, we have √ √ E2 = ε0.62986142, E1 = ε(0.1604619 + 0.398838i), √ √ E4 = ε(−0.91337 + 0.52039i), E3 = ε(0.6055876 + 0.98356i), √ √ E6 = ε(0.28760 + 0.73738i), E5 = ε(−0.91337 − 0.5203899i), √ E7 = − ε0.864976, and 1 (0) = w2000
√
ε(0.12817 + 0.35871i), √ 1 (0) = − ε0.73892, w1100 √ 1 w1010 (0) = ε(0.549281 + .938794i), √ 1 (0) = ε(−0.75798 + 0.52173i), w1001 √ 1 (0) = ε(−0.75798 − 0.52173i), w0110 √ 1 (0) = ε(0.26012 + 0.71879i), w0020 √ 1 (0) = − ε0.75592, w0011
1 w2000 (−τ ) =
√
ε(0.13001 + 0.354975i), √ 1 w1100 (−τ ) = − ε0.724738, √ 1 w1010 (−τ ) = ε(0.21885 + 1.06856i), √ 1 w1001 (−τ ) = ε(−0.57858 + 0.76051i), √ 1 w0110 (−τ ) = ε(−0.57858 − 0.76051i), √ 1 w0020 (−τ ) = ε(−0.2191 + 0.74435i), √ 1 w0011 (−τ ) = − ε0.7701.
Consequently, we have 1 = (0.017687 − 0.102537i)ε, g2100 2 = (−0.11286 − 0.13893i)ε, g1110
1 g1011 = (−0.1316795 − 0.285369i)ε, 2 g0021 = (−0.128979 − 0.130854i)ε,
which lead to B11 = (−0.00728 − 0.103909i)ε,
B12 = (−0.15409 − 0.318484i)ε,
B21 = (−0.08765 − 0.17090i)ε,
B22 = (−0.11857 − 0.1477i)ε.
Moreover, we obtain b11 = −0.156617 − 0.06702i, b21 = −0.15513 − 0.04030i,
b12 = −0.03789 + 1.3248i, b22 = −1.13648 + 1.6247i.
December 9, 2009
11:8
02508
S. Ma et al.
3750
10
x(t)
1
y(t)
9 8
0.5
7 6
0
5 −0.5
4
−1
3 2
−1.5 −2 700
1 705
710
715
0 700
720
705
710
715
720
t
t
5.5
z(t) 5 4.5 2.8
4
2.6
z(t)
3.5
2.4
1
3
2.2
2.5
2 10
0
2 0
200
Fig. 4.
400
600
800
t
8
−1
6
y(t)
1000
4
2
0
−2
x(t)
−2
The time series solution of system (4) and phase portraits at τ = 3.1 and x∗ = 1.8207.
Now, we discuss the restriction on phase portraits near the double Hopf bifurcation point (x∗ , τ ) and rewrite system (26) as r1 = r1 (µ1 + ar12 + br22 ), r2 = r2 (µ2 + cr12 + dr22 ), xε + Re(bi2 )ˆ τε , and where µi = Re(bi1 )ˆ Re(B11 ) = −1, |B11 |
b=
Re(B12 ) = −0.813467, |B22 |
Re(B21 ) = −0.8414566, |B11 |
d=
Re(B22 ) = −1. |B22 |
a= c=
Clearly, one can see that no Hopf bifurcation for the limit cycle occurs in this case. Therefore, the unique attractor for system (4) is asymptotically stable, as illustrated in Fig. 4.
4. Conclusion The dynamical behavior of the Rose–Hindmarsh model with time delay was intensively investigated. By applying the generalized Sturm criterion, which is a powerful tool for real polynomials with constant
December 9, 2009
11:8
02508
Dynamics and Double Hopf Bifurcations of the Rose–Hindmarsh Model
coefficients, classifications of the imaginary roots of the characteristic equation were presented. It was proved that the Rose–Hindmarsh system was absolutely stable for any time delay if corresponding parameters lie in the regime of II(b) and IV. Otherwise, Hopf bifurcation curves were drawn in the (x∗ , τ ) plane by using the software DDE-Biftool. The double Hopf bifurcation points were acquired at the intersection points of two Hopf bifurcation curves. The normal form and universal unfolding at double Hopf bifurcation points were studied by the theory of center manifold theory. Two examples were demonstrated and computational results of normal form at double Hopf bifurcation points reveal that the associated unique attractor was asymptotically stable.
References Buono, P. L. & B´elair, J. [2003] “Restrictions and unfolding of double Hopf bifurcation in functional differential equations,” J. Diff. Eqs. 189, 234–266. Campell, S. A. & B´elair, J. [1995a] “Analytical and symbolically-assisted investigation of Hopf bifurcations in delay-differential equations,” Canad. Appl. Math. Quart. 3, 137–154. Campell, S. A., B´elair, J., Ohira, T. & Milton, J. [1995b] “Limit cycles, tori, and complex dynamics in a second-order differential equations with delayed negative feedback,” J. Dyn. Diff. Eqs. 7, 213–236. Elphick, C., Tiraopegui, E., Brachet, M. E., Coullet, P. & Iooss, G. [1987] “A simple global characterization for normal forms of singular vector fields,” Phys. D 29, 95–127. Grebogi, C., Ott, E., Romeiras, F. & Yorke, J. A. [1986] “Critical exponent of chaotic transients in nonlinear dynamical systems,” Phys. Rev. Lett. 57, 1284–1287. Grebogi, C., Ott, E., Romeiras, F. & Yorke, J. A. [1987] “Critical exponents for crisis-induced intermittency,” Phys. Rev. A 36, 5365–5380. Guckenheimer, J. M. & Holmes, P. J. [1983] Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields (Springer, NY). Hindmarsh, J. L. & Rose, R. M. [1982] “A model of the nerve impulse using two first-order differential equations,” Nature 296, 162–164.
3751
Hindmarsh, J. L. & Rose, R. M. [1984] “A model of neuronal bursting using three coupled first-order differential equations,” Proc. Roy. Soc. Lond. B 221, 87–102. Holden, A. V. & Fan, Y. S. [1992a] “From simple to simple bursting oscillatory behavior via chaos in the Rose–Hindmarsh model for neuronal activity,” Chaos Solit. Fract. 2, 221–236. Holden, A. V. & Fan, Y. S. [1992b] “From simple to complex oscillatory behavior via intermittent chaos in the Rose–Hindmarsh model for neuronal activity,” Chaos Solit. Fract. 2, 346–369. Ma, S. Q., Lu, Q. S. & Hogan, S. J. [2007] “Double Hopf bifurcation for Stuart-Landau system with nonlinear delay feedback and delay dependant parameters,” Adv. Compl. Syst. 10, 1–27. Ma, S. Q., Feng, Z. & Lu, Q. S. [2008a] “A two-patch ecological system with nonlinear transfer rate and noise effect,” Dyn. PDEs 5, 281–298. Ma, S. Q., Lu, Q. S. & Feng, Z. [2008b] “Double Hopf bifurcation for van der Pol-Duffing oscillator with parametric delay feedback control,” J. Math. Anal. Appl. 338, 993–1007. Pomeau, Y. & Manneville, P. “Intermittent transition to turbulence in dissipative dynamic system,” Commun. Math. Phys. 74, 189–197. Shi, X., Lu, Q. S. & Chen, G. [2007] “Anti-phase synchronization of inhibitorily coupled neurons,” Int. J. Bifurcation and Chaos 17, 4355–4364. Wang, Z. & Hu, H. Y. [1999] “Delay independent stability of retarded dynamic systems of multiple degrees of freedom,” J. Sound Vibr. 226, 57–81. Wang, Z. & Hu, H. Y. [2000] “Stability switches of timedelayed dynamic systems with unknown parameters,” J. Sound Vibr. 233, 215–233. Xu, J. & Lu, Q. S. [1999] “Hopf bifurcation of time-delay Li´enard equations,” Int. J. Bifurcation and Chaos 9, 939–951. Yang, L., Hou, X. R. & Zeng, Z. B. [1996a] “A complete discrimination system for polynomial,” Sci. in China (Series E) 26, 424–441. Yang, L., Hou, X. R. & Zeng, Z. B. [1996b] Nonlinear Algebraic Equation Systems and Automated Theory (in Chinese) (Shanghai Scientific and Technological Education Press, Shanghai). Zhao, X. Q. & Zou, X. [2001] “Threshold dynamics in a delayed SIS epidemic model,” J. Math. Anal. Appl. 257, 282–291.
