- Email: [email protected]
The dynamics of a Beddington-type system with impulsive control strategy
Chaos, Solitons and Fractals 29 (2006) 1229–1239 www.elsevier.com/locate/chaos
The dynamics of a Beddington-type system with impulsive control strategy q Zhenqing Li
a,*
, Weiming Wang b, Hailing Wang
c
a
Laboratory of Quantitative Vegetation Ecology, Institute of Botany, The Chinese Academy of Sciences, Beijing 100093, China b Software Engineering Institute of East China Normal University, Shanghai 200062, China c School of Applied Mathematics, University of Electronic Science and Technology of China, Chengdu, Sichuan 610054, China Accepted 24 August 2005
Abstract In this paper, by using the theories and methods of ecology and ordinary differential equation, a prey–predator system with Beddington-type functional response and impulsive control strategy is established. Conditions for the system to be extinct are given by using the theories of impulsive equation and small amplitude perturbation skills. It is proved that the system is permanent via the method of comparison involving multiple Liapunov functions. Furthermore, by using the method of numerical simulation, the influence of the impulsive control strategy on the inherent oscillation are investigated, which shows rich dynamics, such as period doubling bifurcation, crises, symmetry-breaking pitchfork bifurcations, chaotic bands, quasi-periodic oscillation, narrow periodic window, wide periodic window, period-halving bifurcation, etc. That will be useful for study of the dynamic complexity of ecosystems. 2005 Elsevier Ltd. All rights reserved.
1. Introduction Many evolution processes are characterized by the fact at certain moments of time when they experience a change of state abruptly. These processes are subject to short-term perturbations whose duration is negligible in comparison with the duration of the process. Consequently, it is natural to assume that these perturbations act instantaneously, that is, in the form of impulse. It is well known, for example, that many biological phenomena involving thresholds, bursting rhythm models in medicine and biology, optimal control models in economics, pharmacokinetics and frequency modulated systems do exhibit impulsive effects. Thus impulsive differential equations, differential equations involving impulse effects, appear as a natural description of observed evolution phenomena of several real world problems [1]. System with impulse at fixed time is one of the simplest population impulsive differential equations. It may be described as follows:
q
This research was supported by the National Natural Science Foundation of China (Granted No. 30270249) and the Knowledge Innovation Project of the Chinese Academy of Sciences (Granted No. KSCX2-SW-107). * Corresponding author. E-mail addresses: [email protected] (Z. Li), [email protected] (W. Wang). 0960-0779/$ - see front matter 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2005.08.195
1230
Z. Li et al. / Chaos, Solitons and Fractals 29 (2006) 1229–1239
x_ ¼ f ðt; xÞ; Dxðtk Þ ¼ bk xðtk Þ;
t 6¼ tk ; t ¼ tk ; k ¼ 0; 1; 2 . . . ;
ð1:1Þ
þ xðtk Þ, the moments tk (k = 0, 1, 2, . . .) are called mowhere (t, x) 2 D R · Rn, Dxðtk Þ ¼ xðtþ k Þ xðt k Þ, xðt k Þ ¼ limt!tþ k ments of impulsive effect, and satisfy with 0 = t1 < t2 < . . . < tk < , limk!1 = 1, the initial condition is xðtþ k Þ ¼ x0 , the solution x(t) of the system (1.1) is left continuous at tk (k = 0, 1, 2, . . .), that is, xðt k Þ ¼ limh!0þ xðtk hÞ ¼ xðtk Þ. And the bk (k = 0, 1, 2, . . .) are constants, when bk > 0, they can be considered as the birth or release amount of the species, while bk < 0 as the death or capture amount of the species. The field of research of chaotic impulsive differential equations about biological control seems to be a new growing interesting area in the recent years, which some scholars have paid attention to [2–8]. In our paper, we consider a impulsive differential equation—a predator–prey model with a Beddington-type functional response [9], which have already received considerable attention [10–16], by introducing a constant periodic releasing predator and the death of the prey, what follows is the corresponding equation: 8 9 dxðtÞ xðtÞ cxðtÞyðtÞ > > > ¼ rxðtÞ 1 ;> > > dt k a þ xðtÞ þ byðtÞ = > > t 6¼ ðn þ l 1ÞT ; t 6¼ nT ; > > > > dyðtÞ kcxðtÞyðtÞ > > > ; ¼ lyðtÞ; > < dt a þ xðtÞ þ byðtÞ ð1:2Þ DxðtÞ ¼ d1 xðtÞ; > > > t ¼ ðn þ l 1ÞT ; > > DyðtÞ ¼ d2 yðtÞ; > > > > > DxðtÞ ¼ 0; > > : t ¼ nT ; DyðtÞ ¼ p;
where 0 < l 6 1, Dx(t) = x(t+) x(t), Dy(t) = y(t+) y(t), 0 < d1 6 1, (0 < d2 6 1) represents the fraction of prey (predator) which dies at t = (n + l 1)T, T is the period of the impulsive effect, p > 0 is the release amount of predator at t = nT. The rest of this paper is organized as follows: in Section 2, some notations and lemmas are given, and by using the Floquet theory of impulsive equation and small amplitude perturbation skills, we prove the global stability of preyeradication periodic solution and give the condition of permanence of system (1.2). In Section 3, the results of computer simulation are shown, and moreover, these results are discussed briefly. Finally, conclusion and remarks are given in Section 4.
2. Mathematical analysis Let R+ = [0, 1), R2þ ¼ fx 2 R2 : x P 0g, X ¼ int R2þ , N be the set of all non-negative integers. The map f = (f1, f2)T is defined by the right hand of the first two equations of system (1.2). Let V : Rþ R2þ ! Rþ , and V is said to belong to class V0 if (1) V is continuous at ððn 1ÞT ; ðn þ l 1ÞT R2þ [ ððn þ l 1ÞT ; nT R2þ , for each limðt;yÞ!ððnþl1ÞT Þþ ;xÞ V ðt; yÞ ¼ V ððn þ l 1ÞT þ ; xÞ and limðt;yÞ!ðnT þ ;xÞ V ðt; yÞ ¼ V ðnT þ ; xÞ exist. (2) V is locally Lipschitzian in x.
x 2 R2þ ,
n 2 N,
Definition 2.1. V 2 V0, then for ðt; xÞ 2 ððn 1ÞT ; ðn þ l 1ÞT R2þ [ ððn þ l 1ÞT ; nT R2þ , the upper right derivative of V(t, x) with respect to the impulsive differential system (1.2) is defined as 1 Dþ V ðt; xÞ ¼ limþ sup ½V ðt þ h; x þ hf ðt; xÞÞ V ðt; xÞ. h h!0 Definition 2.2. System (1.2) is said to be permanent if there exist positive constants m, M and T0, such that each positive solution (x(t), y(t)) of the system (1.2) satisfies with m 6 x(t) 6 M, m 6 y(t) 6 M, for all t > T0. The solution of system (1.2) is a piecewise continuous function X ðtÞ : Rþ ! R2þ , X(t) is continuous on ((n 1)T, (n + l 1)T) [ ((n + l 1)T, nT) (n 2 N, 0 6 l 6 1). Obviously the smoothness properties of f guarantee the global existence and uniqueness of solutions of system (1.2). (See, for details, article [1,17].)
Z. Li et al. / Chaos, Solitons and Fractals 29 (2006) 1229–1239
1231
It is easy to prove the following lemma. Lemma 2.1. Let X(t) be a solution of system (1.2) with X(0+) P 0, then X(t) P 0 for all t P 0. And further X(t) > 0,t P 0 if X(0+) > 0. Lemma 2.2. There exists a constant M such that x(t) 6 M, y(t) 6 M for each solution (x(t), y(t)) of (1.2) with all t large enough. Finally, we give some basic 8 dyðtÞ > > ¼ lyðtÞ; > > > < dt yðtþ Þ ¼ ð1 d2 ÞyðtÞ; > > > yðtþ Þ ¼ yðtÞ þ p; > > : yð0þ Þ ¼ y 0 .
properties about the following subsystem of system (1.2). t 6¼ ðn þ l 1ÞT ; t 6¼ nT ; t ¼ ðn þ l 1ÞT ; t ¼ nT ;
ð2:1Þ
Clearly
y ðtÞ
8 p expðlðt ðn 1ÞT ÞÞ > > < 1 ð1 d Þ expðlT Þ ;
ðn 1ÞT < t 6 ðn þ l 1ÞT ;
> pð1 d2 Þ expðlðt ðn 1ÞT ÞÞ > : ; 1 ð1 d2 Þ expðlT Þ
ðn þ l 1ÞT < t 6 nT ;
2
(n 2 N), and y ð0þ Þ ¼ y ðnT þ Þ ¼
p pð1 d2 Þ expðllT Þ ; y ðlT þ Þ ¼ 1 ð1 d2 Þ expðlT Þ 1 ð1 d2 Þ expðlT Þ
is a positive periodic solution of system (2.1). Since 8 p > ðn1Þ > > ðyð0þ Þ expðltÞ þ y ðtÞ; < ð1 d2 Þ 1 ð1 d2 Þ expðlT Þ yðtÞ ¼ > p > > ð1 d2 Þn yð0þ Þ expðltÞ þ y ðtÞ; : 1 ð1 d2 Þ expðlT Þ
ðn 1ÞT < t 6 ðn þ l 1ÞT ; ðn þ l 1ÞT < t 6 nT
is the solution of system (2.1) with initial value y0 P 0, where t 2 ððn 1ÞT ; ðn þ l 1ÞT Þ [ ððn þ l 1ÞT ; nT Þ;
n 2 N.
Then we can get: Lemma 2.3. Let y*(t) be a positive periodic solution of system (2.1) and every solution y(t) of system (2.1) 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 (0, y*(t)) of system (1.2). Now, we study the stability of prey-eradication periodic solution. Theorem 2.1. Let (x(t), y(t)) be any solution of (1.2), then (0, y*(t)) is said to be globally asymptotically stable if A1 A2 þ A3 A4 1 rT þ < ln ; lb 1 d1 where
bp expðð1 lÞlT Þ ; A1 ¼ c ln a þ expðlT Þ 1 þ d2 bpð1 d2 Þ ; A3 ¼ c ln a þ expðlT Þ 1 þ d2
bp expðlT Þ A2 ¼ c ln a þ ; expðlT Þ 1 þ d2 bpð1 d2 Þ expðð1 lÞlT Þ A4 ¼ c ln a þ . expðlT Þ 1 þ d2
Proof. The local stability of periodic solution (0, y*(t)) may be determined by considering the behavior of small amplitude perturbations of the solution. Define
1232
Z. Li et al. / Chaos, Solitons and Fractals 29 (2006) 1229–1239
xðtÞ ¼ uðtÞ;
yðtÞ ¼ y ðtÞ þ vðtÞ.
Then uðtÞ vðtÞ
¼ UðtÞ
uð0Þ ; vð0Þ
where U(t) satisfies 0 cy r dU B a þ by ¼B @ cky dt a þ by
0 l
0 6 t < T;
1 C CUðtÞ; A
and U(0) = I, the identity matrix. The linearization of system (1.2) becomes 0 1 d1 uððn þ l 1ÞT Þ uððn þ l 1ÞT þ Þ ¼ 0 1 d2 vððn þ l 1ÞT Þ vððn þ l 1ÞT þ Þ and
uðnT þ Þ vðnT þ Þ
¼
1 0 0 1
uðnT Þ vðnT Þ
.
The stability of the periodic solution (0, y*(t)) is determined by the eigenvalues of 1 0 0 1 d1 UðT Þ; H¼ 0 1 0 1 d2 which are l1 ¼ ð1 d2 Þ expðlT Þ < 1;
l2 ¼ ð1 d1 Þ exp
Z 0
T
r
cy dt ; a þ by
according to Floquet theory of impulsive differential equation [17], (0, y*(t)) is locally asymptotically stable if ju2j < 1, þA3 A4 1 that is to say, rT þ A1 A2lb < lnð1d Þ. In the following, we will prove the global attractivity. Choose a e > 0 such 1 that Z T cðy eÞ dt < 1. r n,ð1 d1 Þ exp a þ bðy eÞ 0 Note that dyðtÞ P lyðtÞ, from Lemma 2.3 and comparison theorem of impulsive equation (see [1, Theorem 3.1.1]), we dt have yðtÞ > y ðtÞ e
ð2:2Þ
for all t large enough. For simplification, we may assume that (2.2) holds for all t P 0. From (1.2) and (2.2), we can obtain 8 < dxðtÞ 6 r cðy eÞ xðtÞ; t 6¼ ðn þ l 1ÞT ; dt a þ bðy eÞ ð2:3Þ : þ t ¼ ðn þ l 1Þt; xðt Þ ¼ ð1 d1 ÞxðtÞ; which leads to
! cðy Þ dt a þ bðy eÞ ðnþl1ÞT ! Z ðnþlÞT cðy Þ ¼ xððn þ l 1ÞT Þð1 d1 Þ exp dt ¼ xððn þ l 1ÞT Þn. r a þ bðy eÞ ðnþl1ÞT
xððn þ lÞT Þ 6 xððn þ l 1ÞT þ Þ exp
Z
ðnþlÞT
r
ð2:4Þ
Hence, x((n + l)T) 6 x(lT)nn and x((n + l)T) ! 0 as n ! 1. Therefore x(t) ! 0 as n ! 1 since 0 < x(t) 6 x((n + l 1)T)(1 d1)exp(rT) for (n + l 1)T 6 (n + l)T. 0 0 Next, we prove that y(t) ! y*(t) as t ! 1. For 0 < < al kc , there must exists a T > 0 such that 0 < x(t) < , t P T . Without loss of generality, we may assume that 0 < x(t) < for all t P 0, then from system (1.2) we get
Z. Li et al. / Chaos, Solitons and Fractals 29 (2006) 1229–1239
lyðtÞ 6
dyðtÞ 6 dt
l þ
1233
kc yðtÞ. a
From Lemma 2.3 and comparison theorem we have z1(t) 6 y(t) 6 z2(t) and z1(t) ! y*(t), z2(t) ! z*(t) as t ! 1, where z1(t) and z2(t) are solutions of 8 dz1 ðtÞ > > t 6¼ ðn þ l 1ÞT ; t 6¼ nT ; > dt ¼ lz1 ðtÞ; > < þ z1 ðt Þ ¼ ð1 d1 Þz1 ðtÞ; t ¼ ðn þ l 1ÞT ; > > z ðtþ Þ ¼ z1 ðtÞ þ p; t ¼ nT ; > > : 1 þ z1 ð0 Þ ¼ yð0þ Þ; and
8 dz2 ðtÞ kc > > ¼ l þ z2 ðtÞ; > > a > < dt z2 ðtþ Þ ¼ ð1 d1 Þz2 ðtÞ; > > > z2 ðtþ Þ ¼ z2 ðtÞ þ p; > > : z2 ð0þ Þ ¼ yð0þ Þ;
t 6¼ ðn þ l 1ÞT ; t 6¼ nT ; t ¼ ðn þ l 1ÞT ; t ¼ nT ;
respectively, so 8 ðt ðn 1ÞT Þ p exp l þ kce > a > > ; > < 1 ð1 d2 Þ exp l þ kce T a z ðtÞ ¼ > > ðt ðn 1ÞT Þ pð1 d2 Þ exp l þ kc > a > ; : 1 ð1 d2 Þ exp l þ kc T a
ðn 1ÞT < t 6 ðn þ l 1ÞT ; ðn þ l 1ÞT < t 6 nT .
Then, for any e1 > 0 there exists a Te > 0 such that y ðtÞ e1 < yðtÞ < z ðtÞ þ e;
t > T~ .
Let e ! 0, we have y ðtÞ e1 < yðtÞ < y ðtÞ þ e1 for t large enough, which implies y(t) ! y*(t) as t ! 1. This completes the proof.
h
Then, we investigate the permanence of system (1.2). Theorem 2.2. System (1.2) is permanent if A1 A2 þ A3 A4 1 > ln ; rT þ 1 d1 lb where Ai (i = 1, 2, 3, 4) are the same definitions as in Theorem 2.1. Proof. Suppose X(t) = (x(t), y(t)) is any solution of (1.2) with X(0) > 0. From Lemmas 2.2, we assume that x(t) 6 M, 2 Þ expðlT Þ y(t) 6 M and M > k, t P 0. From (2.2), we have y(t) > y*(t) for all t large enough and some yðtÞ P pð1d 1expðlT Þ ,m2 for t large enough. Thus we only need to find an m1 > 0 such that x(t) > m1 for t large enough. We will prove this in the following two steps. First, let 0 < m3 < al , e1 > 0 be small enough such that kc ! rm3 T ce1 T A1 A2 þ A3 A4 þ g,ð1 d1 Þ exp rT > 1; 3 k a b l kcm a where
! 3 bp exp ð1 lÞðl kcm ÞT a ; A1 ¼ c ln a þ 3 1 þ d2 exp l kcm a ! bpð1 d2 Þ ; A3 ¼ c ln a þ 3 1 þ d2 exp l kcm a
! 3 T bp exp l kcm a ; A2 ¼ c ln a þ kcm3 exp l a 1 þ d2 0 kcm 1 ðð1lÞ l a 3 T Þ bpð1 d Þe 2 A. A4 ¼ c ln @a þ 3 T 1 þ d2 exp l kcm a
1234
Z. Li et al. / Chaos, Solitons and Fractals 29 (2006) 1229–1239
Easy to prove that x(t) < m3 cannot hold for all t P 0. Otherwise, 8 dyðtÞ kcxðtÞ kcm3 > > 6 yðtÞð lÞ 6 yðtÞð lÞ; t 6¼ ðn þ l 1ÞT ; t 6¼ nT ; < dt a a þ t ¼ ðn þ l 1ÞT ; > > yðt Þ ¼ ð1 d2 ÞyðtÞ; : yðtþ Þ ¼ yðtÞ þ p; t ¼ nT .
ð2:5Þ
Then we have yðtÞ 6 zðtÞ and zðtÞ ! ~zðtÞ (t ! 1), where z(t) is the solution of 8 kcm3 > dzðtÞ > ¼ l zðtÞ; t 6¼ ðn þ l 1ÞT ; t 6¼ nT ; > > > a < dt t ¼ ðn þ l 1ÞT ; zðtþ Þ ¼ ð1 d2 ÞzðtÞ; > > þ > zðt Þ ¼ zðtÞ þ p; t ¼ nT ; > > : þ zð0 Þ ¼ y 0 and
8 3 > ðt ðn 1ÞT Þ p exp l þ kcm > a > > < 1 ð1 d Þ exp l þ kcm3 T ; 2 a ~zðtÞ ¼ 3 > ðt ðn 1ÞT Þ pð1 d Þ exp l þ kcm > 2 a > > ; : 3 T 1 ð1 d2 Þ exp l þ kcm a
ð2:6Þ
ðn 1ÞT < t 6 ðn þ l 1ÞT ; ðn þ l 1ÞT < t 6 nT .
Therefore, there exists a T1 > 0 such that yðtÞ 6 zðtÞ < ~zðtÞ þ e1 and
8 rm3 c < dxðtÞ P r ð~zðtÞ þ e1 Þ xðtÞ; t 6¼ ðn þ l 1ÞT ; a k : dtþ t ¼ ðn þ l 1ÞT xðt Þ ¼ ð1 d1 ÞxðtÞ;
ð2:7Þ
for t > T1. Let N1 2 N and (N1 + l 1)T P T1. Integrating (2.7) on ((n + l 1)T, (n + l)T], n > N1, we can get: þ
xððn þ lÞT Þ P xððn þ l 1ÞT Þ exp
Z
ðnþlÞT
ðnþl1ÞT
¼ xððn þ l 1ÞT Þð1 d1 Þ exp
Z
rm3 c r ð~zðtÞ þ e1 Þ dt k a
ðnþlÞT
ðnþl1ÞT
!
rm3 c ð~zðtÞ þ e1 Þ dt r a k
! ¼ xððn þ l 1ÞT Þg.
ð2:8Þ
Then x((N1 + n + l)T) P x((N1 + l)T)gn ! 1 as n ! 1, which is a contradiction to the boundedness of x(t). Hence there exists a t1 > 0 such that x(t1) P m3. Second, if x(t) P m3 for all t P t1, then our aim is obtained. Hence we only need to consider those solutions which leave the region R ¼ fX ðtÞ 2 R2þ : xðtÞ < m3 g and reenter it again. Let t ¼ inf tPt1 fxðtÞ < m3 g, there are two possible cases for t*. Case I. t* = (n1 + l 1)T, n1 2 N. Then x(t) P m3 for t 2 [t1, t*] and (1 d1)m3 6 x(t*+) = (1 d1)x(t*) < m3. Select n2, n3 2 N such that e1 ln M þp ; ðn2 1ÞT > kcm3 l þ a ð1 d1 Þn2 expððn2 g1 T Þgn3 > ð1 d1 Þn2 expððn2 þ 1Þg1 T Þgn3 > 1; where g1 ¼ ðr rmk 3 ac MÞ < 0. Let T ¼ n2 T þ n3 T . We claim there must be a t2 2 ðt ; t þ T such that x(t2) > m3. Otherwise, consider (2.6) with z(t*+) = y(t*+), and we have 8 ! > p > k ðnðn1 þ1ÞÞ þ > > ð1 d2 Þ exp l þ cma 3 ðt n1 T Þ þ z ðtÞ; ðn 1ÞT < t 6 ðn þ l 1ÞT ; ðzðn1 T Þ > 3 < 1 ð1 d2 Þexp l þ kcm T a ! zðtÞ ¼ > > p > ð1 d Þn1 zðn T þ Þ 3 > exp l þ kcm ðn þ l 1ÞT < t 6 nT ; ðt n1 T Þ þ z ðtÞ; > 2 1 a : 3 T 1 ð1 d2 Þexp l þ kcm a
Z. Li et al. / Chaos, Solitons and Fractals 29 (2006) 1229–1239
1235
and n1 + 1 6 n 6 n1 + n2 + n3. Then ! kcm3 ðt n1 T Þ < e1 jzðtÞ ~zðtÞj < ðM þ pÞ exp l þ a and yðtÞ 6 zðtÞ 6 ~zðtÞ þ e1 ; ðn1 þ n2 1ÞT 6 t 6 t þ T , which implies that (2.7) holds for t þ n2 T 6 t 6 t þ T . As in the first step, we have xðt þ n2 T Þ P xðt þ n2 T Þgn3 . The first equation of system (1.2) gives 8 rm3 c < dxðtÞ P xðtÞðr MÞ; t 6¼ ðn þ l 1ÞT ; dt a k : þ xðt Þ ¼ ð1 d1 ÞxðtÞ; t ¼ ðn þ l 1ÞT .
ð2:9Þ
Integrating (2.9) on [t*, t* + n2T], we have xðt þ n2 T Þ P m3 ð1 d1 Þn2 expðn2 g1 T Þ. Thus xðt þ T Þ P m3 ð1 d1 Þn2 expðn2 g1 T Þgn3 > m3 , and this is a contradiction. Let t ¼ inf tPt fxðtÞ > m3 g, and then x(t) 6 m3, xðtÞ ¼ m3 . For t 2 ðt ; tÞ, we have xðtÞ P m3 ð1 d1 Þn2 þn3 expððn2 þ n3 Þg1 T Þ,m01 . For t > t, the same arguments can be continued since xðtÞ P m3 . Case II. t* 5 (n1 + l 1)T, n1 2 N. Then x(t) P m3 for t 2 [t1, t*] and x(t*) = m3, and suppose t 2 ððn01 þ l 1ÞT ; ðn01 þ lÞT Þ, n01 2 N. There are two possible cases for t 2 ðt ; ðn01 þ lÞT Þ. Case IIa. x(t) 6 m3 for all t 2 ðt ; ðn01 þ lÞT Þ. Similar to Case I, we can prove that there must be a t02 2 ½ðn01 þ lÞT ; 0 ðn1 þ lÞT þ T such that xðt02 Þ > m3 . Here we omit it. Let ~t ¼ inf tPt fxðtÞ > m3 g, then xð~tÞ 6 m3 . For t 2 ðt ; ~tÞ, we have xðtÞ P m3 ð1 d1 Þn2 þn3 expððn2 þ n3 þ 1Þg1 T Þ. Let m1 ,m3 ð1 d1 Þn2 þn3 expððn2 þ n3 þ 1Þg1 T Þ < m01 , so x(t) P m1 for t 2 ðt ; ~tÞ. For t > ~t, the same arguments can be continued since xð~tÞ P m3 . Case IIb. There exists a t 2 ðt ; ðn01 þ lÞT Þ such that x(t) > m3. Let ~t ¼ inf tPt fxðtÞ > m3 g, then xð~tÞ 6 m3 for t 2 ðt ; ~tÞ and xð~tÞ ¼ m3 . For t 2 ðt ; ~tÞ, (2.9) holds true. Integrating (2.9) on ðt ; ~tÞ, we have xðtÞ P xðt Þ expðg1 ðt t ÞÞ P m3 expðg1 T Þ > m1 . Since x(t) P m3 for t > ~t, the same arguments can be continued. Hence x(t) P m1 for all t P t1. The proof is completed. h
3. Computer simulation To study the dynamics of a Beddington-type system with impulsive control strategy, the solution of the system (1.2) with initial conditions in the first quadrant is obtained numerically for biologically feasible range of parametric value and the bifurcation diagram provides a summary of essential dynamical behavior of the system. As we know, the corresponding continuous system (1.2) cannot be solved explicitly. So we have to study the system (1.2) numerically integrated and research the long-term behavior of the solution in Maple. Maple, an international mathematical software, is an exchanged computer algebra system with great ability of symbolic evaluation, numerical calculation, coping with graphics, etc. Its powerful functions library and unique interior programming language provide science calculation and programming with friendly platform [18–20]. Firstly, we consider a set of parameters as the following: r = 3, k = 57, c = 0.76, k = 0.72, l = 0.4, a = 9, b = 0.005, d1 = 0.2, d2 = 0.0001, T = 8 and l = 0.09, and give some figures to help us with analysis. From Theorem 2.1 we know that the prey-eradication periodic solution (0, y*(t)) is globally asymptotically stable provided that p > pmax. A typical prey-eradication periodic solution of the system (1.2) is shown in Fig. 1, where we observe how the variable y(t) oscillates in a stable cycle. In contrast, the prey x(t) rapidly decreases to zero and pmax = 116.5544814. If the amount p of releasing species is smaller than pmax, the prey eradication solution becomes
1236
Z. Li et al. / Chaos, Solitons and Fractals 29 (2006) 1229–1239
unstable and undergoes a transcritical bifurcation, then the prey and predator can coexist on a stable positive periodic solution when p < pmax, and system (1.2) can be permanent which follows from Theorem 2.2. If the parameter is further increasing, the system will exhibit a wide variety of dynamical behavior.
Fig. 1. The dynamics of system (1.2) with impulsive control strategy with r = 3, k = 57, c = 0.76, k = 0.72, l = 0.4, a = 9, b = 0.005, d1 = 0.2, d2 = 0.0001, T = 8 and l = 0.09. when p > pmax = 116.5544814, (0, y*(t)) is globally asymptotically stable, that is to say, the prey will go extinction: (a) time series of the predator population y(t), (b) time series of the prey population, x(t) ! 0, as t ! +1 when p > pmax = 116.5544814, (c) phase portrait of system (1.2) of p = 117.0.
Fig. 2. Bifurcation diagrams of the system (1.2) shows the effects of p with initial values X(0) = (75, 7.2). (a) Prey population x and (b) predator population y, with p over [0, 120].
Fig. 3. The magnified parts of Fig. 2(a). (a) p over [0, 8], (b) p over [20, 30], (c) p over [68, 80], respectively.
Z. Li et al. / Chaos, Solitons and Fractals 29 (2006) 1229–1239
1237
Fig. 4. The magnified parts of Fig. 2(b). (a) p over [0, 8], (b) p over [20, 30], (c) p over [68, 80], respectively.
Fig. 5. Period doubling bifurcation leads to chaos of system (1.2) with initial values X(0) = (75, 7.2). (a) Phase portrait of T-period solution for p = 1.5, (b) phase portrait of 2T-period solution for p = 4.0, (c) phase portrait of 4T-period solution for p = 5.5, (d) phase portrait of 8T-period solution for p = 7.0 and (e) phase portrait of the solution enter chaotic region for p = 7.6.
Fig. 6. Symmetry-breaking pitchfork bifurcations of the system (1.2) with initial values X(0) = (75, 7.2). (a) Phase portrait of period solution for p = 8, (b) phase portrait of period solution for p = 10, (c) phase portrait of period solution for p = 18.
1238
Z. Li et al. / Chaos, Solitons and Fractals 29 (2006) 1229–1239
In Fig. 2, we have displayed typical bifurcation diagrams for the prey population x and the predator population y of the system (1.2) as p increasing from 0 to 120 with initial values X(0) = (75, 7.2), Figs. 3 and 4 are the magnified parts of Fig. 2(a) and (b), and the windows of periodic behavior are more visible. As p increases, the resulting bifurcation diagrams clearly show that system (1.2) has rich dynamics including period doubling bifurcation, crises (the phenomenon of ‘‘crisis’’ in which chaotic attractors can suddenly appear or disappear, or change size discontinuously as a parameter smoothly varies), symmetry-breaking pitchfork bifurcations, chaotic bands, quasi-periodic oscillation, narrow periodic window, wide periodic window, period-halving bifurcation, etc. Fig. 2 shows that adding p will cause the availability of abundance of chaotic regions. As the bifurcation diagrams shown, when 0 < p < 2.197, the system (1.2) has exactly a T-periodic solution, when p > 2.197, there is a period doubling leading chaos (Fig. 5). When p beyond 7.602, symmetry-breaking pitchfork bifurcations of the system (1.2) with initial values X(0) = (75, 7.2), appear (Fig. 6), which have three pitchforks. Periodic windows are intermittently scattered. The abundance of chaos increases as p increase. In course of this variety, the phenomenon of ÔcrisisÕ occurs (Fig. 7), when the p is slightly increased beyond p 25.023, the chaos suddenly disappears, and from p 41.103 to p 51.592 there are a series of crises, which constitute several types of crises, we can see that from the magnified parts of bifurcation diagrams. After these chaotic areas, that is to say, when p > 79.363, the solution of the system (1.2)
Fig. 7. Crises are shown: when the p is slightly increased beyond p 25.023, the chaos suddenly disappears, and from p 41.103 to p 51.592 there are a series of crises exist, (a) p = 23, (b) p = 28, (c) p = 50, (d) p = 60, from (a) to (b) there is a crisis that the chaos suddenly disappears, and from (c) to (d) there is a crisis that the chaos suddenly appears.
Fig. 8. Period-halving bifurcation from chaos to cycles, (a) chaos when p = 70, (b) 2T-period solution when p = 90, (c) T-period solution when p = 108.
Z. Li et al. / Chaos, Solitons and Fractals 29 (2006) 1229–1239
1239
undergoes a cascade of period-halving bifurcation from chaos to cycles (Fig. 8). Periodic halving is the flip bifurcation in the opposite direction.
4. Conclusion and remarks In this paper, we have focused on the dynamics of a Beddington-type system with impulsive control strategy. Conditions for the system to be extinct are given by using the Floquet theory of impulsive equation and small amplitude perturbation skills. It is proved that the system is permanent via the method of comparison involving multiple Liapunov functions. Using the method of numerical simulation, the influence of the impulsive control strategy on the inherent oscillation are investigated, which shows the rich dynamics, such as period doubling bifurcation, crises, symmetrybreaking pitchfork bifurcations, chaotic bands, quasi-periodic oscillation, narrow periodic window, wide periodic window, period-halving bifurcation, and that such bifurcations offer considerable potential for fascinating mathematical analysis. On the other hand, since the complexity of the ecosystem, we cannot get their analytical solution. Then, we have to solve the approximate solution of them. And now, the method of computer simulation is one of the best methods. There are some international mathematical software, such as Matlab, Maple, Mathematica, etc. They all have powerful function library and can provide science calculation and programming with friendly platform. In Maple, there is a package DEtools which helps us work with differential equations. The results indicate that computer simulation is a useful method for study of the complex dynamic systems.
References [1] Lakshmikantham V, Bainov DD, Simeonov PC. Theory of impulsive differential equations. Singapore: World Scientific; 1989. [2] Tang S, Xiao Y, Clancy D. New modelling approach concerning integrated disease control and cost-effectivity. Nonlinear Anal 2005;63:439–71. [3] Lu Z, Chi X, Chen L. Impulsive control strategies in biological control of pesticide. J Theor Popul Biol 2003;64:39–47. [4] Liu B, Zhang Y, Chen L. Dynamic complexities of a Holling I predator–prey model concerning periodic biological and chemical control. Chaos, Solitons & Fractals 2004;22:123–34. [5] Sun J, Qiao F, Wu Q. Impulsive control of a financial model. J Phys Lett A 2005;335:282–8. [6] Zhang Y, Xiu Z, Chen L. Dynamics complexity of a two-prey one-predator system with impulsive effect. Chaos, Solitons & Fractals 2005;26:131–9. [7] Ballinger G, Liu X. Permanence of population growth models with impulsive effects. Math Comput Model 1997;26:59–72. [8] Lenci S, Rega G. Periodic solutions and bifurcations in an impact inverted pendulum under impulsive excitation. Chaos, Solitons & Fractals 2000;11:2453–72. [9] Beddington JR. Mutual interference between parasites or predators and its effect on searching efficiency. J Animal Ecol 1975;44: 331–40. [10] Cantrell RS, Cosner C. On the dynamics of predator–prey models with the Beddington–DeAngelis functional response. J Math Anal Appl 2001;257:206–22. [11] Hwang TW. Global analysis of the predator–prey system with Beddington–DeAngelis functional response. J Math Anal Appl 2003;281:395–401. [12] Hwang TW. Uniqueness of limit cycles of the predator–prey system with Beddington–DeAngelis functional response. J Math Anal Appl 2004;290:113–22. [13] Liu Z, Yuan R. Stability and bifurcation in a delayed predator–prey system with Beddington–DeAngelis functional response. J Math Anal Appl 2004;296:521–37. [14] Qui Z, Yu J, Zou Y. The asymptotic behavior of a chemostat model with the Beddington–DeAngelis functional response. Math Bios 2004;187:175–87. [15] Dimitrov DT, Kojouharov VH. Complete mathematical analysis of predator–prey models with linear prey growth and Beddington–DeAngelis functional response. J Appl Math Comp 2005;162:523–38. [16] Huisman G, Deboer RJ. A formal derivation of the Beddington functional response. J Thero Biol 1997;185:389–400. [17] Bainov DD, Simeonov PS. Impulsive differential equations: asymptotic properties of the solutions. Singapore: World Scientific; 1993. [18] Wang W. Computer algebra system and symbolic computation. Lanzhou, China: Gansu Sci-Tech Press; 2004 [in Chinese]. [19] Wang W, Lian X. A new algorithm for symbolic integral with application. Appl Math Comput 2005;162:949–68. [20] Wang W, Lin C. A new algorithm for integral of trigonometric functions with mechanization. Appl Math Comput 2005;164: 71–82.
The dynamics of a Beddington-type system with impulsive control strategy q Zhenqing Li
a,*
, Weiming Wang b, Hailing Wang
c
a
Laboratory of Quantitative Vegetation Ecology, Institute of Botany, The Chinese Academy of Sciences, Beijing 100093, China b Software Engineering Institute of East China Normal University, Shanghai 200062, China c School of Applied Mathematics, University of Electronic Science and Technology of China, Chengdu, Sichuan 610054, China Accepted 24 August 2005
Abstract In this paper, by using the theories and methods of ecology and ordinary differential equation, a prey–predator system with Beddington-type functional response and impulsive control strategy is established. Conditions for the system to be extinct are given by using the theories of impulsive equation and small amplitude perturbation skills. It is proved that the system is permanent via the method of comparison involving multiple Liapunov functions. Furthermore, by using the method of numerical simulation, the influence of the impulsive control strategy on the inherent oscillation are investigated, which shows rich dynamics, such as period doubling bifurcation, crises, symmetry-breaking pitchfork bifurcations, chaotic bands, quasi-periodic oscillation, narrow periodic window, wide periodic window, period-halving bifurcation, etc. That will be useful for study of the dynamic complexity of ecosystems. 2005 Elsevier Ltd. All rights reserved.
1. Introduction Many evolution processes are characterized by the fact at certain moments of time when they experience a change of state abruptly. These processes are subject to short-term perturbations whose duration is negligible in comparison with the duration of the process. Consequently, it is natural to assume that these perturbations act instantaneously, that is, in the form of impulse. It is well known, for example, that many biological phenomena involving thresholds, bursting rhythm models in medicine and biology, optimal control models in economics, pharmacokinetics and frequency modulated systems do exhibit impulsive effects. Thus impulsive differential equations, differential equations involving impulse effects, appear as a natural description of observed evolution phenomena of several real world problems [1]. System with impulse at fixed time is one of the simplest population impulsive differential equations. It may be described as follows:
q
This research was supported by the National Natural Science Foundation of China (Granted No. 30270249) and the Knowledge Innovation Project of the Chinese Academy of Sciences (Granted No. KSCX2-SW-107). * Corresponding author. E-mail addresses: [email protected] (Z. Li), [email protected] (W. Wang). 0960-0779/$ - see front matter 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2005.08.195
1230
Z. Li et al. / Chaos, Solitons and Fractals 29 (2006) 1229–1239
x_ ¼ f ðt; xÞ; Dxðtk Þ ¼ bk xðtk Þ;
t 6¼ tk ; t ¼ tk ; k ¼ 0; 1; 2 . . . ;
ð1:1Þ
þ xðtk Þ, the moments tk (k = 0, 1, 2, . . .) are called mowhere (t, x) 2 D R · Rn, Dxðtk Þ ¼ xðtþ k Þ xðt k Þ, xðt k Þ ¼ limt!tþ k ments of impulsive effect, and satisfy with 0 = t1 < t2 < . . . < tk < , limk!1 = 1, the initial condition is xðtþ k Þ ¼ x0 , the solution x(t) of the system (1.1) is left continuous at tk (k = 0, 1, 2, . . .), that is, xðt k Þ ¼ limh!0þ xðtk hÞ ¼ xðtk Þ. And the bk (k = 0, 1, 2, . . .) are constants, when bk > 0, they can be considered as the birth or release amount of the species, while bk < 0 as the death or capture amount of the species. The field of research of chaotic impulsive differential equations about biological control seems to be a new growing interesting area in the recent years, which some scholars have paid attention to [2–8]. In our paper, we consider a impulsive differential equation—a predator–prey model with a Beddington-type functional response [9], which have already received considerable attention [10–16], by introducing a constant periodic releasing predator and the death of the prey, what follows is the corresponding equation: 8 9 dxðtÞ xðtÞ cxðtÞyðtÞ > > > ¼ rxðtÞ 1 ;> > > dt k a þ xðtÞ þ byðtÞ = > > t 6¼ ðn þ l 1ÞT ; t 6¼ nT ; > > > > dyðtÞ kcxðtÞyðtÞ > > > ; ¼ lyðtÞ; > < dt a þ xðtÞ þ byðtÞ ð1:2Þ DxðtÞ ¼ d1 xðtÞ; > > > t ¼ ðn þ l 1ÞT ; > > DyðtÞ ¼ d2 yðtÞ; > > > > > DxðtÞ ¼ 0; > > : t ¼ nT ; DyðtÞ ¼ p;
where 0 < l 6 1, Dx(t) = x(t+) x(t), Dy(t) = y(t+) y(t), 0 < d1 6 1, (0 < d2 6 1) represents the fraction of prey (predator) which dies at t = (n + l 1)T, T is the period of the impulsive effect, p > 0 is the release amount of predator at t = nT. The rest of this paper is organized as follows: in Section 2, some notations and lemmas are given, and by using the Floquet theory of impulsive equation and small amplitude perturbation skills, we prove the global stability of preyeradication periodic solution and give the condition of permanence of system (1.2). In Section 3, the results of computer simulation are shown, and moreover, these results are discussed briefly. Finally, conclusion and remarks are given in Section 4.
2. Mathematical analysis Let R+ = [0, 1), R2þ ¼ fx 2 R2 : x P 0g, X ¼ int R2þ , N be the set of all non-negative integers. The map f = (f1, f2)T is defined by the right hand of the first two equations of system (1.2). Let V : Rþ R2þ ! Rþ , and V is said to belong to class V0 if (1) V is continuous at ððn 1ÞT ; ðn þ l 1ÞT R2þ [ ððn þ l 1ÞT ; nT R2þ , for each limðt;yÞ!ððnþl1ÞT Þþ ;xÞ V ðt; yÞ ¼ V ððn þ l 1ÞT þ ; xÞ and limðt;yÞ!ðnT þ ;xÞ V ðt; yÞ ¼ V ðnT þ ; xÞ exist. (2) V is locally Lipschitzian in x.
x 2 R2þ ,
n 2 N,
Definition 2.1. V 2 V0, then for ðt; xÞ 2 ððn 1ÞT ; ðn þ l 1ÞT R2þ [ ððn þ l 1ÞT ; nT R2þ , the upper right derivative of V(t, x) with respect to the impulsive differential system (1.2) is defined as 1 Dþ V ðt; xÞ ¼ limþ sup ½V ðt þ h; x þ hf ðt; xÞÞ V ðt; xÞ. h h!0 Definition 2.2. System (1.2) is said to be permanent if there exist positive constants m, M and T0, such that each positive solution (x(t), y(t)) of the system (1.2) satisfies with m 6 x(t) 6 M, m 6 y(t) 6 M, for all t > T0. The solution of system (1.2) is a piecewise continuous function X ðtÞ : Rþ ! R2þ , X(t) is continuous on ((n 1)T, (n + l 1)T) [ ((n + l 1)T, nT) (n 2 N, 0 6 l 6 1). Obviously the smoothness properties of f guarantee the global existence and uniqueness of solutions of system (1.2). (See, for details, article [1,17].)
Z. Li et al. / Chaos, Solitons and Fractals 29 (2006) 1229–1239
1231
It is easy to prove the following lemma. Lemma 2.1. Let X(t) be a solution of system (1.2) with X(0+) P 0, then X(t) P 0 for all t P 0. And further X(t) > 0,t P 0 if X(0+) > 0. Lemma 2.2. There exists a constant M such that x(t) 6 M, y(t) 6 M for each solution (x(t), y(t)) of (1.2) with all t large enough. Finally, we give some basic 8 dyðtÞ > > ¼ lyðtÞ; > > > < dt yðtþ Þ ¼ ð1 d2 ÞyðtÞ; > > > yðtþ Þ ¼ yðtÞ þ p; > > : yð0þ Þ ¼ y 0 .
properties about the following subsystem of system (1.2). t 6¼ ðn þ l 1ÞT ; t 6¼ nT ; t ¼ ðn þ l 1ÞT ; t ¼ nT ;
ð2:1Þ
Clearly
y ðtÞ
8 p expðlðt ðn 1ÞT ÞÞ > > < 1 ð1 d Þ expðlT Þ ;
ðn 1ÞT < t 6 ðn þ l 1ÞT ;
> pð1 d2 Þ expðlðt ðn 1ÞT ÞÞ > : ; 1 ð1 d2 Þ expðlT Þ
ðn þ l 1ÞT < t 6 nT ;
2
(n 2 N), and y ð0þ Þ ¼ y ðnT þ Þ ¼
p pð1 d2 Þ expðllT Þ ; y ðlT þ Þ ¼ 1 ð1 d2 Þ expðlT Þ 1 ð1 d2 Þ expðlT Þ
is a positive periodic solution of system (2.1). Since 8 p > ðn1Þ > > ðyð0þ Þ expðltÞ þ y ðtÞ; < ð1 d2 Þ 1 ð1 d2 Þ expðlT Þ yðtÞ ¼ > p > > ð1 d2 Þn yð0þ Þ expðltÞ þ y ðtÞ; : 1 ð1 d2 Þ expðlT Þ
ðn 1ÞT < t 6 ðn þ l 1ÞT ; ðn þ l 1ÞT < t 6 nT
is the solution of system (2.1) with initial value y0 P 0, where t 2 ððn 1ÞT ; ðn þ l 1ÞT Þ [ ððn þ l 1ÞT ; nT Þ;
n 2 N.
Then we can get: Lemma 2.3. Let y*(t) be a positive periodic solution of system (2.1) and every solution y(t) of system (2.1) 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 (0, y*(t)) of system (1.2). Now, we study the stability of prey-eradication periodic solution. Theorem 2.1. Let (x(t), y(t)) be any solution of (1.2), then (0, y*(t)) is said to be globally asymptotically stable if A1 A2 þ A3 A4 1 rT þ < ln ; lb 1 d1 where
bp expðð1 lÞlT Þ ; A1 ¼ c ln a þ expðlT Þ 1 þ d2 bpð1 d2 Þ ; A3 ¼ c ln a þ expðlT Þ 1 þ d2
bp expðlT Þ A2 ¼ c ln a þ ; expðlT Þ 1 þ d2 bpð1 d2 Þ expðð1 lÞlT Þ A4 ¼ c ln a þ . expðlT Þ 1 þ d2
Proof. The local stability of periodic solution (0, y*(t)) may be determined by considering the behavior of small amplitude perturbations of the solution. Define
1232
Z. Li et al. / Chaos, Solitons and Fractals 29 (2006) 1229–1239
xðtÞ ¼ uðtÞ;
yðtÞ ¼ y ðtÞ þ vðtÞ.
Then uðtÞ vðtÞ
¼ UðtÞ
uð0Þ ; vð0Þ
where U(t) satisfies 0 cy r dU B a þ by ¼B @ cky dt a þ by
0 l
0 6 t < T;
1 C CUðtÞ; A
and U(0) = I, the identity matrix. The linearization of system (1.2) becomes 0 1 d1 uððn þ l 1ÞT Þ uððn þ l 1ÞT þ Þ ¼ 0 1 d2 vððn þ l 1ÞT Þ vððn þ l 1ÞT þ Þ and
uðnT þ Þ vðnT þ Þ
¼
1 0 0 1
uðnT Þ vðnT Þ
.
The stability of the periodic solution (0, y*(t)) is determined by the eigenvalues of 1 0 0 1 d1 UðT Þ; H¼ 0 1 0 1 d2 which are l1 ¼ ð1 d2 Þ expðlT Þ < 1;
l2 ¼ ð1 d1 Þ exp
Z 0
T
r
cy dt ; a þ by
according to Floquet theory of impulsive differential equation [17], (0, y*(t)) is locally asymptotically stable if ju2j < 1, þA3 A4 1 that is to say, rT þ A1 A2lb < lnð1d Þ. In the following, we will prove the global attractivity. Choose a e > 0 such 1 that Z T cðy eÞ dt < 1. r n,ð1 d1 Þ exp a þ bðy eÞ 0 Note that dyðtÞ P lyðtÞ, from Lemma 2.3 and comparison theorem of impulsive equation (see [1, Theorem 3.1.1]), we dt have yðtÞ > y ðtÞ e
ð2:2Þ
for all t large enough. For simplification, we may assume that (2.2) holds for all t P 0. From (1.2) and (2.2), we can obtain 8 < dxðtÞ 6 r cðy eÞ xðtÞ; t 6¼ ðn þ l 1ÞT ; dt a þ bðy eÞ ð2:3Þ : þ t ¼ ðn þ l 1Þt; xðt Þ ¼ ð1 d1 ÞxðtÞ; which leads to
! cðy Þ dt a þ bðy eÞ ðnþl1ÞT ! Z ðnþlÞT cðy Þ ¼ xððn þ l 1ÞT Þð1 d1 Þ exp dt ¼ xððn þ l 1ÞT Þn. r a þ bðy eÞ ðnþl1ÞT
xððn þ lÞT Þ 6 xððn þ l 1ÞT þ Þ exp
Z
ðnþlÞT
r
ð2:4Þ
Hence, x((n + l)T) 6 x(lT)nn and x((n + l)T) ! 0 as n ! 1. Therefore x(t) ! 0 as n ! 1 since 0 < x(t) 6 x((n + l 1)T)(1 d1)exp(rT) for (n + l 1)T 6 (n + l)T. 0 0 Next, we prove that y(t) ! y*(t) as t ! 1. For 0 < < al kc , there must exists a T > 0 such that 0 < x(t) < , t P T . Without loss of generality, we may assume that 0 < x(t) < for all t P 0, then from system (1.2) we get
Z. Li et al. / Chaos, Solitons and Fractals 29 (2006) 1229–1239
lyðtÞ 6
dyðtÞ 6 dt
l þ
1233
kc yðtÞ. a
From Lemma 2.3 and comparison theorem we have z1(t) 6 y(t) 6 z2(t) and z1(t) ! y*(t), z2(t) ! z*(t) as t ! 1, where z1(t) and z2(t) are solutions of 8 dz1 ðtÞ > > t 6¼ ðn þ l 1ÞT ; t 6¼ nT ; > dt ¼ lz1 ðtÞ; > < þ z1 ðt Þ ¼ ð1 d1 Þz1 ðtÞ; t ¼ ðn þ l 1ÞT ; > > z ðtþ Þ ¼ z1 ðtÞ þ p; t ¼ nT ; > > : 1 þ z1 ð0 Þ ¼ yð0þ Þ; and
8 dz2 ðtÞ kc > > ¼ l þ z2 ðtÞ; > > a > < dt z2 ðtþ Þ ¼ ð1 d1 Þz2 ðtÞ; > > > z2 ðtþ Þ ¼ z2 ðtÞ þ p; > > : z2 ð0þ Þ ¼ yð0þ Þ;
t 6¼ ðn þ l 1ÞT ; t 6¼ nT ; t ¼ ðn þ l 1ÞT ; t ¼ nT ;
respectively, so 8 ðt ðn 1ÞT Þ p exp l þ kce > a > > ; > < 1 ð1 d2 Þ exp l þ kce T a z ðtÞ ¼ > > ðt ðn 1ÞT Þ pð1 d2 Þ exp l þ kc > a > ; : 1 ð1 d2 Þ exp l þ kc T a
ðn 1ÞT < t 6 ðn þ l 1ÞT ; ðn þ l 1ÞT < t 6 nT .
Then, for any e1 > 0 there exists a Te > 0 such that y ðtÞ e1 < yðtÞ < z ðtÞ þ e;
t > T~ .
Let e ! 0, we have y ðtÞ e1 < yðtÞ < y ðtÞ þ e1 for t large enough, which implies y(t) ! y*(t) as t ! 1. This completes the proof.
h
Then, we investigate the permanence of system (1.2). Theorem 2.2. System (1.2) is permanent if A1 A2 þ A3 A4 1 > ln ; rT þ 1 d1 lb where Ai (i = 1, 2, 3, 4) are the same definitions as in Theorem 2.1. Proof. Suppose X(t) = (x(t), y(t)) is any solution of (1.2) with X(0) > 0. From Lemmas 2.2, we assume that x(t) 6 M, 2 Þ expðlT Þ y(t) 6 M and M > k, t P 0. From (2.2), we have y(t) > y*(t) for all t large enough and some yðtÞ P pð1d 1expðlT Þ ,m2 for t large enough. Thus we only need to find an m1 > 0 such that x(t) > m1 for t large enough. We will prove this in the following two steps. First, let 0 < m3 < al , e1 > 0 be small enough such that kc ! rm3 T ce1 T A1 A2 þ A3 A4 þ g,ð1 d1 Þ exp rT > 1; 3 k a b l kcm a where
! 3 bp exp ð1 lÞðl kcm ÞT a ; A1 ¼ c ln a þ 3 1 þ d2 exp l kcm a ! bpð1 d2 Þ ; A3 ¼ c ln a þ 3 1 þ d2 exp l kcm a
! 3 T bp exp l kcm a ; A2 ¼ c ln a þ kcm3 exp l a 1 þ d2 0 kcm 1 ðð1lÞ l a 3 T Þ bpð1 d Þe 2 A. A4 ¼ c ln @a þ 3 T 1 þ d2 exp l kcm a
1234
Z. Li et al. / Chaos, Solitons and Fractals 29 (2006) 1229–1239
Easy to prove that x(t) < m3 cannot hold for all t P 0. Otherwise, 8 dyðtÞ kcxðtÞ kcm3 > > 6 yðtÞð lÞ 6 yðtÞð lÞ; t 6¼ ðn þ l 1ÞT ; t 6¼ nT ; < dt a a þ t ¼ ðn þ l 1ÞT ; > > yðt Þ ¼ ð1 d2 ÞyðtÞ; : yðtþ Þ ¼ yðtÞ þ p; t ¼ nT .
ð2:5Þ
Then we have yðtÞ 6 zðtÞ and zðtÞ ! ~zðtÞ (t ! 1), where z(t) is the solution of 8 kcm3 > dzðtÞ > ¼ l zðtÞ; t 6¼ ðn þ l 1ÞT ; t 6¼ nT ; > > > a < dt t ¼ ðn þ l 1ÞT ; zðtþ Þ ¼ ð1 d2 ÞzðtÞ; > > þ > zðt Þ ¼ zðtÞ þ p; t ¼ nT ; > > : þ zð0 Þ ¼ y 0 and
8 3 > ðt ðn 1ÞT Þ p exp l þ kcm > a > > < 1 ð1 d Þ exp l þ kcm3 T ; 2 a ~zðtÞ ¼ 3 > ðt ðn 1ÞT Þ pð1 d Þ exp l þ kcm > 2 a > > ; : 3 T 1 ð1 d2 Þ exp l þ kcm a
ð2:6Þ
ðn 1ÞT < t 6 ðn þ l 1ÞT ; ðn þ l 1ÞT < t 6 nT .
Therefore, there exists a T1 > 0 such that yðtÞ 6 zðtÞ < ~zðtÞ þ e1 and
8 rm3 c < dxðtÞ P r ð~zðtÞ þ e1 Þ xðtÞ; t 6¼ ðn þ l 1ÞT ; a k : dtþ t ¼ ðn þ l 1ÞT xðt Þ ¼ ð1 d1 ÞxðtÞ;
ð2:7Þ
for t > T1. Let N1 2 N and (N1 + l 1)T P T1. Integrating (2.7) on ((n + l 1)T, (n + l)T], n > N1, we can get: þ
xððn þ lÞT Þ P xððn þ l 1ÞT Þ exp
Z
ðnþlÞT
ðnþl1ÞT
¼ xððn þ l 1ÞT Þð1 d1 Þ exp
Z
rm3 c r ð~zðtÞ þ e1 Þ dt k a
ðnþlÞT
ðnþl1ÞT
!
rm3 c ð~zðtÞ þ e1 Þ dt r a k
! ¼ xððn þ l 1ÞT Þg.
ð2:8Þ
Then x((N1 + n + l)T) P x((N1 + l)T)gn ! 1 as n ! 1, which is a contradiction to the boundedness of x(t). Hence there exists a t1 > 0 such that x(t1) P m3. Second, if x(t) P m3 for all t P t1, then our aim is obtained. Hence we only need to consider those solutions which leave the region R ¼ fX ðtÞ 2 R2þ : xðtÞ < m3 g and reenter it again. Let t ¼ inf tPt1 fxðtÞ < m3 g, there are two possible cases for t*. Case I. t* = (n1 + l 1)T, n1 2 N. Then x(t) P m3 for t 2 [t1, t*] and (1 d1)m3 6 x(t*+) = (1 d1)x(t*) < m3. Select n2, n3 2 N such that e1 ln M þp ; ðn2 1ÞT > kcm3 l þ a ð1 d1 Þn2 expððn2 g1 T Þgn3 > ð1 d1 Þn2 expððn2 þ 1Þg1 T Þgn3 > 1; where g1 ¼ ðr rmk 3 ac MÞ < 0. Let T ¼ n2 T þ n3 T . We claim there must be a t2 2 ðt ; t þ T such that x(t2) > m3. Otherwise, consider (2.6) with z(t*+) = y(t*+), and we have 8 ! > p > k ðnðn1 þ1ÞÞ þ > > ð1 d2 Þ exp l þ cma 3 ðt n1 T Þ þ z ðtÞ; ðn 1ÞT < t 6 ðn þ l 1ÞT ; ðzðn1 T Þ > 3 < 1 ð1 d2 Þexp l þ kcm T a ! zðtÞ ¼ > > p > ð1 d Þn1 zðn T þ Þ 3 > exp l þ kcm ðn þ l 1ÞT < t 6 nT ; ðt n1 T Þ þ z ðtÞ; > 2 1 a : 3 T 1 ð1 d2 Þexp l þ kcm a
Z. Li et al. / Chaos, Solitons and Fractals 29 (2006) 1229–1239
1235
and n1 + 1 6 n 6 n1 + n2 + n3. Then ! kcm3 ðt n1 T Þ < e1 jzðtÞ ~zðtÞj < ðM þ pÞ exp l þ a and yðtÞ 6 zðtÞ 6 ~zðtÞ þ e1 ; ðn1 þ n2 1ÞT 6 t 6 t þ T , which implies that (2.7) holds for t þ n2 T 6 t 6 t þ T . As in the first step, we have xðt þ n2 T Þ P xðt þ n2 T Þgn3 . The first equation of system (1.2) gives 8 rm3 c < dxðtÞ P xðtÞðr MÞ; t 6¼ ðn þ l 1ÞT ; dt a k : þ xðt Þ ¼ ð1 d1 ÞxðtÞ; t ¼ ðn þ l 1ÞT .
ð2:9Þ
Integrating (2.9) on [t*, t* + n2T], we have xðt þ n2 T Þ P m3 ð1 d1 Þn2 expðn2 g1 T Þ. Thus xðt þ T Þ P m3 ð1 d1 Þn2 expðn2 g1 T Þgn3 > m3 , and this is a contradiction. Let t ¼ inf tPt fxðtÞ > m3 g, and then x(t) 6 m3, xðtÞ ¼ m3 . For t 2 ðt ; tÞ, we have xðtÞ P m3 ð1 d1 Þn2 þn3 expððn2 þ n3 Þg1 T Þ,m01 . For t > t, the same arguments can be continued since xðtÞ P m3 . Case II. t* 5 (n1 + l 1)T, n1 2 N. Then x(t) P m3 for t 2 [t1, t*] and x(t*) = m3, and suppose t 2 ððn01 þ l 1ÞT ; ðn01 þ lÞT Þ, n01 2 N. There are two possible cases for t 2 ðt ; ðn01 þ lÞT Þ. Case IIa. x(t) 6 m3 for all t 2 ðt ; ðn01 þ lÞT Þ. Similar to Case I, we can prove that there must be a t02 2 ½ðn01 þ lÞT ; 0 ðn1 þ lÞT þ T such that xðt02 Þ > m3 . Here we omit it. Let ~t ¼ inf tPt fxðtÞ > m3 g, then xð~tÞ 6 m3 . For t 2 ðt ; ~tÞ, we have xðtÞ P m3 ð1 d1 Þn2 þn3 expððn2 þ n3 þ 1Þg1 T Þ. Let m1 ,m3 ð1 d1 Þn2 þn3 expððn2 þ n3 þ 1Þg1 T Þ < m01 , so x(t) P m1 for t 2 ðt ; ~tÞ. For t > ~t, the same arguments can be continued since xð~tÞ P m3 . Case IIb. There exists a t 2 ðt ; ðn01 þ lÞT Þ such that x(t) > m3. Let ~t ¼ inf tPt fxðtÞ > m3 g, then xð~tÞ 6 m3 for t 2 ðt ; ~tÞ and xð~tÞ ¼ m3 . For t 2 ðt ; ~tÞ, (2.9) holds true. Integrating (2.9) on ðt ; ~tÞ, we have xðtÞ P xðt Þ expðg1 ðt t ÞÞ P m3 expðg1 T Þ > m1 . Since x(t) P m3 for t > ~t, the same arguments can be continued. Hence x(t) P m1 for all t P t1. The proof is completed. h
3. Computer simulation To study the dynamics of a Beddington-type system with impulsive control strategy, the solution of the system (1.2) with initial conditions in the first quadrant is obtained numerically for biologically feasible range of parametric value and the bifurcation diagram provides a summary of essential dynamical behavior of the system. As we know, the corresponding continuous system (1.2) cannot be solved explicitly. So we have to study the system (1.2) numerically integrated and research the long-term behavior of the solution in Maple. Maple, an international mathematical software, is an exchanged computer algebra system with great ability of symbolic evaluation, numerical calculation, coping with graphics, etc. Its powerful functions library and unique interior programming language provide science calculation and programming with friendly platform [18–20]. Firstly, we consider a set of parameters as the following: r = 3, k = 57, c = 0.76, k = 0.72, l = 0.4, a = 9, b = 0.005, d1 = 0.2, d2 = 0.0001, T = 8 and l = 0.09, and give some figures to help us with analysis. From Theorem 2.1 we know that the prey-eradication periodic solution (0, y*(t)) is globally asymptotically stable provided that p > pmax. A typical prey-eradication periodic solution of the system (1.2) is shown in Fig. 1, where we observe how the variable y(t) oscillates in a stable cycle. In contrast, the prey x(t) rapidly decreases to zero and pmax = 116.5544814. If the amount p of releasing species is smaller than pmax, the prey eradication solution becomes
1236
Z. Li et al. / Chaos, Solitons and Fractals 29 (2006) 1229–1239
unstable and undergoes a transcritical bifurcation, then the prey and predator can coexist on a stable positive periodic solution when p < pmax, and system (1.2) can be permanent which follows from Theorem 2.2. If the parameter is further increasing, the system will exhibit a wide variety of dynamical behavior.
Fig. 1. The dynamics of system (1.2) with impulsive control strategy with r = 3, k = 57, c = 0.76, k = 0.72, l = 0.4, a = 9, b = 0.005, d1 = 0.2, d2 = 0.0001, T = 8 and l = 0.09. when p > pmax = 116.5544814, (0, y*(t)) is globally asymptotically stable, that is to say, the prey will go extinction: (a) time series of the predator population y(t), (b) time series of the prey population, x(t) ! 0, as t ! +1 when p > pmax = 116.5544814, (c) phase portrait of system (1.2) of p = 117.0.
Fig. 2. Bifurcation diagrams of the system (1.2) shows the effects of p with initial values X(0) = (75, 7.2). (a) Prey population x and (b) predator population y, with p over [0, 120].
Fig. 3. The magnified parts of Fig. 2(a). (a) p over [0, 8], (b) p over [20, 30], (c) p over [68, 80], respectively.
Z. Li et al. / Chaos, Solitons and Fractals 29 (2006) 1229–1239
1237
Fig. 4. The magnified parts of Fig. 2(b). (a) p over [0, 8], (b) p over [20, 30], (c) p over [68, 80], respectively.
Fig. 5. Period doubling bifurcation leads to chaos of system (1.2) with initial values X(0) = (75, 7.2). (a) Phase portrait of T-period solution for p = 1.5, (b) phase portrait of 2T-period solution for p = 4.0, (c) phase portrait of 4T-period solution for p = 5.5, (d) phase portrait of 8T-period solution for p = 7.0 and (e) phase portrait of the solution enter chaotic region for p = 7.6.
Fig. 6. Symmetry-breaking pitchfork bifurcations of the system (1.2) with initial values X(0) = (75, 7.2). (a) Phase portrait of period solution for p = 8, (b) phase portrait of period solution for p = 10, (c) phase portrait of period solution for p = 18.
1238
Z. Li et al. / Chaos, Solitons and Fractals 29 (2006) 1229–1239
In Fig. 2, we have displayed typical bifurcation diagrams for the prey population x and the predator population y of the system (1.2) as p increasing from 0 to 120 with initial values X(0) = (75, 7.2), Figs. 3 and 4 are the magnified parts of Fig. 2(a) and (b), and the windows of periodic behavior are more visible. As p increases, the resulting bifurcation diagrams clearly show that system (1.2) has rich dynamics including period doubling bifurcation, crises (the phenomenon of ‘‘crisis’’ in which chaotic attractors can suddenly appear or disappear, or change size discontinuously as a parameter smoothly varies), symmetry-breaking pitchfork bifurcations, chaotic bands, quasi-periodic oscillation, narrow periodic window, wide periodic window, period-halving bifurcation, etc. Fig. 2 shows that adding p will cause the availability of abundance of chaotic regions. As the bifurcation diagrams shown, when 0 < p < 2.197, the system (1.2) has exactly a T-periodic solution, when p > 2.197, there is a period doubling leading chaos (Fig. 5). When p beyond 7.602, symmetry-breaking pitchfork bifurcations of the system (1.2) with initial values X(0) = (75, 7.2), appear (Fig. 6), which have three pitchforks. Periodic windows are intermittently scattered. The abundance of chaos increases as p increase. In course of this variety, the phenomenon of ÔcrisisÕ occurs (Fig. 7), when the p is slightly increased beyond p 25.023, the chaos suddenly disappears, and from p 41.103 to p 51.592 there are a series of crises, which constitute several types of crises, we can see that from the magnified parts of bifurcation diagrams. After these chaotic areas, that is to say, when p > 79.363, the solution of the system (1.2)
Fig. 7. Crises are shown: when the p is slightly increased beyond p 25.023, the chaos suddenly disappears, and from p 41.103 to p 51.592 there are a series of crises exist, (a) p = 23, (b) p = 28, (c) p = 50, (d) p = 60, from (a) to (b) there is a crisis that the chaos suddenly disappears, and from (c) to (d) there is a crisis that the chaos suddenly appears.
Fig. 8. Period-halving bifurcation from chaos to cycles, (a) chaos when p = 70, (b) 2T-period solution when p = 90, (c) T-period solution when p = 108.
Z. Li et al. / Chaos, Solitons and Fractals 29 (2006) 1229–1239
1239
undergoes a cascade of period-halving bifurcation from chaos to cycles (Fig. 8). Periodic halving is the flip bifurcation in the opposite direction.
4. Conclusion and remarks In this paper, we have focused on the dynamics of a Beddington-type system with impulsive control strategy. Conditions for the system to be extinct are given by using the Floquet theory of impulsive equation and small amplitude perturbation skills. It is proved that the system is permanent via the method of comparison involving multiple Liapunov functions. Using the method of numerical simulation, the influence of the impulsive control strategy on the inherent oscillation are investigated, which shows the rich dynamics, such as period doubling bifurcation, crises, symmetrybreaking pitchfork bifurcations, chaotic bands, quasi-periodic oscillation, narrow periodic window, wide periodic window, period-halving bifurcation, and that such bifurcations offer considerable potential for fascinating mathematical analysis. On the other hand, since the complexity of the ecosystem, we cannot get their analytical solution. Then, we have to solve the approximate solution of them. And now, the method of computer simulation is one of the best methods. There are some international mathematical software, such as Matlab, Maple, Mathematica, etc. They all have powerful function library and can provide science calculation and programming with friendly platform. In Maple, there is a package DEtools which helps us work with differential equations. The results indicate that computer simulation is a useful method for study of the complex dynamic systems.
References [1] Lakshmikantham V, Bainov DD, Simeonov PC. Theory of impulsive differential equations. Singapore: World Scientific; 1989. [2] Tang S, Xiao Y, Clancy D. New modelling approach concerning integrated disease control and cost-effectivity. Nonlinear Anal 2005;63:439–71. [3] Lu Z, Chi X, Chen L. Impulsive control strategies in biological control of pesticide. J Theor Popul Biol 2003;64:39–47. [4] Liu B, Zhang Y, Chen L. Dynamic complexities of a Holling I predator–prey model concerning periodic biological and chemical control. Chaos, Solitons & Fractals 2004;22:123–34. [5] Sun J, Qiao F, Wu Q. Impulsive control of a financial model. J Phys Lett A 2005;335:282–8. [6] Zhang Y, Xiu Z, Chen L. Dynamics complexity of a two-prey one-predator system with impulsive effect. Chaos, Solitons & Fractals 2005;26:131–9. [7] Ballinger G, Liu X. Permanence of population growth models with impulsive effects. Math Comput Model 1997;26:59–72. [8] Lenci S, Rega G. Periodic solutions and bifurcations in an impact inverted pendulum under impulsive excitation. Chaos, Solitons & Fractals 2000;11:2453–72. [9] Beddington JR. Mutual interference between parasites or predators and its effect on searching efficiency. J Animal Ecol 1975;44: 331–40. [10] Cantrell RS, Cosner C. On the dynamics of predator–prey models with the Beddington–DeAngelis functional response. J Math Anal Appl 2001;257:206–22. [11] Hwang TW. Global analysis of the predator–prey system with Beddington–DeAngelis functional response. J Math Anal Appl 2003;281:395–401. [12] Hwang TW. Uniqueness of limit cycles of the predator–prey system with Beddington–DeAngelis functional response. J Math Anal Appl 2004;290:113–22. [13] Liu Z, Yuan R. Stability and bifurcation in a delayed predator–prey system with Beddington–DeAngelis functional response. J Math Anal Appl 2004;296:521–37. [14] Qui Z, Yu J, Zou Y. The asymptotic behavior of a chemostat model with the Beddington–DeAngelis functional response. Math Bios 2004;187:175–87. [15] Dimitrov DT, Kojouharov VH. Complete mathematical analysis of predator–prey models with linear prey growth and Beddington–DeAngelis functional response. J Appl Math Comp 2005;162:523–38. [16] Huisman G, Deboer RJ. A formal derivation of the Beddington functional response. J Thero Biol 1997;185:389–400. [17] Bainov DD, Simeonov PS. Impulsive differential equations: asymptotic properties of the solutions. Singapore: World Scientific; 1993. [18] Wang W. Computer algebra system and symbolic computation. Lanzhou, China: Gansu Sci-Tech Press; 2004 [in Chinese]. [19] Wang W, Lian X. A new algorithm for symbolic integral with application. Appl Math Comput 2005;162:949–68. [20] Wang W, Lin C. A new algorithm for integral of trigonometric functions with mechanization. Appl Math Comput 2005;164: 71–82.