Electronic Journal of Qualitative Theory of Differential Equations 2010, No. 63, 1-22; http://www.math.u-szeged.hu/ejqtde/

Bifurcation analysis of Rössler system with multiple delayed feedback^{∗}

Meihong Xu^{a}, Yuan Wei^{b}, Junjie Wei^{a} ^{†}

aDepartment of Mathematics, Harbin Institute of Technology Harbin, 150001, P.R.China b Beijing Institute of Strength and Environment, Beijing, P.R. China

Abstract

In this paper, regarding the delay as parameter, we investigate the effect of delay on the dynamics of a Rössler system with multiple delayed feedback proposed by Ghosh and Chowdhury. At first we consider the stability of equilibrium and the existence of Hopf bifurcations. Then an explicit algorithm for determining the direction and the stability of the bifurcating periodic solutions is derived by using the normal form theory and center manifold argument. Finally, we give a numerical simulation example which indicates that chaotic oscillation is converted into a stable steady state or a stable periodic orbit when the delay passes through certain critical values. Keywords: Rössler system; delayed feedback control; Hopf bifurcation

1 Introduction

The study of chaotic systems has increasingly gained interest of many researchers since the pioneering work of Lorenz [8]. For a quite long period of time, people thought that chaos was neither predictable nor controllable. Recently the trend of analyzing and understanding chaos has been extended to controlling and utilizing chaos. The main goal of chaos control was to eliminate chaotic behavior and to stabilize the chaotic system at one of the system’s equilibrium points. More specially, when it is useful, we want to generate chaos intentionally. Until now, many advanced theories and methodologies have been developed for controlling chaos. Many scientists have more concerns with delayed control (see Guan, Chen and Peng [5], Shu et al.[16], Zhang and Su [21]). The existing control method can be classified, mainly, into two categories. The first one, the OGY method developed by Ott, Grebogi and Yorke [9] in 1990s of the last century has completely changed the chaos research topic. The second one, proposed by Pyragas [10, 11], using time-delayed controlling forces. Compared with the first one, it is much simpler and more convenient on controlling chaos in continuous dynamics system. Here, we mainly study the Rössler system with delayed controlling method developed by Pyragas. Rössler system is described by the following three-dimensional smooth autonomous system (see Rössler [13])

∗Supported by the National Natural Science Foundation of China (No10771045) and the Program of Excellent Team in HIT.

†Corresponding author. E-mail address: weijj@hit.edu.cn

\begin{matrix}ẋ(t) = −y(t) − z(t), \\ ẏ(t) = x(t) + β_{1}y(t), \\ ż(t) = β_{2} + z(t) (x(t) − γ) , \end{matrix}  (1)

which is chaotic when β_{1} = β_{2} = 0.2, γ = 5.7.

Rössler system is a quite simple set of differential equations with chaos to simplify the Lorenz model of turbulence that contains just one (second order) nonlinearity in one variable. Due to its simplicity, the Rössler system has become a standard one to issue the effectiveness of the chaos control strategy. Recently, many literatures adopted controlling strategy for the Rössler system. In the last years there are many studies on Rössler system. For example, Pyragas [10] stabilized unstable periodic orbits of a Rössler system to a desired periodic orbit by self-controlling feedback. Tao et al. [18] used the speed feedback control such that the controlled Rössler system will gradually converge to unsteadily equilibrium point. Tian et al. [19] used a nonlinear open-plus-closed-loop (NOPCL) control to Rössler system. Yang et. al. [20] presented an impulsive control to achieve periodic motions for the Rössler system. Moreover there are extensive study, for example Amhed et al. [1], Chang et al. [2], Chen et al. [3], Rasussen et al. [12]. Recently, Ghosh et al. [4] have proposed the multiple delayed system in the following form:

\begin{matrix}ẋ(t) = −y(t) − z(t) + α_{1}x(t − τ_{1}) + α_{2}x(t − τ_{2}), \\ ẏ(t) = x(t) + β_{1}y(t), \\ ż(t) = β_{2} + z(t) (x(t) − γ) , \end{matrix}  (2)

where α_{i}, β_{i} (i = 1, 2) and γ are all positive constants. They studied the system (2) in numerical simulations mainly. The purpose of the present paper is to investigate system (2) analytically and numerically. Our analytical results show that the stability changes as the delays vary. Meanwhile, our numerical simulations indicate that chaotic oscillation is converted into a stable steady state or a stable periodic orbit when the delay passes through certain critical values. This shows that the chaos property changes as the delay varies.

The rest of the paper is organized as follows. In Section 2, we study the distribution of the eigenvalues by using the result due to Ruan and Wei [14, 15] on the analysis of distribution of the zeros of exponential polynomial. Hence the stability and existence of Hopf bifurcations are obtained. In Section 3, the direction and stability of the Hopf bifurcation are determined by using the center manifold and normal forms theory. Some numerical simulations are carried out for supporting the analysis results in Section 4. Conclusions and discussions are given in Section 5.

2 Analysis of stability and bifurcation

In this section, we shall study the stability of the interior equilibrium and the existence of local Hopf bifurcations. For convenience, denote

A = 1 + (α_{1} + α_{2}) β_{1}.

Proposition 2.1. (i) If γ^{2}A < 4β_{1}β_{2}, then system (2) has no equilibrium, and if γ^{2}A = 4β_{1}β_{2}, the system has only one equilibrium given by

(\frac{γ}{2} , \frac{−γ}{2β_{1}} , \frac{2β_{2}}{γ} ) .

(ii) If

(H_{0}) γ^{2} > \frac{4β_{1}β_{2}}{A}

holds, then the system (2) has two equilibria (x_{0}, y_{0}, z_{0}) and (x_{1}, y_{1}, z_{1}), where

x_{0} = −β_{1}X_{+}, y_{0} = X_{+}, z_{0} = \frac{β_{2}}{β_{1}X_{+} + γ} , x_{1} = −β_{1}X_{−}, y_{1} = X_{−}, z_{1} = \frac{β_{2}}{β_{1}X_{−} + γ} ,

and

X_{±} = \frac{1}{2β_{1}} [−γ ± \sqrt{γ^{2} − \frac{4β_{1}β_{2}}{A}}].

Thorough out this paper, we always assume that (H_{0}) is satisfied and only consider the dynamics of system (2) near the equilibrium (x_{0}, y_{0}, z_{0}).

Let u_{1} = x − x_{0}, u_{2} = y − y_{0} and u_{3} = z − z_{0}. Then system (2) can be written in the following form

