- Email: [email protected]
The study of predator–prey system with defensive ability of prey and impulsive perturbations on the predator
Chaos, Solitons and Fractals 23 (2005) 631–643 www.elsevier.com/locate/chaos
The study of predator–prey system with defensive ability of prey and impulsive perturbations on the predator Shuwen Zhang a
a,b,* ,
Lingzhen Dong c, Lansun Chen
a
Department of Applied Mathematics, Dalian University of Technology, Liaoning, Dalian 116024, PR China b Institute of Biomathematics, Anshan Normal University, Liaoning Anshan 114005, PR China c Department of Mathematics, Taiyuan University of Technology, Shanxi, Taiyuan 030024, PR China Accepted 5 May 2004
Abstract Predator–prey system with non-monotonic functional response and impulsive perturbations on the predator is established. By using Floquet theorem and small amplitude perturbation skills, a locally asymptotically stable preyeradication periodic solution is obtained when the impulsive period is less than the critical value. Otherwise, if the impulsive period is larger than the critical value, the system is permanent. Further, using numerical simulation method the influences of the impulsive perturbations on the inherent oscillation are investigated. With the increasing of the impulsive value, the system displays a series of complex phenomena, which include (1) quasi-periodic oscillating, (2) period-doubling, (3) period-halfing, (4) non-unique dynamics (meaning that several attractors coexist), (5) attractor crisis and (6) chaotic bands with periodic windows. 2004 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 attached per unit time per predator as the prey density changes. Holling [1] gave three different kinds of functional response, which are monotonic in the first quadrant. But some experiments and observations indicate that the non-monotonic response occur at a level: when the nutrient concentrations reaches a high level an inhibitory effect on the specific growth rate may occur. To model such an inhibitory effect, Andrews [2] suggested a function mx pðxÞ ¼ ; a þ bx þ x2 called the Monod–Haldane function, or Holling type-IV function. Sokol and Howell [3] proposed a simplified Monod– Haldane function of the form: mx ; pðxÞ ¼ a þ x2 which describes the phenomenon of group defense whereby predation is decreased, or even prevented altogether, due to the increased ability of the prey to better defend or disguise themselves when their numbers are large enough. An example of this phenomenon is introduced by Tener [4].
*
Corresponding author. E-mail addresses: [email protected] (S. Zhang), [email protected] (L. Chen).
0960-0779/$ - see front matter 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2004.05.044
632
S. Zhang et al. / Chaos, Solitons and Fractals 23 (2005) 631–643
To study the predator–prey interaction when the prey exhibits group defense, Ruan and Xiao [5] investigated the following model 8 > < x0 ðtÞ ¼ rxðtÞ 1 xðtÞ axðtÞyðtÞ ; 2 k 1 þx ðtÞ ð1:1Þ > : y 0 ðtÞ ¼ yðtÞ lxðtÞ2 d ; a1 þx ðtÞ where xðtÞ and yðtÞ are functions of time representing population densities of prey and predator, all parameters are positive constants, k is the carrying capacity of the prey and d is the death of the predator, l is the maximum predation rate, and a1 is the so-called half-saturation constant. For system (1.1), some useful results in [5] for us will be shown in Section 2. With the model (1.1), let us take into account the possible effects of the human’s exploitative activities. In the papers [6,7], the authors assumed the human’s activities occur continuously. However, in practical situations, human activities always happen in a short time or instantaneous. For example, biological control is the using of a specially chosen living organism to control a particular pest, which a component of an integrated pest management strategy [8]. The chosen organism might be a predator [9–11], the natural enemy can be stocked in fixed moments to eradiate the pest or regulate it to densities below the threshold for economical damage. One of the first successful cases of biological control in greenhouses was the use of the parasitiod Encarsia formosa against the greenhouse whitefly Trialeurodes vaporariorum on tomaroes and cucumbers [10,12]. Probably, the reasonable model which describes this phenomena might be a predator–prey differential equation with impulsive perturbation. Besides, fishery management for another important example. There is a spate of interest in bioeconomic analysis of exploitation or renewable resources fishery management. In China, in order to make fish resource develop persistently, the harvest of such recourses have been restricted in a number of ways. It includes area of take, time of year, and the take of specific species. Time-area closures are used extensively to control human activity. These closures may be temporary or permanent. This implies that the population may be harvested and removed from the system, or stocked and added into the system continuously or instantaneously. Once human activities occur in regular pulses, studying predator–prey system with impulsive perturbations will be more significant. Although impulsive perturbations make the differential system more intractable, some impulsive systems have been recently studied in population dynamics in relation with: impulsive birth [13,14], impulsive vaccination [15,16], chemotherapeutic treatment of disease [17,18] and population ecology [19]. In this paper, with the idea of group defence and impulsive perturbations, we will study the following predator–prey system with periodic constant impulsive immigration or stock of predator. 9 8 xðtÞyðtÞ > > = > x0 ðtÞ ¼ rxðtÞ 1 xðtÞ ; 2 > k a1 þx ðtÞ > > t 6¼ nT ; > < 0 > ; ; y ðtÞ ¼ yðtÞ d þ a1lxðtÞ 2 þx ðtÞ ð1:2Þ >
> þ þ > > Þ ¼ xðnT Þ; yðnT Þ ¼ yðnT Þ þ s; xðnT > > t ¼ nT ; : X ð0þ Þ ¼ x0 ¼ ðx0 ; y 0 ÞT ; where T is the period of the impulsive immigration or stock of the predator, s > 0 is the size of immigration or stock of the predator. We obtain the conditions for extinction of prey and permanence of the system. Further, using numerical simulation method, we analyze the complexities of system with the increasing of the size of the impulsive immigration or stock, which includes (1) quasi-periodic oscillating, (2) period-doubling, (3) period-halfing, (4) nonunique dynamics (meaning that several attractor coexist), (5) attractor crisis and(6) chaotic bands with periodic windows. This paper is arranged as follows. In Section 3, some notations and lemmas are given. In Section 4, using the Floquet theory of impulsive equation and small amplitude perturbation skills, we prove the local stability of prey-eradication periodic solution and give the condition of permanence of system (1.2). In Section 5, the results of numerical analysis are shown, moreover, these results are discussed briefly.
2. Formulation of the model System (1.1) was studied completely in [5], and it is shown that the system exhibits complex dynamics behavior. Their main results are summarized in the following Lemmas.
S. Zhang et al. / Chaos, Solitons and Fractals 23 (2005) 631–643
633
System (1.1) has at most four equilibria, (0,0), ðk; 0Þ, and two interior equilibria ðx1 ; y2 Þ, ðx2 ; y2 Þ, where pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi l l2 4ad 2 x1 ; y1 ¼ r 1 ða þ x21 Þ; x1 ¼ 2d k pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi l þ l2 4ad 2 x2 ; y2 ¼ r 1 ða þ x22 Þ: x2 ¼ 2d k Lemma 2.1. If 4ad 2 < l2 < 163 ad 2 and x1 < k < x2 , then system (1.1) has a globally asymptotically stable positive equilibrium ðx1 ; y1 Þ. pffiffiffiffiffiffiffiffiffiffiffiffi 2l l2 4ad 2 , then system (1.1) has three equilibria: two hyperbolic Lemma 2.2. If l2 > 163 ad 2 and x2 > k > x3 , where x3 ¼ 2d saddles (0,0) and ðk; 0Þ and a positive unstable focus (or node) ðx1 ; y1 Þ. Moreover, system (1.1) has a unique limit cycle, which is stable. pffiffiffiffiffiffiffiffiffiffiffiffi pffiffi 2l l2 4ad 2 , then system (1.1) has four equilibriums: two Lemma 2.3. If 163 ad 2 < l2 < 18þ23 6 ad 2 and x3 < k < ld , where x3 ¼ 2d hyperbolic saddles (0,0) and ðx2 ; y2 Þ, a hyperbolic stable node ðk; 0Þ, and an unstable focus (or node) ðx1 ; y1 Þ. Moreover, there exists a unique limit cycle, which is stable. Remark 2.1. When 18þ23 very complicated.
pffiffi 6
ad 2 6 l2 and k < x3 , the dynamics of system (1.1) in the interior of the first quadrant could be
For simplicity, set x1 ¼ x, x2 ¼ a11 y. We transform system (1.2) into 9 8 0 1 ðtÞx2 ðtÞ x1 ðtÞ ¼ x1 ðtÞða bx1 ðtÞÞ x1þex > 2 ðtÞ ; = > > 1 > t 6¼ nT ; > < x0 ðtÞ ¼ x ðtÞ d þ mx1 ðtÞ ; ; 2 2 2 1þex1 ðtÞ > x ðnT þ Þ ¼ x ðnT Þ; x ðnT þ Þ ¼ x ðnT Þ þ p; ) > 1 1 2 2 > > > t ¼ nT ; : xð0þ Þ ¼ x0 ¼ ðx01 ; x02 ÞT ;
ð2:1Þ
where a ¼ r, b ¼ kr , m ¼ al1 , e ¼ a11 , p ¼ a11 s. We will consider the dynamics of system (2.1) and study the influences of the impulsive perturbation on the inherent oscillation by numerical methods, and we get some interesting results. 3. Preliminaries In this section, we will give some definitions, notations and lemmas which will be useful for our main results. Let Rþ ¼ ½0; 1Þ, R2þ ¼ fx 2 R2 jx P 0g. Denote f ¼ ðf1 ; f2 Þ the map defined by the right hand of the first two equations of system (2.1), and N the set of all non-negative integers. Let V : Rþ R2þ ! Rþ , then V is said to belong to class V0 if (1) V is continuous in ðnT ; ðn þ 1ÞT R2þ and for each x 2 R2þ ; n 2 N , limðt;yÞ!ðnT þ ;xÞ ¼ V ðnT þ ; xÞ exists. (2) V is locally Lipschitzian in x. Definition 3.1. Let V 2 V0 , then for ðt; xÞ 2 ðnT ; ðn þ 1ÞT R2þ , the upper right derivative of V ðt; xÞ with respect to the impulsive differential system (2.1) is defined as 1 Dþ V ðt; xÞ ¼ limþ sup ½V ðt þ h; x þ hf ðt; xÞÞ V ðt; xÞ: h!0 h Definition 3.2. System (2.1) is said to be permanent if there exist two positive constants m; M and T0 , such that each positive solution ðx1 ðtÞ; x2 ðtÞÞ of the system (2.1) satisfies m 6 xi ðtÞ 6 M, for all t > T0 ; i ¼ 1; 2. The solution of system (2.1) is a piecewise continuous function x : Rþ ! R2þ , xðtÞ is continuous on ðnT ; ðn þ 1ÞT , n 2 N and xðnT þ Þ ¼ limt!nT þ xðtÞ exists. The smoothness of f guarantee the global existence and uniqueness of the solution of system (2.1). For the details, it is referred to the books [20,21].
634
S. Zhang et al. / Chaos, Solitons and Fractals 23 (2005) 631–643
The following lemma is obvious. Lemma 3.1. Let xðtÞ be a solution of system (2.1) with xð0þ Þ P 0, then xðtÞ P 0 for all t P 0. And further xðtÞ > 0, for all t > 0 if xð0þ Þ > 0. And we will use the following important comparison theorem on impulsive differential equation [20]: Lemma 3.2. Suppose V 2 V0 . Assume that þ D V ðt; xÞ 6 gðt; V ðt; xÞÞ; t 6¼ nT ; V ðt; xðtþ ÞÞ 6 wn ðV ðt; xÞÞ; t ¼ nT ;
ð3:1Þ
where g : Rþ Rþ ! R is continuous in ðnT ; ðn þ 1ÞT Rþ and for u 2 Rþ ; n 2 N , limðt;yÞ!ðnT þ ;uÞ gðt; yÞ ¼ gðnT þ ; uÞ exists, wn : Rþ ! Rþ is non-decreasing. Let rðtÞ be the maximal solution of the scalar impulsive differential equation 8 0 < u ðtÞ ¼ gðt; uðtÞÞ; t 6¼ nT ; ð3:2Þ uðtþ Þ ¼ wn ðuðtÞÞ; t ¼ nT ; : uð0þ Þ ¼ u0 ; existing on ½0; 1Þ. Then V ð0þ ; x0 Þ 6 u0 implies that V ðt; xðtÞÞ 6 rðtÞ, t P 0, where xðtÞ is any solution of (2.1). Finally, we give some basic properties about the following subsystem of system (2.1). 8 0 t 6¼ nT ; < y ðtÞ ¼ dyðtÞ; yðtþ Þ ¼ yðtÞ þ p; t ¼ nT ; : yð0þ Þ ¼ y0 :
ð3:3Þ
Clearly y ðtÞ ¼
p expðdðt nT ÞÞ p ; t 2 ðnT ; ðn þ 1ÞT ; n 2 N ; y ð0þ Þ ¼ 1 expðdT Þ 1 expðdT Þ
is a positive periodic solution of system (3.3). Since p expðdtÞ þ y ðtÞ yðtÞ ¼ yð0þ Þ 1 expðdT Þ is the solution of system (3.3) with initial value y0 P 0, where t 2 ðnT ; ðn þ 1ÞT ; n 2 N , we get Lemma 3.3. For a positive periodic solution y ðtÞ of system (3.3) and every solution yðtÞ of system (3.3) with y0 P 0, we have jyðtÞ y ðtÞj ! 0, when t ! 1. Therefore, we obtain the complete expression for the prey-eradication periodic solution of system (2.1) p expðdðt nT ÞÞ : ð0; x2 ðtÞÞ ¼ 0; 1 expðdT Þ for t 2 ðnT ; ðn þ 1ÞT . 4. Extinction and permanence Firstly, we study the stability of prey-eradication periodic solution. Theorem 4.1. Let ðx1 ðtÞ; x2 ðtÞÞ be any solution of (2.1), then ð0; x2 ðtÞÞ is locally asymptotically stable provided that T < adp . Proof. The local stability of periodic solution ð0; x2 ðtÞÞ may be determined by considering the behavior of small amplitude perturbations of the solution. Define x1 ðtÞ ¼ uðtÞ; there may be written
x2 ðtÞ ¼ x2 ðtÞ þ vðtÞ:
S. Zhang et al. / Chaos, Solitons and Fractals 23 (2005) 631–643
uðtÞ vðtÞ
uð0Þ ; ¼ UðtÞ vð0Þ
635
06t < T;
where UðtÞ satisfies dU ¼ dt
a x2 ðtÞ 0 UðtÞ; mx2 ðtÞ d
and Uð0Þ ¼ I, the identity matrix. The linearization of the third and fourth equations of system (2.1) becomes
uðnT þ Þ vðnT þ Þ
¼
1 0 0 1
uðnT Þ : vðnT Þ
Hence, if both eigenvalues of 1 0 M¼ UðT Þ; 0 1 have absolute values less than one, then the periodic solution ð0; x2 ðtÞÞ is locally stable. Since all eigenvalues of M are Z T l2 ¼ exp ða x2 ðtÞÞdt ; l1 ¼ expðdT Þ < 1; 0 p . ad
According to Floquet theory of impulsive differential equation, the prey-eradication j l2 j< 1 if and only if T < solution ð0; x2 ðtÞÞ is locally stable. This completes the proof. h Remark. If the period of pulses T is more than Tmax ¼ adp , the prey-eradication solution becomes unstable. Theorem 4.2. There exists a constant M > 0, such that xi ðtÞ 6 M, i ¼ 1; 2 for each solution xðtÞ ¼ ðx1 ðtÞ; x2 ðtÞÞ of system (2.1) with all t large enough. Proof. Define V ðt; xÞ as V ðt; xÞ ¼ mx1 ðtÞ þ x2 ðtÞ: It is clear that V 2 V0 . We calculate the upper right derivative of V ðt; xÞ along a solution of system (2.1) and get the following impulsive differential equation þ D V ðtÞ þ LV ðtÞ ¼ mða þ LÞx1 bmx21 þ ðL dÞx2 ; t 6¼ nT ; ð4:1Þ t ¼ nT : V ðtþ Þ ¼ V ðtÞ þ p; Let 0 < L < d, then mða þ LÞx1 bmx21 þ ðL dÞx2 is bounded. Select L0 and L1 such that þ D V ðtÞ 6 L0 V ðtÞ þ L1 ; t 6¼ nT ; V ðtþ Þ ¼ V ðtÞ þ p; t ¼ nT : where L0 ; L1 are two positive constant. According to Lemma 3.2, we have L1 P ð1 expðnL0 T ÞÞ L1 expðL0 T Þ expðL0 ðt nT ÞÞ þ : expðL0 tÞ þ V ðtÞ 6 V ð0þ Þ L0 L0 expðL0 T Þ 1 where t 2 ðnT ; ðn þ 1ÞT . Hence lim V ðtÞ 6
t!1
L1 p expðL0 T Þ : þ L0 expðL0 T Þ 1
Therefore V ðt; xÞ is ultimately bounded. We obtain that each positive solution of system (2.1) is uniformly ultimately bounded. This completes the proof. h In the following, we investigate the permanence of system (2.1).
636
S. Zhang et al. / Chaos, Solitons and Fractals 23 (2005) 631–643
Theorem 4.3. System (2.1) is permanent if T > adp . Proof. Suppose xðtÞ is a solution of system (2.1) with x0 > 0. From Theorem 4.2 we may assume xi ðtÞ 6 M, i ¼ 1; 2 and p expðdT Þ M > ab, t P 0. Let m2 ¼ 1expðdT e2 , e2 > 0. According to Lemmas 3.2 and 3.3, we have x2 ðtÞ > m2 , for all t large Þ enough. In the following, we want to find m1 > 0 such that x1 ðtÞ P m1 for all t large enough. We will do it in the following two steps for convenience. mm3 Step 1. Since T > adp , we can select m3 > 0, e1 > 0 small enough such that m3 < minfab, p1ffieg, d ¼ 1þem 2 < d, 3 p r ¼ aT bm3 T dd e1 T > 0. We will prove there exists t1 2 ð0; 1Þ such that x1 ðt1 Þ P m3 . Otherwise, according to the above assumption, we get x02 ðtÞ 6 x2 ðtÞðd þ dÞ: ÞÞ By Lemmas 3.2 and 3.3, we get x2 ðtÞ 6 yðtÞ and yðtÞ ! yðtÞ, where yðtÞ ¼ p expððdþdÞðtnT , t 2 ðnT ; ðn þ 1ÞT , and yðtÞ is 1expððdþdÞT Þ the solution of the following equation 8 0 > < y ðtÞ ¼ yðtÞðd þ dÞ; t 6¼ nT ; yðtþ Þ ¼ yðtÞ þ p; t ¼ nT ; ð4:2Þ > : yð0þ Þ ¼ x02 > 0:
Therefore there exists a T1 > 0 such that x2 ðtÞ 6 yðtÞ 6 yðtÞ þ e1 ; x01 ðtÞ P x1 ðtÞða bm3 ðyðtÞ þ e1 ÞÞ:
ð4:3Þ
Let N1 2 N and N1 T P T1 . Integrating (4.3) on ðnT ; ðn þ 1ÞT ðn P N1 Þ, we have Z ðnþ1ÞT ða bm3 ðyðtÞ þ e1 ÞÞdt ¼ x1 ðnT Þ expðrÞ: x1 ððn þ 1ÞT Þ P x1 ðnT Þ exp nT
Then x1 ððN1 þ kÞT Þ P x1 ðN1 T Þ expðkrÞ ! 1; k ! 1; which is a contradiction to the boundedness of x1 ðtÞ. Step 2. If x1 ðtÞ P m3 ; t P t1 , then our aim is obtained. Otherwise, if x1 ðtÞ < m3 for some t > t1 , we may set t ¼ inf t P t1 fx1 ðtÞ < m3 g. We have x1 ðtÞ P m3 , for t 2 ½t; t Þ. It is easy to see x1 ðt Þ ¼ m3 since x1 ðtÞ is continuous. Suppose t 2 ½n1 T ; ðn1 þ 1ÞT Þ; n1 2 N . Select n2 ; n3 2 N such that e1 ðd þ dÞ; n2 T > T2 ¼ ln M þp expððn2 þ 1Þr1 T Þ expðn3 rÞ > 1; where r1 ¼ a bm3 M < 0. Set T ¼ n2 T þ n3 T . We claim that there must exist a t0 2 ððn1 þ 1ÞT ; ðn1 þ 1ÞT þ T such that x1 ðt0 Þ P m3 . Otherwise, x1 ðtÞ < m3 ; t 2 ððn1 þ 1ÞT ; ðn1 þ 1ÞT þ T . Considering (4.2) with yððn1 þ 1ÞT þ Þ ¼ x2 ððn1 þ 1ÞT þ Þ, we have p yðtÞ ¼ yððn1 þ 1ÞT þ Þ expððd þ dÞðt ðn1 þ 1ÞT ÞÞ þ yðtÞ 1 expððd þ dÞT Þ t 2 ðnT ; ðn þ 1ÞT ; n1 þ 1 6 n < n1 þ 1 þ n2 þ n3 . Then j yðtÞ yðtÞ j< ðM þ pÞ expððd dÞn2 T Þ < e1 ; x2 ðtÞ 6 yðtÞ 6 yðtÞ þ e1 ; for ðn1 þ n2 þ 1ÞT 6 t 6 ðn1 þ 1ÞT þ T , which implies (4.3) holds on ½ðn1 þ 1 þ n2 ÞT ; ðn1 þ 1ÞT þ T . As in step 1,we have x1 ððn1 þ 1 þ n2 þ n3 ÞT Þ P x1 ððn1 þ 1 þ n2 ÞT Þ expðn3 rÞ: There are two possible cases for t 2 ðt ; ðn1 þ 1ÞT . Case a. If x1 ðtÞ < m3 for t 2 ðt ; ðn1 þ 1ÞT , then x1 ðtÞ < m3 for all t 2 ðt ; ðn1 þ 1 þ n2 ÞT . System (2.1) gives x01 ðtÞ P x1 ðtÞða bm3 MÞ ¼ r1 x1 ðtÞ:
ð4:4Þ
S. Zhang et al. / Chaos, Solitons and Fractals 23 (2005) 631–643
637
Integrating (4.4) on ½t ; ðn1 þ 1 þ n2 ÞT yields x1 ððn1 þ 1 þ n2 ÞT Þ P m3 expðr1 ðn2 þ 1ÞT Þ; Then x1 ððn1 þ 1 þ n2 þ n3 ÞT Þ P m3 expðr1 ðn2 þ 1ÞT Þ expðn3 rÞ > m3 ; which is a contradiction. Let t ¼ inf fx1 ðtÞ P m3 g, then x1 ðtÞ ¼ m3 and (4.4) holds for t 2 ½t ; tÞ. Integrating (4.4) on tPt ½t ; tÞ yields x1 ðtÞ P x1 ðt Þ expðr1 ðt t ÞÞ P m3 expðr1 ð1 þ n2 þ n3 ÞT Þ: Let m3 expðr1 ð1 þ n2 þ n3 ÞT Þ ¼ m1 . For t > t the same argument can be continued since x1 ðtÞ P m3 ; t P t1 , hence x1 ðtÞ P m1 for all t > t1 . Case b. There exists a t0 2 ðt ; ðn1 þ 1ÞT such that x1 ðt0 Þ P m3 . Let t ¼ inf fx1 ðtÞ P m3 g, then x1 ðtÞ < m3 for t 2 ½t ; tÞ t>t and xðtÞ ¼ m3 . For t 2 ½t ; tÞ (4.4) holds. Integrating (4.4) on ½t ; tÞ, we have x1 ðtÞ P x1 ðt Þ expðr1 ðt t ÞÞ P m3 expðr1 T Þ > m1 : This process can be continued since x1 ðtÞ P m3 , and we have x1 ðtÞ P m1 for t > t1 . Thus in both cases we get x1 ðtÞ P m1 for t P t1 .This completes the proof. h
5. Numerical analysis In this section we will study the influence of impulsive perturbation p on inherent oscillation. For (A) a ¼ 10, b ¼ 0:25, d ¼ 0:1, m ¼ 0:007, e ¼ 0:001, p ¼ 0, T ¼ 8, by Lemma 2.1 we know system (2.1) has a globally asymptotically stable positive equilibrium. For (B) a ¼ 6, b ¼ 196 , d ¼ 0:1, m ¼ 0:025, e ¼ 0:01, p ¼ 0, T ¼ 8, by Lemma 2.2 we know system (2.1) has an unstable positive equilibrium and a unique stable limit cycle. For (C) a ¼ 5, b ¼ 225 , d ¼ 0:2, m ¼ 0:05, e ¼ 0:01, p ¼ 0, T ¼ 8, by Lemma 2.3 we know system (2.1) has two positive equilibrium and a unique stable limit cycle. For (D) a ¼ 10, b ¼ 6, d ¼ 0:1, m ¼ 0:95, e ¼ 0:001, p ¼ 0, T ¼ 8, by Remark 2.1 we know the dynamics of system (2.1) in the interior of the first quadrant could be very complicated. Since the corresponding continuous system (2.1) (p ¼ 0) cannot be solved explicitly and system (2.1) cannot be rewritten as equivalent difference equations, it is difficult to study them analytically. However, the influence of p may be documented by stroboscopically sampling one of the variables over a range of p values. Stroboscopic map is a special case of the Poincare map for periodically forced system or periodically pulsed system. Fixing points of the stroboscopic map correspond to periodic solutions of system (2.1) having the same period as the pulsing term; periodic points of period k about stroboscopic map correspond to entrained periodic solutions of system (2.1) having exactly k times the period of the pulsing; invariant circles correspond to quasi-periodic solutions of system (2.1); system (2.1) possibly appears chaotic (strange) attractors. For (A), the system (2.1) incarnates T -periodic solution with p increasing from 0.001 to 7. When p > 7, x1 will go extinct. Complexity doesn’t occur with respect to system (2.1) in this case. For (B) and (C), from bifurcation diagrams (Figs. 1 and 2) or called final state diagrams in [22], we can easily see that the dynamical behavior of these two cases is very complicated, which includes quasi-periodic oscillating, many chaotic bands, narrow periodic windows, wide periodic windows, and cries (the phenomenon of ‘‘crisis’’ in chaotic attractors can suddenly appear or disappear, or change size discontinuously as a parameter smoothly varies, was first extensively analogized by Grebogi et al., (see [23]). For (B), when p is sufficiently small (p < p1 0:025), system (2.1) experiences quasi-periodic oscillating (Fig. 3(a) limit cycle), and if p tends to zero, the annular region (Fig. 3(b)) will shrink to limit cycle. However, when p > p1 quasiperiodic oscillating is destroyed and 2T -periodic solution occurs (Fig. 4(a)) and is stable if p < p2 0:09. When p > p2 , it is unstable (chaotic crises happen) and it comes in chaotic area (Fig. 5) with periodic windows. When p > p3 0:8, the chaos suddenly disappear and appear T -periodic solution. Whereas the parameter p further increases, 2T -periodic solution occurs at p4 ¼ 2:8 (Fig. 6) again. For (C), when p is sufficiently small (p < p1 0:085), system (2.1) experiences quasi-periodic oscillating (Fig. 7(a) limit cycle), if p tends to zero, the annular region (Fig. 7(b)) will shrink to limit cycle. However, when p > p1 quasiperiodic oscillating will be destroyed and it comes in chaotic area (Fig. 8(a)) with periodic windows (Fig. 9). When
638
S. Zhang et al. / Chaos, Solitons and Fractals 23 (2005) 631–643
Fig. 1. Bifurcation diagrams of system (2.1) with a ¼ 6, b ¼ 196 , d ¼ 0:1, m ¼ 0:025, e ¼ 0:01, T ¼ 8. (a) and (c) x1 are plotted for p over [0.001,5], [0,0.8], (b) and (d) x2 are plotted for p over [0.001,5], [0,0.8]. Show the effect of parameter p on the dynamical behavior.
Fig. 2. Bifurcation diagrams of system (2.1) with a ¼ 5, b ¼ 225 , d ¼ 0:2, m ¼ 0:05, e ¼ 0:01, T ¼ 8. (a) and (c) x1 are plotted for p over [0.001,8], [0,0.7], (b) and (d) x2 are plotted for p over [0.001,8], [0,0.7]. Show the effect of parameter p on the dynamical behavior.
p ¼ p2 0:71 the chaos suddenly disappear and appear 2T -periodic solution (Fig. 10(a)). When p ¼ p3 2:17 the chaos abruptly appears (Fig. 8(b)) and p ¼ p4 2:25 the chaos abruptly disappears.
S. Zhang et al. / Chaos, Solitons and Fractals 23 (2005) 631–643
639
Fig. 3. Phase portraits of system (2.1). The system (2.1) are plotted with initial value x0 ¼ ð1; 2:4Þ, (a) p ¼ 0, (b) p ¼ 0:02.
Fig. 4. Periodic windows. (a) Phase portrait of 2T -periodic solution for p ¼ 0:06, (b) phase portrait of 5T -periodic solution for p ¼ 0:66.
Fig. 5. Strange attractor. (a) Phase portrait of system (2.1) for p ¼ 0:3, (b) phase portrait of system (2.1) for p ¼ 0:76:
Fig. 6. Periodic solution. (a) Phase portrait of T -periodic solution for p ¼ 1:6, (b) phase portrait of 2T periodic solution for p ¼ 2:8.
640
S. Zhang et al. / Chaos, Solitons and Fractals 23 (2005) 631–643
Fig. 7. Phase portraits of system (2.1). The system are ploted with initial value x0 ¼ ð1; 2:4Þ, (a) p ¼ 0, (b) p ¼ 0:05.
Fig. 8. Strange attractor. (a) Phase portrait of system (2.1) for p ¼ 0:61, (b) phase portrait of system (2.1) for p ¼ 2:18.
Fig. 9. Periodic windows. (a) Phase portrait of 3T -periodic solution for p ¼ 0:3, (b) phase portrait of 4T -periodic solution for p ¼ 0:47.
Fig. 10. Periodic solution. (a) Phase portrait of 2T -periodic solution for p ¼ 1:2, (b) phase portrait of 2T -periodic solution for p ¼ 2:48.
S. Zhang et al. / Chaos, Solitons and Fractals 23 (2005) 631–643
641
For (D), the bifurcation diagrams (Fig. 11) clearly show that: with p increasing from 1 to 8, the system (2.1) experiences the process of periodic doubling cascade fi chaos fi periodic halfing, which is characterized by (1) periodic doubling, (2) period halfing, (3) non-unique dynamics, (4) periodic windows. When p < p1 3:63, system (2.1) incarnates a 2T -periodic solution of the period of the impulsive perturbation and it is stable. when p > p1 , it is unstable and there is a cascade of period doubling bifurcations (Fig. 12) leading to chaos (Fig. 13). which is followed by a cascade of periodic halfing bifurcations from chaos to T -periodic solution (Fig. 14). A
Fig. 11. Bifurcation diagrams of system (2.1) with a ¼ 10, b ¼ 6, d ¼ 0:1, m ¼ 0:95, e ¼ 0:001, T ¼ 8. (a, a1) x1 are plotted for p over [1,8], [6.4,8], (b, b1) x2 are plotted for p over [1,8], [6.4,8]. Show the effect of parameter p on the dynamical behavior.
Fig. 12. Period doubling cascade. (a) Phase portrait of T -periodic solution for p ¼ 3, (b) phase portrait of 2T -periodic solution for p ¼ 4.
Fig. 13. A strange attractor. (a) Phase portrait of system (2.1) of p ¼ 6:6, (b) and (c) time series of x1 and x2 .
642
S. Zhang et al. / Chaos, Solitons and Fractals 23 (2005) 631–643
Fig. 14. Period halfing cascade. (a) Phase portrait of 8T -period solution for p ¼ 7:09, (b) phase portrait of 4T -period solution for p ¼ 7:12.
Fig. 15. Coexistence of 2T -period solution with a stranger attractor when p ¼ 7:32. (a) Solution with x0 ¼ ð1:2; 3Þ will finally tend to a 2T -period solution, (b) solution with x0 ¼ ð0:2; 7:5Þ will tend to a strange attractor.
Fig. 16. (a) Phase portrait and (b), (c) time series of solution with x0 ¼ ð0:5; 7Þ when p > adT ¼ 8. x1 ðtÞ ! 0 as t ! 0.
typical chaos oscillation is captured when p ¼ 5:71. This periodic-doubling route to chaos is the hallmark of the logistic and Ricker maps [24,25] and has been studied extensively by Mathematicians [26,27]. For the predator–prey system, chaotic behaviors are usually obtained by continuous system with periodic forcing [28,29]. Periodic halfing is the flip bifurcation in the opposite direction, which is also observed in [30,31]. However, when p ¼ 7:32, it appears that attractor is non-unique[32]: different attractors can co-exist, obviously, which one of the attractors is reached depends on the initial values. For example in Fig. 15 2T -periodic solution and a strange attractor co-exist. If p > adT ¼ 8, by Theorem 4.1 we know that the system will go extinct with x1 ðtÞ ! 0 as t ! 0 (Fig. 16).
6. Conclusion In this paper, we have investigated predator–prey system with defensive ability of prey and impulsive perturbations on the predator. Using Floquet theorem and small amplitude perturbation skills, we have proved that prey-eradication periodic solution ð0; x2 ðtÞÞ is locally asymptotically stable when the impulsive period T < Tmax ¼ adp (critical value). Otherwise, if the impulsive period T > Tmax , the prey-eradication periodic solution becomes unstable and the system is permanent. That is, we have established conditions guaranteeing the system to be permanent and driving the prey to be extinct.
S. Zhang et al. / Chaos, Solitons and Fractals 23 (2005) 631–643
643
Numerical analysis indicates that the complex dynamics of system (2.1) depends on the values of impulsive perturbations p and all parameters. By choosing impulsive perturbations p as bifurcation parameter, we have obtained bifurcation diagrams (Figs. 1, 2, 11). Figs. 1 and 2 have shown that there exists complexity for system (2.1) including quasi-periodic oscillating, many chaotic bands, narrow periodic windows, wide periodic windows and cries. Moreover, bifurcation diagram Fig. 11 displays richer structure, which consists of periodic doubling cascade, chaotic bands with periodic windows, periodic halfing cascade and non-unique dynamics. All these results show that dynamical behavior of system (2.1) becomes more complex under periodically impulsive perturbations.
Acknowledgements This Work is Supported by the National Natural Science Foundation of China (10171106).
References [1] Holling CS. The functional response of predator to prey density and its role in mimicry and population regulation. Mem Ent Sec Can 1965;45:1–60. [2] Andrews JF. A mathematical model for the continuous culture of microorganisms utilizing inhibitory substrates. Biotechnol Bioeng 1968;10:707–23. [3] Sugie J, Howell JA. Kinetics of phenol oxidation by washed cell. Biotechnol Bioeng 1980;23:2039–49. [4] Tener JS. Muskoxen. Ottawa: Queen’s Printer; 1965. [5] Ruan S, Xiao D. Global analysis in a predator–prey system with nonmonotonic functional response. SIAM J Appl Math 2001;61:1445–72. [6] Brauer F, Soudack AC. Coexistence properties of some predator–prey systems under constant rate harvesting and stocking. J Math Biol 1981;12:101–14. [7] Brauer F, Soudack AC. Constant-rate stocking of predator–prey systems. J Math Biol 1981;11:1–14. [8] Croft BA. Arthropod biological control agents and pesticides. New York: John Wiley Sons; 1990. [9] Debach P, Rosen D. Biological control by natural enemies. 2nd ed. Cambridge: Cambridge University Press; 1991. [10] Van lenteren JC. Measures of success in biological of anthrop pods by augmentation of natural enemies. In: Wratten S, Gurr G, editors. Measures of success in biological control. Dordrdcht: Kluwer Academic Publishers; 2000. p. 77–89. [11] Luff ML. The potential of predators for pest control. Agri Ecos Environ 1983;10:159–81. [12] Van Lenteren JC, Woets J. Biological and integrated pest control in greenhouses. Ann Rev Ent 1988;33:239–50. [13] Roberts MG, Kao RR. The dynamics of an infectious disease in a population with birth pulses. Math Biosci 1998;149:23–36. [14] Tang SY, Chen LS. Density-dependent birth rate,birth pulse and their population dynamic consequences. J Math Biol 2002;44:185–99. [15] Shulgin B, Stone L, Agur Z. Pulse vaccination strategy in the SIR epidemic model. Bull Math Biol 1998;60:1–26. [16] D’Onofrio A. Stability properties of pulse vaccination strategy in SEIR epidemic model. Math Comput Modell 1997;26:59–72. [17] Panetta JC. A mathematical model of periodically pulsed chemotherapy: tumor recurrence and metastasis in a competitive environment. Bull Math Biol 1996;58:425–47. [18] Lakmeche A, Arino O. Bifurcation of non trivial periodic solutions of impulsive differential equations arising chemotherapeutic treatment. Dyn Continuous Discrete Impulsive Syst 2000;7:165–87. [19] Ballinger G, Liu X. Permanence of population growth models with impulsive effects. Math Comput Modell 1997;26:59–72. [20] Lakshmikantham V, Bainov DD, Simeonov PC. Theory of impulsive differential equations. Singapore: World Scientific; 1989. [21] Bainov DD, Simeonov DD. Impulsive differential equations:periodic solutions and applications. England: Longman; 1993. [22] Davies B. Exploring chaos, theory and experiment. Reading, MA: Perseus Books; 1999. [23] Grebogi C, Ott E, Yorke JA. Crises,sudden changes in chaotic attractors and chaotic transients. Physick D 1983;7:181–200. [24] May RM. Biological population with nonoverlapping generations: stable points, stable cycles, and chaos. Science 1974;186:645–57. [25] May RM, Oster GF. Bifurcations and dynamic complexity in simple ecological models. Am Nature 1976;110:573–99. [26] Collet P, Eeckmann JP. Iterated maps of the inter val as dynamical systems. Boston: Birkhauser; 1980. [27] Venkatesan A, Parthasarathy S, Lakshmanan M. Occurrence of multiple period-doubling bifurcation route to chaos in periodically pulsed chaotic dynamical systems. Chaos, Solitons & Fractals 2003;18:891–8. [28] Vandermeer J, Stone L, Blasius B. Categories of chaos and fractal basin boundaries in forced predator–prey models. Chaos, Solitons & Fractals 2001;12:265–76. [29] Gakkhar S, Naji RK. Chaos in seasonally perturbed ratio-dependent prey–predator system. Chaos, Solitons & Fractals 2000;11:107–18. [30] Neubert MG, Caswell H. Density-dependent vital rates and their population dynamic consequences. J Math Biol 2000;41:103–21. [31] Wikan A. From chaos to chaos. An analysis of a discrete age-structured prey–predator model. J Math Biol 2001;43:471–500. [32] Tang SY, Chen LS. Chaos in functional response host–parasitoid ecosystem models. Chaos, Solitons & Fractals 2002;13:875–84.
The study of predator–prey system with defensive ability of prey and impulsive perturbations on the predator Shuwen Zhang a
a,b,* ,
Lingzhen Dong c, Lansun Chen
a
Department of Applied Mathematics, Dalian University of Technology, Liaoning, Dalian 116024, PR China b Institute of Biomathematics, Anshan Normal University, Liaoning Anshan 114005, PR China c Department of Mathematics, Taiyuan University of Technology, Shanxi, Taiyuan 030024, PR China Accepted 5 May 2004
Abstract Predator–prey system with non-monotonic functional response and impulsive perturbations on the predator is established. By using Floquet theorem and small amplitude perturbation skills, a locally asymptotically stable preyeradication periodic solution is obtained when the impulsive period is less than the critical value. Otherwise, if the impulsive period is larger than the critical value, the system is permanent. Further, using numerical simulation method the influences of the impulsive perturbations on the inherent oscillation are investigated. With the increasing of the impulsive value, the system displays a series of complex phenomena, which include (1) quasi-periodic oscillating, (2) period-doubling, (3) period-halfing, (4) non-unique dynamics (meaning that several attractors coexist), (5) attractor crisis and (6) chaotic bands with periodic windows. 2004 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 attached per unit time per predator as the prey density changes. Holling [1] gave three different kinds of functional response, which are monotonic in the first quadrant. But some experiments and observations indicate that the non-monotonic response occur at a level: when the nutrient concentrations reaches a high level an inhibitory effect on the specific growth rate may occur. To model such an inhibitory effect, Andrews [2] suggested a function mx pðxÞ ¼ ; a þ bx þ x2 called the Monod–Haldane function, or Holling type-IV function. Sokol and Howell [3] proposed a simplified Monod– Haldane function of the form: mx ; pðxÞ ¼ a þ x2 which describes the phenomenon of group defense whereby predation is decreased, or even prevented altogether, due to the increased ability of the prey to better defend or disguise themselves when their numbers are large enough. An example of this phenomenon is introduced by Tener [4].
*
Corresponding author. E-mail addresses: [email protected] (S. Zhang), [email protected] (L. Chen).
0960-0779/$ - see front matter 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2004.05.044
632
S. Zhang et al. / Chaos, Solitons and Fractals 23 (2005) 631–643
To study the predator–prey interaction when the prey exhibits group defense, Ruan and Xiao [5] investigated the following model 8 > < x0 ðtÞ ¼ rxðtÞ 1 xðtÞ axðtÞyðtÞ ; 2 k 1 þx ðtÞ ð1:1Þ > : y 0 ðtÞ ¼ yðtÞ lxðtÞ2 d ; a1 þx ðtÞ where xðtÞ and yðtÞ are functions of time representing population densities of prey and predator, all parameters are positive constants, k is the carrying capacity of the prey and d is the death of the predator, l is the maximum predation rate, and a1 is the so-called half-saturation constant. For system (1.1), some useful results in [5] for us will be shown in Section 2. With the model (1.1), let us take into account the possible effects of the human’s exploitative activities. In the papers [6,7], the authors assumed the human’s activities occur continuously. However, in practical situations, human activities always happen in a short time or instantaneous. For example, biological control is the using of a specially chosen living organism to control a particular pest, which a component of an integrated pest management strategy [8]. The chosen organism might be a predator [9–11], the natural enemy can be stocked in fixed moments to eradiate the pest or regulate it to densities below the threshold for economical damage. One of the first successful cases of biological control in greenhouses was the use of the parasitiod Encarsia formosa against the greenhouse whitefly Trialeurodes vaporariorum on tomaroes and cucumbers [10,12]. Probably, the reasonable model which describes this phenomena might be a predator–prey differential equation with impulsive perturbation. Besides, fishery management for another important example. There is a spate of interest in bioeconomic analysis of exploitation or renewable resources fishery management. In China, in order to make fish resource develop persistently, the harvest of such recourses have been restricted in a number of ways. It includes area of take, time of year, and the take of specific species. Time-area closures are used extensively to control human activity. These closures may be temporary or permanent. This implies that the population may be harvested and removed from the system, or stocked and added into the system continuously or instantaneously. Once human activities occur in regular pulses, studying predator–prey system with impulsive perturbations will be more significant. Although impulsive perturbations make the differential system more intractable, some impulsive systems have been recently studied in population dynamics in relation with: impulsive birth [13,14], impulsive vaccination [15,16], chemotherapeutic treatment of disease [17,18] and population ecology [19]. In this paper, with the idea of group defence and impulsive perturbations, we will study the following predator–prey system with periodic constant impulsive immigration or stock of predator. 9 8 xðtÞyðtÞ > > = > x0 ðtÞ ¼ rxðtÞ 1 xðtÞ ; 2 > k a1 þx ðtÞ > > t 6¼ nT ; > < 0 > ; ; y ðtÞ ¼ yðtÞ d þ a1lxðtÞ 2 þx ðtÞ ð1:2Þ >
> þ þ > > Þ ¼ xðnT Þ; yðnT Þ ¼ yðnT Þ þ s; xðnT > > t ¼ nT ; : X ð0þ Þ ¼ x0 ¼ ðx0 ; y 0 ÞT ; where T is the period of the impulsive immigration or stock of the predator, s > 0 is the size of immigration or stock of the predator. We obtain the conditions for extinction of prey and permanence of the system. Further, using numerical simulation method, we analyze the complexities of system with the increasing of the size of the impulsive immigration or stock, which includes (1) quasi-periodic oscillating, (2) period-doubling, (3) period-halfing, (4) nonunique dynamics (meaning that several attractor coexist), (5) attractor crisis and(6) chaotic bands with periodic windows. This paper is arranged as follows. In Section 3, some notations and lemmas are given. In Section 4, using the Floquet theory of impulsive equation and small amplitude perturbation skills, we prove the local stability of prey-eradication periodic solution and give the condition of permanence of system (1.2). In Section 5, the results of numerical analysis are shown, moreover, these results are discussed briefly.
2. Formulation of the model System (1.1) was studied completely in [5], and it is shown that the system exhibits complex dynamics behavior. Their main results are summarized in the following Lemmas.
S. Zhang et al. / Chaos, Solitons and Fractals 23 (2005) 631–643
633
System (1.1) has at most four equilibria, (0,0), ðk; 0Þ, and two interior equilibria ðx1 ; y2 Þ, ðx2 ; y2 Þ, where pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi l l2 4ad 2 x1 ; y1 ¼ r 1 ða þ x21 Þ; x1 ¼ 2d k pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi l þ l2 4ad 2 x2 ; y2 ¼ r 1 ða þ x22 Þ: x2 ¼ 2d k Lemma 2.1. If 4ad 2 < l2 < 163 ad 2 and x1 < k < x2 , then system (1.1) has a globally asymptotically stable positive equilibrium ðx1 ; y1 Þ. pffiffiffiffiffiffiffiffiffiffiffiffi 2l l2 4ad 2 , then system (1.1) has three equilibria: two hyperbolic Lemma 2.2. If l2 > 163 ad 2 and x2 > k > x3 , where x3 ¼ 2d saddles (0,0) and ðk; 0Þ and a positive unstable focus (or node) ðx1 ; y1 Þ. Moreover, system (1.1) has a unique limit cycle, which is stable. pffiffiffiffiffiffiffiffiffiffiffiffi pffiffi 2l l2 4ad 2 , then system (1.1) has four equilibriums: two Lemma 2.3. If 163 ad 2 < l2 < 18þ23 6 ad 2 and x3 < k < ld , where x3 ¼ 2d hyperbolic saddles (0,0) and ðx2 ; y2 Þ, a hyperbolic stable node ðk; 0Þ, and an unstable focus (or node) ðx1 ; y1 Þ. Moreover, there exists a unique limit cycle, which is stable. Remark 2.1. When 18þ23 very complicated.
pffiffi 6
ad 2 6 l2 and k < x3 , the dynamics of system (1.1) in the interior of the first quadrant could be
For simplicity, set x1 ¼ x, x2 ¼ a11 y. We transform system (1.2) into 9 8 0 1 ðtÞx2 ðtÞ x1 ðtÞ ¼ x1 ðtÞða bx1 ðtÞÞ x1þex > 2 ðtÞ ; = > > 1 > t 6¼ nT ; > < x0 ðtÞ ¼ x ðtÞ d þ mx1 ðtÞ ; ; 2 2 2 1þex1 ðtÞ > x ðnT þ Þ ¼ x ðnT Þ; x ðnT þ Þ ¼ x ðnT Þ þ p; ) > 1 1 2 2 > > > t ¼ nT ; : xð0þ Þ ¼ x0 ¼ ðx01 ; x02 ÞT ;
ð2:1Þ
where a ¼ r, b ¼ kr , m ¼ al1 , e ¼ a11 , p ¼ a11 s. We will consider the dynamics of system (2.1) and study the influences of the impulsive perturbation on the inherent oscillation by numerical methods, and we get some interesting results. 3. Preliminaries In this section, we will give some definitions, notations and lemmas which will be useful for our main results. Let Rþ ¼ ½0; 1Þ, R2þ ¼ fx 2 R2 jx P 0g. Denote f ¼ ðf1 ; f2 Þ the map defined by the right hand of the first two equations of system (2.1), and N the set of all non-negative integers. Let V : Rþ R2þ ! Rþ , then V is said to belong to class V0 if (1) V is continuous in ðnT ; ðn þ 1ÞT R2þ and for each x 2 R2þ ; n 2 N , limðt;yÞ!ðnT þ ;xÞ ¼ V ðnT þ ; xÞ exists. (2) V is locally Lipschitzian in x. Definition 3.1. Let V 2 V0 , then for ðt; xÞ 2 ðnT ; ðn þ 1ÞT R2þ , the upper right derivative of V ðt; xÞ with respect to the impulsive differential system (2.1) is defined as 1 Dþ V ðt; xÞ ¼ limþ sup ½V ðt þ h; x þ hf ðt; xÞÞ V ðt; xÞ: h!0 h Definition 3.2. System (2.1) is said to be permanent if there exist two positive constants m; M and T0 , such that each positive solution ðx1 ðtÞ; x2 ðtÞÞ of the system (2.1) satisfies m 6 xi ðtÞ 6 M, for all t > T0 ; i ¼ 1; 2. The solution of system (2.1) is a piecewise continuous function x : Rþ ! R2þ , xðtÞ is continuous on ðnT ; ðn þ 1ÞT , n 2 N and xðnT þ Þ ¼ limt!nT þ xðtÞ exists. The smoothness of f guarantee the global existence and uniqueness of the solution of system (2.1). For the details, it is referred to the books [20,21].
634
S. Zhang et al. / Chaos, Solitons and Fractals 23 (2005) 631–643
The following lemma is obvious. Lemma 3.1. Let xðtÞ be a solution of system (2.1) with xð0þ Þ P 0, then xðtÞ P 0 for all t P 0. And further xðtÞ > 0, for all t > 0 if xð0þ Þ > 0. And we will use the following important comparison theorem on impulsive differential equation [20]: Lemma 3.2. Suppose V 2 V0 . Assume that þ D V ðt; xÞ 6 gðt; V ðt; xÞÞ; t 6¼ nT ; V ðt; xðtþ ÞÞ 6 wn ðV ðt; xÞÞ; t ¼ nT ;
ð3:1Þ
where g : Rþ Rþ ! R is continuous in ðnT ; ðn þ 1ÞT Rþ and for u 2 Rþ ; n 2 N , limðt;yÞ!ðnT þ ;uÞ gðt; yÞ ¼ gðnT þ ; uÞ exists, wn : Rþ ! Rþ is non-decreasing. Let rðtÞ be the maximal solution of the scalar impulsive differential equation 8 0 < u ðtÞ ¼ gðt; uðtÞÞ; t 6¼ nT ; ð3:2Þ uðtþ Þ ¼ wn ðuðtÞÞ; t ¼ nT ; : uð0þ Þ ¼ u0 ; existing on ½0; 1Þ. Then V ð0þ ; x0 Þ 6 u0 implies that V ðt; xðtÞÞ 6 rðtÞ, t P 0, where xðtÞ is any solution of (2.1). Finally, we give some basic properties about the following subsystem of system (2.1). 8 0 t 6¼ nT ; < y ðtÞ ¼ dyðtÞ; yðtþ Þ ¼ yðtÞ þ p; t ¼ nT ; : yð0þ Þ ¼ y0 :
ð3:3Þ
Clearly y ðtÞ ¼
p expðdðt nT ÞÞ p ; t 2 ðnT ; ðn þ 1ÞT ; n 2 N ; y ð0þ Þ ¼ 1 expðdT Þ 1 expðdT Þ
is a positive periodic solution of system (3.3). Since p expðdtÞ þ y ðtÞ yðtÞ ¼ yð0þ Þ 1 expðdT Þ is the solution of system (3.3) with initial value y0 P 0, where t 2 ðnT ; ðn þ 1ÞT ; n 2 N , we get Lemma 3.3. For a positive periodic solution y ðtÞ of system (3.3) and every solution yðtÞ of system (3.3) with y0 P 0, we have jyðtÞ y ðtÞj ! 0, when t ! 1. Therefore, we obtain the complete expression for the prey-eradication periodic solution of system (2.1) p expðdðt nT ÞÞ : ð0; x2 ðtÞÞ ¼ 0; 1 expðdT Þ for t 2 ðnT ; ðn þ 1ÞT . 4. Extinction and permanence Firstly, we study the stability of prey-eradication periodic solution. Theorem 4.1. Let ðx1 ðtÞ; x2 ðtÞÞ be any solution of (2.1), then ð0; x2 ðtÞÞ is locally asymptotically stable provided that T < adp . Proof. The local stability of periodic solution ð0; x2 ðtÞÞ may be determined by considering the behavior of small amplitude perturbations of the solution. Define x1 ðtÞ ¼ uðtÞ; there may be written
x2 ðtÞ ¼ x2 ðtÞ þ vðtÞ:
S. Zhang et al. / Chaos, Solitons and Fractals 23 (2005) 631–643
uðtÞ vðtÞ
uð0Þ ; ¼ UðtÞ vð0Þ
635
06t < T;
where UðtÞ satisfies dU ¼ dt
a x2 ðtÞ 0 UðtÞ; mx2 ðtÞ d
and Uð0Þ ¼ I, the identity matrix. The linearization of the third and fourth equations of system (2.1) becomes
uðnT þ Þ vðnT þ Þ
¼
1 0 0 1
uðnT Þ : vðnT Þ
Hence, if both eigenvalues of 1 0 M¼ UðT Þ; 0 1 have absolute values less than one, then the periodic solution ð0; x2 ðtÞÞ is locally stable. Since all eigenvalues of M are Z T l2 ¼ exp ða x2 ðtÞÞdt ; l1 ¼ expðdT Þ < 1; 0 p . ad
According to Floquet theory of impulsive differential equation, the prey-eradication j l2 j< 1 if and only if T < solution ð0; x2 ðtÞÞ is locally stable. This completes the proof. h Remark. If the period of pulses T is more than Tmax ¼ adp , the prey-eradication solution becomes unstable. Theorem 4.2. There exists a constant M > 0, such that xi ðtÞ 6 M, i ¼ 1; 2 for each solution xðtÞ ¼ ðx1 ðtÞ; x2 ðtÞÞ of system (2.1) with all t large enough. Proof. Define V ðt; xÞ as V ðt; xÞ ¼ mx1 ðtÞ þ x2 ðtÞ: It is clear that V 2 V0 . We calculate the upper right derivative of V ðt; xÞ along a solution of system (2.1) and get the following impulsive differential equation þ D V ðtÞ þ LV ðtÞ ¼ mða þ LÞx1 bmx21 þ ðL dÞx2 ; t 6¼ nT ; ð4:1Þ t ¼ nT : V ðtþ Þ ¼ V ðtÞ þ p; Let 0 < L < d, then mða þ LÞx1 bmx21 þ ðL dÞx2 is bounded. Select L0 and L1 such that þ D V ðtÞ 6 L0 V ðtÞ þ L1 ; t 6¼ nT ; V ðtþ Þ ¼ V ðtÞ þ p; t ¼ nT : where L0 ; L1 are two positive constant. According to Lemma 3.2, we have L1 P ð1 expðnL0 T ÞÞ L1 expðL0 T Þ expðL0 ðt nT ÞÞ þ : expðL0 tÞ þ V ðtÞ 6 V ð0þ Þ L0 L0 expðL0 T Þ 1 where t 2 ðnT ; ðn þ 1ÞT . Hence lim V ðtÞ 6
t!1
L1 p expðL0 T Þ : þ L0 expðL0 T Þ 1
Therefore V ðt; xÞ is ultimately bounded. We obtain that each positive solution of system (2.1) is uniformly ultimately bounded. This completes the proof. h In the following, we investigate the permanence of system (2.1).
636
S. Zhang et al. / Chaos, Solitons and Fractals 23 (2005) 631–643
Theorem 4.3. System (2.1) is permanent if T > adp . Proof. Suppose xðtÞ is a solution of system (2.1) with x0 > 0. From Theorem 4.2 we may assume xi ðtÞ 6 M, i ¼ 1; 2 and p expðdT Þ M > ab, t P 0. Let m2 ¼ 1expðdT e2 , e2 > 0. According to Lemmas 3.2 and 3.3, we have x2 ðtÞ > m2 , for all t large Þ enough. In the following, we want to find m1 > 0 such that x1 ðtÞ P m1 for all t large enough. We will do it in the following two steps for convenience. mm3 Step 1. Since T > adp , we can select m3 > 0, e1 > 0 small enough such that m3 < minfab, p1ffieg, d ¼ 1þem 2 < d, 3 p r ¼ aT bm3 T dd e1 T > 0. We will prove there exists t1 2 ð0; 1Þ such that x1 ðt1 Þ P m3 . Otherwise, according to the above assumption, we get x02 ðtÞ 6 x2 ðtÞðd þ dÞ: ÞÞ By Lemmas 3.2 and 3.3, we get x2 ðtÞ 6 yðtÞ and yðtÞ ! yðtÞ, where yðtÞ ¼ p expððdþdÞðtnT , t 2 ðnT ; ðn þ 1ÞT , and yðtÞ is 1expððdþdÞT Þ the solution of the following equation 8 0 > < y ðtÞ ¼ yðtÞðd þ dÞ; t 6¼ nT ; yðtþ Þ ¼ yðtÞ þ p; t ¼ nT ; ð4:2Þ > : yð0þ Þ ¼ x02 > 0:
Therefore there exists a T1 > 0 such that x2 ðtÞ 6 yðtÞ 6 yðtÞ þ e1 ; x01 ðtÞ P x1 ðtÞða bm3 ðyðtÞ þ e1 ÞÞ:
ð4:3Þ
Let N1 2 N and N1 T P T1 . Integrating (4.3) on ðnT ; ðn þ 1ÞT ðn P N1 Þ, we have Z ðnþ1ÞT ða bm3 ðyðtÞ þ e1 ÞÞdt ¼ x1 ðnT Þ expðrÞ: x1 ððn þ 1ÞT Þ P x1 ðnT Þ exp nT
Then x1 ððN1 þ kÞT Þ P x1 ðN1 T Þ expðkrÞ ! 1; k ! 1; which is a contradiction to the boundedness of x1 ðtÞ. Step 2. If x1 ðtÞ P m3 ; t P t1 , then our aim is obtained. Otherwise, if x1 ðtÞ < m3 for some t > t1 , we may set t ¼ inf t P t1 fx1 ðtÞ < m3 g. We have x1 ðtÞ P m3 , for t 2 ½t; t Þ. It is easy to see x1 ðt Þ ¼ m3 since x1 ðtÞ is continuous. Suppose t 2 ½n1 T ; ðn1 þ 1ÞT Þ; n1 2 N . Select n2 ; n3 2 N such that e1 ðd þ dÞ; n2 T > T2 ¼ ln M þp expððn2 þ 1Þr1 T Þ expðn3 rÞ > 1; where r1 ¼ a bm3 M < 0. Set T ¼ n2 T þ n3 T . We claim that there must exist a t0 2 ððn1 þ 1ÞT ; ðn1 þ 1ÞT þ T such that x1 ðt0 Þ P m3 . Otherwise, x1 ðtÞ < m3 ; t 2 ððn1 þ 1ÞT ; ðn1 þ 1ÞT þ T . Considering (4.2) with yððn1 þ 1ÞT þ Þ ¼ x2 ððn1 þ 1ÞT þ Þ, we have p yðtÞ ¼ yððn1 þ 1ÞT þ Þ expððd þ dÞðt ðn1 þ 1ÞT ÞÞ þ yðtÞ 1 expððd þ dÞT Þ t 2 ðnT ; ðn þ 1ÞT ; n1 þ 1 6 n < n1 þ 1 þ n2 þ n3 . Then j yðtÞ yðtÞ j< ðM þ pÞ expððd dÞn2 T Þ < e1 ; x2 ðtÞ 6 yðtÞ 6 yðtÞ þ e1 ; for ðn1 þ n2 þ 1ÞT 6 t 6 ðn1 þ 1ÞT þ T , which implies (4.3) holds on ½ðn1 þ 1 þ n2 ÞT ; ðn1 þ 1ÞT þ T . As in step 1,we have x1 ððn1 þ 1 þ n2 þ n3 ÞT Þ P x1 ððn1 þ 1 þ n2 ÞT Þ expðn3 rÞ: There are two possible cases for t 2 ðt ; ðn1 þ 1ÞT . Case a. If x1 ðtÞ < m3 for t 2 ðt ; ðn1 þ 1ÞT , then x1 ðtÞ < m3 for all t 2 ðt ; ðn1 þ 1 þ n2 ÞT . System (2.1) gives x01 ðtÞ P x1 ðtÞða bm3 MÞ ¼ r1 x1 ðtÞ:
ð4:4Þ
S. Zhang et al. / Chaos, Solitons and Fractals 23 (2005) 631–643
637
Integrating (4.4) on ½t ; ðn1 þ 1 þ n2 ÞT yields x1 ððn1 þ 1 þ n2 ÞT Þ P m3 expðr1 ðn2 þ 1ÞT Þ; Then x1 ððn1 þ 1 þ n2 þ n3 ÞT Þ P m3 expðr1 ðn2 þ 1ÞT Þ expðn3 rÞ > m3 ; which is a contradiction. Let t ¼ inf fx1 ðtÞ P m3 g, then x1 ðtÞ ¼ m3 and (4.4) holds for t 2 ½t ; tÞ. Integrating (4.4) on tPt ½t ; tÞ yields x1 ðtÞ P x1 ðt Þ expðr1 ðt t ÞÞ P m3 expðr1 ð1 þ n2 þ n3 ÞT Þ: Let m3 expðr1 ð1 þ n2 þ n3 ÞT Þ ¼ m1 . For t > t the same argument can be continued since x1 ðtÞ P m3 ; t P t1 , hence x1 ðtÞ P m1 for all t > t1 . Case b. There exists a t0 2 ðt ; ðn1 þ 1ÞT such that x1 ðt0 Þ P m3 . Let t ¼ inf fx1 ðtÞ P m3 g, then x1 ðtÞ < m3 for t 2 ½t ; tÞ t>t and xðtÞ ¼ m3 . For t 2 ½t ; tÞ (4.4) holds. Integrating (4.4) on ½t ; tÞ, we have x1 ðtÞ P x1 ðt Þ expðr1 ðt t ÞÞ P m3 expðr1 T Þ > m1 : This process can be continued since x1 ðtÞ P m3 , and we have x1 ðtÞ P m1 for t > t1 . Thus in both cases we get x1 ðtÞ P m1 for t P t1 .This completes the proof. h
5. Numerical analysis In this section we will study the influence of impulsive perturbation p on inherent oscillation. For (A) a ¼ 10, b ¼ 0:25, d ¼ 0:1, m ¼ 0:007, e ¼ 0:001, p ¼ 0, T ¼ 8, by Lemma 2.1 we know system (2.1) has a globally asymptotically stable positive equilibrium. For (B) a ¼ 6, b ¼ 196 , d ¼ 0:1, m ¼ 0:025, e ¼ 0:01, p ¼ 0, T ¼ 8, by Lemma 2.2 we know system (2.1) has an unstable positive equilibrium and a unique stable limit cycle. For (C) a ¼ 5, b ¼ 225 , d ¼ 0:2, m ¼ 0:05, e ¼ 0:01, p ¼ 0, T ¼ 8, by Lemma 2.3 we know system (2.1) has two positive equilibrium and a unique stable limit cycle. For (D) a ¼ 10, b ¼ 6, d ¼ 0:1, m ¼ 0:95, e ¼ 0:001, p ¼ 0, T ¼ 8, by Remark 2.1 we know the dynamics of system (2.1) in the interior of the first quadrant could be very complicated. Since the corresponding continuous system (2.1) (p ¼ 0) cannot be solved explicitly and system (2.1) cannot be rewritten as equivalent difference equations, it is difficult to study them analytically. However, the influence of p may be documented by stroboscopically sampling one of the variables over a range of p values. Stroboscopic map is a special case of the Poincare map for periodically forced system or periodically pulsed system. Fixing points of the stroboscopic map correspond to periodic solutions of system (2.1) having the same period as the pulsing term; periodic points of period k about stroboscopic map correspond to entrained periodic solutions of system (2.1) having exactly k times the period of the pulsing; invariant circles correspond to quasi-periodic solutions of system (2.1); system (2.1) possibly appears chaotic (strange) attractors. For (A), the system (2.1) incarnates T -periodic solution with p increasing from 0.001 to 7. When p > 7, x1 will go extinct. Complexity doesn’t occur with respect to system (2.1) in this case. For (B) and (C), from bifurcation diagrams (Figs. 1 and 2) or called final state diagrams in [22], we can easily see that the dynamical behavior of these two cases is very complicated, which includes quasi-periodic oscillating, many chaotic bands, narrow periodic windows, wide periodic windows, and cries (the phenomenon of ‘‘crisis’’ in chaotic attractors can suddenly appear or disappear, or change size discontinuously as a parameter smoothly varies, was first extensively analogized by Grebogi et al., (see [23]). For (B), when p is sufficiently small (p < p1 0:025), system (2.1) experiences quasi-periodic oscillating (Fig. 3(a) limit cycle), and if p tends to zero, the annular region (Fig. 3(b)) will shrink to limit cycle. However, when p > p1 quasiperiodic oscillating is destroyed and 2T -periodic solution occurs (Fig. 4(a)) and is stable if p < p2 0:09. When p > p2 , it is unstable (chaotic crises happen) and it comes in chaotic area (Fig. 5) with periodic windows. When p > p3 0:8, the chaos suddenly disappear and appear T -periodic solution. Whereas the parameter p further increases, 2T -periodic solution occurs at p4 ¼ 2:8 (Fig. 6) again. For (C), when p is sufficiently small (p < p1 0:085), system (2.1) experiences quasi-periodic oscillating (Fig. 7(a) limit cycle), if p tends to zero, the annular region (Fig. 7(b)) will shrink to limit cycle. However, when p > p1 quasiperiodic oscillating will be destroyed and it comes in chaotic area (Fig. 8(a)) with periodic windows (Fig. 9). When
638
S. Zhang et al. / Chaos, Solitons and Fractals 23 (2005) 631–643
Fig. 1. Bifurcation diagrams of system (2.1) with a ¼ 6, b ¼ 196 , d ¼ 0:1, m ¼ 0:025, e ¼ 0:01, T ¼ 8. (a) and (c) x1 are plotted for p over [0.001,5], [0,0.8], (b) and (d) x2 are plotted for p over [0.001,5], [0,0.8]. Show the effect of parameter p on the dynamical behavior.
Fig. 2. Bifurcation diagrams of system (2.1) with a ¼ 5, b ¼ 225 , d ¼ 0:2, m ¼ 0:05, e ¼ 0:01, T ¼ 8. (a) and (c) x1 are plotted for p over [0.001,8], [0,0.7], (b) and (d) x2 are plotted for p over [0.001,8], [0,0.7]. Show the effect of parameter p on the dynamical behavior.
p ¼ p2 0:71 the chaos suddenly disappear and appear 2T -periodic solution (Fig. 10(a)). When p ¼ p3 2:17 the chaos abruptly appears (Fig. 8(b)) and p ¼ p4 2:25 the chaos abruptly disappears.
S. Zhang et al. / Chaos, Solitons and Fractals 23 (2005) 631–643
639
Fig. 3. Phase portraits of system (2.1). The system (2.1) are plotted with initial value x0 ¼ ð1; 2:4Þ, (a) p ¼ 0, (b) p ¼ 0:02.
Fig. 4. Periodic windows. (a) Phase portrait of 2T -periodic solution for p ¼ 0:06, (b) phase portrait of 5T -periodic solution for p ¼ 0:66.
Fig. 5. Strange attractor. (a) Phase portrait of system (2.1) for p ¼ 0:3, (b) phase portrait of system (2.1) for p ¼ 0:76:
Fig. 6. Periodic solution. (a) Phase portrait of T -periodic solution for p ¼ 1:6, (b) phase portrait of 2T periodic solution for p ¼ 2:8.
640
S. Zhang et al. / Chaos, Solitons and Fractals 23 (2005) 631–643
Fig. 7. Phase portraits of system (2.1). The system are ploted with initial value x0 ¼ ð1; 2:4Þ, (a) p ¼ 0, (b) p ¼ 0:05.
Fig. 8. Strange attractor. (a) Phase portrait of system (2.1) for p ¼ 0:61, (b) phase portrait of system (2.1) for p ¼ 2:18.
Fig. 9. Periodic windows. (a) Phase portrait of 3T -periodic solution for p ¼ 0:3, (b) phase portrait of 4T -periodic solution for p ¼ 0:47.
Fig. 10. Periodic solution. (a) Phase portrait of 2T -periodic solution for p ¼ 1:2, (b) phase portrait of 2T -periodic solution for p ¼ 2:48.
S. Zhang et al. / Chaos, Solitons and Fractals 23 (2005) 631–643
641
For (D), the bifurcation diagrams (Fig. 11) clearly show that: with p increasing from 1 to 8, the system (2.1) experiences the process of periodic doubling cascade fi chaos fi periodic halfing, which is characterized by (1) periodic doubling, (2) period halfing, (3) non-unique dynamics, (4) periodic windows. When p < p1 3:63, system (2.1) incarnates a 2T -periodic solution of the period of the impulsive perturbation and it is stable. when p > p1 , it is unstable and there is a cascade of period doubling bifurcations (Fig. 12) leading to chaos (Fig. 13). which is followed by a cascade of periodic halfing bifurcations from chaos to T -periodic solution (Fig. 14). A
Fig. 11. Bifurcation diagrams of system (2.1) with a ¼ 10, b ¼ 6, d ¼ 0:1, m ¼ 0:95, e ¼ 0:001, T ¼ 8. (a, a1) x1 are plotted for p over [1,8], [6.4,8], (b, b1) x2 are plotted for p over [1,8], [6.4,8]. Show the effect of parameter p on the dynamical behavior.
Fig. 12. Period doubling cascade. (a) Phase portrait of T -periodic solution for p ¼ 3, (b) phase portrait of 2T -periodic solution for p ¼ 4.
Fig. 13. A strange attractor. (a) Phase portrait of system (2.1) of p ¼ 6:6, (b) and (c) time series of x1 and x2 .
642
S. Zhang et al. / Chaos, Solitons and Fractals 23 (2005) 631–643
Fig. 14. Period halfing cascade. (a) Phase portrait of 8T -period solution for p ¼ 7:09, (b) phase portrait of 4T -period solution for p ¼ 7:12.
Fig. 15. Coexistence of 2T -period solution with a stranger attractor when p ¼ 7:32. (a) Solution with x0 ¼ ð1:2; 3Þ will finally tend to a 2T -period solution, (b) solution with x0 ¼ ð0:2; 7:5Þ will tend to a strange attractor.
Fig. 16. (a) Phase portrait and (b), (c) time series of solution with x0 ¼ ð0:5; 7Þ when p > adT ¼ 8. x1 ðtÞ ! 0 as t ! 0.
typical chaos oscillation is captured when p ¼ 5:71. This periodic-doubling route to chaos is the hallmark of the logistic and Ricker maps [24,25] and has been studied extensively by Mathematicians [26,27]. For the predator–prey system, chaotic behaviors are usually obtained by continuous system with periodic forcing [28,29]. Periodic halfing is the flip bifurcation in the opposite direction, which is also observed in [30,31]. However, when p ¼ 7:32, it appears that attractor is non-unique[32]: different attractors can co-exist, obviously, which one of the attractors is reached depends on the initial values. For example in Fig. 15 2T -periodic solution and a strange attractor co-exist. If p > adT ¼ 8, by Theorem 4.1 we know that the system will go extinct with x1 ðtÞ ! 0 as t ! 0 (Fig. 16).
6. Conclusion In this paper, we have investigated predator–prey system with defensive ability of prey and impulsive perturbations on the predator. Using Floquet theorem and small amplitude perturbation skills, we have proved that prey-eradication periodic solution ð0; x2 ðtÞÞ is locally asymptotically stable when the impulsive period T < Tmax ¼ adp (critical value). Otherwise, if the impulsive period T > Tmax , the prey-eradication periodic solution becomes unstable and the system is permanent. That is, we have established conditions guaranteeing the system to be permanent and driving the prey to be extinct.
S. Zhang et al. / Chaos, Solitons and Fractals 23 (2005) 631–643
643
Numerical analysis indicates that the complex dynamics of system (2.1) depends on the values of impulsive perturbations p and all parameters. By choosing impulsive perturbations p as bifurcation parameter, we have obtained bifurcation diagrams (Figs. 1, 2, 11). Figs. 1 and 2 have shown that there exists complexity for system (2.1) including quasi-periodic oscillating, many chaotic bands, narrow periodic windows, wide periodic windows and cries. Moreover, bifurcation diagram Fig. 11 displays richer structure, which consists of periodic doubling cascade, chaotic bands with periodic windows, periodic halfing cascade and non-unique dynamics. All these results show that dynamical behavior of system (2.1) becomes more complex under periodically impulsive perturbations.
Acknowledgements This Work is Supported by the National Natural Science Foundation of China (10171106).
References [1] Holling CS. The functional response of predator to prey density and its role in mimicry and population regulation. Mem Ent Sec Can 1965;45:1–60. [2] Andrews JF. A mathematical model for the continuous culture of microorganisms utilizing inhibitory substrates. Biotechnol Bioeng 1968;10:707–23. [3] Sugie J, Howell JA. Kinetics of phenol oxidation by washed cell. Biotechnol Bioeng 1980;23:2039–49. [4] Tener JS. Muskoxen. Ottawa: Queen’s Printer; 1965. [5] Ruan S, Xiao D. Global analysis in a predator–prey system with nonmonotonic functional response. SIAM J Appl Math 2001;61:1445–72. [6] Brauer F, Soudack AC. Coexistence properties of some predator–prey systems under constant rate harvesting and stocking. J Math Biol 1981;12:101–14. [7] Brauer F, Soudack AC. Constant-rate stocking of predator–prey systems. J Math Biol 1981;11:1–14. [8] Croft BA. Arthropod biological control agents and pesticides. New York: John Wiley Sons; 1990. [9] Debach P, Rosen D. Biological control by natural enemies. 2nd ed. Cambridge: Cambridge University Press; 1991. [10] Van lenteren JC. Measures of success in biological of anthrop pods by augmentation of natural enemies. In: Wratten S, Gurr G, editors. Measures of success in biological control. Dordrdcht: Kluwer Academic Publishers; 2000. p. 77–89. [11] Luff ML. The potential of predators for pest control. Agri Ecos Environ 1983;10:159–81. [12] Van Lenteren JC, Woets J. Biological and integrated pest control in greenhouses. Ann Rev Ent 1988;33:239–50. [13] Roberts MG, Kao RR. The dynamics of an infectious disease in a population with birth pulses. Math Biosci 1998;149:23–36. [14] Tang SY, Chen LS. Density-dependent birth rate,birth pulse and their population dynamic consequences. J Math Biol 2002;44:185–99. [15] Shulgin B, Stone L, Agur Z. Pulse vaccination strategy in the SIR epidemic model. Bull Math Biol 1998;60:1–26. [16] D’Onofrio A. Stability properties of pulse vaccination strategy in SEIR epidemic model. Math Comput Modell 1997;26:59–72. [17] Panetta JC. A mathematical model of periodically pulsed chemotherapy: tumor recurrence and metastasis in a competitive environment. Bull Math Biol 1996;58:425–47. [18] Lakmeche A, Arino O. Bifurcation of non trivial periodic solutions of impulsive differential equations arising chemotherapeutic treatment. Dyn Continuous Discrete Impulsive Syst 2000;7:165–87. [19] Ballinger G, Liu X. Permanence of population growth models with impulsive effects. Math Comput Modell 1997;26:59–72. [20] Lakshmikantham V, Bainov DD, Simeonov PC. Theory of impulsive differential equations. Singapore: World Scientific; 1989. [21] Bainov DD, Simeonov DD. Impulsive differential equations:periodic solutions and applications. England: Longman; 1993. [22] Davies B. Exploring chaos, theory and experiment. Reading, MA: Perseus Books; 1999. [23] Grebogi C, Ott E, Yorke JA. Crises,sudden changes in chaotic attractors and chaotic transients. Physick D 1983;7:181–200. [24] May RM. Biological population with nonoverlapping generations: stable points, stable cycles, and chaos. Science 1974;186:645–57. [25] May RM, Oster GF. Bifurcations and dynamic complexity in simple ecological models. Am Nature 1976;110:573–99. [26] Collet P, Eeckmann JP. Iterated maps of the inter val as dynamical systems. Boston: Birkhauser; 1980. [27] Venkatesan A, Parthasarathy S, Lakshmanan M. Occurrence of multiple period-doubling bifurcation route to chaos in periodically pulsed chaotic dynamical systems. Chaos, Solitons & Fractals 2003;18:891–8. [28] Vandermeer J, Stone L, Blasius B. Categories of chaos and fractal basin boundaries in forced predator–prey models. Chaos, Solitons & Fractals 2001;12:265–76. [29] Gakkhar S, Naji RK. Chaos in seasonally perturbed ratio-dependent prey–predator system. Chaos, Solitons & Fractals 2000;11:107–18. [30] Neubert MG, Caswell H. Density-dependent vital rates and their population dynamic consequences. J Math Biol 2000;41:103–21. [31] Wikan A. From chaos to chaos. An analysis of a discrete age-structured prey–predator model. J Math Biol 2001;43:471–500. [32] Tang SY, Chen LS. Chaos in functional response host–parasitoid ecosystem models. Chaos, Solitons & Fractals 2002;13:875–84.