- Email: [email protected]
Bifurcation and chaos in a ratio-dependent predator–prey system with time delay
Chaos, Solitons and Fractals 39 (2009) 1883–1895 www.elsevier.com/locate/chaos
Bifurcation and chaos in a ratio-dependent predator–prey system with time delay Qintao Gan a
a,*
, Rui Xu
a,b
, Pinghua Yang
a
Institute of Applied Mathematics, Shijiazhuang Mechanical Engineering College, Shijiazhuang 050003, PR China b Department of Applied Mathematics, Xi’an Jiaotong University, Xi’an 710049, PR China Accepted 20 June 2007
Abstract In this paper, a ratio-dependent predator–prey model with time delay is investigated. We first consider the local stability of a positive equilibrium and the existence of Hopf bifurcations. By using the normal form theory and center manifold reduction, we derive explicit formulae which determine the stability, direction and other properties of bifurcating periodic solutions. Finally, we consider the effect of impulses on the dynamics of the above time-delayed population model. Numerical simulations show that the system with constant periodic impulsive perturbations admits rich complex dynamic, such as periodic doubling cascade and chaos. 2007 Elsevier Ltd. All rights reserved.
1. Introduction In population dynamics, a functional response of the predator to the prey density refers to the change in the density of prey per unit time per predator as a function of the prey density. Standard Lotka–Volterra type models, on which a large body of existing predator–prey theory is built, assume that the per capita rate of predation depends solely on the prey numbers only (it is usually called prey-dependent functional response). Predator–prey systems with prey-dependent response have been studied extensively and the dynamics of such systems are now very well understood (see, for example, Kuang [20] and the references cited therein). Recently, the traditional prey-dependent predator–prey models have been challenged by several biologists (see, for example, Arditi and Ginzburg [1], Arditi et al. [2] and Gutierrez [13]) based on the fact that functional and numerical responses over typical ecological timescales ought to depend on the densities of both prey and predators (most likely and simply on their ratio), especially when predators have to search for food (and therefore have to share or compete for food). Such a functional response is called a ratio-dependent response function. These hypotheses are strongly supported by numerous field and laboratory experiment and observations [1–4,15]. Based on the Michaelis–Menten or Holling type II function, Arditi and Ginzburg [1] proposed a ratio-dependent function of the form
*
Corresponding author. E-mail address: [email protected] (Q. Gan).
0960-0779/$ - see front matter 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2007.06.122
1884
Q. Gan et al. / Chaos, Solitons and Fractals 39 (2009) 1883–1895
P
x cx=y cx ¼ ¼ ; y m þ x=y my þ x
and the following ratio-dependent predator–prey model: x_ ¼ xða bxÞ cxy=ðmy þ xÞ; y_ ¼ yðd þ fx=ðmy þ xÞÞ:
ð1:1Þ
Here, x(t) and y(t) represent population densities of the prey and the predator at time t, respectively; parameters a; b; c; d; f ; m are positive constants in which a/b is the carrying capacity of the prey, d > 0 is the death rate of the predator, and a; c; m; and f represent the prey intrinsic growth rate, capturing rate, half saturation constant and conversion rate, respectively. The ratio-dependent predator–prey model (1.1) has been studied by several researchers recently and very rich dynamics have been observed (see, for example, [7,8,11,15,17,18,21,39,40]). On the other hand, time delays of one type or another have been incorporated into biological models by many researchers, we refer to the monographs of Gopalsamy [17], Kuang [21] and references cited therein for general delayed biological systems and for studies on delayed predator–prey systems. In general, delay differential equations exhibit much more complicated dynamics than ordinary differential equations since a time delay could cause a stable equilibrium to become unstable and cause the populations to fluctuate (see, for example, [9,19,24–26,30–35,37,38,41]). In this paper, we are concerned with the following ratio-dependent predator–prey system with time delay xðt sÞ a1 xðtÞyðtÞ x_ ðtÞ ¼ r1 xðtÞ 1 ; k myðtÞ þ xðtÞ ð1:2Þ a2 xðtÞ y_ ðtÞ ¼ yðtÞ r2 þ ; myðtÞ þ xðtÞ where x(t), y(t) denote the density of the prey and the predator populations at time t, respectively; the parameter s P 0 is a time delay due to regulation of the prey. Parameters a1, a2, r1, r2, k and m are positive constants whose ecological meanings were mentioned above for model (1.1). The initial conditions for system (1.2) take the form xðhÞ ¼ /ðhÞ;
yðhÞ ¼ wðhÞ;
/ðhÞ P 0;
wðhÞ P 0;
h 2 ½s; 0;
/ð0Þ > 0;
wð0Þ > 0;
ð1:3Þ
Cð½s; 0; R2þ0 Þ,
the Banach space of continuous functions mapping the interval [s, 0] into R2þ0 , where ð/ðhÞ; wðhÞÞ 2 2 where Rþ0 ¼ fðx1 ; x2 Þ : x1 P 0; x2 P 0g. It is well-known by the fundamental theory of functional differential equations [14], system (1.2) has a unique solution (x(t), y(t)) satisfying initial conditions (1.3). It is easy to show that all solutions of system (1.2) with initial conditions (1.3) are defined on [0, + 1) and remain positive for all t P 0. We note that impulsive differential equations are suitable for the mathematical simulation of evolutionary process in which the parameters undergo relatively long periods of smooth variation followed by a short-term rapid change (i.e., jumps) in their values. Recently, equations with impulsive effect have been found in almost every domain of applied science. Numerous examples are given in Bainovs and his collaborator’s book [5,23]. Some impulsive differential equations have recently been introduced in population dynamics in relation to impulsive birth [28,36], impulsive vaccination [10,29], chemotherapeutic treatment of disease [22,27], and population ecology [6]. Motivated by the work above, in this paper, we further discuss the effect of impulses on the dynamics of system (1.2). To this end, we discuss the following impulsive equations 8 9 a1 xðtÞyðtÞ > > > = x_ ðtÞ ¼ r1 xðtÞ 1 xðtsÞ ; > k myðtÞþxðtÞ > > t–nT ; > < a2 xðtÞ > ; y_ ðtÞ ¼ yðtÞ r2 þ myðtÞþxðtÞ ; ð1:4Þ > > þ > > DxðtÞ ¼ xðt Þ xðt Þ ¼ 0; > > t ¼ nT : : DyðtÞ ¼ yðtþ Þ yðt Þ ¼ fyðtþ Þ; The organization of this paper is as follows. In Section 2, we discuss the local stability of a positive equilibrium of system (1.2) and the existence of Hopf bifurcations. In Section 3, we study the stability of the bifurcating periodic solutions and the direction of the Hopf bifurcation in system (1.2) by using the normal form theory and center manifold reduction. In Section 4, we numerically study the effect of impulsive perturbation on the dynamics of system (1.2). Numerical simulations show that system (1.4) admits rich complex dynamics such as periodic doubling cascade and chaos.
Q. Gan et al. / Chaos, Solitons and Fractals 39 (2009) 1883–1895
1885
2. Stability of a positive equilibrium and the existence of Hopf bifurcations In this section, we discuss the local stability of a positive equilibrium and the existence of Hopf bifurcations in system (1.2). It is easy to show that system (1.2) has a unique positive equilibrium E* = (x*, y*) if the following holds: (H1) a2 > r2 and ma2r1 > a1(a2 r2), where x ¼ k Set
ka1 ða2 r2 Þ ; ma2 r1
y ¼
kða2 r2 Þ ka1 ða2 r2 Þ2 : m2 a2 r1 r2 mr2
xðt sÞ a1 xðtÞyðtÞ f ð1Þ ¼ r1 xðtÞ 1 ; k myðtÞ þ xðtÞ a2 xðtÞ f ð2Þ ¼ yðtÞ r2 þ ; myðtÞ þ xðtÞ a1 r2 ða2 r2 Þ a1 r22 ; p2 ¼ ; 2 a22 ma2 r2 ðr2 a2 Þ r 1 x p4 ¼ ; p5 ¼ : k a2 p1 ¼
p3 ¼
ða2 r2 Þ2 ; ma2
ð2:1Þ
Linearizing system (1.2) at E*, we derive that x_ ðtÞ ¼ p1 xðtÞ þ p2 yðtÞ þ p5 xðt sÞ; _ ¼ p3 xðtÞ þ p4 yðtÞ; yðtÞ
ð2:2Þ
where pi (i = 1, 2, 3, 4, 5) are defined in (2.1). Hence, the characteristic equation of system (1.2) at the positive equilibrium E* takes the form k2 ðp1 þ p4 Þk p5 ðk p4 Þeks ¼ 0:
ð2:3Þ
It is well-known that the zero steady state of system (2.2) is asymptotically stable if all roots of Eq. (2.3) have negative real parts, and is unstable if Eq. (2.3) has a root with positive real part. In the sequel, we shall investigate the distribution of roots of Eq. (2.3). If ix(x > 0) is a root of Eq. (2.3), then x2 iðp1 þ p4 Þx p5 ðix p4 Þeixs ¼ 0: Separating the real and imaginary parts, we obtain x2 ¼ p5 x sin xs p4 p5 cos xs; ðp1 þ p4 Þx ¼ p5 x cos xs þ p4 p5 sin xs:
ð2:4Þ
It follows from (2.4) that x4 þ ½ðp1 þ p4 Þ2 p25 x2 p24 p25 ¼ 0:
ð2:5Þ
Clearly, Eq. (2.5) has only one positive root x0 defined by " ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi!#12 r 2 1 2 2 2 2 2 2 p5 ðp1 þ p4 Þ þ 4p4 p5 : x0 ¼ p5 ðp1 þ p4 Þ þ 2 Define sj ¼
1 x0 ðx20 þ ðp1 þ p4 Þp4 Þ þ 2jp ; arcsin x0 p5 ðp24 þ x20 Þ
j ¼ 0; 1; . . . :
Then (sj, x0) solves Eq. (2.4). This means that when s = sj, Eq. (2.3) has a pair of purely imaginary roots ± ix0. Now let us consider the behavior of roots of Eq. (2.3) near sj. Denote kðsÞ ¼ aðsÞ þ ixðsÞ
ð2:6Þ
ð2:7Þ
1886
Q. Gan et al. / Chaos, Solitons and Fractals 39 (2009) 1883–1895
the root of Eq. (2.3) such that aðsj Þ ¼ 0;
xðsj Þ ¼ x0 :
Substituting k(s) into Eq. (2.3) and differentiating both sides of it with respect to s, we have 1 dkðsÞ ð2k ðp1 þ p4 ÞÞeks p5 s þ : ¼ p5 kðk p4 Þ p5 kðk p4 Þ k ds We therefore derive that 1 dk ð2k ðp1 þ p4 ÞÞeks p5 s ¼ Re þ Re Re p5 kðk p4 Þ p5 kðk p4 Þ k s¼sj ds s¼sj s¼sj ) (
1 ¼ Re ðp1 þ p4 Þ cos x0 sj 2x0 sin x0 sj :þi½2x0 cos x0 sj ðp1 þ p4 Þ sin x0 sj p5 x20 þ ip4 p5 x0 p5 þ Re p5 x20 þ ip4 p5 x0 1 ¼ fp5 x20 ½ðp1 þ p4 Þ cos x0 sj þ 2x0 sin x0 sj þp4 p5 x0 ½2x0 cos x0 sj ðp1 þ p4 Þ sin x0 sj p25 x20 C 1 ¼ fðp1 þ p4 Þx0 ½p5 x0 cos x0 sj þ p4 p5 sin x0 sj 2x20 ½p5 x0 sin x0 sj p4 p5 cos x0 sj p25 x20 C o x2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x2 n ¼ 0 ðp1 þ p4 Þ2 þ 2x20 p25 ¼ 0 ðp25 ðp1 þ p4 Þ2 Þ2 þ 4p24 p25 ; C C where we have used Eqs. (2.4) and (2.6), and C ¼ p25 x40 þ p24 p25 x20 > 0. Hence, it follows that ( ) ( ) 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 dk dk x ¼ sign Re ¼ sign 0 ðp25 ðp1 þ p4 Þ2 Þ2 þ 4p24 p25 > 0: sign Re C ds s¼sj ds s¼sj
ð2:8Þ
Therefore, when the delay s near sj is increased, the root of Eq. (2.3) crosses the imaginary axis from left to right. In addition, note that when s = 0, Eq. (2.3) has roots with negative real parts only if (H2) p1 + p4 + p5 < 0 and p4p5 > 0. Thus, by the well-known Rouche´ theorem we obtain the following results on the distribution of roots of Eq. (2.3). Lemma 2.1. Let sj (j = 0,1, . . .) be defined as in (2.7). Then all roots of Eq. (2.3) have negative real parts for all s 2 [0, s0). However, Eq. (2.3) has at least one root with positive real part when s > s0, and Eq. (2.3) has a pair of purely imaginary root ± ix0 when s = s0. More detail, for s 2 (sj,sj + 1] (j = 0, 1, . . .), Eq. (2.3) has 2(j + 1) roots with positive real parts. Moreover, all roots of Eq. (2.3) with s = sj (j = 0, 1, . . .) have negative real parts except ±ix0. Applying Lemma 2.1, the transversality condition (2.8) and Theorem 11.1 developed in [14], we obtain the following results. Theorem 2.1. Let (H1) and (H2) hold. let x0 and sj (j = 0, 1, . . .) be defined as in (2.6) and (2.7), respectively. (i) The positive equilibrium E* of system (1.2) is asymptotically stable for all s 2 [0, s0) and unstable for s > s0. (ii) System (1.2) undergoes a Hopf Bifurcation at the positive equilibrium E* when s = sj (j = 0, 1, . . .).
3. Direction of Hopf bifurcation In Section 2, we have shown that system (1.2) admits a series of periodic solutions bifurcating from the positive equilibrium E * at the critical values of s. In this section, we derive explicit formulae to determine the properties of the Hopf bifurcation at critical values sj by using the normal form theory and center manifold reduction (see, for example, Hassard et al. [16]). Without loss of generality, denote the critical values sj by ~s, and set s ¼ ~s þ l. Then l = 0 is a Hopf bifurcation value of system (1.2). Thus, we can work in the phase space C ¼ Cð½~s; 0; R2 Þ.
Q. Gan et al. / Chaos, Solitons and Fractals 39 (2009) 1883–1895
Let u1 ðtÞ ¼ xðtÞ x ;
1887
u2 ðtÞ ¼ yðtÞ y :
Then system (1.2) is transformed into 8 P ð1Þ i j 1 l > > < u_ 1 ðtÞ ¼ p1 u1 ðtÞ þ p2 u2 ðtÞ þ p5 u2 ðt ~sÞ þ iþjþlP2 i!j!l! fijl u1 ðtÞu2 ðtÞu1 ðt ~sÞ; P 1 ð2Þ i > > f u1 ðtÞuj2 ðtÞ; : u_ 2 ðtÞ ¼ p3 u1 ðtÞ þ p4 u2 ðtÞ þ i!j! ij
ð3:1Þ
iþjP2
where ð1Þ fijl
¼ ; oxi oy j oxðt ~sÞl ðx ;y ;x Þ oiþjþl f ð1Þ
ð2Þ
fij ¼
oiþj f ð2Þ ; oxi oy j ðx ;y Þ
i; j; l P 0;
here, f (1) and f (2) are defined in (2.1). For the simplicity of notations, we rewrite (3.1) as _ ¼ Ll ðut Þ þ f ðl; ut Þ; uðtÞ T
ð3:2Þ 2
where u(t) = (u1(t),u2(t)) 2 R , ut(h) 2 C is defined by ut(h) = u(t + h), and Ll:C ! R, f:R · C 2 R are given by p5 0 p1 p2 Ll / ¼ /ð0Þ þ /ð~sÞ; ð3:3Þ p3 p4 0 0 and
1 0 P ð1Þ 1 f /i1 ð0Þ/j2 ð0Þ/l1 ð~sÞ i!j!l! ijl C B iþjþlP2 C f ðl; /Þ ¼ B P 1 ð2Þ i A @ j f / ð0Þ/ ð0Þ ij 2 1 i!j!
ð3:4Þ
iþjP2
respectively. By Riesz representation theorem, there exists a function g(h,l) of bounded variation for h 2 ½~s; 0 such that Z 0 Ll / ¼ dgðh; 0Þ/ðhÞ for / 2 C: ð3:5Þ ~s
In fact, we can choose p1 p2 p5 gðh; lÞ ¼ dðhÞ 0 p3 p4
0 dðh þ ~sÞ; 0
ð3:6Þ
where d is the Dirac delta function. For / 2 C 1 ð½~s; 0; R2 Þ, define ( d/ðhÞ ; h 2 ½~s; 0Þ; AðlÞ/ ¼ Rdh 0 dgðl; sÞ/ðsÞ; h ¼ 0 ~s and
RðlÞ/ ¼
0;
h 2 ½~s; 0Þ;
f ðl; /Þ;
h ¼ 0:
Then system (3.2) is equivalent to u_ t ¼ AðlÞut þ RðlÞut ;
ð3:7Þ
where xt(h) = x(t + h) for h 2 ½~s; 0. For w 2 C 1 ð½0; ~s; ðR2 Þ Þ, define ( dwðsÞ ; s 2 ð0; ~s; A wðsÞ ¼ R 0 ds dgT ðt; 0ÞwðtÞ; s ¼ 0 ~s and a bilinear inner product hwðsÞ; /ðhÞi ¼ wð0Þ/ð0Þ
Z
0
~s
Z
h
n¼0
hÞdgðhÞ/ðnÞdn; wðn
ð3:8Þ
1888
Q. Gan et al. / Chaos, Solitons and Fractals 39 (2009) 1883–1895
where g(h) = g(h,0). Then A(0) and A* are adjoint operators. By discussions in Section 2 and foregoing assumption, we know that ±ix0 are eigenvalues of A(0). Thus, they are also eigenvalues of A*. We first need to compute the eigenvector of A(0) and A* corresponding to ix0 and ix0, respectively. Suppose that q(h) = (1,q)T eix0 h is the eigenvector of A(0) corresponding to ix0. Then A(0)q(h) = ix0q(h). It follows from the definition of A(0), (3.5) and (3.6) that 0 p1 þ p5 eix0 ~s ix0 p2 qð0Þ ¼ : 0 p3 p4 ix0 We therefore derive that qð0Þ ¼ ð1; qÞT ¼ 1;
p3 ix0 p4
T :
On the other hand, suppose that q*(s) = D(1,r)eix0 s is the eigenvector of A* corresponding to ix0. From the definition of A*, (3.5) and (3.6) we have 0 p1 þ p5 eix0 ~s þ ix0 p3 ðq ð0ÞÞT ¼ ; 0 p2 p4 þ ix0 which yields
q ð0Þ ¼ Dð1; rÞ ¼ D 1;
p2 : p4 þ ix0
In order to assure hq*(s), q(h)i = 1, we need to determine the value of D. From (3.7), we have Z 0Z h Þð1; qÞT ÞeiðnhÞx0 dgðhÞð1; qÞT einx0 dn ð1; r hq ðsÞ; qðhÞi ¼ D ð1; r r ¼ D 1 þ q
Z
~s
n¼0
0 ix0 h
Þhe ð1; r
T
dgðhÞð1; qÞ
¼ Df1 þ q r þ p5~seix0 ~s g:
~s
Thus, we can choose 1 D¼ r þ p5~seix0 ~s 1þq such that hq ðsÞ; qðhÞi ¼ 1; hq ðsÞ; qðhÞi ¼ 0: In the remainder of this section, we use the same notations as in Hassard et al. [16]. We first compute the coordinates to describe the center manifold C0 at l = 0. Let ut be the solution of Eq. (3.2) with l = 0. Define zðtÞ ¼ hq ; ut i; W ðt; hÞ ¼ ut ðhÞ 2RefzðtÞqðhÞg: ð3:9Þ On the center manifold C0 we have W ðt; hÞ ¼ W ðzðtÞ; zðtÞ; hÞ; where z2 z2 þ W 11 ðhÞzz þ W 02 ðhÞ þ ; 2 2 q . Note that W is real if ut is real. we z and z are local coordinates for center manifold C0 in the direction of q* and consider only real solutions. For the solution ut 2 C0 of (3.2), since l = 0, we have W ðz; z; hÞ ¼ W 20 ðhÞ
z_ ¼ ix0 z þ hq ðhÞ; f ð0; W ðz; z; hÞ þ 2RefzqðhÞgÞi ¼ ix0 z þ q ðhÞf ð0; W ðz; z; hÞ þ 2RefzqðhÞgÞ def
q ð0Þf0 ðz; zÞ: ¼ ix0 z þ q ð0Þf ð0; W ðz; z; 0Þ þ 2Refzqð0ÞgÞ ¼ ix0 z þ
ð3:10Þ
We rewrite (3.10) as z_ ¼ ix0 z þ gðz; zÞ with gðz; zÞ ¼ q ð0Þf0 ðz; zÞ ¼ g20
z2 z2 z2z þ g11 zz þ g02 þ g21 þ : 2 2 2
ð3:11Þ
Q. Gan et al. / Chaos, Solitons and Fractals 39 (2009) 1883–1895
1889
Noting that ut ðhÞ ¼ ðu1t ðhÞ; u2t ðhÞÞ ¼ W ðt; hÞ þ zqðhÞ þ z qðhÞ and q(h) = (1,q)T eix0 h , we have ð1Þ
u1t ð0Þ ¼ z þ z þ W 20 ð0Þ
z2 z2 ð1Þ ð1Þ þ W 11 ð0Þzz þ W 02 ð0Þ þ ; 2 2 ð1Þ
u1t ð~sÞ ¼ eix0 ~s z þ eix0 ~sz þ W 20 ð~sÞ ð2Þ
z þ W 20 ð0Þ u2t ð0Þ ¼ qz þ q
z2 z2 ð1Þ ð1Þ þ W 11 ð~sÞzz þ W 02 ð~sÞ þ ; 2 2
z2 z2 ð2Þ ð2Þ þ W 11 ð0Þzz þ W 02 ð0Þ þ : 2 2
Thus, it follows from (3.4) and (3.11) that gðz; zÞ ¼ q ð0Þf0 ðz; zÞ 0 P B iþjþlP2
ÞB ¼ Dð1; r @
1 ð1Þ 1 f ui ð0Þuj2t ð0Þul1t ð~sÞ i!j!l! ijl 1t
P
iþjP2
1 ð2Þ i f u1t ð0Þuj2t ð0Þ i!j! ij
C C A
1 ð1Þ 1 ð1Þ 1 1 ð1Þ ð1Þ f20ð2Þ þ f02ð2Þ q2 þ f11ð2Þ q z2 ¼ D f110 q þ f101 eix0 ~s þ f020 q2 þ f200 z2 þ r 2 2 2 2 ð1Þ ð1Þ ð1Þ ð1Þ þ 2f110 Refqg þ 2f101 Refeix0 ~s g þ f020 q q þ f200 zz 1 ð1Þ 2 1 ð1Þ 2 ð2Þ ð1Þ ð1Þ ix0 ~s Þ zz þ f110 þ f101 þ f200 z f20ð2Þ þ f02ð2Þ q þr q þ f11 ðq þ q e þ f020 q q 2 2 1 1 1 ð2Þ ð1Þ 1 ð1Þ ð1Þ 2 þ f11ð2Þ q z2 þ f110 f20ð2Þ þ f02ð2Þ q þr q þ W 20 ð0Þ þ qW 11 ð0Þ W ð0Þ 2 2 2 20 2 1 ð1Þ ð2Þ ð1Þ 1 ð1Þ ð1Þ ð1Þ þW 11 ð0Þ þ f101 W 20 ð0Þeix0 ~s þ W 20 ð~sÞ þ qeix0 ~s W 11 ð0Þ þ W 11 ð~sÞ 2 2 1 1 ð1Þ ð1Þ ð2Þ ð1Þ ð1Þ ð1Þ 1 W ð2Þ ð0Þ þ f W ð0Þ þ ð0Þ þ f030 q2 q W þ f020 qW 11 ð0Þ þ q 20 200 11 2 2 20 2 1 ð1Þ 1 1 ð1Þ ð1Þ 1 ð1Þ z2z þ r f20ð2Þ W ð1Þ q þ f210 q þ q þ f300 þ f120 q2 þ q 11 ð0Þ þ W 20 ð0Þ 2 2 2 2 1 1 ð2Þ 1 ð2Þ ð2Þ ð2Þ ð2Þ ð1Þ ð1Þ 2 W ð2Þ ð0Þ þ f W ð0Þ þ ð0Þ þ qW ð0Þ þ ð0Þ z z þ : W q W þf02 qW 11 ð0Þ þ q 20 11 11 11 20 2 2 20 2 Comparing the coefficients in (3.11), we get ð1Þ ð1Þ ð1Þ ð1Þ f20ð2Þ þ f02ð2Þ q2 þ 2f11ð2Þ q D; g20 ¼ 2f110 q þ 2f101 eix0 ~s þ f020 q2 þ f200 þ r ð1Þ ð1Þ ð1Þ ð1Þ ð2Þ Þ D; f20ð2Þ þ f02ð2Þ q g11 ¼ 2f110 Refqg þ 2f101 Refeix0 ~s g þ f020 q q þ f200 þ r q þ f11 ðq þ q ð1Þ ð1Þ ix0 ~s ð1Þ 2 ð1Þ þ 2f101 þ f200 2 þ 2f11ð2Þ q D; f20ð2Þ þ f02ð2Þ q e þ f020 q þr g02 ¼ 2f110 q ð1Þ ð1Þ ð2Þ ð1Þ ð2Þ q þ W 20 ð0Þ þ 2qW 11 ð0Þ þ 2W 11 ð0Þ g21 ¼ f110 W 20 ð0Þ ð1Þ ð1Þ ð1Þ ð1Þ ð1Þ þf101 W 20 ð0Þeix0 ~s þ W 20 ð~sÞ þ 2qeix0 ~s W 11 ð0Þ þ 2W 11 ð~sÞ ð1Þ ð2Þ ð1Þ ð1Þ ð1Þ ð1Þ 2 W ð2Þ þ f020 2qW 11 ð0Þ þ q 20 ð0Þ þ f200 2W 11 ð0Þ þ W 20 ð0Þ þ f030 q q
ð1Þ ð1Þ ð1Þ ð2Þ ð1Þ ð1Þ Þ þ r f20 2W 11 ð0Þ þ W 20 ð0Þ q þ f210 ð2q þ q þ f300 þ f120 q2 þ 2q ð2Þ ð2Þ ð2Þ ð2Þ ð2Þ ð1Þ W ð2Þ W ð1Þ 2W 11 ð0Þ þ W 20 ð0Þ þ 2qW 11 ð0Þ þ q D: þf02 2qW 11 ð0Þ þ q 20 ð0Þ þ f11 20 ð0Þ
ð3:12Þ
1890
Q. Gan et al. / Chaos, Solitons and Fractals 39 (2009) 1883–1895
We now compute W20(h) and W11(h). It follows from (3.7) and (3.9) that W_ ¼ u_ t z_ q z_ q AW 2Refq ð0Þf0 qðhÞg; ¼ AW 2Refq ð0Þf0 qð0Þg þ f0 ;
h 2 ð0; ~s; h ¼ 0;
ð3:13Þ
def
¼ AW þ H ðz; z; hÞ;
where H ðz; z; hÞ ¼ H 20 ðhÞ
z2 z2 þ H 11 ðhÞzz þ H 02 ðhÞ þ : 2 2
ð3:14Þ
On the other hand, on C0 near the origin W_ ¼ W z z_ þ W zz_ :
ð3:15Þ
We derive from (3.13)–(3.15) that ðA 2ix0 ÞW 20 ðhÞ ¼ H 20 ðhÞ;
AW 11 ðhÞ ¼ H 11 ðhÞ; :
ð3:16Þ
It follows from (3.11) and (3.13) that for h 2 ½~s; 0Þ, qðhÞ: H ðz; z; hÞ ¼ q ð0Þf0 qðhÞ q ð0Þf 0 qðhÞ ¼ gqðhÞ g
ð3:17Þ
Comparing the coefficients in (3.14) gives that for h 2 ½~s; 0Þ, H 20 ðhÞ ¼ g20 qðhÞ g02 qðhÞ
ð3:18Þ
H 11 ðhÞ ¼ g11 qðhÞ g11 qðhÞ:
ð3:19Þ
and
We derive from (3.16) and (3.18) and the definition of A that W_ 20 ðhÞ ¼ 2ix0 W 20 ðhÞ þ g20 qðhÞ þ g02 qðhÞ: Noting that q(h) = q(0)eix0 h , it follows that ig20 i g02 qð0Þeix0 h þ E1 e2ix0 h ; qð0Þeix0 h þ x0 3x0 ð1Þ ð2Þ where E1 ¼ E1 ; E1 2 R2 is a constant vector. W 20 ðhÞ ¼
ð3:20Þ
Similarly, from (3.16) and (3.19), we can obtain ig i g11 W 11 ðhÞ ¼ 11 qð0Þeix0 h þ qð0Þeix0 h þ E2 ; x0 x0 ð1Þ ð2Þ where E2 ¼ E2 ; E2 2 R2 is also a constant vector.
ð3:21Þ
In what follows, we seek appropriate E1 and E2. From the definition of A and (3.16), we obtain Z 0 dgðhÞW 20 ðhÞ ¼ 2ix0 W 20 ð0Þ H 20 ð0Þ
ð3:22Þ
~s
and Z
0
dgðhÞW 11 ðhÞ ¼ H 11 ð0Þ;
ð3:23Þ
~s
where g(h) = g(0,h). From (3.13), it follows that ð1Þ
H 20 ð0Þ ¼ g20 qð0Þ g02 qð0Þ þ
ð1Þ
ð1Þ
ð1Þ
2f110 q þ 2f101 eix0 ~s þ f020 q2 þ f200 ð2Þ
ð2Þ
! ð3:24Þ
ð2Þ
f20 þ f02 q2 þ 2f11 q
and ð1Þ
H 11 ð0Þ ¼ g11 qð0Þ g11 qð0Þ þ
ð1Þ
ð1Þ
ð1Þ
2f110 Refqg þ 2f101 Refeix0 ~s g þ f020 q q þ f200 ð2Þ
ð2Þ
ð2Þ
Þ f20 þ f02 q q þ f11 ðq þ q
! :
ð3:25Þ
Q. Gan et al. / Chaos, Solitons and Fractals 39 (2009) 1883–1895
1891
Substituting (3.20) and (3.24) into (3.22) and noticing that Z 0 eix0 h dgðhÞ qð0Þ ¼ 0 ix0 I ~s
and Z ix0 I
0
eix0 h dgðhÞ qð0Þ ¼ 0;
~s
we obtain Z 2ix0 I
0
e
2ix0 h
dgðhÞ E1 ¼
ð1Þ
ð2Þ
p2 E1 ¼ 2ix0 p4
It follows that ð1Þ E1
ð1Þ ð1Þ ð1Þ ð1Þ 1 2f110 q þ 2f101 eix0 ~s þ f020 q2 þ f200 ¼ ð2Þ ð2Þ ð2Þ A f þ f q2 þ 2f q 20
and ð2Þ E1
02
11
1 2ix0 p1 p5 e2ix0 ~s ¼ A p3
ð1Þ
ð1Þ
ð2Þ
! ;
ð2Þ
f20 þ f02 q2 þ 2f11 q
~s
which leads to 2ix0 p1 p5 e2ix0 ~s p3
ð1Þ
2f110 q þ 2f101 eix0 ~s þ f020 q2 þ f200
ð1Þ
ð1Þ
ð1Þ
ð1Þ
2f110 q þ 2f101 eix0 ~s þ f020 q2 þ f200 ð2Þ
ð2Þ
ð2Þ
f20 þ f02 q2 þ 2f11 q
! :
2ix0 p4 p2
ð1Þ ð1Þ ð1Þ ð1Þ 2f110 q þ 2f101 eix0 ~s þ f020 q2 þ f200 ; ð2Þ ð2Þ ð2Þ f þ f q2 þ 2f q 20
02
11
where 2ix0 p1 p5 e2ix0 ~s A ¼ p3
p2 : 2ix0 p4
Similarly, substituting (3.21) and (3.25) into (3.23), we get ! ð1Þ ð1Þ ð1Þ ð1Þ p1 p5 p2 q þ f200 2f110 Refqg þ 2f101 Refeix0 ~s g þ f020 q E2 ¼ ð2Þ ð2Þ ð2Þ p3 p4 Þ f20 þ f02 q q þ f11 ðq þ q and hence, ð1Þ E2
ð2Þ
E2
ð1Þ ð1Þ ð1Þ ð1Þ ix ~s q þ f200 p2 1 2f110 Refqg þ 2f101 Refe 0 g þ f020 q ¼ ; ð2Þ p4 p5 f ð2Þ þ f ð2Þ q Þ p4 20 02 q þ f11 ðq þ q ð1Þ ð1Þ ð1Þ ð1Þ ix ~s q þ f200 1 p1 p5 2f110 Refqg þ 2f101 Refe 0 g þ f020 q ¼ : ð2Þ ð2Þ ð2Þ p4 p5 p Þ f þ f q q þ f ðq þ q 4
20
02
11
Thus, we can determine W20(h) and W11(h) from (3.20) and (3.21). Furthermore, we can determine g21. Therefore, each gij in (3.12) is determined by the parameters and delay in system (3.1). Thus, we can compute the following values: ! i jg02 j2 g 2 þ 21 ; g11 g20 2jg11 j c1 ð0Þ ¼ 3 2 2x0 Refc1 ð0Þg ; b ¼ 2Refc1 ð0Þg; Refk0 ð~sÞg 2 Imfc1 ð0Þg þ l2 Imfk0 ð~sÞg ; T2 ¼ x0 l2 ¼
ð3:26Þ
which determine the quantities of bifurcating periodic solutions in the center manifold at the critical value ~s, i.e., l2 determines the direction of the Hopf bifurcation: if l2 > 0 (l2 < 0), then the Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solutions exist for s > ~sðs < ~sÞ; b2 determines the stability of the bifurcating periodic
1892
Q. Gan et al. / Chaos, Solitons and Fractals 39 (2009) 1883–1895
solutions: the bifurcating periodic solutions are stable (unstable) if b2 < 0(b2 > 0); and T2 determines the period of the bifurcating periodic solutions: the period increase (decrease) if T2 > 0(T2 < 0). From what has been discussed above, we see that if the values of r1, r2, a1, a2, k, m and s are given, we may determine the stability and direction of periodic solutions bifurcating from the positive equilibrium E* at the critical point sj. In the following, we give an example to illustrate the result above. Example 1. In system (1.2), we choose r1 = 2, r2 = 0.4, a1 = a2 = 1, k = 20, m = 2, then system (1.2) has a positive 2jp equilibrium E* = (17, 12.75). It is easy to show that sj 0:8768 þ 1:7124 . By Theorem 2.1, we see that the positive * equilibrium E is stable when s < s0 (see, Fig. 1), and system (1.2) undergoes a Hopf bifurcation at sj. Further, we can determine the stability and direction of periodic solutions bifurcating from the positive equilibrium E* at the critical point sj. For example, when s = s0 0.8768, c1(0) 0.0014 0.0018i. It follows from (3.26) that l2 > 0 and b < 0. Therefore, the bifurcation takes place when s crosses s0 to the right (s > s0), and the corresponding periodic orbits are orbitally asymptotically stable, as depicted in Fig. 2.
4. Chaotic behavior in system (1.4) In this section, we numerically study the chaotic behavior in system (1.4). In system (1.4), we set r1 = 2, r2 = 0.4, a1 = a2 = 1, k = 20, m = 2, f = 0.8,T = 2, 0.9 6 s 6 1.2. The influences of s may be documented by stroboscopically sampling some of the variables over a range of s values. We numerically integrate system (1.4) for 500 pulsing cycles at each value of s. For each s we plot the last 200 measures of the prey x and the predator y. Since we sample at forcing period, the T periodic solutions appear as fixed points, the 2T periodic solutions appear as two cycles, and so forth. The resulting bifurcation diagrams (see, Fig. 3) clearly show that: with the increasing of s from 0.9 to 1.2, system (1.4) experiences process of cycles ! periodic doubling cascade (see, Fig. 4) ! chaos (see, Fig. 5). When s is small ðs < s1 0:958Þ, the T periodic solution is stable. From Fig. 3, we note that if s > s1 0:958, the T-periodic solution is unstable and the periodic doubling cascade occurs. When s > s2 1:09, there is a chaotic band. A typical chaotic oscillation is captured when s = 1.1.
30
predator y
prey x
25 20 15 10 5 0
13
13
12
12
11
11
10
10 predator y
35
9 8 7
50
100
150 t
200
250
300
8 7
6
6
5
5 4
4 0
9
3
0
50
100
150 t
200
250
3
300
0
5
10
15 20 prey x
25
30
35
Fig. 1. The temporal solution found by numerical integration of system (1.2) with r1 = 2, r2 = 0.4, a1 = a2 = 1, k = 20, m = 2, s = 0.8 < s0 0.8768, and initial values x = 2.5, y = 4.
40 35
14
14
12
12
10
10
predator y
prey x
25 20 15
predator y
30
8
8
6
6
4
4
10 5 0
0
50
100
150 t
200
250
300
2
0
50
100
150 t
200
250
300
2
0
5
10
15
20 25 prey x
30
35
40
Fig. 2. The temporal solution found by numerical integration of system (1.2) with r1 = 2, r2 = 0.4, a1 = a2 = 1, k = 20, m = 2, s = 0.9 > s0 0.8768, and initial values x = 2.5, y = 4.
Q. Gan et al. / Chaos, Solitons and Fractals 39 (2009) 1883–1895
1893
9.4
5
9.2
4.8
9
4.6
8.8
4.4 predator y
predator y
Fig. 3. Bifurcation diagrams of system (1.4) showing the effect of s with r1 = 2, r2 = 0.4, a1 = a2 = 1, k = 20, m = 2, f = 0.8, T = 2, 0.9 6 s 6 1.2 and initial values x = 2.5, y = 4.
8.6 8.4 8.2
4 3.8
8
3.6
7.8
3.4
7.6
3.2
7.4 17.4 17.45 17.5 17.55 17.6 17.65 17.7 17.75 17.8 prey x
3
5
0
5
10
15
20 25 prey x
30
35
40
45
5.5 5
predator y
4.5 predator y
4.2
4 3.5
4.5 4 3.5
3 2.5
3
0
10
20
30
40
prey x
50
2.5
0
10
20
30
40
50
prey x
Fig. 4. Period doubling cascade: (a) phase portrait of T-periodic solution for s = 0.92; (b) phase portrait of 2T-periodic solution for s = 1; (c) phase portrait of 4T-periodic solution for s = 1.05; (d) phase portrait of 8T-periodic solution for s = 1.08.
6
60
6
5.5
50
5.5
4
20
3.5
10
3 2.5
30
0
10
20
30 prey x
40
50
60
0 600 650 700 750 800 850 900 950 1000 t
5 predator y
40 prey x
predator y
5 4.5
4.5 4 3.5 3 2.5 600 650 700 750 800 850 900 950 1000 t
Fig. 5. Chaos of system (1.4) for s = 1.1: (a) time series of prey x; (b) time series of predator y; (c) phase portrait.
Acknowledgements This work was supported by the National Natural Science Foundation of China (No. 10671209).
1894
Q. Gan et al. / Chaos, Solitons and Fractals 39 (2009) 1883–1895
References [1] Arditi R, Ginzburg LR. Coupling in predator–prey dynamics: ratio-dependence. J Theor Biol 1989;139:311–26. [2] Arditi R, Ginzburg LR, Akcakaya HR. Variation in plankton densities among lakes: a case for ratio-dependent models. Am Nat 1991;138:287–1296. [3] Arditi R, Perrin N, Saiah H. Functional response and heterogeneities: an experiment test with cladocerans. OIKOS 1991;60:9–75. [4] Arditi R, Saiah H. Empirical evidence of the role of heterogeneity in ratio-dependent consumption. Ecology 1992;73:544–1551. [5] Bainov DD, Simeonov DD. Impulsive differential equations: periodic solutions and applications. Pitman monographs and surveys in pure and applied mathematics, 66. Harlow: Longman Scientific & Technical; 1993. [6] Ballinger G, Liu X. Permanence of population growth models with impulsive effects. Math Comput Modell 1997;26:9–72. [7] Beretta E, Kuang Y. Global analyses in some delayed ratio-dependent predator–prey systems. Nonlinear Anal 1998;32:81–408. [8] Berezovskaya F, Karev G, Arditi R. Parametric analysis of the ratio-dependent predator–prey model. J Math Biol 2001;43:21–246. [9] Chen Y, Yu J, Sun C. Stability and Hopf bifurcation analysis in a three-level food chain system with delay. Chaos, Solitons & Fractals 2007;31:83–694. [10] D’Onofrio A. Stability properties of pulse vaccination strategy in SEIR epidemic model. Math Biosci 2002;179:7–72. [11] Fan M, Wang K. Periodicity in a delayed ratio-dependent predator–prey system. J Math Anal Appl 2001;262:179–90. [13] Gutierrez AP. The physiological basis of ratio-dependent predator–prey theory: a metabolic pool model of Nicholson’s blowflies as an example. Ecology 1992;73:1552–63. [14] Hale J, Lunel S. Introduction to functional differential equations. New York: Springer-Verlag; 1993. [15] Hanski I. The functional response of predator: worries about scale. TREE 1991;6:41–142. [16] Hassard B, Kazarinoff D, Wan Y. Theory and applications of Hopf bifurcation. Cambridge: Cambridge University Press; 1981. [17] Hsu SB, Hwang TW, Kuang Y. Global analysis of the Michaelis–Menten ratio-dependent predator–prey system. J Math Biol 2001;42:489–506. [18] Hsu SB, Hwang TW, Kuang Y. Rich dynamics of a ratio-dependent one prey two predator model. J Math Biol 2001;43:377–96. [19] Jiang W, Wei J. Bifurcation analysis in a limit cycle oscillator with delayed feedback. Chaos, Solitons & Fractals 2005;23:817–31. [20] Kuang Y. Delay differential equations with applications in population dynamics. New York: Academic Press; 1993. [21] Kuang Y, Beretta E. Global qualitative analysis of a ratio-dependent predator–prey system. J Math Biol 1998;36:389–406. [22] Lakmeche A, Arino O. Bifurcation of nontrivial periodic solutions of impulsive differential equations arising chemotherapeutic treatment. Dyn Continuous, Discrete Impulsive Syst 2000;7:265–87. [23] Lakshmikantham V, Bainov DD, Simeonov PS. Theory of impulsive differential equations. Series in modern applied mathematics, vol. 6. New Jersey: World Scientific; 1989. [24] Li S, Liao X, Li C. Hopf bifurcations in a Volterra prey–predator model with strong kernel. Chaos, Solitons & Fractals 2004;22:713–22. [25] Liu Z, Yuan R. Stability and Hopf bifurcations in a harvested one-prey-two-predator model with delays. Chaos, Solitons & Fractals 2006;27:1395–407. [26] Meng X, Wei J. Stability and Hopf bifurcations of mutual system with time delay. Chaos, Solitons & Fractals 2004;21:729–40. [27] Panetta JC. A mathematical model of periodically pulsed chemotherapy: tumor recurrence and metastasis in a competitive environment. Bull Math Biol 1996;58:425–47. [28] Roberts MG, Kao RR. The dynamics of an infectious disease in a population with birth pulses. Math Biosci 1998;149:23–36. [29] Shulgin B, Stone L, Agur Z. Pulse vaccination strategy in the SIR epidemic model. Bull Math Biol 1998;60:1123–48. [30] Song Y, Han M, Peng Y. Stability and Hopf bifurcations in a competitive Lotka–Volterra system with two delays. Chaos, Solitons & Fractals 2004;22:1139–48. [31] Song Y, Han M, Wei J. Stability and Hopf bifurcation analysis on a simplified BAM neural network with delays. Physica D 2005;200:185–204. [32] Song Y, Wei J. Bifurcation analysis for Chen’s system with delayed feedback and its application to control of chaos. Chaos, Solitons & Fractals 2004;22:75–91. [33] Sun C, Lin Y, Han M. Stability and Hopf bifurcations for an epidemic disease model with delay. Chaos, Solitons & Fractals 2006;30:204–16. [34] Sun C, Lin Y, Han M. Analysis of stability and Hopf bifurcations for an epidemic disease model with delay. Chaos, Solitons & Fractals 2006;30:204–16. [35] Sun C, Han M, Lin Y. Analysis of stability and Hopf bifurcations for a delayed logistic equation. Chaos, Solitons & Fractals 2007;31:672–82. [36] Tang S, Chen L. Density-dependent birth rate, birth pulses and their population dynamic consequences. J Math Biol 2002;44:185–99. [37] Wang F, Zhang S, Chen L, Sun L. Bifurcation and complexity of Monod type predator–prey system in a pulsed chemostat. Chaos, Solitons & Fractals 2006;27:447–58. [38] Wang L, Zou X. Hopf bifurcation in bidirectional associative memory neural network with delays: analysis and computation. J Comput Appl Math 2004;167:73–90. [39] Xu R, Chen L. Persistence and stability for a two-species ratio-dependent predator–prey system with time delay in a two-patch environment. Comput Math Appl 2000;40:577–88.
Q. Gan et al. / Chaos, Solitons and Fractals 39 (2009) 1883–1895
1895
[40] Xu R, Chen L. Persistence and global stability for n-species ratio-dependent predator–prey system with time delays. J Math Anal Appl 2002;275:27–43. [41] Zhou L, Tang Y, Hussein S. Stability and Hopf bifurcations for a delay competition diffusion system. Chaos, Solitons & Fractals 2002;14:1201–25.
Bifurcation and chaos in a ratio-dependent predator–prey system with time delay Qintao Gan a
a,*
, Rui Xu
a,b
, Pinghua Yang
a
Institute of Applied Mathematics, Shijiazhuang Mechanical Engineering College, Shijiazhuang 050003, PR China b Department of Applied Mathematics, Xi’an Jiaotong University, Xi’an 710049, PR China Accepted 20 June 2007
Abstract In this paper, a ratio-dependent predator–prey model with time delay is investigated. We first consider the local stability of a positive equilibrium and the existence of Hopf bifurcations. By using the normal form theory and center manifold reduction, we derive explicit formulae which determine the stability, direction and other properties of bifurcating periodic solutions. Finally, we consider the effect of impulses on the dynamics of the above time-delayed population model. Numerical simulations show that the system with constant periodic impulsive perturbations admits rich complex dynamic, such as periodic doubling cascade and chaos. 2007 Elsevier Ltd. All rights reserved.
1. Introduction In population dynamics, a functional response of the predator to the prey density refers to the change in the density of prey per unit time per predator as a function of the prey density. Standard Lotka–Volterra type models, on which a large body of existing predator–prey theory is built, assume that the per capita rate of predation depends solely on the prey numbers only (it is usually called prey-dependent functional response). Predator–prey systems with prey-dependent response have been studied extensively and the dynamics of such systems are now very well understood (see, for example, Kuang [20] and the references cited therein). Recently, the traditional prey-dependent predator–prey models have been challenged by several biologists (see, for example, Arditi and Ginzburg [1], Arditi et al. [2] and Gutierrez [13]) based on the fact that functional and numerical responses over typical ecological timescales ought to depend on the densities of both prey and predators (most likely and simply on their ratio), especially when predators have to search for food (and therefore have to share or compete for food). Such a functional response is called a ratio-dependent response function. These hypotheses are strongly supported by numerous field and laboratory experiment and observations [1–4,15]. Based on the Michaelis–Menten or Holling type II function, Arditi and Ginzburg [1] proposed a ratio-dependent function of the form
*
Corresponding author. E-mail address: [email protected] (Q. Gan).
0960-0779/$ - see front matter 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2007.06.122
1884
Q. Gan et al. / Chaos, Solitons and Fractals 39 (2009) 1883–1895
P
x cx=y cx ¼ ¼ ; y m þ x=y my þ x
and the following ratio-dependent predator–prey model: x_ ¼ xða bxÞ cxy=ðmy þ xÞ; y_ ¼ yðd þ fx=ðmy þ xÞÞ:
ð1:1Þ
Here, x(t) and y(t) represent population densities of the prey and the predator at time t, respectively; parameters a; b; c; d; f ; m are positive constants in which a/b is the carrying capacity of the prey, d > 0 is the death rate of the predator, and a; c; m; and f represent the prey intrinsic growth rate, capturing rate, half saturation constant and conversion rate, respectively. The ratio-dependent predator–prey model (1.1) has been studied by several researchers recently and very rich dynamics have been observed (see, for example, [7,8,11,15,17,18,21,39,40]). On the other hand, time delays of one type or another have been incorporated into biological models by many researchers, we refer to the monographs of Gopalsamy [17], Kuang [21] and references cited therein for general delayed biological systems and for studies on delayed predator–prey systems. In general, delay differential equations exhibit much more complicated dynamics than ordinary differential equations since a time delay could cause a stable equilibrium to become unstable and cause the populations to fluctuate (see, for example, [9,19,24–26,30–35,37,38,41]). In this paper, we are concerned with the following ratio-dependent predator–prey system with time delay xðt sÞ a1 xðtÞyðtÞ x_ ðtÞ ¼ r1 xðtÞ 1 ; k myðtÞ þ xðtÞ ð1:2Þ a2 xðtÞ y_ ðtÞ ¼ yðtÞ r2 þ ; myðtÞ þ xðtÞ where x(t), y(t) denote the density of the prey and the predator populations at time t, respectively; the parameter s P 0 is a time delay due to regulation of the prey. Parameters a1, a2, r1, r2, k and m are positive constants whose ecological meanings were mentioned above for model (1.1). The initial conditions for system (1.2) take the form xðhÞ ¼ /ðhÞ;
yðhÞ ¼ wðhÞ;
/ðhÞ P 0;
wðhÞ P 0;
h 2 ½s; 0;
/ð0Þ > 0;
wð0Þ > 0;
ð1:3Þ
Cð½s; 0; R2þ0 Þ,
the Banach space of continuous functions mapping the interval [s, 0] into R2þ0 , where ð/ðhÞ; wðhÞÞ 2 2 where Rþ0 ¼ fðx1 ; x2 Þ : x1 P 0; x2 P 0g. It is well-known by the fundamental theory of functional differential equations [14], system (1.2) has a unique solution (x(t), y(t)) satisfying initial conditions (1.3). It is easy to show that all solutions of system (1.2) with initial conditions (1.3) are defined on [0, + 1) and remain positive for all t P 0. We note that impulsive differential equations are suitable for the mathematical simulation of evolutionary process in which the parameters undergo relatively long periods of smooth variation followed by a short-term rapid change (i.e., jumps) in their values. Recently, equations with impulsive effect have been found in almost every domain of applied science. Numerous examples are given in Bainovs and his collaborator’s book [5,23]. Some impulsive differential equations have recently been introduced in population dynamics in relation to impulsive birth [28,36], impulsive vaccination [10,29], chemotherapeutic treatment of disease [22,27], and population ecology [6]. Motivated by the work above, in this paper, we further discuss the effect of impulses on the dynamics of system (1.2). To this end, we discuss the following impulsive equations 8 9 a1 xðtÞyðtÞ > > > = x_ ðtÞ ¼ r1 xðtÞ 1 xðtsÞ ; > k myðtÞþxðtÞ > > t–nT ; > < a2 xðtÞ > ; y_ ðtÞ ¼ yðtÞ r2 þ myðtÞþxðtÞ ; ð1:4Þ > > þ > > DxðtÞ ¼ xðt Þ xðt Þ ¼ 0; > > t ¼ nT : : DyðtÞ ¼ yðtþ Þ yðt Þ ¼ fyðtþ Þ; The organization of this paper is as follows. In Section 2, we discuss the local stability of a positive equilibrium of system (1.2) and the existence of Hopf bifurcations. In Section 3, we study the stability of the bifurcating periodic solutions and the direction of the Hopf bifurcation in system (1.2) by using the normal form theory and center manifold reduction. In Section 4, we numerically study the effect of impulsive perturbation on the dynamics of system (1.2). Numerical simulations show that system (1.4) admits rich complex dynamics such as periodic doubling cascade and chaos.
Q. Gan et al. / Chaos, Solitons and Fractals 39 (2009) 1883–1895
1885
2. Stability of a positive equilibrium and the existence of Hopf bifurcations In this section, we discuss the local stability of a positive equilibrium and the existence of Hopf bifurcations in system (1.2). It is easy to show that system (1.2) has a unique positive equilibrium E* = (x*, y*) if the following holds: (H1) a2 > r2 and ma2r1 > a1(a2 r2), where x ¼ k Set
ka1 ða2 r2 Þ ; ma2 r1
y ¼
kða2 r2 Þ ka1 ða2 r2 Þ2 : m2 a2 r1 r2 mr2
xðt sÞ a1 xðtÞyðtÞ f ð1Þ ¼ r1 xðtÞ 1 ; k myðtÞ þ xðtÞ a2 xðtÞ f ð2Þ ¼ yðtÞ r2 þ ; myðtÞ þ xðtÞ a1 r2 ða2 r2 Þ a1 r22 ; p2 ¼ ; 2 a22 ma2 r2 ðr2 a2 Þ r 1 x p4 ¼ ; p5 ¼ : k a2 p1 ¼
p3 ¼
ða2 r2 Þ2 ; ma2
ð2:1Þ
Linearizing system (1.2) at E*, we derive that x_ ðtÞ ¼ p1 xðtÞ þ p2 yðtÞ þ p5 xðt sÞ; _ ¼ p3 xðtÞ þ p4 yðtÞ; yðtÞ
ð2:2Þ
where pi (i = 1, 2, 3, 4, 5) are defined in (2.1). Hence, the characteristic equation of system (1.2) at the positive equilibrium E* takes the form k2 ðp1 þ p4 Þk p5 ðk p4 Þeks ¼ 0:
ð2:3Þ
It is well-known that the zero steady state of system (2.2) is asymptotically stable if all roots of Eq. (2.3) have negative real parts, and is unstable if Eq. (2.3) has a root with positive real part. In the sequel, we shall investigate the distribution of roots of Eq. (2.3). If ix(x > 0) is a root of Eq. (2.3), then x2 iðp1 þ p4 Þx p5 ðix p4 Þeixs ¼ 0: Separating the real and imaginary parts, we obtain x2 ¼ p5 x sin xs p4 p5 cos xs; ðp1 þ p4 Þx ¼ p5 x cos xs þ p4 p5 sin xs:
ð2:4Þ
It follows from (2.4) that x4 þ ½ðp1 þ p4 Þ2 p25 x2 p24 p25 ¼ 0:
ð2:5Þ
Clearly, Eq. (2.5) has only one positive root x0 defined by " ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi!#12 r 2 1 2 2 2 2 2 2 p5 ðp1 þ p4 Þ þ 4p4 p5 : x0 ¼ p5 ðp1 þ p4 Þ þ 2 Define sj ¼
1 x0 ðx20 þ ðp1 þ p4 Þp4 Þ þ 2jp ; arcsin x0 p5 ðp24 þ x20 Þ
j ¼ 0; 1; . . . :
Then (sj, x0) solves Eq. (2.4). This means that when s = sj, Eq. (2.3) has a pair of purely imaginary roots ± ix0. Now let us consider the behavior of roots of Eq. (2.3) near sj. Denote kðsÞ ¼ aðsÞ þ ixðsÞ
ð2:6Þ
ð2:7Þ
1886
Q. Gan et al. / Chaos, Solitons and Fractals 39 (2009) 1883–1895
the root of Eq. (2.3) such that aðsj Þ ¼ 0;
xðsj Þ ¼ x0 :
Substituting k(s) into Eq. (2.3) and differentiating both sides of it with respect to s, we have 1 dkðsÞ ð2k ðp1 þ p4 ÞÞeks p5 s þ : ¼ p5 kðk p4 Þ p5 kðk p4 Þ k ds We therefore derive that 1 dk ð2k ðp1 þ p4 ÞÞeks p5 s ¼ Re þ Re Re p5 kðk p4 Þ p5 kðk p4 Þ k s¼sj ds s¼sj s¼sj ) (
1 ¼ Re ðp1 þ p4 Þ cos x0 sj 2x0 sin x0 sj :þi½2x0 cos x0 sj ðp1 þ p4 Þ sin x0 sj p5 x20 þ ip4 p5 x0 p5 þ Re p5 x20 þ ip4 p5 x0 1 ¼ fp5 x20 ½ðp1 þ p4 Þ cos x0 sj þ 2x0 sin x0 sj þp4 p5 x0 ½2x0 cos x0 sj ðp1 þ p4 Þ sin x0 sj p25 x20 C 1 ¼ fðp1 þ p4 Þx0 ½p5 x0 cos x0 sj þ p4 p5 sin x0 sj 2x20 ½p5 x0 sin x0 sj p4 p5 cos x0 sj p25 x20 C o x2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x2 n ¼ 0 ðp1 þ p4 Þ2 þ 2x20 p25 ¼ 0 ðp25 ðp1 þ p4 Þ2 Þ2 þ 4p24 p25 ; C C where we have used Eqs. (2.4) and (2.6), and C ¼ p25 x40 þ p24 p25 x20 > 0. Hence, it follows that ( ) ( ) 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 dk dk x ¼ sign Re ¼ sign 0 ðp25 ðp1 þ p4 Þ2 Þ2 þ 4p24 p25 > 0: sign Re C ds s¼sj ds s¼sj
ð2:8Þ
Therefore, when the delay s near sj is increased, the root of Eq. (2.3) crosses the imaginary axis from left to right. In addition, note that when s = 0, Eq. (2.3) has roots with negative real parts only if (H2) p1 + p4 + p5 < 0 and p4p5 > 0. Thus, by the well-known Rouche´ theorem we obtain the following results on the distribution of roots of Eq. (2.3). Lemma 2.1. Let sj (j = 0,1, . . .) be defined as in (2.7). Then all roots of Eq. (2.3) have negative real parts for all s 2 [0, s0). However, Eq. (2.3) has at least one root with positive real part when s > s0, and Eq. (2.3) has a pair of purely imaginary root ± ix0 when s = s0. More detail, for s 2 (sj,sj + 1] (j = 0, 1, . . .), Eq. (2.3) has 2(j + 1) roots with positive real parts. Moreover, all roots of Eq. (2.3) with s = sj (j = 0, 1, . . .) have negative real parts except ±ix0. Applying Lemma 2.1, the transversality condition (2.8) and Theorem 11.1 developed in [14], we obtain the following results. Theorem 2.1. Let (H1) and (H2) hold. let x0 and sj (j = 0, 1, . . .) be defined as in (2.6) and (2.7), respectively. (i) The positive equilibrium E* of system (1.2) is asymptotically stable for all s 2 [0, s0) and unstable for s > s0. (ii) System (1.2) undergoes a Hopf Bifurcation at the positive equilibrium E* when s = sj (j = 0, 1, . . .).
3. Direction of Hopf bifurcation In Section 2, we have shown that system (1.2) admits a series of periodic solutions bifurcating from the positive equilibrium E * at the critical values of s. In this section, we derive explicit formulae to determine the properties of the Hopf bifurcation at critical values sj by using the normal form theory and center manifold reduction (see, for example, Hassard et al. [16]). Without loss of generality, denote the critical values sj by ~s, and set s ¼ ~s þ l. Then l = 0 is a Hopf bifurcation value of system (1.2). Thus, we can work in the phase space C ¼ Cð½~s; 0; R2 Þ.
Q. Gan et al. / Chaos, Solitons and Fractals 39 (2009) 1883–1895
Let u1 ðtÞ ¼ xðtÞ x ;
1887
u2 ðtÞ ¼ yðtÞ y :
Then system (1.2) is transformed into 8 P ð1Þ i j 1 l > > < u_ 1 ðtÞ ¼ p1 u1 ðtÞ þ p2 u2 ðtÞ þ p5 u2 ðt ~sÞ þ iþjþlP2 i!j!l! fijl u1 ðtÞu2 ðtÞu1 ðt ~sÞ; P 1 ð2Þ i > > f u1 ðtÞuj2 ðtÞ; : u_ 2 ðtÞ ¼ p3 u1 ðtÞ þ p4 u2 ðtÞ þ i!j! ij
ð3:1Þ
iþjP2
where ð1Þ fijl
¼ ; oxi oy j oxðt ~sÞl ðx ;y ;x Þ oiþjþl f ð1Þ
ð2Þ
fij ¼
oiþj f ð2Þ ; oxi oy j ðx ;y Þ
i; j; l P 0;
here, f (1) and f (2) are defined in (2.1). For the simplicity of notations, we rewrite (3.1) as _ ¼ Ll ðut Þ þ f ðl; ut Þ; uðtÞ T
ð3:2Þ 2
where u(t) = (u1(t),u2(t)) 2 R , ut(h) 2 C is defined by ut(h) = u(t + h), and Ll:C ! R, f:R · C 2 R are given by p5 0 p1 p2 Ll / ¼ /ð0Þ þ /ð~sÞ; ð3:3Þ p3 p4 0 0 and
1 0 P ð1Þ 1 f /i1 ð0Þ/j2 ð0Þ/l1 ð~sÞ i!j!l! ijl C B iþjþlP2 C f ðl; /Þ ¼ B P 1 ð2Þ i A @ j f / ð0Þ/ ð0Þ ij 2 1 i!j!
ð3:4Þ
iþjP2
respectively. By Riesz representation theorem, there exists a function g(h,l) of bounded variation for h 2 ½~s; 0 such that Z 0 Ll / ¼ dgðh; 0Þ/ðhÞ for / 2 C: ð3:5Þ ~s
In fact, we can choose p1 p2 p5 gðh; lÞ ¼ dðhÞ 0 p3 p4
0 dðh þ ~sÞ; 0
ð3:6Þ
where d is the Dirac delta function. For / 2 C 1 ð½~s; 0; R2 Þ, define ( d/ðhÞ ; h 2 ½~s; 0Þ; AðlÞ/ ¼ Rdh 0 dgðl; sÞ/ðsÞ; h ¼ 0 ~s and
RðlÞ/ ¼
0;
h 2 ½~s; 0Þ;
f ðl; /Þ;
h ¼ 0:
Then system (3.2) is equivalent to u_ t ¼ AðlÞut þ RðlÞut ;
ð3:7Þ
where xt(h) = x(t + h) for h 2 ½~s; 0. For w 2 C 1 ð½0; ~s; ðR2 Þ Þ, define ( dwðsÞ ; s 2 ð0; ~s; A wðsÞ ¼ R 0 ds dgT ðt; 0ÞwðtÞ; s ¼ 0 ~s and a bilinear inner product hwðsÞ; /ðhÞi ¼ wð0Þ/ð0Þ
Z
0
~s
Z
h
n¼0
hÞdgðhÞ/ðnÞdn; wðn
ð3:8Þ
1888
Q. Gan et al. / Chaos, Solitons and Fractals 39 (2009) 1883–1895
where g(h) = g(h,0). Then A(0) and A* are adjoint operators. By discussions in Section 2 and foregoing assumption, we know that ±ix0 are eigenvalues of A(0). Thus, they are also eigenvalues of A*. We first need to compute the eigenvector of A(0) and A* corresponding to ix0 and ix0, respectively. Suppose that q(h) = (1,q)T eix0 h is the eigenvector of A(0) corresponding to ix0. Then A(0)q(h) = ix0q(h). It follows from the definition of A(0), (3.5) and (3.6) that 0 p1 þ p5 eix0 ~s ix0 p2 qð0Þ ¼ : 0 p3 p4 ix0 We therefore derive that qð0Þ ¼ ð1; qÞT ¼ 1;
p3 ix0 p4
T :
On the other hand, suppose that q*(s) = D(1,r)eix0 s is the eigenvector of A* corresponding to ix0. From the definition of A*, (3.5) and (3.6) we have 0 p1 þ p5 eix0 ~s þ ix0 p3 ðq ð0ÞÞT ¼ ; 0 p2 p4 þ ix0 which yields
q ð0Þ ¼ Dð1; rÞ ¼ D 1;
p2 : p4 þ ix0
In order to assure hq*(s), q(h)i = 1, we need to determine the value of D. From (3.7), we have Z 0Z h Þð1; qÞT ÞeiðnhÞx0 dgðhÞð1; qÞT einx0 dn ð1; r hq ðsÞ; qðhÞi ¼ D ð1; r r ¼ D 1 þ q
Z
~s
n¼0
0 ix0 h
Þhe ð1; r
T
dgðhÞð1; qÞ
¼ Df1 þ q r þ p5~seix0 ~s g:
~s
Thus, we can choose 1 D¼ r þ p5~seix0 ~s 1þq such that hq ðsÞ; qðhÞi ¼ 1; hq ðsÞ; qðhÞi ¼ 0: In the remainder of this section, we use the same notations as in Hassard et al. [16]. We first compute the coordinates to describe the center manifold C0 at l = 0. Let ut be the solution of Eq. (3.2) with l = 0. Define zðtÞ ¼ hq ; ut i; W ðt; hÞ ¼ ut ðhÞ 2RefzðtÞqðhÞg: ð3:9Þ On the center manifold C0 we have W ðt; hÞ ¼ W ðzðtÞ; zðtÞ; hÞ; where z2 z2 þ W 11 ðhÞzz þ W 02 ðhÞ þ ; 2 2 q . Note that W is real if ut is real. we z and z are local coordinates for center manifold C0 in the direction of q* and consider only real solutions. For the solution ut 2 C0 of (3.2), since l = 0, we have W ðz; z; hÞ ¼ W 20 ðhÞ
z_ ¼ ix0 z þ hq ðhÞ; f ð0; W ðz; z; hÞ þ 2RefzqðhÞgÞi ¼ ix0 z þ q ðhÞf ð0; W ðz; z; hÞ þ 2RefzqðhÞgÞ def
q ð0Þf0 ðz; zÞ: ¼ ix0 z þ q ð0Þf ð0; W ðz; z; 0Þ þ 2Refzqð0ÞgÞ ¼ ix0 z þ
ð3:10Þ
We rewrite (3.10) as z_ ¼ ix0 z þ gðz; zÞ with gðz; zÞ ¼ q ð0Þf0 ðz; zÞ ¼ g20
z2 z2 z2z þ g11 zz þ g02 þ g21 þ : 2 2 2
ð3:11Þ
Q. Gan et al. / Chaos, Solitons and Fractals 39 (2009) 1883–1895
1889
Noting that ut ðhÞ ¼ ðu1t ðhÞ; u2t ðhÞÞ ¼ W ðt; hÞ þ zqðhÞ þ z qðhÞ and q(h) = (1,q)T eix0 h , we have ð1Þ
u1t ð0Þ ¼ z þ z þ W 20 ð0Þ
z2 z2 ð1Þ ð1Þ þ W 11 ð0Þzz þ W 02 ð0Þ þ ; 2 2 ð1Þ
u1t ð~sÞ ¼ eix0 ~s z þ eix0 ~sz þ W 20 ð~sÞ ð2Þ
z þ W 20 ð0Þ u2t ð0Þ ¼ qz þ q
z2 z2 ð1Þ ð1Þ þ W 11 ð~sÞzz þ W 02 ð~sÞ þ ; 2 2
z2 z2 ð2Þ ð2Þ þ W 11 ð0Þzz þ W 02 ð0Þ þ : 2 2
Thus, it follows from (3.4) and (3.11) that gðz; zÞ ¼ q ð0Þf0 ðz; zÞ 0 P B iþjþlP2
ÞB ¼ Dð1; r @
1 ð1Þ 1 f ui ð0Þuj2t ð0Þul1t ð~sÞ i!j!l! ijl 1t
P
iþjP2
1 ð2Þ i f u1t ð0Þuj2t ð0Þ i!j! ij
C C A
1 ð1Þ 1 ð1Þ 1 1 ð1Þ ð1Þ f20ð2Þ þ f02ð2Þ q2 þ f11ð2Þ q z2 ¼ D f110 q þ f101 eix0 ~s þ f020 q2 þ f200 z2 þ r 2 2 2 2 ð1Þ ð1Þ ð1Þ ð1Þ þ 2f110 Refqg þ 2f101 Refeix0 ~s g þ f020 q q þ f200 zz 1 ð1Þ 2 1 ð1Þ 2 ð2Þ ð1Þ ð1Þ ix0 ~s Þ zz þ f110 þ f101 þ f200 z f20ð2Þ þ f02ð2Þ q þr q þ f11 ðq þ q e þ f020 q q 2 2 1 1 1 ð2Þ ð1Þ 1 ð1Þ ð1Þ 2 þ f11ð2Þ q z2 þ f110 f20ð2Þ þ f02ð2Þ q þr q þ W 20 ð0Þ þ qW 11 ð0Þ W ð0Þ 2 2 2 20 2 1 ð1Þ ð2Þ ð1Þ 1 ð1Þ ð1Þ ð1Þ þW 11 ð0Þ þ f101 W 20 ð0Þeix0 ~s þ W 20 ð~sÞ þ qeix0 ~s W 11 ð0Þ þ W 11 ð~sÞ 2 2 1 1 ð1Þ ð1Þ ð2Þ ð1Þ ð1Þ ð1Þ 1 W ð2Þ ð0Þ þ f W ð0Þ þ ð0Þ þ f030 q2 q W þ f020 qW 11 ð0Þ þ q 20 200 11 2 2 20 2 1 ð1Þ 1 1 ð1Þ ð1Þ 1 ð1Þ z2z þ r f20ð2Þ W ð1Þ q þ f210 q þ q þ f300 þ f120 q2 þ q 11 ð0Þ þ W 20 ð0Þ 2 2 2 2 1 1 ð2Þ 1 ð2Þ ð2Þ ð2Þ ð2Þ ð1Þ ð1Þ 2 W ð2Þ ð0Þ þ f W ð0Þ þ ð0Þ þ qW ð0Þ þ ð0Þ z z þ : W q W þf02 qW 11 ð0Þ þ q 20 11 11 11 20 2 2 20 2 Comparing the coefficients in (3.11), we get ð1Þ ð1Þ ð1Þ ð1Þ f20ð2Þ þ f02ð2Þ q2 þ 2f11ð2Þ q D; g20 ¼ 2f110 q þ 2f101 eix0 ~s þ f020 q2 þ f200 þ r ð1Þ ð1Þ ð1Þ ð1Þ ð2Þ Þ D; f20ð2Þ þ f02ð2Þ q g11 ¼ 2f110 Refqg þ 2f101 Refeix0 ~s g þ f020 q q þ f200 þ r q þ f11 ðq þ q ð1Þ ð1Þ ix0 ~s ð1Þ 2 ð1Þ þ 2f101 þ f200 2 þ 2f11ð2Þ q D; f20ð2Þ þ f02ð2Þ q e þ f020 q þr g02 ¼ 2f110 q ð1Þ ð1Þ ð2Þ ð1Þ ð2Þ q þ W 20 ð0Þ þ 2qW 11 ð0Þ þ 2W 11 ð0Þ g21 ¼ f110 W 20 ð0Þ ð1Þ ð1Þ ð1Þ ð1Þ ð1Þ þf101 W 20 ð0Þeix0 ~s þ W 20 ð~sÞ þ 2qeix0 ~s W 11 ð0Þ þ 2W 11 ð~sÞ ð1Þ ð2Þ ð1Þ ð1Þ ð1Þ ð1Þ 2 W ð2Þ þ f020 2qW 11 ð0Þ þ q 20 ð0Þ þ f200 2W 11 ð0Þ þ W 20 ð0Þ þ f030 q q
ð1Þ ð1Þ ð1Þ ð2Þ ð1Þ ð1Þ Þ þ r f20 2W 11 ð0Þ þ W 20 ð0Þ q þ f210 ð2q þ q þ f300 þ f120 q2 þ 2q ð2Þ ð2Þ ð2Þ ð2Þ ð2Þ ð1Þ W ð2Þ W ð1Þ 2W 11 ð0Þ þ W 20 ð0Þ þ 2qW 11 ð0Þ þ q D: þf02 2qW 11 ð0Þ þ q 20 ð0Þ þ f11 20 ð0Þ
ð3:12Þ
1890
Q. Gan et al. / Chaos, Solitons and Fractals 39 (2009) 1883–1895
We now compute W20(h) and W11(h). It follows from (3.7) and (3.9) that W_ ¼ u_ t z_ q z_ q AW 2Refq ð0Þf0 qðhÞg; ¼ AW 2Refq ð0Þf0 qð0Þg þ f0 ;
h 2 ð0; ~s; h ¼ 0;
ð3:13Þ
def
¼ AW þ H ðz; z; hÞ;
where H ðz; z; hÞ ¼ H 20 ðhÞ
z2 z2 þ H 11 ðhÞzz þ H 02 ðhÞ þ : 2 2
ð3:14Þ
On the other hand, on C0 near the origin W_ ¼ W z z_ þ W zz_ :
ð3:15Þ
We derive from (3.13)–(3.15) that ðA 2ix0 ÞW 20 ðhÞ ¼ H 20 ðhÞ;
AW 11 ðhÞ ¼ H 11 ðhÞ; :
ð3:16Þ
It follows from (3.11) and (3.13) that for h 2 ½~s; 0Þ, qðhÞ: H ðz; z; hÞ ¼ q ð0Þf0 qðhÞ q ð0Þf 0 qðhÞ ¼ gqðhÞ g
ð3:17Þ
Comparing the coefficients in (3.14) gives that for h 2 ½~s; 0Þ, H 20 ðhÞ ¼ g20 qðhÞ g02 qðhÞ
ð3:18Þ
H 11 ðhÞ ¼ g11 qðhÞ g11 qðhÞ:
ð3:19Þ
and
We derive from (3.16) and (3.18) and the definition of A that W_ 20 ðhÞ ¼ 2ix0 W 20 ðhÞ þ g20 qðhÞ þ g02 qðhÞ: Noting that q(h) = q(0)eix0 h , it follows that ig20 i g02 qð0Þeix0 h þ E1 e2ix0 h ; qð0Þeix0 h þ x0 3x0 ð1Þ ð2Þ where E1 ¼ E1 ; E1 2 R2 is a constant vector. W 20 ðhÞ ¼
ð3:20Þ
Similarly, from (3.16) and (3.19), we can obtain ig i g11 W 11 ðhÞ ¼ 11 qð0Þeix0 h þ qð0Þeix0 h þ E2 ; x0 x0 ð1Þ ð2Þ where E2 ¼ E2 ; E2 2 R2 is also a constant vector.
ð3:21Þ
In what follows, we seek appropriate E1 and E2. From the definition of A and (3.16), we obtain Z 0 dgðhÞW 20 ðhÞ ¼ 2ix0 W 20 ð0Þ H 20 ð0Þ
ð3:22Þ
~s
and Z
0
dgðhÞW 11 ðhÞ ¼ H 11 ð0Þ;
ð3:23Þ
~s
where g(h) = g(0,h). From (3.13), it follows that ð1Þ
H 20 ð0Þ ¼ g20 qð0Þ g02 qð0Þ þ
ð1Þ
ð1Þ
ð1Þ
2f110 q þ 2f101 eix0 ~s þ f020 q2 þ f200 ð2Þ
ð2Þ
! ð3:24Þ
ð2Þ
f20 þ f02 q2 þ 2f11 q
and ð1Þ
H 11 ð0Þ ¼ g11 qð0Þ g11 qð0Þ þ
ð1Þ
ð1Þ
ð1Þ
2f110 Refqg þ 2f101 Refeix0 ~s g þ f020 q q þ f200 ð2Þ
ð2Þ
ð2Þ
Þ f20 þ f02 q q þ f11 ðq þ q
! :
ð3:25Þ
Q. Gan et al. / Chaos, Solitons and Fractals 39 (2009) 1883–1895
1891
Substituting (3.20) and (3.24) into (3.22) and noticing that Z 0 eix0 h dgðhÞ qð0Þ ¼ 0 ix0 I ~s
and Z ix0 I
0
eix0 h dgðhÞ qð0Þ ¼ 0;
~s
we obtain Z 2ix0 I
0
e
2ix0 h
dgðhÞ E1 ¼
ð1Þ
ð2Þ
p2 E1 ¼ 2ix0 p4
It follows that ð1Þ E1
ð1Þ ð1Þ ð1Þ ð1Þ 1 2f110 q þ 2f101 eix0 ~s þ f020 q2 þ f200 ¼ ð2Þ ð2Þ ð2Þ A f þ f q2 þ 2f q 20
and ð2Þ E1
02
11
1 2ix0 p1 p5 e2ix0 ~s ¼ A p3
ð1Þ
ð1Þ
ð2Þ
! ;
ð2Þ
f20 þ f02 q2 þ 2f11 q
~s
which leads to 2ix0 p1 p5 e2ix0 ~s p3
ð1Þ
2f110 q þ 2f101 eix0 ~s þ f020 q2 þ f200
ð1Þ
ð1Þ
ð1Þ
ð1Þ
2f110 q þ 2f101 eix0 ~s þ f020 q2 þ f200 ð2Þ
ð2Þ
ð2Þ
f20 þ f02 q2 þ 2f11 q
! :
2ix0 p4 p2
ð1Þ ð1Þ ð1Þ ð1Þ 2f110 q þ 2f101 eix0 ~s þ f020 q2 þ f200 ; ð2Þ ð2Þ ð2Þ f þ f q2 þ 2f q 20
02
11
where 2ix0 p1 p5 e2ix0 ~s A ¼ p3
p2 : 2ix0 p4
Similarly, substituting (3.21) and (3.25) into (3.23), we get ! ð1Þ ð1Þ ð1Þ ð1Þ p1 p5 p2 q þ f200 2f110 Refqg þ 2f101 Refeix0 ~s g þ f020 q E2 ¼ ð2Þ ð2Þ ð2Þ p3 p4 Þ f20 þ f02 q q þ f11 ðq þ q and hence, ð1Þ E2
ð2Þ
E2
ð1Þ ð1Þ ð1Þ ð1Þ ix ~s q þ f200 p2 1 2f110 Refqg þ 2f101 Refe 0 g þ f020 q ¼ ; ð2Þ p4 p5 f ð2Þ þ f ð2Þ q Þ p4 20 02 q þ f11 ðq þ q ð1Þ ð1Þ ð1Þ ð1Þ ix ~s q þ f200 1 p1 p5 2f110 Refqg þ 2f101 Refe 0 g þ f020 q ¼ : ð2Þ ð2Þ ð2Þ p4 p5 p Þ f þ f q q þ f ðq þ q 4
20
02
11
Thus, we can determine W20(h) and W11(h) from (3.20) and (3.21). Furthermore, we can determine g21. Therefore, each gij in (3.12) is determined by the parameters and delay in system (3.1). Thus, we can compute the following values: ! i jg02 j2 g 2 þ 21 ; g11 g20 2jg11 j c1 ð0Þ ¼ 3 2 2x0 Refc1 ð0Þg ; b ¼ 2Refc1 ð0Þg; Refk0 ð~sÞg 2 Imfc1 ð0Þg þ l2 Imfk0 ð~sÞg ; T2 ¼ x0 l2 ¼
ð3:26Þ
which determine the quantities of bifurcating periodic solutions in the center manifold at the critical value ~s, i.e., l2 determines the direction of the Hopf bifurcation: if l2 > 0 (l2 < 0), then the Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solutions exist for s > ~sðs < ~sÞ; b2 determines the stability of the bifurcating periodic
1892
Q. Gan et al. / Chaos, Solitons and Fractals 39 (2009) 1883–1895
solutions: the bifurcating periodic solutions are stable (unstable) if b2 < 0(b2 > 0); and T2 determines the period of the bifurcating periodic solutions: the period increase (decrease) if T2 > 0(T2 < 0). From what has been discussed above, we see that if the values of r1, r2, a1, a2, k, m and s are given, we may determine the stability and direction of periodic solutions bifurcating from the positive equilibrium E* at the critical point sj. In the following, we give an example to illustrate the result above. Example 1. In system (1.2), we choose r1 = 2, r2 = 0.4, a1 = a2 = 1, k = 20, m = 2, then system (1.2) has a positive 2jp equilibrium E* = (17, 12.75). It is easy to show that sj 0:8768 þ 1:7124 . By Theorem 2.1, we see that the positive * equilibrium E is stable when s < s0 (see, Fig. 1), and system (1.2) undergoes a Hopf bifurcation at sj. Further, we can determine the stability and direction of periodic solutions bifurcating from the positive equilibrium E* at the critical point sj. For example, when s = s0 0.8768, c1(0) 0.0014 0.0018i. It follows from (3.26) that l2 > 0 and b < 0. Therefore, the bifurcation takes place when s crosses s0 to the right (s > s0), and the corresponding periodic orbits are orbitally asymptotically stable, as depicted in Fig. 2.
4. Chaotic behavior in system (1.4) In this section, we numerically study the chaotic behavior in system (1.4). In system (1.4), we set r1 = 2, r2 = 0.4, a1 = a2 = 1, k = 20, m = 2, f = 0.8,T = 2, 0.9 6 s 6 1.2. The influences of s may be documented by stroboscopically sampling some of the variables over a range of s values. We numerically integrate system (1.4) for 500 pulsing cycles at each value of s. For each s we plot the last 200 measures of the prey x and the predator y. Since we sample at forcing period, the T periodic solutions appear as fixed points, the 2T periodic solutions appear as two cycles, and so forth. The resulting bifurcation diagrams (see, Fig. 3) clearly show that: with the increasing of s from 0.9 to 1.2, system (1.4) experiences process of cycles ! periodic doubling cascade (see, Fig. 4) ! chaos (see, Fig. 5). When s is small ðs < s1 0:958Þ, the T periodic solution is stable. From Fig. 3, we note that if s > s1 0:958, the T-periodic solution is unstable and the periodic doubling cascade occurs. When s > s2 1:09, there is a chaotic band. A typical chaotic oscillation is captured when s = 1.1.
30
predator y
prey x
25 20 15 10 5 0
13
13
12
12
11
11
10
10 predator y
35
9 8 7
50
100
150 t
200
250
300
8 7
6
6
5
5 4
4 0
9
3
0
50
100
150 t
200
250
3
300
0
5
10
15 20 prey x
25
30
35
Fig. 1. The temporal solution found by numerical integration of system (1.2) with r1 = 2, r2 = 0.4, a1 = a2 = 1, k = 20, m = 2, s = 0.8 < s0 0.8768, and initial values x = 2.5, y = 4.
40 35
14
14
12
12
10
10
predator y
prey x
25 20 15
predator y
30
8
8
6
6
4
4
10 5 0
0
50
100
150 t
200
250
300
2
0
50
100
150 t
200
250
300
2
0
5
10
15
20 25 prey x
30
35
40
Fig. 2. The temporal solution found by numerical integration of system (1.2) with r1 = 2, r2 = 0.4, a1 = a2 = 1, k = 20, m = 2, s = 0.9 > s0 0.8768, and initial values x = 2.5, y = 4.
Q. Gan et al. / Chaos, Solitons and Fractals 39 (2009) 1883–1895
1893
9.4
5
9.2
4.8
9
4.6
8.8
4.4 predator y
predator y
Fig. 3. Bifurcation diagrams of system (1.4) showing the effect of s with r1 = 2, r2 = 0.4, a1 = a2 = 1, k = 20, m = 2, f = 0.8, T = 2, 0.9 6 s 6 1.2 and initial values x = 2.5, y = 4.
8.6 8.4 8.2
4 3.8
8
3.6
7.8
3.4
7.6
3.2
7.4 17.4 17.45 17.5 17.55 17.6 17.65 17.7 17.75 17.8 prey x
3
5
0
5
10
15
20 25 prey x
30
35
40
45
5.5 5
predator y
4.5 predator y
4.2
4 3.5
4.5 4 3.5
3 2.5
3
0
10
20
30
40
prey x
50
2.5
0
10
20
30
40
50
prey x
Fig. 4. Period doubling cascade: (a) phase portrait of T-periodic solution for s = 0.92; (b) phase portrait of 2T-periodic solution for s = 1; (c) phase portrait of 4T-periodic solution for s = 1.05; (d) phase portrait of 8T-periodic solution for s = 1.08.
6
60
6
5.5
50
5.5
4
20
3.5
10
3 2.5
30
0
10
20
30 prey x
40
50
60
0 600 650 700 750 800 850 900 950 1000 t
5 predator y
40 prey x
predator y
5 4.5
4.5 4 3.5 3 2.5 600 650 700 750 800 850 900 950 1000 t
Fig. 5. Chaos of system (1.4) for s = 1.1: (a) time series of prey x; (b) time series of predator y; (c) phase portrait.
Acknowledgements This work was supported by the National Natural Science Foundation of China (No. 10671209).
1894
Q. Gan et al. / Chaos, Solitons and Fractals 39 (2009) 1883–1895
References [1] Arditi R, Ginzburg LR. Coupling in predator–prey dynamics: ratio-dependence. J Theor Biol 1989;139:311–26. [2] Arditi R, Ginzburg LR, Akcakaya HR. Variation in plankton densities among lakes: a case for ratio-dependent models. Am Nat 1991;138:287–1296. [3] Arditi R, Perrin N, Saiah H. Functional response and heterogeneities: an experiment test with cladocerans. OIKOS 1991;60:9–75. [4] Arditi R, Saiah H. Empirical evidence of the role of heterogeneity in ratio-dependent consumption. Ecology 1992;73:544–1551. [5] Bainov DD, Simeonov DD. Impulsive differential equations: periodic solutions and applications. Pitman monographs and surveys in pure and applied mathematics, 66. Harlow: Longman Scientific & Technical; 1993. [6] Ballinger G, Liu X. Permanence of population growth models with impulsive effects. Math Comput Modell 1997;26:9–72. [7] Beretta E, Kuang Y. Global analyses in some delayed ratio-dependent predator–prey systems. Nonlinear Anal 1998;32:81–408. [8] Berezovskaya F, Karev G, Arditi R. Parametric analysis of the ratio-dependent predator–prey model. J Math Biol 2001;43:21–246. [9] Chen Y, Yu J, Sun C. Stability and Hopf bifurcation analysis in a three-level food chain system with delay. Chaos, Solitons & Fractals 2007;31:83–694. [10] D’Onofrio A. Stability properties of pulse vaccination strategy in SEIR epidemic model. Math Biosci 2002;179:7–72. [11] Fan M, Wang K. Periodicity in a delayed ratio-dependent predator–prey system. J Math Anal Appl 2001;262:179–90. [13] Gutierrez AP. The physiological basis of ratio-dependent predator–prey theory: a metabolic pool model of Nicholson’s blowflies as an example. Ecology 1992;73:1552–63. [14] Hale J, Lunel S. Introduction to functional differential equations. New York: Springer-Verlag; 1993. [15] Hanski I. The functional response of predator: worries about scale. TREE 1991;6:41–142. [16] Hassard B, Kazarinoff D, Wan Y. Theory and applications of Hopf bifurcation. Cambridge: Cambridge University Press; 1981. [17] Hsu SB, Hwang TW, Kuang Y. Global analysis of the Michaelis–Menten ratio-dependent predator–prey system. J Math Biol 2001;42:489–506. [18] Hsu SB, Hwang TW, Kuang Y. Rich dynamics of a ratio-dependent one prey two predator model. J Math Biol 2001;43:377–96. [19] Jiang W, Wei J. Bifurcation analysis in a limit cycle oscillator with delayed feedback. Chaos, Solitons & Fractals 2005;23:817–31. [20] Kuang Y. Delay differential equations with applications in population dynamics. New York: Academic Press; 1993. [21] Kuang Y, Beretta E. Global qualitative analysis of a ratio-dependent predator–prey system. J Math Biol 1998;36:389–406. [22] Lakmeche A, Arino O. Bifurcation of nontrivial periodic solutions of impulsive differential equations arising chemotherapeutic treatment. Dyn Continuous, Discrete Impulsive Syst 2000;7:265–87. [23] Lakshmikantham V, Bainov DD, Simeonov PS. Theory of impulsive differential equations. Series in modern applied mathematics, vol. 6. New Jersey: World Scientific; 1989. [24] Li S, Liao X, Li C. Hopf bifurcations in a Volterra prey–predator model with strong kernel. Chaos, Solitons & Fractals 2004;22:713–22. [25] Liu Z, Yuan R. Stability and Hopf bifurcations in a harvested one-prey-two-predator model with delays. Chaos, Solitons & Fractals 2006;27:1395–407. [26] Meng X, Wei J. Stability and Hopf bifurcations of mutual system with time delay. Chaos, Solitons & Fractals 2004;21:729–40. [27] Panetta JC. A mathematical model of periodically pulsed chemotherapy: tumor recurrence and metastasis in a competitive environment. Bull Math Biol 1996;58:425–47. [28] Roberts MG, Kao RR. The dynamics of an infectious disease in a population with birth pulses. Math Biosci 1998;149:23–36. [29] Shulgin B, Stone L, Agur Z. Pulse vaccination strategy in the SIR epidemic model. Bull Math Biol 1998;60:1123–48. [30] Song Y, Han M, Peng Y. Stability and Hopf bifurcations in a competitive Lotka–Volterra system with two delays. Chaos, Solitons & Fractals 2004;22:1139–48. [31] Song Y, Han M, Wei J. Stability and Hopf bifurcation analysis on a simplified BAM neural network with delays. Physica D 2005;200:185–204. [32] Song Y, Wei J. Bifurcation analysis for Chen’s system with delayed feedback and its application to control of chaos. Chaos, Solitons & Fractals 2004;22:75–91. [33] Sun C, Lin Y, Han M. Stability and Hopf bifurcations for an epidemic disease model with delay. Chaos, Solitons & Fractals 2006;30:204–16. [34] Sun C, Lin Y, Han M. Analysis of stability and Hopf bifurcations for an epidemic disease model with delay. Chaos, Solitons & Fractals 2006;30:204–16. [35] Sun C, Han M, Lin Y. Analysis of stability and Hopf bifurcations for a delayed logistic equation. Chaos, Solitons & Fractals 2007;31:672–82. [36] Tang S, Chen L. Density-dependent birth rate, birth pulses and their population dynamic consequences. J Math Biol 2002;44:185–99. [37] Wang F, Zhang S, Chen L, Sun L. Bifurcation and complexity of Monod type predator–prey system in a pulsed chemostat. Chaos, Solitons & Fractals 2006;27:447–58. [38] Wang L, Zou X. Hopf bifurcation in bidirectional associative memory neural network with delays: analysis and computation. J Comput Appl Math 2004;167:73–90. [39] Xu R, Chen L. Persistence and stability for a two-species ratio-dependent predator–prey system with time delay in a two-patch environment. Comput Math Appl 2000;40:577–88.
Q. Gan et al. / Chaos, Solitons and Fractals 39 (2009) 1883–1895
1895
[40] Xu R, Chen L. Persistence and global stability for n-species ratio-dependent predator–prey system with time delays. J Math Anal Appl 2002;275:27–43. [41] Zhou L, Tang Y, Hussein S. Stability and Hopf bifurcations for a delay competition diffusion system. Chaos, Solitons & Fractals 2002;14:1201–25.