\begin{matrix}u̇_{1}(t) = −u_{2}(t) − u_{3}(t) + α_{1}u_{1}(t − τ_{1}) + α_{2}u_{1}(t − τ_{2}), \\ u̇_{2}(t) = u_{1}(t) + β_{1}u_{2}(t), \\ u̇_{3}(t) = u_{1}(t)u_{3}(t) + u_{3}(t)(x_{0} − γ) + z_{0}u_{1}(t). \end{matrix}  (3)

The characteristic equation of system (3) at the equilibrium (0, 0, 0) is

λ^{3} +a_{2}λ^{2} +a_{1}λ+a_{0} −α_{1}(λ^{2} +a_{2}λ+b_{0})e^{−λτ_{1}} −α_{2}(λ^{2} +a_{2}λ+b_{0})e^{−λτ_{2}} = 0,  (4)

where

a_{0} = γ − x_{0} − β_{1}z_{0},
a_{1} = −β_{1}γ + β_{1}x_{0} + 1 + z_{0},
a_{2} = γ − x_{0} − β_{1},
b_{0} = −β_{1}γ + β_{1}x_{0}.

Now we employ the method due to Ruan and Wei [14, 15] to investigate the distribution of roots of Eq.(4).

When τ_{1} = 0 and τ_{2} = 0, Eq.(4) becomes

λ^{3} + (a_{2} − α_{1} − α_{2})λ^{2} + (a_{1} − α_{1}a_{2} − α_{2}a_{2})λ + a_{0} − α_{1}b_{0} − α_{2}b_{0} = 0. (5)

For convenience, we make the following hypothesis:

(H_{1}) \{ \begin{matrix}a_{0} − α_{1}b_{0} − α_{2}b_{0} > 0, a_{2} − α_{1} − α_{2} > 0, \\ (a_{2} − α_{1} − α_{2})(a_{1} − α_{1}a_{2} − α_{2}a_{2}) − a_{0} + α_{1}b_{0} + α_{2}b_{0} > 0. \end{matrix}

By Routh-Hurwitz criterion we know that if (H_{1}) holds, then all roots of Eq.(5) have negative real parts. Let iν(τ_{1})(ν > 0) be a root of Eq.(4) with τ_{2} = 0. Then it follows that

\{ \begin{matrix}−a_{2}ν^{2} + a_{0} + α_{2}ν^{2} − α_{2}b_{0} = (α_{1}b_{0} − α_{1}ν^{2}) cos ντ_{1} + α_{1}a_{2}ν sin ντ_{1}, \\ −ν^{3} + a_{1}ν − α_{2}a_{2}ν = α_{1}a_{2}ν cos ντ_{1} − (α_{1}b_{0} − α_{1}ν^{2}) sin ντ_{1}, \end{matrix}  (6)

which leads to

ν^{6} + pν^{4} + qν^{2} + r = 0,  (7)

where

p = (α_{2} − a_{2})^{2} + 2(α_{2}a_{2} − a_{1}) − α_{1}^{2},
q = 2(α_{2} − a_{2})(a_{0} − α_{2}b_{0}) − (α_{2}a_{2} − a_{1})^{2} + 2b_{0}α_{1}^{2} − α_{1}^{2}a_{2}^{2},
r = (a_{0} − α_{2}b_{0})^{2} − α_{1}^{2}b_{0}^{2}.

If we set z = ν^{2}, then Eq.(7) becomes

z^{3} + pz^{2} + qz + r = 0.  (8)

Denote

h(z) = z^{3} + pz^{2} + qz + r.

Hence

\frac{dh(z)}{dz} = 3z^{2} + 2pz + q.

Set

△ = p^{2} − 3q.

When r ≥ 0 and △ > 0, the following equation

3z^{2} + 2pz + q = 0

has two real roots

z_{1}^{∗} = \frac{−p + \sqrt{△}}{3} and z_{2}^{∗} = \frac{−p − \sqrt{△}}{3}.

Applying the Claim 3 and Lemma 2.1 of Ruan and Wei [14], we have the following conclusions. Lemma 2.1. For Eq.(8), we have following conclusions:

(i) If r < 0, then Eq.(8) has at least one positive root. (ii) If r ≥ 0 and △ ≤ 0, then Eq.(8) has no positive root. (iii) If r ≥ 0 and △ > 0, then Eq.(8) has one positive roots if and only if z_{1}^{∗} > 0

and h(z_{1}^{∗}) ≤ 0. Without loss of generality, we assume that Eq.(8) has three positive roots,

denoted by x_{1}, x_{2} and x_{3}, respectively. Then Eq.(7) has three positive roots ν_{k} = \sqrt{x_{k}}, k = 1, 2, 3. Substituting ν_{k} into (6) gives

cos ν_{k}τ_{1} = \frac{(b_{0} − ν_{k}^{2})[(α_{2} − a_{2})ν_{k}^{2} + a_{0} − α_{2}b_{0}] + a_{2}ν_{k}[(a_{1} − α_{2}a_{2})ν_{k} − ν_{k}^{3}]}{α_{1}(b_{0} − ν_{k}^{2})^{2} + α_{1}a_{2}^{2}ν_{k}^{2}}.

Let

τ_{1k}^{(j)} = \frac{1}{ν_{k}} [arccos \frac{(b_{0} − ν_{k}^{2})[(α_{2} − a_{2})ν_{k}^{2} + a_{0} − α_{2}b_{0}] + a_{2}ν_{k}[(a_{1} − α_{2}a_{2})ν_{k} − ν_{k}^{3}]}{α_{1}(b_{0} − ν_{k}^{2})^{2} + α_{1}a_{2}^{2}ν_{k}^{2}}+2jπ],

EJQTDE, 2010 No. 63, p. 5 where k = 1, 2, 3; j = 0, 1, ⋯, and τ_{1k}^{(0)}ν_{k} ∈ (0, 2π] is determined by the sign of sin(ν_{k}τ_{1}^{(0)}). These show that (ν_{k}, τ_{1k}^{(j)}) is a root of Eq.(6). This shows that ±iν_{k} is a pair of purely imaginary roots of Eq.(4) with τ_{2} = 0 when τ_{1} = τ_{1k}^{(j)}.

Define

τ_{1}^{0} = min_{1≤k≤3,j≥0}{τ_{1k}^{(j)}}.

Let λ(τ_{1}) = α(τ_{1}) + iν(τ_{1}) be the root of Eq.(4) with τ_{2} = 0 satisfying α(τ_{1k}^{(j)}) = 0, ν(τ_{1k}^{(j)}) = ν_{k}. Then the following transversality condition holds.

Lemma 2.2. Suppose that z_{k} = ν_{k}^{2}, and h′(z_{k}) ≠ 0. Then \frac{d(Reλ(τ_{1k}^{(j)}))}{dτ_{1}} ≠ 0 and the sign of \frac{d(Reλ(τ_{1k}^{(j)}))}{dτ_{1}} is coincident with that of h′(z_{k}).

Proof. Substituting λ(τ_{1}) into Eq(4) with τ_{2} = 0 and differentiating both sides with respect to τ_{1}, it follows that

(\frac{dλ}{dτ_{1}} )^{−1} = −\frac{[3λ^{2} + 2(a_{2} − α_{2})λ + a_{1} − α_{2}a_{2}]e^{λτ_{1}}}{α_{1}λ(λ^{2} + a_{2}λ + b_{0})} + \frac{2λ + a_{2}}{λ(λ^{2} + a_{2}λ + b_{0})} + \frac{τ_{1}}{λ} .

And hence,

[\frac{d(Reλ(τ_{1}))}{dτ_{1}}]_{τ_{1}=τ_{1k}^{(j)}}^{−1} = Re [\frac{(3λ^{2} + 2(a_{2} − α_{2})λ + a_{1} − α_{2}a_{2})e^{λτ_{1}}}{−α_{1}λ(λ^{2} + a_{2}λ + b_{0})} ]_{τ_{1}=τ_{1k}^{(j)}}

+ Re [\frac{2λ + a_{2}}{λ(λ^{2} + a_{2}λ + b_{0})}]_{τ_{1}=τ_{1k}^{(j)}} + Re(\frac{τ_{1}}{λ})_{τ_{1}=τ_{1k}^{(j)}}

= \frac{1}{Λ}\{3ν_{k}^{6} + 2[(α_{2} − a_{2})^{2} + 2(α_{2}a_{2} − a_{2}) − α_{1}^{2}]ν_{k}^{4}

+ [(a_{1} − α_{2}a_{2})^{2} + 2(α_{2} − a_{2})(a_{0} − α_{2}b_{0}) − α_{1}^{2}a_{2}^{2} + 2b_{0}α_{1}^{2}]ν_{k}^{2}}

= \frac{1}{Λ}(3ν_{k}^{6} + 2pν_{k}^{4} + qν_{k}^{2})

= \frac{z_{k}}{Λ} h′(z_{k}),

where

Λ = α_{1}^{2}[a_{2}^{2}ν_{k}^{4} + (ν_{k}^{2} − b_{0})^{2}ν_{k}^{2}].

Notice that Λ > 0 and z_{k} > 0, we conclude that

sign[\frac{d(Reλ(τ_{1}))}{dτ_{1}}]_{τ_{1}=τ_{1k}^{(j)}}^{−1} = sign\{h′(z_{k})\}.

This completes the proof.

By the Lemmas 2.1, 2.2 and applying the Hopf bifurcation theorem for functional differential equations (see Hale [6], Chapter 11, Theorem 1.1), we can conclude the existence of Hopf bifurcation as stated in the following theorem.

Theorem 2.3. Suppose that (H_{1}) is satisfied and τ_{2} = 0.

(i) If r > 0 and △ = p^{2} − 3q ≤ 0, then all roots of Eq.(4) have negative real parts for all τ_{1} ≥ 0, and hence the equilibrium (x_{0}, y_{0}, z_{0}) of system (2) is asymptotically stable for all τ_{1} ≥ 0.

(ii) If either r < 0 or r ≥ 0 and △ > 0, z_{1}^{∗} > 0 and h(z_{1}^{∗}) ≤ 0 hold, then h(z) has at least a positive root z_{k}, all roots of Eq.(4) have negative real parts for τ_{1} ∈ [0, τ_{1}^{0}). Hence the equilibrium (x_{0}, y_{0}, z_{0})of system (2) is asymptotically stable when τ_{1} ∈ [0, τ_{1}^{0}).

(iii) If all conditions as stated in (ii) and h′(z_{k}) ≠ 0 hold, then system (2) undergoes a Hopf bifurcation at the equilibrium (x_{0}, y_{0}, z_{0}), when τ_{1} = τ_{1k}^{(j)}(j = 0, 1, 2, ⋯).

We have known that (H_{1}) ensures that all roots of Eq.(5) have negative real parts. Now we consider (H_{1}) is not satisfied. For convenience, denote

a = (a_{2} − α_{1} − α_{2}), b = (a_{1} − α_{1}a_{2} − α_{2}a_{2}), c = a_{0} − α_{1}b_{0} − α_{2}b_{0}.

Then Eq.(5) becomes

λ^{3} + aλ^{2} + bλ + c = 0.

Let λ = X − a/3. Then it reduces to

X^{3} + p_{1}X + q_{1} = 0,  (9)

where p_{1} = b − \frac{1}{3}a^{2} and q_{1} = \frac{2}{27}a^{3} − \frac{1}{3}ab + c.

Denote

Δ_{1} = (\frac{p_{1}}{3} )^{3} + (\frac{q_{1}}{2} )^{2}, ε = \frac{−1}{2} + \frac{\sqrt{3}}{2}i,

α = \sqrt[3]{−\frac{q_{1}}{2} + \sqrt{Δ_{1}}} and β = \sqrt[3]{−\frac{q_{1}}{2} − \sqrt{Δ_{1}}}.

Then from Cardano’s formula for the third degree algebra equation we have the followings.

Proposition 2.2. If Δ_{1} > 0, then Eq.(9) has a real root α+β and a pair of conjugate complex roots −\frac{α+β}{2} ± i\frac{\sqrt{3}}{2}(α − β), that is, that Eq.(5) has a real root α + β − \frac{a}{3} and a pair of conjugate complex roots −(\frac{α+β}{2} + \frac{a}{3}) ± i\frac{\sqrt{3}}{2}(α − β).

We make the following assumptions:

(H_{2}) Δ_{1} > 0, α + β − \frac{a}{3} < 0, \frac{α + β}{2} + \frac{a}{3} < 0, α − β ≠ 0.

Theorem 2.4. Suppose that (H_{2}) is satisfied and τ_{2} = 0.

(i) If r > 0 and △ = p^{2} − 3q ≤ 0, then at least one of the roots of Eq.(4) has positive real parts for all τ_{1} ≥ 0, and hence the equilibrium (x_{0}, y_{0}, z_{0}) of system (2) is unstable for all τ_{1} ≥ 0.

(ii) If either r < 0 or r ≥ 0 and △ > 0, z_{1}^{∗} > 0 and h(z_{1}^{∗}) ≤ 0 hold, then h(z) has at least a positive root z_{k}, at least one of the roots of Eq.(4) has positive real parts for τ_{1} ∈ [0, τ_{1}^{0}). Hence the equilibrium (x_{0}, y_{0}, z_{0}) of system (2) is unstable when τ_{1} ∈ [0, τ_{1}^{0}). Additionally, if \frac{d(Reλ(τ_{1}))}{dτ_{1}}|_{τ_{1}=τ_{1}^{0}} < 0, then the equilibrium (x_{0}, y_{0}, z_{0})of system (2) is asymptotically stable when τ_{1} ∈ (τ_{1}^{0}, τ_{1}^{1}), where τ_{1}^{1} is the second critical value.

(iii) If all conditions as stated in (ii) and h′(z_{k}) ≠ 0 hold, then system (2) undergoes a Hopf bifurcation at the equilibrium (x_{0}, y_{0}, z_{0}), when τ_{1} = τ_{1k}^{(j)}(j = 0, 1, 2, ⋯).

From the discussions above, we know that there may exist stability switches as τ_{1} varies for system (2) with τ_{2} = 0. So denote I as stable interval of τ_{1}, that is the equilibrium (x_{0}, y_{0}, z_{0}) of Eq.(2) is asymptotically stable when τ_{1} ∈ I and τ_{2} = 0. Let τ_{1} ∈ I, and λ = iω(τ_{2}) (ω > 0) be a root of Eq.(4). Then we obtain

\{ \begin{matrix}−a_{2}ω^{2} + a_{0} − α_{1}(b_{0} − ω^{2}) cos ωτ_{1} − α_{1}a_{2}ω sin ωτ_{1} = α_{2}(b_{0} − ω^{2}) cos ωτ_{2} + α_{2}a_{2}ω sin ωτ_{2}, \\ a_{1}ω − ω^{3} − α_{1}a_{2}ω cos ωτ_{1} + α_{1}(b_{0} − ω^{2}) sin ωτ_{1} = α_{2}a_{2}ω cos ωτ_{2} − α_{2}(b_{0} − ω^{2}) sin ωτ_{2}. \end{matrix}

Then we have

\begin{matrix}ω^{6} + (a_{2}^{2} − 2a_{1} − α_{2}^{2} + α_{1}^{2})ω^{4} + [a_{1}^{2} − 2a_{0}a_{2} + (2b_{0} − a_{2}^{2})(α_{2}^{2} − α_{1}^{2})]ω^{2} + a_{0}^{2} \\ −b_{0}^{2}(α_{2}^{2} − α_{1}^{2}) − 2[α_{1}(a_{0} − a_{2}ω^{2})(b_{0} − ω^{2}) + α_{1}a_{2}ω(a_{1}ω − ω^{3})] cos ωτ_{1} \\ −2[α_{1}a_{2}ω(a_{0} − a_{2}ω^{2}) − α_{1}(a_{1}ω − ω^{3})(b_{0} − ω^{2})] sin ωτ_{1} = 0, \end{matrix} (10)

and

cos ωτ_{2} = \frac{1}{α_{2}[a_{2}^{2}ω^{2}+(b_{0}−ω^{2})^{2}]}\{(b_{0} − ω^{2})[−a_{2}ω^{2} + a_{0} − α_{1}(b_{0} − ω^{2}) cos ωτ_{1} − α_{1}a_{2}ω sin ωτ_{1}]

+ a_{2}ω[a_{1}ω − ω^{3} − α_{1}a_{2}ω cos ωτ_{1} + α_{1}(b_{0} − ω^{2}) sin ωτ_{1}]}. (11)

We know that Eq.(10) has at most finite positive roots. Without loss of generality, we assume that Eq.(10) has N positive roots, denoted by ω_{k} (k = 1, 2, ⋯, N). According to (11), define

τ_{2k}^{(j)} = \frac{1}{ω_{k}} [± arccos(A) + 2jπ], k = 1, 2, ⋯, N; j = 0, 1, 2, ⋯,

where A is the value of right hand of (11) with ω = ω_{k} and the sign of arc-cosine function are determined by sin ωτ_{2}. Then ±iω_{k} is a pair of purely imaginary roots of Eq.(4) with τ_{2k}^{(j)}. Denote

τ_{2}^{0} = τ_{2k_{0}}^{(0)} = min_{k∈{1,2,··· ,N}}{τ_{2k}^{(0)}}, ω_{0} = ω_{k_{0}}.

Let λ(τ_{2}) = α(τ_{2}) + iω(τ_{2}) be the root of Eq.(4) satisfying α(τ_{2k}^{(j)}) = 0, ω(τ_{2k}^{(j)}) = ω_{k}. By computation, we obtain

α′(τ_{2}^{0}) = Λ_{1}{RS − T U − 2α_{2}^{2}ω_{0}^{2} + 2b_{0}α_{2}^{2}ω_{0}^{2} − a_{2}^{2}α_{2}^{2}ω_{0}^{2} + (RBω_{0} − AT ω_{0}) sin ω_{0}τ_{2}^{0} + (−QB − P A)ω_{0} cos ω_{0}(τ_{1} + τ_{2}^{0}) + (−AQ + P B)ω_{0} sin ω_{0}(τ_{1} + τ_{2}^{0})

− (ARω_{0} + BT ω_{0}) cos ω_{0}τ_{2}^{0} + (SP + UQ) cos ω_{0}τ_{1} + (QS − P U) sin ω_{0}τ_{1}}^{−1},

where

P = α_{1}(−a_{2} + b_{0}τ_{1} − τ_{1}ω_{0}^{2}), Q = α_{1}(−2ω_{0} + a_{2}τ_{1}ω_{0}), R = −3ω_{0}^{2} + a_{1},
Λ_{1} = a_{2}^{2}α_{2}^{2}ω_{0}^{4} + (b_{0}^{2} − ω_{0}^{2})^{2}α_{2}^{2}ω_{0}^{2}, A = α_{1}a_{2}ω_{0}, B = α_{1}(b_{0} − ω_{0}),
S = a_{1}ω_{0}^{2} − ω_{0}^{4}, T = 2a_{2}ω_{0}, U = −a_{2}ω_{0}^{3} + a_{0}ω_{0}.

Summarizing the discussions above, we have the following conclusions.

Theorem 2.5. Suppose that either (H_{1}) or (H_{2}) is satisfied, τ_{1} ∈ I and Eq.(10) has positive roots. Then all roots of Eq.(4) have negative real parts for τ_{2} ∈ [0, τ_{2}^{0}). Furthermore, the equilibrium (x_{0}, y_{0}, z_{0}) of system (2) is asymptotically stable when τ_{2} ∈ [0, τ_{2}^{0}). Additionally, if α′(τ_{2}^{0}) ≠ 0, then system (2) undergoes a Hopf bifurcation at the equilibrium (x_{0}, y_{0}, z_{0}) when τ_{2} = τ_{2}^{0}.

3 Stability and direction of the Hopf bifurcation

In the previous section, we obtained conditions for Hopf bifurcation to occur when τ_{2} = τ_{2}^{0}. In this section we study the direction of the Hopf bifurcation and the stability of the bifurcating periodic solutions when τ_{2} = τ_{2}^{0}, using techniques from normal form and center manifold theory (see e.g. Hassard et al.[7]).

We assume τ_{1} > τ_{2}^{0} and denote τ_{2} = τ_{2}^{0} + μ. Then the system (3) can be written as an FDE in C = C([−τ_{1}, 0], R^{3}) as

u̇(t) = L_{μ}(u_{t}) + F (μ, u_{t}),  (12)

where u = (u_{1}, u_{2}, u_{3})^{T}, u_{t}(θ) = u(t + θ) ∈ C, and L_{μ} : C → R^{3}, F : R × C → R^{3} are given, respectively, by

L_{μ}φ = Aφ(0) + B_{1}φ(−τ_{1}) + B_{2}φ(−(τ_{2}^{0} + μ)),  (13)

where

A = ( \begin{matrix}0 & −1 & −1 \\ 1 & β_{1} & 0 \\ z_{0} & 0 & x_{0} − γ \end{matrix} ) ,

B_{1} = ( \begin{matrix}α_{1} & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix} ) , B_{2} = ( \begin{matrix}α_{2} & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix} ) ,

φ(t) = (φ_{1}(t), φ_{2}(t), φ_{3}(t))^{T} ,

and

F (μ, φ) = ( \begin{matrix}0 \\ 0 \\ φ_{1}(0)φ_{3}(0) \end{matrix} ) .

From the discussion in Section 2, we know that system (2) undergoes a Hopf bifurcation at (0, 0, 0) when μ = 0, and the associated characteristic equation of system (2) with μ = 0 has a pair of simple imaginary roots ±iω_{0}.

By the Riesz Representation Theorem, there exists a function η(θ, μ) of bounded variation for θ ∈ [−τ_{1}, 0], such that

L_{μ}φ = ∫_{−τ_{1}}^{0} dη(θ, μ)φ(θ), for φ ∈ C.  (14)

In fact, we can choose

η(θ, μ) = \begin{cases}A + B_{1} + B_{2}, θ = 0, \\ B_{1} + B_{2},  θ ∈ (−τ_{2}, 0), \\ B_{1},  θ ∈ (−τ_{1}, −τ_{2}], \\ 0,  θ = −τ_{1}. \end{cases}

For φ ∈ C^{1}([−τ_{1}, 0], R^{3}), we set

A(μ)φ = \{ \begin{matrix}\frac{dφ(θ)}{dθ}, θ ∈ [−τ_{1}, 0), \\ ∫_{−τ_{1}}^{0} dη(ξ, μ)φ(ξ), θ = 0, \end{matrix}

and

R(μ)φ = \{ \begin{matrix}0,  θ ∈ [−τ_{1}, 0), \\ f(μ, φ),  θ = 0. \end{matrix}

Then system (12) can be rewritten as

u̇_{t} = A(μ)u_{t} + R(μ)u_{t},  (15)

where

u_{t}(θ) = u(t + θ)

for θ ∈ [−τ_{1}, 0]. For ψ ∈ C^{1}([0, τ_{1}], (R^{3})^{∗}), define

A^{∗}ψ(s) = \{ \begin{matrix}−\frac{dψ(s)}{ds}, s ∈ (0, τ_{1}], \\ ∫_{−τ_{1}}^{0} ψ(−t)dη(t, 0), s = 0. \end{matrix}

For φ ∈ C^{1}([−τ_{1}, 0], R^{3}) and ψ ∈ C^{1}([0, τ_{1}], (R^{3})^{∗}), using the bilinear form

< ψ, φ >= ψ̄(0)φ(0) − ∫_{−τ_{1}}^{0} ∫_{ξ=0}^{θ} ψ̄(ξ − θ)dη(θ)φ(ξ)dξ,  (16)

where η(θ) = η(θ, 0), we know that A^{∗} and A = A(0) are adjoint operators. By the discussion in Section 2, we know that ±iω_{0} are eigenvalues of A(0). Thus they are eigenvalues of A^{∗}.

We know that the eigenvector of A(0) corresponding to the eigenvalue iω_{0} denoted by q(θ) satisfies A(0)q(θ) = iω_{0}q(θ). By the definition of A(0) we obtain

q(θ) = (1, \frac{1}{iω_{0} − β_{1}} , \frac{z_{0}}{iω_{0} − x_{0} + γ})e^{iω_{0}θ}  (17)

Similarly, it can be verified that

q^{∗}(s) = D(1, \frac{1}{iω_{0} + β_{1}} , \frac{1}{iω_{0} + x_{0} − γ})e^{iω_{0}s}

is a eigenvector of A^{∗} corresponding to −iω_{0}, where D is a constant such that < q^{∗}(s), q(θ) >= 1. Denote

a = \frac{1}{iω_{0} − β_{1}} , b = \frac{z_{0}}{iω_{0} − x_{0} + γ} , a^{∗} = \frac{1}{iω_{0} + β_{1}} , and b^{∗} = \frac{1}{iω_{0} + x_{0} − γ} .

By (16) we have

< q^{∗}(s), q(θ) >= \overline{D}(1, \overline{a^{∗}}, \overline{b^{∗}})(1, a, b)^{T} − ∫_{−τ_{1}}^{0} ∫_{ξ=0}^{θ}(1, \overline{a^{∗}}, \overline{b^{∗}})e^{−i(ξ−θ)ω_{0}}dη(θ)(1, a, b)^{T} e^{iξω_{0}}dξ = \overline{D}[1 + a\overline{a^{∗}} + b\overline{b^{∗}} − ∫_{−τ_{1}}^{0}(1, \overline{a^{∗}}, \overline{b^{∗}})θe^{iθω_{0}}dη(θ)(1, a, b)^{T}]

= \overline{D}[1 + a\overline{a^{∗}} + b\overline{b^{∗}} + α_{1}τ_{1}e^{−iω_{0}τ_{1}} + α_{2}τ_{2}^{0}e^{−iω_{0}τ_{2}^{0}}].

Hence, we can choose

D = \frac{1}{1 + a^{∗}ā + b^{∗}b̄ + α_{1}τ_{1}e^{iω_{0}τ_{1}} + α_{2}τ_{2}^{0}e^{iω_{0}τ_{2}^{0}}}

so that < q^{∗}, q >= 1. Clearly, < q^{∗}, q̄ >= 0.

Following the algorithms in Hassard et al. [1981] to describe the center manifold C_{0} at μ = 0. Let u_{t} be the solution of Eq.(15) when μ = 0. Define

z(t) =< q^{∗}(s), u_{t}(θ) >, W (t, θ) = u_{t}(θ) − 2Re{z(t)q(θ)}.  (18)

On the center manifold C_{0} we have

W (t, θ) = W (z(t), \overline{z}(t), θ),

where

W (z, \overline{z}, θ) = W_{20}(θ)\frac{z^{2}}{2} + W_{11}(θ)z\overline{z} + W_{02}(θ)\frac{\overline{z}^{2}}{2} + ⋯,

z and \overline{z} are local coordinates for center manifold C_{0} in the direction of q^{∗} and \overline{q^{∗}}. Note that W is real if u_{t} is real. We consider only real solution.

For solution u_{t} ∈ C_{0} of (13), since μ = 0,

ż(t) = iω_{0}z+ < q^{∗}(θ), f(0, W (z, \overline{z}, θ) + 2Re\{zq(θ)\}) > = iω_{0}z + \overline{q^{∗}}f(0, W (z, \overline{z}, 0) + 2Re\{zq(0)\})

=^{def} iω_{0}z + \overline{q^{∗}}(0)f_{0}(z, \overline{z}).

We rewrite this equation as

ż(t) = iω_{0}z(t) + g(z, \overline{z}),  (19)

where

g(z, \overline{z}) = \overline{q^{∗}}(0)f_{0}(z, \overline{z}) = g_{20} \frac{z^{2}}{2} + g_{11}z\overline{z} + g_{02} \frac{\overline{z}^{2}}{2} + g_{21} \frac{z^{2}\overline{z}}{2} + · · · .  (20)

By (18), we have

u_{t}(θ) = (u_{1t}, u_{2t}, u_{3t}) = W (t, θ) + zq(θ) + \overline{zq}(θ),

and then

\begin{matrix}u_{1t}(0) = z + \overline{z} + W_{20}^{(1)}(0)\frac{z^{2}}{2} + W_{11}^{(1)}(0)z\overline{z} + W_{02}^{(1)}(0)\frac{\overline{z}^{2}}{2} + ⋯, \\ u_{3t}(0) = bz + \overline{bz} + W_{20}^{(3)}(0)\frac{z^{2}}{2} + W_{11}^{(3)}(0)z\overline{z} + W_{02}^{(3)}(0)\frac{\overline{z}^{2}}{2} + · · · . \end{matrix}

From(20), we have

g(z, \overline{z}) = q̄^{∗}(0)f_{0}(z, \overline{z}) = D̄b̄^{∗}u_{1t}(0)u_{3t}(0)

= D̄b̄^{∗}(z + \overline{z} + W_{20}^{(1)}(0)\frac{z^{2}}{2} + W_{11}^{(1)}(0)z\overline{z} + W_{02}^{(1)}(0)\frac{\overline{z}^{2}}{2} + ⋯)(bz + \overline{bz} + W_{20}^{(3)}(0)\frac{z^{2}}{2}

+ W_{11}^{(3)}(0)z\overline{z} + W_{02}^{(3)}(0)\frac{\overline{z}^{2}}{2} + ⋯).

(21)

Comparing the coefficients with (20), we have

g_{20} = 2D̄b̄^{∗}b,
g_{11} = D̄b̄^{∗}(\overline{b} + b),
g_{02} = 2D̄b̄^{∗}\overline{b},
g_{21} = 2\overline{D}b̄^{∗}(\frac{1}{2}W_{20}^{(1)}(0)\overline{b} + \frac{1}{2}W_{20}^{(3)}(0) + W_{11}^{(3)}(0) + bW_{11}^{(1)}(0)).

We still need to compute W_{20}(θ) and W_{11}(θ). From (15) and (18), we have

\begin{matrix}Ẇ = u̇_{t} − żq − \overline{zq} = \{ \begin{matrix}AW − 2Re\{\overline{q^{∗}}(0)f_{0}q(θ)\},  θ ∈ [−τ_{1}^{∗}, 0), \\ AW − 2Re\{\overline{q^{∗}}(0)f_{0}q(0)\} + f_{0}, θ = 0, \end{matrix} \\ =^{def} AW + H(z, \overline{z}, θ), \end{matrix}  (22)

where

H(z, \overline{z}, θ) = H_{20}(θ)\frac{z^{2}}{2} + H_{11}(θ)z\overline{z} + H_{02}(θ)\frac{\overline{z}^{2}}{2} + ⋯,  (23)
A(0)W (t, θ) = \frac{1}{2}A(0)W_{20}(θ)z^{2} + A(0)W_{11}(θ)z\overline{z} + ⋯,  (24)

and

Ẇ = W_{z}ż + W\overline{z} ˙\overline{z} = (W_{20}(θ)z + W_{11}(θ)\overline{z} + W_{02}(θ)\frac{\overline{z}^{2}}{2} + ⋯)(iω_{0}z + g(z, \overline{z}))
+ (W_{02}(θ)\overline{z} + W_{11}(θ)z + ⋯)(iω_{0}\overline{z} + \overline{g}(z, \overline{z}))
= 2iω_{0}W_{20}(θ)z\overline{z} + ⋯,
A(0)W (t, θ) − Ẇ = (A(0) − 2iω_{0})W_{20}(θ)\frac{\overline{z}^{2}}{2} + A(0)W_{11}()z\overline{z} + · · · .

Hence

(A(0)−2iω_{0})W_{20}(θ)\frac{\overline{z}^{2}}{2} +A(0)W_{11}(θ)z\overline{z}+⋯ = −H_{20}(θ)\frac{z^{2}}{2} −H_{11}(θ)z\overline{z}−H_{02}(θ)\frac{\overline{z^{2}}}{2} +· · · .

(25)

Comparing the coefficients, we obtain

(A − 2iω_{0}I)W_{20}(θ) = −H_{20}(θ), AW_{11}(θ) = −H_{11}(θ), · · · .  (26)

By (22), we know that for θ ∈ [−τ_{1}, 0),

H(z, \overline{z}, θ) = −\overline{q^{∗}}(0)f_{0}q(θ) − q^{∗}(0)\overline{f_{0}}\overline{q}(θ) = −gq(θ) − \overline{gq}(θ).  (27)

Comparing the coefficients with (23) gives that

\begin{matrix}H_{20}(θ) = −g_{20}q(θ) − \overline{g}_{02}\overline{q}(θ), \\ H_{11}(θ) = −g_{11}q(θ) − \overline{g}_{11}\overline{q}(θ). \end{matrix}  (28)

From (26), (27) and the definition of A, we obtain

Ẇ_{20} = 2iω_{0}W_{20}(θ) + g_{20}q(θ) + \overline{g}_{02}\overline{q}(θ).

Solving for W_{20}(θ), we obtain

W_{20}(θ) = \frac{ig_{20}}{ω_{0}} q(0)e^{iω_{0}θ} + \frac{i\overline{g}_{02}}{3ω_{0}} \overline{q}(0)e^{−iω_{0}θ} + E_{1}e^{2iω_{0}θ},

and similarly

W_{11}(θ) = −\frac{ig_{11}}{ω_{0}} q(0)e^{iω_{0}θ} + \frac{i\overline{g}_{11}}{ω_{0}} \overline{q}(0)e^{−iω_{0}θ} + E_{2},

where E_{1} and E_{2} are both three-dimensional vectors, and can be determined by setting θ = 0 in H. In fact, since

H(z, \overline{z}, 0) = −2Re\{\overline{q^{∗}}(0)f_{0}q(0)\} + f_{0},

we have

H_{20}(0) = −g_{20}q(0) − \overline{g}_{02}\overline{q}(0) + f_{z^{2}},

and

H_{11}(0) = −g_{11}q(0) − \overline{g}_{11}\overline{q}(0) + f_{z\overline{z}},

where

f_{0} = f_{z^{2}} \frac{z^{2}}{2} + f_{z\overline{z}}z\overline{z} + f\overline{z}^{2} \frac{\overline{z}^{2}}{2} + · · · .

Hence, combining with the definition of A, we obtain

∫_{−τ_{1}}^{0} dη(θ)W_{20}(θ) = 2iω_{0}W_{20}(0) + g_{20}q(0) + \overline{g}_{02}\overline{q}(0) − f_{z^{2}},

and

∫_{−τ_{1}}^{0} dη(θ)W_{11}(θ) = g_{11}q(0) + \overline{g}_{11}\overline{q}(0) − f_{z\overline{z}}.

Notice that

(iω_{0}I − ∫_{−τ_{1}}^{0} e^{iω_{0}θ}dη(θ)) q(0) = 0, (−iω_{0}I − ∫_{−τ_{1}}^{0} e^{−iω_{0}θ}dη(θ)) \overline{q}(0) = 0

we have

(2iω_{0}I − ∫_{−τ_{1}}^{0} e^{2iω_{0}θ}dη(θ)) E_{1} = f_{z^{2}}.

Similarly, we have

− (∫_{−τ_{1}}^{0} dη(θ)) E_{2} = f_{z\overline{z}}.

Hence, we get

( \begin{matrix}2iω_{0} − α_{1}e^{−2iω_{0}τ_{1}} − α_{2}e^{−2iω_{0}τ_{2}^{0}} & 1 & 1 \\ −1 & 2iω_{0} − β_{1} & 0 \\ −z_{0} & 0 & 2iω_{0} − x_{0} + γ \end{matrix} ) E_{1} = ( \begin{matrix}0 \\ 0 \\ b \end{matrix} )

and

( \begin{matrix}−α_{1} − α_{2} & 1 & 1 \\ −1 & −β_{1} & 0 \\ −z_{0} & 0 & −x_{0} + γ \end{matrix} ) E_{2} = ( \begin{matrix}0 \\ 0 \\ b + \overline{b} \end{matrix} ) .

Then g_{21} can be expressed by the parameters.

Based on the above analysis, we can see that each g_{ij} can be determined by the parameters. Thus we can compute the following quantities:

C_{1}(0) = \frac{i}{2ω_{0}} (g_{11}g_{20} − 2|g_{11}|^{2} − \frac{|g_{02}|^{2}}{3} ) + \frac{g_{21}}{2} ,

μ_{2} = −\frac{Re(C_{1}(0))}{Re(λ′(τ_{2}^{0}))},

β_{2} = 2Re(C_{1}(0)),

T_{2} = −\frac{Im(C_{1}(0)) + μ_{2}Im(λ′(τ_{2}^{0}))}{ω_{0}}.

Theorem 3.1. For system (2)

(i) μ_{2} determines the direction of the Hopf bifurcation: if μ_{2} > 0 (resp.μ_{2} < 0), then the bifurcating periodic solutions exist for τ_{2} in a right hand side neighborhood (τ_{2}^{0}, τ_{2}^{0} + ϵ) (resp. left hand side neighborhood (τ_{2}^{0} − ϵ, τ_{2}^{0})) of the bifurcation value τ_{2}^{0}.

(ii) β_{2} determines the stability of bifurcating periodic solutions: the bifurcating periodic solutions are orbitally asymptotically stable (resp. unstable) if β_{2} < 0 (resp.β_{2} > 0);

(iii) T_{2} determines the period of the bifurcating periodic solutions: the period increases (resp. decreases) if T_{2} > 0 (resp. T_{2} < 0).

4 Numerical examples

In this section, we shall use MATLAB to perform some numerical simulations on system (2).

In the following, we choose a set of parameters as follows::

(a) α_{1} = 0.2, α_{2} = −0.2, β_{1} = 0.2, β_{2} = 0.2, γ = 5.7.

With these parameters, one can find that (H_{2}) is satisfied. When τ_{2} = 0, by a direct computation we obtain that Eq.(7) has two positive roots ν_{1} ≐ 0.8327 and ν_{2} ≐ 1.1506. Substituting them and the data (a) into Eq.(6) gives, respectively,

τ_{11}^{(j)} ≐ 2.2945 + 7.5456j (j = 0, 1, 2, ⋯), τ_{12}^{(j)} ≐ −1.1036 + 5.4608j (j = 1, 2, ⋯).

Eq.(4) has pure imaginary roots when τ_{1} = τ_{11}^{(j)} or τ_{1} = τ_{12}^{(j)}. Further α′(τ_{11}^{(j)}) < 0 and α′(τ_{12}^{(j)}) > 0. By Theorem 2.4 we know that the stability switches exist, that is the equilibrium is unstable when τ_{1} ∈ [0, 2.2945), and asymptotically stable when τ_{1} ∈ (2.2945, 4.3572). The results are illustrated in Fig.1-Fig.4.

Fig.1. The equilibrium is unstable, and chaos exists for system (2) with the data (a) and τ_{1} = τ_{2} = 0.

Fig.2. The equilibrium is unstable and chaos phenomenon still exists for system (2) with the data (a), τ_{1} ∈ [0, 2.2945) and τ_{2} = 0, where τ_{1} = 1.

Fig.3. The equilibrium becomes stable and the chaos phenomenon disappears for system (2) with the data (a), τ_{1} ∈ (2.2945, 4.3572) and τ_{2} = 0, where τ_{1} = 3.

Fig.4. The equilibrium is unstable, and a bifurcating periodic solution appears for system (2) with the data (a) and τ_{1} > 4.3572 and close to 4.3572, and τ_{2} = 0, where τ_{1} = 5.

Let τ_{1} = 3.0 ∈ (2.2945, 4.3572), we obtain τ_{2}^{0} ≐ 1.6295. By the Theorem 2.4 we know that (x_{0}, y_{0}, z_{0}) is asymptotically stable for τ_{1} = 3 and τ_{2} ∈ [0, 1.6295). Furthermore by direct computation using the algorithm derived in Section 3, we have C_{1}(0) ≐ −0.0003 + 0.0002i, β_{2} ≐ −0.0006 < 0, and μ_{2} > 0. We know that, at τ_{2}^{0} ≐ 1.6295, the bifurcating periodic solution is orbitally asymptotically stable, and the direction of the Hopf bifurcation is forward, which are illustrated in Fig.5-Fig.6. On the other hand, the numerical simulations show that the bifurcating periodic solutions disappear when the delay τ_{2} is far away τ_{2}^{0} = 1.6295, and chaos occurs again. This is shown in Fig.7.

Fig.5. The equilibrium is asymptotically stable for system (2) with the data (a), τ_{1} ∈ (2.2945, 4.3572) and τ_{2} ∈ [0, τ_{2}^{0}), where τ_{1} = 3, and τ_{2} = 1.

Fig.6. The equilibrium is unstable, and a bifurcating periodic solution appears for system (2) with the data (a), τ_{1} ∈ (2.2945, 4.3572) and τ_{2} > τ_{2}^{0} is close to τ_{2}^{0}, where τ_{1} = 3 and τ_{2} = 2 > 1.6295.

Fig.7. Chaos occurs again for system (2) with the data (a), τ_{1} ∈ (2.2945, 4.3572) and τ_{2} > τ_{2}^{0} increasing further, where τ_{1} = 3, τ_{2} = 3.5 > 1.6295.

5 Conclusion

Bifurcation in Rössler system with single delay has been observed by many researchers. However, there are few papers on the bifurcation of Rössler system with multiple delays.

In this paper we have analyzed the Rössler system with multiple delays on two different conditions. We find out that there are stability switches for the interior equilibrium when τ_{1} varies in the case of τ_{2} = 0. Then for τ_{1} in a stability interval, regarding the delay τ_{2} as parameter, we show that there exists a first critical value of τ_{2} at which the interior equilibrium loses its stability and the Hopf bifurcation occurs. We also investigate the direction of the Hopf bifurcation and the stability of the bifurcating periodic solutions, by using the center manifold theory and normal form method.

Our theoretical results and numerical simulations show that, for a Rössler system with chaos phenomena, the chaos oscillation can be controlled by delays. For example, the multiple delayed Rössler system we studied possess chaos oscillation when τ_{1} = τ_{2} = 0. The chaos disappears when the delays increase, and the stability of the equilibrium is lost at same time, and the periodic solutions occur from Hopf bifurcation. As the delays increasing further, the numerical simulations show that the periodic solution disappears and the chaos oscillation appears again. Acknowledgement We would like to express our gratitude to the referee for his or her valuable comments and suggestions that led to a truly significant improvement of the manuscript.

References

[1] Ahmed, E., El-Sayed, A. M. A., El-Saka, H. A. A., On some Routh-Hurwitz conditions for fractional order differential equations and their applications in Lorenz, Rössler, Chua and Chen systems, Phys. Lett. A. 2006, 358: 1-4.

[2] Chang, J-F., Hung, M-L. and Yang, Y-S. et al., Controlling chaos of the family of Rössler system using sliding mode control, Chaos, Solitons and Fractals. 2006, 37: 609-622.

[3] Chen, Q., Zhang, Y. and Hong, Y., Generation and control of striped attractors of Rössler system with feedback, Chaos, Solitons and Fractals. 2007, 34: 693-703.

[4] Ghosh, D., Roy Chowdhury, A. and Saha, P., Multiple delay Rössler systembifurcation and chaos control, Chaos, Solitions and Fractals. 2008, 35: 472-485.

[5] Guan, X., Chen, C. and Peng, H. et al., Time-delay feedback control of timedelay chaotic system, Int. J. Bifur. Chaos. 2003, 13: 193-205.

[6] Hale J., Theory of Functional Differential Equations. Springer, New York. (1977).

[7] Hassard, B., Kazarinoff, N. and Wan, Y., Theory and Application of Hopf Bifurcation, Cambridge: Cambridge University Press, 1981.

[8] Lorenz, E., Deterministic non-periodic flows, J. Atmos. Sci. 1963, 20: 130-41.

[9] Ott, E., Grebogi, C. and Yorke, J., Controlling chaos, Phys. Rev. Lett. 1990, 64:1196-1199.

[10] Pyragas, K., Continuous control of chaos by self-controlling feedback, Phys. Lett. A 1992, 170: 421-428.

[11] Pyragas, K., Experimental control of chaos by delayed self-controlling feedback, Phys. Lett. A. 1993, 180: 99-102.

[12] Rasmussen, J., Mosekilde, E. and Reick, C., Bifurcation in two coupled Rössler system, Mathematics and Computers in Simulation. 1996, 40: 247-270.

[13] Rössler, E., An equation for continuous chaos, Phys. Lett. A. 1976, 57: 397-398.

[14] Ruan, S. and Wei, J., On the zeros of a third degree exponential polynomial with applications to a delayed model for the control of testosterone secretion, IMA. J. Math. Appl. Med. Biol. 2001, 18: 41-52.

[15] Ruan, S. and Wei, J., On the zeros of transcendental functions with applications to stability of delay differential equations with two delays, Dynamics of Continuous, Discrete and Impulsive Systems. 2003, 10: 863-874.

[16] Shu, Y., Tan, P. and Li, C., Control of n-dimensional continuous-time system with delay, Phys. Lett. A. 2004, 323: 251-259.

[17] Starkov, K. E. and Starkov Jr. K. K., Localization of periodic orbits of the Rössler system under variation of its parameters, Chaos, Solitons and Fractals. 2007, 33: 1445-1449.

[18] Tao, C., Yang, C. and Luo, Y. et al., Speed feedback control of chaotic system Chaos, Solitons and Fractals. 2005, 23: 259-263.

[19] Tian, Y., Tade, M. and Tang, J., Nonlinear open-plus-closed-loop (NOPCL) control of dynamic systems, Chaos, Solitons and Fractals. 2000, 11: 1029-1035.

[20] Yang, T., Yang, C. and Yang, L., Control of Rössler system to periodic motions using impulsive control methods, Phys. Lett. A. 1997, 232: 356-361.

[21] Zhang, Y. and Sun, J., Delay-dependent stability criteria for coupled chaotic systems via unidirectional linear feedback approach, Chaos, Solitons and Fractals. 2004, 22: 199-205.

(Received February 10, 2010)