- Email: [email protected]
The dynamic complexity of a three-species Beddington-type food chain with impulsive control strategy
Chaos, Solitons and Fractals 32 (2007) 1772–1785 www.elsevier.com/locate/chaos
The dynamic complexity of a three-species Beddington-type food chain with impulsive control strategy q Weiming Wang a
a,b
, Hailing Wang c, Zhenqing Li
a,*
Laboratory of Quantitative Vegetation Ecology, Institute of Botany, The Chinese Academy of Sciences, Beijing 100093, China b School of Mathematics and Information Science, Wenzhou University, Wenzhou Zhejiang 325035, China c School of Applied Mathematics, University of Electronic Science and Technology of China, Chengdu Sichuan 610054, China Accepted 1 December 2005
Communicated by Prof. L. Marck-Crnjac
Abstract In this paper, by using theories and methods of ecology and ordinary differential equation, the dynamics complexity of 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 Floquet theory of impulsive equation and small amplitude perturbation skills. Furthermore, by using the method of numerical simulation with the international software Maple, the influence of the impulsive perturbations on the inherent oscillation is investigated, which shows rich dynamics, such as quasi-periodic oscillation, narrow periodic window, wide periodic window, chaotic bands, period doubling bifurcation, symmetry-breaking pitchfork bifurcation, period-halving bifurcation and crises, etc. The numerical results indicate that computer simulation is a useful method for studying the complex dynamic systems. Ó 2006 Elsevier Ltd. All rights reserved.
1. Introduction The recent 40 years have seen an explosion of interest in the study of nonlinear dynamical systems. More recently, dynamical systems have benefited from an infusion of interests and techniques from a variety of fields [1]. It is currently very much in vogue to study chaotic behavior, and much of it is prompted by new and potential application in different fields [2–24]. In the recent decade, simple multi-species systems comprising of a three-species food chain were discussed by a number of researchers, and there is also a large body of literature on models with perturbations, particularly in prey–predator systems [3,6–24].
q
This research was supported by the National Natural Science Foundation of China (Grant Nos. 30270249 and 10371090) and the Knowledge Innovation Project of the Chinese Academy of Sciences (Grant No. KSCX2-SW-107). * Corresponding author. E-mail addresses: [email protected] (W. Wang), [email protected] (H. Wang), [email protected] (Z. Li). 0960-0779/$ - see front matter Ó 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2005.12.025
W. Wang et al. / Chaos, Solitons and Fractals 32 (2007) 1772–1785
1773
A well-known model of such systems is the predator–prey model with a Beddington-type functional response [25–32], which can be described by the following differential equations: 8 dxðtÞ c1 xðtÞyðtÞ > > > dt ¼ xðtÞða bxðtÞÞ a þ xðtÞ þ b yðtÞ ; > > 1 1 > > > < dyðtÞ k 1 c1 xðtÞyðtÞ c2 yðtÞzðtÞ d 1 yðtÞ; ¼ ð1:1Þ > dt a1 þ xðtÞ þ b1 yðtÞ a2 þ yðtÞ þ b2 zðtÞ > > > > > > dzðtÞ k 2 c2 yðtÞzðtÞ > : d 2 zðtÞ. ¼ dt a2 þ yðtÞ þ b2 zðtÞ In our paper, we consider an impulsive differential equation—a predator–prey model with a Beddington-type functional response, by introducing a proportion periodic impulsive catching or poisoning for the prey populations, and a constant periodic releasing for the predator at different fixed moment, respectively, what follows is the corresponding equations: 8 9 dxðtÞ c1 xðtÞyðtÞ > > > > ¼ xðtÞða bxðtÞÞ > > > dt > a1 þ xðtÞ þ b1 yðtÞ > > > > > > > > > dyðtÞ = > k c xðtÞyðtÞ c yðtÞzðtÞ > 1 1 2 > d yðtÞ ¼ t 6¼ ðn þ l 1ÞT ; t 6¼ nT ; > 1 > > > dt a1 þ xðtÞ þ b1 yðtÞ a2 þ yðtÞ þ b2 zðtÞ > > > > > > > > > > > dzðtÞ k 2 c2 yðtÞzðtÞ > > > > ; ¼ d zðtÞ 2 > > dt a þ yðtÞ þ b zðtÞ > 2 2 < 9 ð1:2Þ DxðtÞ ¼ d1 xðtÞ > > = > > > > t ¼ ðn þ l 1ÞT ; DyðtÞ ¼ d2 yðtÞ > > > > ; > > > DzðtÞ ¼ d3 zðtÞ > > > > 9 > > > DxðtÞ ¼ 0 > > = > > > > DyðtÞ ¼ 0 t ¼ nT ; > > > > ; : DzðtÞ ¼ p where x(t), y(t) and z(t) are functions of time representing population density of the prey, mid-level and top predator, respectively, all parameters are positive constants, the parameter ki (i = 1, 2) are conversion efficiencies, di (i = 1, 2) are the mortality rates of the mid-level predator and the top predator, ci (i = 1, 2) are the maximum numbers of prey that can be eaten by a predator per unit of time, ai (i = 1, 2) are saturation constants and bi (i = 1, 2) scale the impact of the predator interference. T is the period of the impulsive effect and n 2 N. Observe the simple relation of these three species: z only prey on y and y prey on x and nutrient recycling is not accounted for. This simple relation produces the socalled simple food chain system with impulsive perturbations. With model (1.2), we can take into account the possible exterior effects under which the population density changes very rapidly. An impulsive increase of the predator population density is possible by artificial breeding of the species or release some species (p > 0). The rest of this paper is organized as follows: In Section 2, some notations and theorems are given. In Section 3, the results of computer simulation are shown, and moreover, these results are discussed briefly. Finally, conclusions and remarks are given in Section 4.
2. Mathematical analysis In this section, we will give some definitions, notations and some lemmas which will be useful for our main results. Let R+ = [0, 1), R3þ ¼ fx 2 R3 : x P 0g; X ¼ int R3þ , N be the set of all non-negative integers. The map f = (f1, f2, f3)T defined by the system (1.1). Let V : Rþ R3þ ! Rþ , then V is said to belong to class V0 if (1) V is continuous in ððn 1ÞT ; ðn þ l lÞT R3þ and ððn þ l 1ÞT ; nT R3þ , and for all x 2 R3þ , n 2 N, limðt;yÞ!ððnþl1ÞT þ ;xÞ V ðt; yÞ ¼ V ððn þ l lÞT þ ; xÞ and limt;y!ðnT þ ;xÞ V ðt; yÞ ¼ V ðnT þ ; xÞ exist; (2) V is locally Lipschitzian in x.
1774
W. Wang et al. / Chaos, Solitons and Fractals 32 (2007) 1772–1785
Definition 2.1. V 2 V0, then for ðt; xÞ 2 ððn 1ÞT ; ðn þ l 1ÞT R3þ and ððn þ l 1ÞT ; nT R3þ , 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 The solution of system (1.2) is a piecewise continuous function X ðtÞ : Rþ ! R3þ ; X ðtÞ is continuous on T, T = ((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). It is easy to prove the following lemmas. Lemma 2.1. Let X(t) be a solution to system (1.2) with X(0+) P 0, then X(t) P 0 for all t P 0. And furthermore, 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, z(t) 6 M for each solution (x(t), y(t), z(t)) of (1.2) with sufficiently large t. Then, we give some basic properties about the following subsystem of system (1.2). 8 dzðtÞ > > ¼ d 2 zðtÞ; t 6¼ ðn þ l 1ÞT ; t 6¼ nT ; > > dt > > < þ zðt Þ ¼ ð1 d3 ÞzðtÞ; t ¼ ðn þ l 1ÞT ; > > zðtþ Þ ¼ zðtÞ þ p; t ¼ nT ; > > > > : þ zð0 Þ ¼ z0 .
ð2:1Þ
Clearly 8 p expðd 2 ðt ðn 1ÞT ÞÞ > > ; ðn 1ÞT < t 6 ðn þ l 1ÞT ; > < 1 ð1 d3 Þ expðd 2 T Þ z ðtÞ ¼ > pð1 d3 Þ expðd 2 ðt ðn 1ÞT ÞÞ > > ; ðn þ l 1ÞT < t 6 nT ; : 1 ð1 d3 Þ expðd 2 T Þ p pð1 d3 Þ expðd 2 lT Þ þ þ þ ; z ðlT Þ ¼ n 2 N ; z ð0 Þ ¼ z ðnT Þ ¼ 1 ð1 d3 Þ expðd 2 T Þ 1 ð1 d3 Þ expðd 2 T Þ is a positive periodic solution of system (1.2). Since 8 p > ðn1Þ þ > > expðd 2 tÞ þ z ðtÞ; ðn 1ÞT < t 6 ðn þ l 1ÞT ; Þ zð0 Þ ð1 d 3 < 1 ð1 d3 Þ expðd 2 T Þ zðtÞ ¼ > p n > þ > expðd 2 tÞ þ z ðtÞ; ðn þ l 1ÞT < t 6 nT Þ zð0 Þ ð1 d : 3 1 ð1 d3 Þ expðd 2 T Þ is the solution of system (1.2) with initial value z0 P 0, where t 2 T. Then we can get: Lemma 2.3. Let z*(t) be a positive periodic solution of system (1.2) and all solution z(t) of system (1.2) with z0 P 0, we have jz(t) z*(t)j ! 0, when t ! 1. Therefore, we obtain the complete expression for prey and mid-level predator eradication periodic solution (0, 0, z*(t)) and ab ; 0; z ðtÞ to system (1.2) for t 2 T. Now, we study the stability of prey and mid-level predator eradication periodic solution to system (1.2). Theorem 2.1. Let (x(t), y(t), z(t)) be any solution to (1.2), then (1) (0, 0, z*(t)) is unstable; (2) ðab ; 0; z ðtÞÞ is locally asymptotically stable provided that
c2 ða1 b þ aÞðA1 A2 þ A3 A4 d 1 d 2 b2 T Þ þ k 1 c1 ab2 d 2 T 1 < ln ; b2 d 2 ða1 b þ aÞ 1 d2
W. Wang et al. / Chaos, Solitons and Fractals 32 (2007) 1772–1785
where
1775
a2 expðd 2 T Þ þ b2 p expðð1 lÞd 2 T Þ a2 þ a2 d3 ; expðd 2 T Þ þ d3 1 ða2 þ b2 Þ expðd 2 T Þ a2 þ a2 d3 ; A2 ¼ ln expðd 2 T Þ þ d3 1 a2 expðd 2 T Þ a2 þ a2 d3 þ b2 p þ b2 pd3 ; A3 ¼ ln expðd 2 T Þ þ d3 1 a2 expðd 2 T Þ þ ð1 d3 Þb2 p expðð1 lÞd 2 T Þ a2 þ a2 d3 . A4 ¼ ln expðd 2 T Þ þ d3 1 A1 ¼ ln
Proof. (1) The local stability of periodic solution (0, 0, z*(t)) is determined by considering the behavior of small amplitude perturbations of the solution. Define xðtÞ ¼ uðtÞ; Then
zðtÞ ¼ wðtÞ þ z ðtÞ.
yðtÞ ¼ vðtÞ;
1 0 1 uð0Þ uðtÞ C B C B @ vðtÞ A ¼ UðtÞ@ vð0Þ A; wð0Þ wðtÞ 0
0 6 t < T;
where U(t) satisfies 0 a 0 c 2 z B 0 d1 dU B a2 þ b2 z ¼B B dt @ k 2 c 2 z 0 a2 þ b2 z
0
1
C 0 C CUðtÞ C A d 2
and U(0) = I, the identity matrix. The linearization of system (1.2) becomes 1 0 1 0 10 uððn þ l 1ÞT Þ 0 0 1 d1 uððn þ l 1ÞT þ Þ B C B C B C 1 d2 0 A@ vððn þ l 1ÞT Þ A @ vððn þ l 1ÞT þ Þ A ¼ @ 0 wððn þ l 1ÞT Þ 0 0 1 d3 wððn þ l 1ÞT þ Þ and 1 10 1 0 uðnT Þ 1 0 0 uðnT þ Þ C CB C B B @ vðnT þ Þ A ¼ @ 0 1 0 A@ vðnT Þ A. þ wðnT Þ 0 0 1 wðnT Þ 0
Hence, if all eigenvalues of 0 1 d1 0 B H¼@ 0 1 d2 0
0
10
1 1 0 0 CB C A@ 0 1 0 AUðT Þ
0 0 1 d3
0 0 1
have absolute less than one, the periodic solution (0, 0, z*(t)) is locally stable. But one of the Floquet multipliers is u1 = exp((1 d1)aT) > 0, according to the Floquet theory, the solution X*(t) is unstable. (2) The local stability of periodic solution ðab ; 0; z ðtÞÞ can be similarly discussed. Define a xðtÞ ¼ uðtÞ þ ; yðtÞ ¼ vðtÞ; zðtÞ ¼ wðtÞ þ z ðtÞ. b then
0
1 0 1 uðtÞ uð0Þ B C B C @ vðtÞ A ¼ UðtÞ@ vð0Þ A; wðtÞ wð0Þ
0 6 t < T;
1776
W. Wang et al. / Chaos, Solitons and Fractals 32 (2007) 1772–1785
where U(t) satisfies 0 c1 ab a B a1 þ ab B B a k 1 c1 b dU B c 2 z ¼B 0 d1 a B dt a1 þ b a2 þ b2 z B @ k 2 c 2 z 0 a2 þ b2 z
0
1
C C C C 0 CUðtÞ C C A d 2
and U(0) = I, the identity matrix. The linearization of system (1.2) becomes 1 10 1 0 0 uððn þ l 1ÞT Þ 1 d1 0 0 uððn þ l 1ÞT þ Þ C CB C B B 1 d2 0 A@ vððn þ l 1ÞT Þ A @ vððn þ l 1ÞT þ Þ A ¼ @ 0 wððn þ l 1ÞT Þ 0 0 1 d3 wððn þ l 1ÞT þ Þ and 0
1 0 10 1 1 0 0 uðnT Þ uðnT þ Þ B C B CB C @ vðnT þ Þ A ¼ @ 0 1 0 A@ vðnT Þ A. 0 0 1 wðnT Þ wðnT þ Þ The stability of the periodic solution ðab ; 0; z ðtÞÞ is determined by the eigenvalues of 1 0 10 1 0 0 0 0 1 d1 C B CB 1 d2 0 A@ 0 1 0 AUðT Þ H¼@ 0 0 0 1 0 0 1 d3 which are u1 ¼ ð1 d1 Þ expðaT Þ < 1; Z T k 1 c1 ab c 2 z d u2 ¼ ð1 d2 Þ exp 1 dt ; a1 þ ab a2 þ b2 z 0 u3 ¼ ð1 d3 Þ expðd 2 T Þ < 1. According to Floquet theory, the solution ðab ; 0; z ðtÞÞ is locally stable if ju2j < 1, that is to say, c2 ða1 b þ aÞðA1 A2 þ A3 A4 d 1 d 2 b2 T Þ þ k 1 c1 ab2 d 2 T 1 < ln . b2 d 2 ða1 b þ aÞ 1 d2 This completes the proof.
h
3. Computer simulation It is well known that the correspondingly continuous system (1.2) cannot be solved explicitly. So we have to study the system (1.2) numerically synthetically and research the long-term behavior of the solution with Maple. Maple, an international mathematical software, is an interactive computer algebra system with great ability of symbolic evaluation, numerical calculation, coping with graphics, etc. Its powerful function library and unique interior programming language provide science calculation and programming with friendly platform [33–38]. We consider three sets of parameters for our analysis as follows. Case A: a = 5, b = 0.5, c1 = 2, c2 = 2, a1 = 0.1, a2 = 0.5, b1 = 1, b2 = 0.1, d1 = 0.1, d2 = 0.9, k1 = 0.7, k2 = 0.5, T = 8, d1 = 0.2, d2 = 0.0001, d3 = 0.03, l = 0.6. Case B: a = 2, b = 0.1, c1 = 2, c2 = 3.5, a1 = 0.1, a2 = 0.7, b1 = 7, b2 = 0.1, d1 = 0.1, d2 = 0.9, k1 = 0.7, k2 = 0.5, T = 16, d1 = 0.2, d2 = 0.0001, d3 = 0.07, l = 0.6. Case C: a = 3, b = 0.25, c1 = 2, c2 = 3.5, a1 = 0.1, a2 = 0.5, b1 = 1, b2 = 0.1, d1 = 0.1, d2 = 0.9, k1 = 0.7, k2 = 0.5, T = 5, d1 = 0.2, d2 = 0.0001, d3 = 0.08, l = 0.68.
W. Wang et al. / Chaos, Solitons and Fractals 32 (2007) 1772–1785
1777
For Case A, in Fig. 1, we have displayed typical bifurcation diagrams for the prey population x(t), the mid-level predator population y(t) and the top predator population z(t) of the system (1.2) as p increasing from 0 to 3.5 with initial values X(0) = (15, 1, 0.4). As the similarity of the bifurcation diagrams, we only give the magnified parts of Fig. 1(a) in Fig. 2, and the windows of periodic behavior are more clearly visible. As p increases, the resulting bifurcation diagrams clearly show that system (1.2) has rich dynamics including: quasi-periodic oscillation, narrow periodic window, wide periodic window, chaotic bands, crises (the phenomenon of ‘‘crisis’’ in which chaotic attractors can suddenly appear or disappear, or change size discontinuously as a parameter smoothly varies), period doubling bifurcation, period-halving bifurcation, symmetry-breaking pitchfork bifurcation, tangent bifurcation, etc. In Fig. 1(a) and Fig. 2, the bifurcation diagrams are drawn as a function of x. As the bifurcation diagrams shown, the dynamics behavior of prey population x may dramatically be affected by small changes in the value of the parameter p. For 0 < p < 1.342, Fig. 2(a) shows that a quasi-periodic solution oscillation, in which there are frequency-lockings, symmetry-breaking pitchfork bifurcation and tangent bifurcation (Fig. 3). When the parameter p is slightly increased beyond p > 1.342, there is a 4T-period solution of system (1.2) appears (Fig. 4), and after this, the parameter p is slightly
Fig. 1. Bifurcation diagrams of system (1.2) with a = 5, b = 0.5, c1 = 2, c2 = 2, a1 = 0.1, a2 = 0.5, b1 = 1, b2 = 0.1, d1 = 0.1, d2 = 0.9, k1 = 0.7, k2 = 0.5, T = 8, d1 = 0.2, d2 = 0.0001, d3 = 0.03, l = 0.6, and with initial values X(0) = (15, 1, 0.4). (a) Prey population x, (b) mid-level predator population y and (c) the top population z, with p over [0, 3.5].
Fig. 2. The magnified parts of Fig. 1(a): (a) p over [0.5, 1.1] and (b) p over [1.2, 2.4], respectively.
Fig. 3. Quasi-periodic solutions oscillation of system (1.2) with initial values X(0) = (15, 1, 0.4), there are frequency-lockings, symmetry-breaking pitchfork bifurcation and tangent bifurcation in it. (a) p = 0.25, (b) p = 0.5, (c) p = 0.8 and (d) p = 1.3.
1778
W. Wang et al. / Chaos, Solitons and Fractals 32 (2007) 1772–1785
Fig. 4. A 4T-periodic solution when the parameter p is slightly increased beyond p > 1.342.
increased beyond p > 2.09, the 4T-period solution becomes unstable and there is a cascade of period doubling bifurcation leads to chaos, which can been seen in Fig. 5, when p increased beyond 2.254, there is a cascade of period-halving bifurcation leads to 4T-period solution (Fig. 6). There is a crisis appears again when p 2.313, leads to a chaos band, and then, the chaos disappears, a T-period solution of system (1.2) appears (Fig. 7). If p > 2.665954781, as Theorem 2.1 shows, the mid-level predator will go extinct y(t) ! 0 as t ! 1 (Fig. 8).
Fig. 5. Period doubling bifurcation leads to chaos of system (1.2) with initial values X(0) = (15, 1, 0.4), (a) phase portrait of 8T-period solution for p = 2.1 and (b) phase portrait of the four attractors for p = 2.18.
Fig. 6. Period-halving bifurcation leads to chaos of system (1.2) with initial values X(0) = (15, 1, 0.4), (a) phase portrait of 4T-period solution for p = 2.223 and (b) phase portrait of the solutions enter chaotic region for p = 2.23.
Fig. 7. Crises are shown: when the parameter p is slightly increased beyond p = 2.313, the chaos suddenly disappears, (a) p = 2.31 and (b) p = 2.5, from (a) to (b) there is a crisis that the chaos suddenly disappears.
W. Wang et al. / Chaos, Solitons and Fractals 32 (2007) 1772–1785
1779
Fig. 8. When p > pmax ¼ 2:665954781; ab ; 0; z ðtÞ is locally asymptotically stable, that is to say, the mid-level predator will go extinction: (a) time series of the prey population x(t), (b) time series of the mid-level predator population y(t), when p > pmax = 2.665954781, y(t) ! 0, as t ! +1, (c) time series of the top predator population z(t) and (d) phase portrait of system (1.2) of p = 2.666.
For Case B, in Fig. 9, we have displayed typical bifurcation diagrams for the prey population x(t), the mid-level predator population y(t) and the top predator population z(t) of the system (1.2) as p increasing from 0 to 4 with initial values X(0) = (15, 1, 0.4). As the similarity of the bifurcation diagrams, we only give the magnified parts of Fig. 9(a) in Fig. 10. As p increases, the resulting bifurcation diagrams clearly show that system (1.2) has rich dynamics including: quasi-periodic oscillation, period doubling bifurcation, symmetry-breaking pitchfork bifurcation, tangent bifurcation, crisis, chaotic bands, narrow periodic window, wide periodic window, period-halving bifurcation, etc. In Fig. 9(a) and Fig. 10, the bifurcation diagrams are drawn as a function of x. As the bifurcation diagrams shown, the dynamics behavior of prey population x may dramatically be affected by small changes in the value of the parameter p. For 0 < p < 0.069, Fig. 9(a) shows that a quasi-periodic solution oscillation, when 0.069 < p < 0.153, T-period solution of system (1.2) is still stable (Fig. 11). When p > 0.153, it becomes unstable and there is a cascade of period doubling bifurcation leads to chaos, it enters a chaotic region (Fig. 12). Unlike the common period doubling route to chaos, here we observe the frequent occurrence of sudden changes of attractors (crises) in the range of [0.254, 0.321], where
Fig. 9. Bifurcation diagrams of system (1.2) with a = 2, b = 0.1, c1 = 2, c2 = 3.5, a1 = 0.1, a2 = 0.7, b1 = 7, b2 = 0.1, d1 = 0.1, d2 = 0.9, k1 = 0.7, k2 = 0.5, T = 16, d1 = 0.2, d2 = 0.0001, d3 = 0.07, l = 0.6, and with initial values X(0) = (15, 1, 0.4). (a) Prey population x, (b) mid-level predator population y and (c) the top population z, with p over [0, 4].
Fig. 10. The magnified parts of Fig. 9(a): (a) p over [0, 0.5], (b) p over [0.5, 1.1] and (c) p over [1.2, 2.4], respectively.
1780
W. Wang et al. / Chaos, Solitons and Fractals 32 (2007) 1772–1785
Fig. 11. Quasi-periodic solutions oscillation of system (1.2) with initial values X(0) = (15, 1, 0.4), (a) p = 0.068 and (b) p = 0.13.
Fig. 12. Period doubling bifurcation leads to chaos of system (1.2) with initial values X(0) = (15, 1, 0.4), (a) phase portrait of 2T-period solution for p = 0.2, (b) phase portrait of 4T-period solution for p = 0.35, (c) phase portrait of 8T-period solution for p = 0.38 and (d) phase portrait of the solutions enter chaotic region for p = 0.42.
Fig. 13. Crises are shown: when the parameter p is slightly increased beyond p = 0.677, the chaos suddenly disappears and multiply attractors appears, (a) p = 0.68, (b) p = 0.75 and (c) p = 0.85.
Fig. 14. Crises are shown: when the parameter p is slightly increased beyond p = 0.923, the chaos suddenly disappears and multiply attractors appears, (a) p = 0.92, (b) p = 0.94 and (c) p = 1.0.
multiple attractors coexist. And the same instance occur on the range of [0.677, 0.858] and [0.923, 0.990] (Figs. 13 and 14). When the parameter p is slightly increased beyond p > 1.239, the chaos disappears, period-halving bifurcation leads
W. Wang et al. / Chaos, Solitons and Fractals 32 (2007) 1772–1785
1781
to T-period solution (Fig. 15). If p > 3.751536111, as Theorem 2.1 shows, the mid-level predator will go extinct y(t) ! 0 as t ! 1 (Fig. 16). For Case C, in Fig. 17, we have displayed typical bifurcation diagrams for the prey population x(t), the mid-level predator population y(t) and the top predator population z(t) of the system (1.2) as p increasing from 0 to 1 with initial values X(0) = (15, 1, 0.4). As the similarity of the bifurcation diagrams, we only give the magnified parts of Fig. 17(a) in Fig. 18. As p increases, the resulting bifurcation diagrams clearly show that system (1.2) has rich dynamics including: quasi-periodic oscillation, period doubling bifurcation, symmetry-breaking pitchfork bifurcation, tangent bifurcation, crisis, chaotic bands, narrow periodic window, wide periodic window, period-halving bifurcation, etc. In Figs. 17(a) and 18, the bifurcation diagrams are drawn as a function of x. As the bifurcation diagrams shown, the dynamics behavior of prey population x may dramatically be affected by small changes in the value of the parameter p.
Fig. 15. Period-halving bifurcation leads to T-period solution, (a) phase portrait of the solutions in chaotic region for p = 1.24, (b) phase portrait of 4T-period solution for p = 1.4, (c) phase portrait of 2T-period solution for p = 1.6 and (d) phase portrait of T-period solution for p = 3.7.
Fig. 16. When p > pmax ¼ 3:751536111; ab ; 0; z ðtÞ is locally asymptotically stable, that is to say, the mid-level predator will go extinction: (a) time series of the prey population x(t), (b) time series of the mid-level predator population y(t), when p > pmax = 3.751536111, y(t) ! 0, as t ! +1, (c) time series of the top predator population z(t) and (d) phase portrait of system (1.2) of p = 3.8.
Fig. 17. Bifurcation diagrams of system (1.2) with a = 3, b = 0.25, c1 = 2, c2 = 3.5, a1 = 0.1, a2 = 0.5, b1 = 1, b2 = 0.1, d1 = 0.1, d2 = 0.9, k1 = 0.7, k2 = 0.5, T = 5, d1 = 0.2, d2 = 0.0001, d3 = 0.08, l = 0.68 and with initial values X(0) = (15,1,0.4). (a) prey population x, (b) mid-level predator population y and (c) the top population z, with p over [0, 1].
1782
W. Wang et al. / Chaos, Solitons and Fractals 32 (2007) 1772–1785
Fig. 18. The magnified parts of Fig. 17(a): (a) p over [0, 0.04], (b) p over [0.04, 0.2], (c) p over [0.2, 0.28] and (d) p over [0.3, 0.6], respectively.
Fig. 19. Quasi-periodic solutions oscillation of system (1.2) with initial values X(0) = (15, 1, 0.4), (a) p = 0.0015 and (b) p = 0.04.
Fig. 20. Chaos regions with frequency-lockings, symmetry-breaking pitchfork bifurcation and tangent bifurcation in it. (a) p = 0.015, (b) p = 0.02, (c) p = 0.035 and (d) p = 0.08.
Fig. 21. Period doubling bifurcation leads to chaos of system (1.2) with initial values X(0) = (15, 1, 0.4), (a) phase portrait of 2T-period solution for p = 0.11, (b) phase portrait of 4T-period solution for p = 0.12, (c) phase portrait of 8T-period solution for p = 0.148 and (d) phase portrait of the solutions enter chaotic region for p = 0.165.
W. Wang et al. / Chaos, Solitons and Fractals 32 (2007) 1772–1785
1783
For 0 < p < 0.0058, Fig. 17(a) shows that a quasi-periodic solution oscillation (Fig. 19), when p > 0.0123, a crisis occurs, that is to say, leads to a chaos region, in which there are frequency-lockings, symmetry-breaking pitchfork bifurcation, tangent bifurcation, narrow periodic window and wide periodic window Fig. 20. When the parameter p is slightly increased beyond p > 0.098, the chaos disappears, period doubling bifurcation leads to another chaos region Fig. 21,
Fig. 22. Period doubling bifurcation leads to chaos of system (1.2) with initial values X(0) = (15, 1, 0.4), (a) phase portrait of 2T-period solution for p = 0.176, (b) phase portrait of 4T-period solution for p = 0.18 and (c) phase portrait of the solutions enter chaotic region for p = 0.19.
Fig. 23. Chaos regions with frequency-lockings, symmetry-breaking pitchfork bifurcation and tangent bifurcation in it, that is to say, there are several crises in the region [0.22, 0.7]. (a) p = 0.222, (b) p = 0.225, (c) p = 0.25, (d) p = 0.4, (e) p = 0.55 and (f) p = 0.7.
Fig. 24. When p > pmax ¼ 0:7183181912; ab ; 0; z ðtÞ is locally asymptotically stable, that is to say, the mid-level predator will go extinction: (a) time series of the prey population x(t), (b) time series of the mid-level predator population y(t), when p > pmax = 0.7183181912, y(t) ! 0, as t ! +1, (c) time series of the top predator population z(t) and (d) phase portrait of system (1.2) of p = 0.72.
1784
W. Wang et al. / Chaos, Solitons and Fractals 32 (2007) 1772–1785
then the chaos disappears again when p > 0.173, another period doubling bifurcation leads to the third chaos region Fig. 22, in the following, when the p > 0.228, it becomes unstable and there is a cascade of period-halving bifurcation leads to the fourth chaos region, and in this course, multiply attractors coexist in the range of [0.239, 0.244], which we can observe in Fig. 18. Then the last crisis appears when p is slightly increased beyond p > 0.450, after this, period-halving bifurcation leads to T-period solution (Fig. 23). If p > 0.7183181912, as Theorem 2.1 shows, the mid-level predator will go extinct y(t) ! 0 as t ! 1 (Fig. 24).
4. Conclusions and remarks We have considered the dynamic behavior of a three-species food chain with Beddington-type functional response and periodic constant impulsive perturbations of top predator. Conditions for extinction of prey and mid-level predator is established and the boundary of the system (1.2) is given. We have observed in the computer simulation that the dynamics of a population may dramatically be affected by small changes in the value of the parameter p, at the same time, we can see from the bifurcation diagrams, with different parameter sets, we have different bifurcation diagrams, that is to say, we have different routes to chaos. On the other hand, for the complexity of the ecosystem, we cannot get its analytical solution. Then, we have to solve the numerical solutions. It is well known that the method of computer simulation is one of the best methods. There are some international mathematical softwares, 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 studying the complex dynamical systems.
References [1] Devaney RL. An introduction to chaotic dynamical systems. Menlo Park, CA: Benjamin/Cummings publishing company; 1986. [2] May RM. Necessity and chance: deterministic chaos in ecology and evolution. Am Math Soc 1995;32:291–308. [3] 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. [4] 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. [5] Gakkhar S, Naji MA. Order and chaos in predator to prey ratio-dependent food chain. Chaos, Solitons & Fractals 2003; 18:229–39. [6] Tang S, Chen L. Chaos in functional response host-parasitoid ecosystem models. Chaos, Solitons & Fractals 2002;13:875–84. [7] Zhang S, Chen L. Chaos in three species food chain system with impulsive perturbations. Chaos, Solitons & Fractals 2005;24: 73–83. [8] Varriale MC, Gomes AA. A study of a three species food chain. Ecol Moel 1998;110:119–33. [9] Kuznetsov YA, Inaldi SR. Remarks on food chains dynamics. Math Bios 1996;134:1–33. [10] McCann K, Yodis P. Bifurcation structure of a three-species food chains model. Theor Popu Biol 1995;48:93–125. [11] Gragnani A, DeFeo O, Rinaldi S. Food chain in the chemostat: relationships between mean yield and complex dynamics. Bull Math Biol 1998;60:703–19. [12] Upadhyay RK, Rai V. Crisis-limited chaotic dynamics in ecological system. Chaos, Solitons & Fractals 2001;12:205–18. [13] Rinaldi S, Casagrandi R, Gragnani A. Reduced order models for the prediction of the time of occurrence of extreme episodes. Chaos, Solitons & Fractals 2001;12:313–20. [14] Letellier C, Aziz-Alaoui MA. Analysis of the dynamics of a realistic ecological model. Chaos, Solitons & Fractals 2002;13:95–107. [15] Lenbury Y. Singular perturbation analysis of a model for a predator–prey system invaded by a parasite. Biosystems 1996;39: 251–62. [16] Ballinger G, Liu X. Permanence of population growth models with impulsive effects. Math Comput Model 1997;26:59–72. [17] Lenci S, Rega G. Periodic solutions and bifurcations in an impact inverted pendulum under impulsive excitation. Chaos, Solitons & Fractals 2000;11:2453–72. [18] Liu X, Chen L. Complex dynamics of Holling type II Lotka–Volterra predator–prey system with impulsive perturbations on the predator. Chaos, Solitons & Fractals 2003;16:311–20. [19] 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. [20] Zhang S, Dong L, Chen L. The study of predator–prey system with defensive ability of prey and impulsive perturbations on the predator. Chaos, Solitons & Fractals 2005;23:631–43. [21] Zhang S, Chen L. A Holling II functional response food chain model with impulsive perturbations. Chaos, Solitons & Fractals 2005;24:1269–78.
W. Wang et al. / Chaos, Solitons and Fractals 32 (2007) 1772–1785
1785
[22] Zhang S, Wang F, Chen L. A food chain model with impulsive perturbations and Holling IV functional response. Chaos, Solitons & Fractals 2005;26:855–66. [23] Zhang Y, Xiu Z, Chen L. Dynamic complexity of a two-prey one-predator system with impulsive effect. Chaos, Solitons & Fractals 2005;26:131–9. [24] Zhang Y, Xiu Z, Chen L. Chaos in a food chain chemostat with pulsed input and washout. Chaos, Solitons & Fractals 2005;26: 159–66. [25] 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. [26] Beddington JR. Mutual interference between parasites or predators and its effect on searching efficiency. J Animal Ecol 1975;44: 331–40. [27] Hwang TW. Global analysis of the predator–prey system with Beddington–DeAngelis functional response. J Math Anal Appl 2003;281:395–401. [28] Hwang TW. Uniqueness of limit cycles of the predator–prey system with Beddington–DeAngelis functional response. J Math Anal Appl 2004;290:113–22. [29] 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. [30] 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. [31] Dimitrov TD, Kojouharov VH. Complete mathematical analysis of predator–prey models with linear prey growth and Beddington–DeAngelis functional response. J Appl Math Comput 2005;162:523–38. [32] Li Z, Wang W, Wang H. The dynamics of a Beddington-type system with impulsive control strategy, Chaos, Solitons & Fractals, doi:10.1016/j.chaos.2005.08.195, in press. [33] Lakshmikantham V, Bainov DD, Simeonov PC. Theory of impulsive differential equations. Singapore: World Scientific; 1989. [34] Bainov DD, Simeonov PS. Impulsive differential equations: asymptotic properties of the solutions. Singapore: World Scientific; 1993. [35] Wang W. Computer algebra system and symbolic computation. Lanzhou, China: Gansu Sci-Tech press; 2004 (in Chinese). [36] Wang W, Lian X. A new algorithm for symbolic integral with application. Appl Math Comput 2005;162:949–68. [37] Wang W, Lin C. A new algorithm for integral of trigonometric functions with mechanization. Appl Math Comput 2005;164: 71–82. [38] Li Z, Wang W. Mechanization for solving SPP by reducing order method. Appl Math Comput 2005;169:1028–37.
The dynamic complexity of a three-species Beddington-type food chain with impulsive control strategy q Weiming Wang a
a,b
, Hailing Wang c, Zhenqing Li
a,*
Laboratory of Quantitative Vegetation Ecology, Institute of Botany, The Chinese Academy of Sciences, Beijing 100093, China b School of Mathematics and Information Science, Wenzhou University, Wenzhou Zhejiang 325035, China c School of Applied Mathematics, University of Electronic Science and Technology of China, Chengdu Sichuan 610054, China Accepted 1 December 2005
Communicated by Prof. L. Marck-Crnjac
Abstract In this paper, by using theories and methods of ecology and ordinary differential equation, the dynamics complexity of 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 Floquet theory of impulsive equation and small amplitude perturbation skills. Furthermore, by using the method of numerical simulation with the international software Maple, the influence of the impulsive perturbations on the inherent oscillation is investigated, which shows rich dynamics, such as quasi-periodic oscillation, narrow periodic window, wide periodic window, chaotic bands, period doubling bifurcation, symmetry-breaking pitchfork bifurcation, period-halving bifurcation and crises, etc. The numerical results indicate that computer simulation is a useful method for studying the complex dynamic systems. Ó 2006 Elsevier Ltd. All rights reserved.
1. Introduction The recent 40 years have seen an explosion of interest in the study of nonlinear dynamical systems. More recently, dynamical systems have benefited from an infusion of interests and techniques from a variety of fields [1]. It is currently very much in vogue to study chaotic behavior, and much of it is prompted by new and potential application in different fields [2–24]. In the recent decade, simple multi-species systems comprising of a three-species food chain were discussed by a number of researchers, and there is also a large body of literature on models with perturbations, particularly in prey–predator systems [3,6–24].
q
This research was supported by the National Natural Science Foundation of China (Grant Nos. 30270249 and 10371090) and the Knowledge Innovation Project of the Chinese Academy of Sciences (Grant No. KSCX2-SW-107). * Corresponding author. E-mail addresses: [email protected] (W. Wang), [email protected] (H. Wang), [email protected] (Z. Li). 0960-0779/$ - see front matter Ó 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2005.12.025
W. Wang et al. / Chaos, Solitons and Fractals 32 (2007) 1772–1785
1773
A well-known model of such systems is the predator–prey model with a Beddington-type functional response [25–32], which can be described by the following differential equations: 8 dxðtÞ c1 xðtÞyðtÞ > > > dt ¼ xðtÞða bxðtÞÞ a þ xðtÞ þ b yðtÞ ; > > 1 1 > > > < dyðtÞ k 1 c1 xðtÞyðtÞ c2 yðtÞzðtÞ d 1 yðtÞ; ¼ ð1:1Þ > dt a1 þ xðtÞ þ b1 yðtÞ a2 þ yðtÞ þ b2 zðtÞ > > > > > > dzðtÞ k 2 c2 yðtÞzðtÞ > : d 2 zðtÞ. ¼ dt a2 þ yðtÞ þ b2 zðtÞ In our paper, we consider an impulsive differential equation—a predator–prey model with a Beddington-type functional response, by introducing a proportion periodic impulsive catching or poisoning for the prey populations, and a constant periodic releasing for the predator at different fixed moment, respectively, what follows is the corresponding equations: 8 9 dxðtÞ c1 xðtÞyðtÞ > > > > ¼ xðtÞða bxðtÞÞ > > > dt > a1 þ xðtÞ þ b1 yðtÞ > > > > > > > > > dyðtÞ = > k c xðtÞyðtÞ c yðtÞzðtÞ > 1 1 2 > d yðtÞ ¼ t 6¼ ðn þ l 1ÞT ; t 6¼ nT ; > 1 > > > dt a1 þ xðtÞ þ b1 yðtÞ a2 þ yðtÞ þ b2 zðtÞ > > > > > > > > > > > dzðtÞ k 2 c2 yðtÞzðtÞ > > > > ; ¼ d zðtÞ 2 > > dt a þ yðtÞ þ b zðtÞ > 2 2 < 9 ð1:2Þ DxðtÞ ¼ d1 xðtÞ > > = > > > > t ¼ ðn þ l 1ÞT ; DyðtÞ ¼ d2 yðtÞ > > > > ; > > > DzðtÞ ¼ d3 zðtÞ > > > > 9 > > > DxðtÞ ¼ 0 > > = > > > > DyðtÞ ¼ 0 t ¼ nT ; > > > > ; : DzðtÞ ¼ p where x(t), y(t) and z(t) are functions of time representing population density of the prey, mid-level and top predator, respectively, all parameters are positive constants, the parameter ki (i = 1, 2) are conversion efficiencies, di (i = 1, 2) are the mortality rates of the mid-level predator and the top predator, ci (i = 1, 2) are the maximum numbers of prey that can be eaten by a predator per unit of time, ai (i = 1, 2) are saturation constants and bi (i = 1, 2) scale the impact of the predator interference. T is the period of the impulsive effect and n 2 N. Observe the simple relation of these three species: z only prey on y and y prey on x and nutrient recycling is not accounted for. This simple relation produces the socalled simple food chain system with impulsive perturbations. With model (1.2), we can take into account the possible exterior effects under which the population density changes very rapidly. An impulsive increase of the predator population density is possible by artificial breeding of the species or release some species (p > 0). The rest of this paper is organized as follows: In Section 2, some notations and theorems are given. In Section 3, the results of computer simulation are shown, and moreover, these results are discussed briefly. Finally, conclusions and remarks are given in Section 4.
2. Mathematical analysis In this section, we will give some definitions, notations and some lemmas which will be useful for our main results. Let R+ = [0, 1), R3þ ¼ fx 2 R3 : x P 0g; X ¼ int R3þ , N be the set of all non-negative integers. The map f = (f1, f2, f3)T defined by the system (1.1). Let V : Rþ R3þ ! Rþ , then V is said to belong to class V0 if (1) V is continuous in ððn 1ÞT ; ðn þ l lÞT R3þ and ððn þ l 1ÞT ; nT R3þ , and for all x 2 R3þ , n 2 N, limðt;yÞ!ððnþl1ÞT þ ;xÞ V ðt; yÞ ¼ V ððn þ l lÞT þ ; xÞ and limt;y!ðnT þ ;xÞ V ðt; yÞ ¼ V ðnT þ ; xÞ exist; (2) V is locally Lipschitzian in x.
1774
W. Wang et al. / Chaos, Solitons and Fractals 32 (2007) 1772–1785
Definition 2.1. V 2 V0, then for ðt; xÞ 2 ððn 1ÞT ; ðn þ l 1ÞT R3þ and ððn þ l 1ÞT ; nT R3þ , 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 The solution of system (1.2) is a piecewise continuous function X ðtÞ : Rþ ! R3þ ; X ðtÞ is continuous on T, T = ((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). It is easy to prove the following lemmas. Lemma 2.1. Let X(t) be a solution to system (1.2) with X(0+) P 0, then X(t) P 0 for all t P 0. And furthermore, 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, z(t) 6 M for each solution (x(t), y(t), z(t)) of (1.2) with sufficiently large t. Then, we give some basic properties about the following subsystem of system (1.2). 8 dzðtÞ > > ¼ d 2 zðtÞ; t 6¼ ðn þ l 1ÞT ; t 6¼ nT ; > > dt > > < þ zðt Þ ¼ ð1 d3 ÞzðtÞ; t ¼ ðn þ l 1ÞT ; > > zðtþ Þ ¼ zðtÞ þ p; t ¼ nT ; > > > > : þ zð0 Þ ¼ z0 .
ð2:1Þ
Clearly 8 p expðd 2 ðt ðn 1ÞT ÞÞ > > ; ðn 1ÞT < t 6 ðn þ l 1ÞT ; > < 1 ð1 d3 Þ expðd 2 T Þ z ðtÞ ¼ > pð1 d3 Þ expðd 2 ðt ðn 1ÞT ÞÞ > > ; ðn þ l 1ÞT < t 6 nT ; : 1 ð1 d3 Þ expðd 2 T Þ p pð1 d3 Þ expðd 2 lT Þ þ þ þ ; z ðlT Þ ¼ n 2 N ; z ð0 Þ ¼ z ðnT Þ ¼ 1 ð1 d3 Þ expðd 2 T Þ 1 ð1 d3 Þ expðd 2 T Þ is a positive periodic solution of system (1.2). Since 8 p > ðn1Þ þ > > expðd 2 tÞ þ z ðtÞ; ðn 1ÞT < t 6 ðn þ l 1ÞT ; Þ zð0 Þ ð1 d 3 < 1 ð1 d3 Þ expðd 2 T Þ zðtÞ ¼ > p n > þ > expðd 2 tÞ þ z ðtÞ; ðn þ l 1ÞT < t 6 nT Þ zð0 Þ ð1 d : 3 1 ð1 d3 Þ expðd 2 T Þ is the solution of system (1.2) with initial value z0 P 0, where t 2 T. Then we can get: Lemma 2.3. Let z*(t) be a positive periodic solution of system (1.2) and all solution z(t) of system (1.2) with z0 P 0, we have jz(t) z*(t)j ! 0, when t ! 1. Therefore, we obtain the complete expression for prey and mid-level predator eradication periodic solution (0, 0, z*(t)) and ab ; 0; z ðtÞ to system (1.2) for t 2 T. Now, we study the stability of prey and mid-level predator eradication periodic solution to system (1.2). Theorem 2.1. Let (x(t), y(t), z(t)) be any solution to (1.2), then (1) (0, 0, z*(t)) is unstable; (2) ðab ; 0; z ðtÞÞ is locally asymptotically stable provided that
c2 ða1 b þ aÞðA1 A2 þ A3 A4 d 1 d 2 b2 T Þ þ k 1 c1 ab2 d 2 T 1 < ln ; b2 d 2 ða1 b þ aÞ 1 d2
W. Wang et al. / Chaos, Solitons and Fractals 32 (2007) 1772–1785
where
1775
a2 expðd 2 T Þ þ b2 p expðð1 lÞd 2 T Þ a2 þ a2 d3 ; expðd 2 T Þ þ d3 1 ða2 þ b2 Þ expðd 2 T Þ a2 þ a2 d3 ; A2 ¼ ln expðd 2 T Þ þ d3 1 a2 expðd 2 T Þ a2 þ a2 d3 þ b2 p þ b2 pd3 ; A3 ¼ ln expðd 2 T Þ þ d3 1 a2 expðd 2 T Þ þ ð1 d3 Þb2 p expðð1 lÞd 2 T Þ a2 þ a2 d3 . A4 ¼ ln expðd 2 T Þ þ d3 1 A1 ¼ ln
Proof. (1) The local stability of periodic solution (0, 0, z*(t)) is determined by considering the behavior of small amplitude perturbations of the solution. Define xðtÞ ¼ uðtÞ; Then
zðtÞ ¼ wðtÞ þ z ðtÞ.
yðtÞ ¼ vðtÞ;
1 0 1 uð0Þ uðtÞ C B C B @ vðtÞ A ¼ UðtÞ@ vð0Þ A; wð0Þ wðtÞ 0
0 6 t < T;
where U(t) satisfies 0 a 0 c 2 z B 0 d1 dU B a2 þ b2 z ¼B B dt @ k 2 c 2 z 0 a2 þ b2 z
0
1
C 0 C CUðtÞ C A d 2
and U(0) = I, the identity matrix. The linearization of system (1.2) becomes 1 0 1 0 10 uððn þ l 1ÞT Þ 0 0 1 d1 uððn þ l 1ÞT þ Þ B C B C B C 1 d2 0 A@ vððn þ l 1ÞT Þ A @ vððn þ l 1ÞT þ Þ A ¼ @ 0 wððn þ l 1ÞT Þ 0 0 1 d3 wððn þ l 1ÞT þ Þ and 1 10 1 0 uðnT Þ 1 0 0 uðnT þ Þ C CB C B B @ vðnT þ Þ A ¼ @ 0 1 0 A@ vðnT Þ A. þ wðnT Þ 0 0 1 wðnT Þ 0
Hence, if all eigenvalues of 0 1 d1 0 B H¼@ 0 1 d2 0
0
10
1 1 0 0 CB C A@ 0 1 0 AUðT Þ
0 0 1 d3
0 0 1
have absolute less than one, the periodic solution (0, 0, z*(t)) is locally stable. But one of the Floquet multipliers is u1 = exp((1 d1)aT) > 0, according to the Floquet theory, the solution X*(t) is unstable. (2) The local stability of periodic solution ðab ; 0; z ðtÞÞ can be similarly discussed. Define a xðtÞ ¼ uðtÞ þ ; yðtÞ ¼ vðtÞ; zðtÞ ¼ wðtÞ þ z ðtÞ. b then
0
1 0 1 uðtÞ uð0Þ B C B C @ vðtÞ A ¼ UðtÞ@ vð0Þ A; wðtÞ wð0Þ
0 6 t < T;
1776
W. Wang et al. / Chaos, Solitons and Fractals 32 (2007) 1772–1785
where U(t) satisfies 0 c1 ab a B a1 þ ab B B a k 1 c1 b dU B c 2 z ¼B 0 d1 a B dt a1 þ b a2 þ b2 z B @ k 2 c 2 z 0 a2 þ b2 z
0
1
C C C C 0 CUðtÞ C C A d 2
and U(0) = I, the identity matrix. The linearization of system (1.2) becomes 1 10 1 0 0 uððn þ l 1ÞT Þ 1 d1 0 0 uððn þ l 1ÞT þ Þ C CB C B B 1 d2 0 A@ vððn þ l 1ÞT Þ A @ vððn þ l 1ÞT þ Þ A ¼ @ 0 wððn þ l 1ÞT Þ 0 0 1 d3 wððn þ l 1ÞT þ Þ and 0
1 0 10 1 1 0 0 uðnT Þ uðnT þ Þ B C B CB C @ vðnT þ Þ A ¼ @ 0 1 0 A@ vðnT Þ A. 0 0 1 wðnT Þ wðnT þ Þ The stability of the periodic solution ðab ; 0; z ðtÞÞ is determined by the eigenvalues of 1 0 10 1 0 0 0 0 1 d1 C B CB 1 d2 0 A@ 0 1 0 AUðT Þ H¼@ 0 0 0 1 0 0 1 d3 which are u1 ¼ ð1 d1 Þ expðaT Þ < 1; Z T k 1 c1 ab c 2 z d u2 ¼ ð1 d2 Þ exp 1 dt ; a1 þ ab a2 þ b2 z 0 u3 ¼ ð1 d3 Þ expðd 2 T Þ < 1. According to Floquet theory, the solution ðab ; 0; z ðtÞÞ is locally stable if ju2j < 1, that is to say, c2 ða1 b þ aÞðA1 A2 þ A3 A4 d 1 d 2 b2 T Þ þ k 1 c1 ab2 d 2 T 1 < ln . b2 d 2 ða1 b þ aÞ 1 d2 This completes the proof.
h
3. Computer simulation It is well known that the correspondingly continuous system (1.2) cannot be solved explicitly. So we have to study the system (1.2) numerically synthetically and research the long-term behavior of the solution with Maple. Maple, an international mathematical software, is an interactive computer algebra system with great ability of symbolic evaluation, numerical calculation, coping with graphics, etc. Its powerful function library and unique interior programming language provide science calculation and programming with friendly platform [33–38]. We consider three sets of parameters for our analysis as follows. Case A: a = 5, b = 0.5, c1 = 2, c2 = 2, a1 = 0.1, a2 = 0.5, b1 = 1, b2 = 0.1, d1 = 0.1, d2 = 0.9, k1 = 0.7, k2 = 0.5, T = 8, d1 = 0.2, d2 = 0.0001, d3 = 0.03, l = 0.6. Case B: a = 2, b = 0.1, c1 = 2, c2 = 3.5, a1 = 0.1, a2 = 0.7, b1 = 7, b2 = 0.1, d1 = 0.1, d2 = 0.9, k1 = 0.7, k2 = 0.5, T = 16, d1 = 0.2, d2 = 0.0001, d3 = 0.07, l = 0.6. Case C: a = 3, b = 0.25, c1 = 2, c2 = 3.5, a1 = 0.1, a2 = 0.5, b1 = 1, b2 = 0.1, d1 = 0.1, d2 = 0.9, k1 = 0.7, k2 = 0.5, T = 5, d1 = 0.2, d2 = 0.0001, d3 = 0.08, l = 0.68.
W. Wang et al. / Chaos, Solitons and Fractals 32 (2007) 1772–1785
1777
For Case A, in Fig. 1, we have displayed typical bifurcation diagrams for the prey population x(t), the mid-level predator population y(t) and the top predator population z(t) of the system (1.2) as p increasing from 0 to 3.5 with initial values X(0) = (15, 1, 0.4). As the similarity of the bifurcation diagrams, we only give the magnified parts of Fig. 1(a) in Fig. 2, and the windows of periodic behavior are more clearly visible. As p increases, the resulting bifurcation diagrams clearly show that system (1.2) has rich dynamics including: quasi-periodic oscillation, narrow periodic window, wide periodic window, chaotic bands, crises (the phenomenon of ‘‘crisis’’ in which chaotic attractors can suddenly appear or disappear, or change size discontinuously as a parameter smoothly varies), period doubling bifurcation, period-halving bifurcation, symmetry-breaking pitchfork bifurcation, tangent bifurcation, etc. In Fig. 1(a) and Fig. 2, the bifurcation diagrams are drawn as a function of x. As the bifurcation diagrams shown, the dynamics behavior of prey population x may dramatically be affected by small changes in the value of the parameter p. For 0 < p < 1.342, Fig. 2(a) shows that a quasi-periodic solution oscillation, in which there are frequency-lockings, symmetry-breaking pitchfork bifurcation and tangent bifurcation (Fig. 3). When the parameter p is slightly increased beyond p > 1.342, there is a 4T-period solution of system (1.2) appears (Fig. 4), and after this, the parameter p is slightly
Fig. 1. Bifurcation diagrams of system (1.2) with a = 5, b = 0.5, c1 = 2, c2 = 2, a1 = 0.1, a2 = 0.5, b1 = 1, b2 = 0.1, d1 = 0.1, d2 = 0.9, k1 = 0.7, k2 = 0.5, T = 8, d1 = 0.2, d2 = 0.0001, d3 = 0.03, l = 0.6, and with initial values X(0) = (15, 1, 0.4). (a) Prey population x, (b) mid-level predator population y and (c) the top population z, with p over [0, 3.5].
Fig. 2. The magnified parts of Fig. 1(a): (a) p over [0.5, 1.1] and (b) p over [1.2, 2.4], respectively.
Fig. 3. Quasi-periodic solutions oscillation of system (1.2) with initial values X(0) = (15, 1, 0.4), there are frequency-lockings, symmetry-breaking pitchfork bifurcation and tangent bifurcation in it. (a) p = 0.25, (b) p = 0.5, (c) p = 0.8 and (d) p = 1.3.
1778
W. Wang et al. / Chaos, Solitons and Fractals 32 (2007) 1772–1785
Fig. 4. A 4T-periodic solution when the parameter p is slightly increased beyond p > 1.342.
increased beyond p > 2.09, the 4T-period solution becomes unstable and there is a cascade of period doubling bifurcation leads to chaos, which can been seen in Fig. 5, when p increased beyond 2.254, there is a cascade of period-halving bifurcation leads to 4T-period solution (Fig. 6). There is a crisis appears again when p 2.313, leads to a chaos band, and then, the chaos disappears, a T-period solution of system (1.2) appears (Fig. 7). If p > 2.665954781, as Theorem 2.1 shows, the mid-level predator will go extinct y(t) ! 0 as t ! 1 (Fig. 8).
Fig. 5. Period doubling bifurcation leads to chaos of system (1.2) with initial values X(0) = (15, 1, 0.4), (a) phase portrait of 8T-period solution for p = 2.1 and (b) phase portrait of the four attractors for p = 2.18.
Fig. 6. Period-halving bifurcation leads to chaos of system (1.2) with initial values X(0) = (15, 1, 0.4), (a) phase portrait of 4T-period solution for p = 2.223 and (b) phase portrait of the solutions enter chaotic region for p = 2.23.
Fig. 7. Crises are shown: when the parameter p is slightly increased beyond p = 2.313, the chaos suddenly disappears, (a) p = 2.31 and (b) p = 2.5, from (a) to (b) there is a crisis that the chaos suddenly disappears.
W. Wang et al. / Chaos, Solitons and Fractals 32 (2007) 1772–1785
1779
Fig. 8. When p > pmax ¼ 2:665954781; ab ; 0; z ðtÞ is locally asymptotically stable, that is to say, the mid-level predator will go extinction: (a) time series of the prey population x(t), (b) time series of the mid-level predator population y(t), when p > pmax = 2.665954781, y(t) ! 0, as t ! +1, (c) time series of the top predator population z(t) and (d) phase portrait of system (1.2) of p = 2.666.
For Case B, in Fig. 9, we have displayed typical bifurcation diagrams for the prey population x(t), the mid-level predator population y(t) and the top predator population z(t) of the system (1.2) as p increasing from 0 to 4 with initial values X(0) = (15, 1, 0.4). As the similarity of the bifurcation diagrams, we only give the magnified parts of Fig. 9(a) in Fig. 10. As p increases, the resulting bifurcation diagrams clearly show that system (1.2) has rich dynamics including: quasi-periodic oscillation, period doubling bifurcation, symmetry-breaking pitchfork bifurcation, tangent bifurcation, crisis, chaotic bands, narrow periodic window, wide periodic window, period-halving bifurcation, etc. In Fig. 9(a) and Fig. 10, the bifurcation diagrams are drawn as a function of x. As the bifurcation diagrams shown, the dynamics behavior of prey population x may dramatically be affected by small changes in the value of the parameter p. For 0 < p < 0.069, Fig. 9(a) shows that a quasi-periodic solution oscillation, when 0.069 < p < 0.153, T-period solution of system (1.2) is still stable (Fig. 11). When p > 0.153, it becomes unstable and there is a cascade of period doubling bifurcation leads to chaos, it enters a chaotic region (Fig. 12). Unlike the common period doubling route to chaos, here we observe the frequent occurrence of sudden changes of attractors (crises) in the range of [0.254, 0.321], where
Fig. 9. Bifurcation diagrams of system (1.2) with a = 2, b = 0.1, c1 = 2, c2 = 3.5, a1 = 0.1, a2 = 0.7, b1 = 7, b2 = 0.1, d1 = 0.1, d2 = 0.9, k1 = 0.7, k2 = 0.5, T = 16, d1 = 0.2, d2 = 0.0001, d3 = 0.07, l = 0.6, and with initial values X(0) = (15, 1, 0.4). (a) Prey population x, (b) mid-level predator population y and (c) the top population z, with p over [0, 4].
Fig. 10. The magnified parts of Fig. 9(a): (a) p over [0, 0.5], (b) p over [0.5, 1.1] and (c) p over [1.2, 2.4], respectively.
1780
W. Wang et al. / Chaos, Solitons and Fractals 32 (2007) 1772–1785
Fig. 11. Quasi-periodic solutions oscillation of system (1.2) with initial values X(0) = (15, 1, 0.4), (a) p = 0.068 and (b) p = 0.13.
Fig. 12. Period doubling bifurcation leads to chaos of system (1.2) with initial values X(0) = (15, 1, 0.4), (a) phase portrait of 2T-period solution for p = 0.2, (b) phase portrait of 4T-period solution for p = 0.35, (c) phase portrait of 8T-period solution for p = 0.38 and (d) phase portrait of the solutions enter chaotic region for p = 0.42.
Fig. 13. Crises are shown: when the parameter p is slightly increased beyond p = 0.677, the chaos suddenly disappears and multiply attractors appears, (a) p = 0.68, (b) p = 0.75 and (c) p = 0.85.
Fig. 14. Crises are shown: when the parameter p is slightly increased beyond p = 0.923, the chaos suddenly disappears and multiply attractors appears, (a) p = 0.92, (b) p = 0.94 and (c) p = 1.0.
multiple attractors coexist. And the same instance occur on the range of [0.677, 0.858] and [0.923, 0.990] (Figs. 13 and 14). When the parameter p is slightly increased beyond p > 1.239, the chaos disappears, period-halving bifurcation leads
W. Wang et al. / Chaos, Solitons and Fractals 32 (2007) 1772–1785
1781
to T-period solution (Fig. 15). If p > 3.751536111, as Theorem 2.1 shows, the mid-level predator will go extinct y(t) ! 0 as t ! 1 (Fig. 16). For Case C, in Fig. 17, we have displayed typical bifurcation diagrams for the prey population x(t), the mid-level predator population y(t) and the top predator population z(t) of the system (1.2) as p increasing from 0 to 1 with initial values X(0) = (15, 1, 0.4). As the similarity of the bifurcation diagrams, we only give the magnified parts of Fig. 17(a) in Fig. 18. As p increases, the resulting bifurcation diagrams clearly show that system (1.2) has rich dynamics including: quasi-periodic oscillation, period doubling bifurcation, symmetry-breaking pitchfork bifurcation, tangent bifurcation, crisis, chaotic bands, narrow periodic window, wide periodic window, period-halving bifurcation, etc. In Figs. 17(a) and 18, the bifurcation diagrams are drawn as a function of x. As the bifurcation diagrams shown, the dynamics behavior of prey population x may dramatically be affected by small changes in the value of the parameter p.
Fig. 15. Period-halving bifurcation leads to T-period solution, (a) phase portrait of the solutions in chaotic region for p = 1.24, (b) phase portrait of 4T-period solution for p = 1.4, (c) phase portrait of 2T-period solution for p = 1.6 and (d) phase portrait of T-period solution for p = 3.7.
Fig. 16. When p > pmax ¼ 3:751536111; ab ; 0; z ðtÞ is locally asymptotically stable, that is to say, the mid-level predator will go extinction: (a) time series of the prey population x(t), (b) time series of the mid-level predator population y(t), when p > pmax = 3.751536111, y(t) ! 0, as t ! +1, (c) time series of the top predator population z(t) and (d) phase portrait of system (1.2) of p = 3.8.
Fig. 17. Bifurcation diagrams of system (1.2) with a = 3, b = 0.25, c1 = 2, c2 = 3.5, a1 = 0.1, a2 = 0.5, b1 = 1, b2 = 0.1, d1 = 0.1, d2 = 0.9, k1 = 0.7, k2 = 0.5, T = 5, d1 = 0.2, d2 = 0.0001, d3 = 0.08, l = 0.68 and with initial values X(0) = (15,1,0.4). (a) prey population x, (b) mid-level predator population y and (c) the top population z, with p over [0, 1].
1782
W. Wang et al. / Chaos, Solitons and Fractals 32 (2007) 1772–1785
Fig. 18. The magnified parts of Fig. 17(a): (a) p over [0, 0.04], (b) p over [0.04, 0.2], (c) p over [0.2, 0.28] and (d) p over [0.3, 0.6], respectively.
Fig. 19. Quasi-periodic solutions oscillation of system (1.2) with initial values X(0) = (15, 1, 0.4), (a) p = 0.0015 and (b) p = 0.04.
Fig. 20. Chaos regions with frequency-lockings, symmetry-breaking pitchfork bifurcation and tangent bifurcation in it. (a) p = 0.015, (b) p = 0.02, (c) p = 0.035 and (d) p = 0.08.
Fig. 21. Period doubling bifurcation leads to chaos of system (1.2) with initial values X(0) = (15, 1, 0.4), (a) phase portrait of 2T-period solution for p = 0.11, (b) phase portrait of 4T-period solution for p = 0.12, (c) phase portrait of 8T-period solution for p = 0.148 and (d) phase portrait of the solutions enter chaotic region for p = 0.165.
W. Wang et al. / Chaos, Solitons and Fractals 32 (2007) 1772–1785
1783
For 0 < p < 0.0058, Fig. 17(a) shows that a quasi-periodic solution oscillation (Fig. 19), when p > 0.0123, a crisis occurs, that is to say, leads to a chaos region, in which there are frequency-lockings, symmetry-breaking pitchfork bifurcation, tangent bifurcation, narrow periodic window and wide periodic window Fig. 20. When the parameter p is slightly increased beyond p > 0.098, the chaos disappears, period doubling bifurcation leads to another chaos region Fig. 21,
Fig. 22. Period doubling bifurcation leads to chaos of system (1.2) with initial values X(0) = (15, 1, 0.4), (a) phase portrait of 2T-period solution for p = 0.176, (b) phase portrait of 4T-period solution for p = 0.18 and (c) phase portrait of the solutions enter chaotic region for p = 0.19.
Fig. 23. Chaos regions with frequency-lockings, symmetry-breaking pitchfork bifurcation and tangent bifurcation in it, that is to say, there are several crises in the region [0.22, 0.7]. (a) p = 0.222, (b) p = 0.225, (c) p = 0.25, (d) p = 0.4, (e) p = 0.55 and (f) p = 0.7.
Fig. 24. When p > pmax ¼ 0:7183181912; ab ; 0; z ðtÞ is locally asymptotically stable, that is to say, the mid-level predator will go extinction: (a) time series of the prey population x(t), (b) time series of the mid-level predator population y(t), when p > pmax = 0.7183181912, y(t) ! 0, as t ! +1, (c) time series of the top predator population z(t) and (d) phase portrait of system (1.2) of p = 0.72.
1784
W. Wang et al. / Chaos, Solitons and Fractals 32 (2007) 1772–1785
then the chaos disappears again when p > 0.173, another period doubling bifurcation leads to the third chaos region Fig. 22, in the following, when the p > 0.228, it becomes unstable and there is a cascade of period-halving bifurcation leads to the fourth chaos region, and in this course, multiply attractors coexist in the range of [0.239, 0.244], which we can observe in Fig. 18. Then the last crisis appears when p is slightly increased beyond p > 0.450, after this, period-halving bifurcation leads to T-period solution (Fig. 23). If p > 0.7183181912, as Theorem 2.1 shows, the mid-level predator will go extinct y(t) ! 0 as t ! 1 (Fig. 24).
4. Conclusions and remarks We have considered the dynamic behavior of a three-species food chain with Beddington-type functional response and periodic constant impulsive perturbations of top predator. Conditions for extinction of prey and mid-level predator is established and the boundary of the system (1.2) is given. We have observed in the computer simulation that the dynamics of a population may dramatically be affected by small changes in the value of the parameter p, at the same time, we can see from the bifurcation diagrams, with different parameter sets, we have different bifurcation diagrams, that is to say, we have different routes to chaos. On the other hand, for the complexity of the ecosystem, we cannot get its analytical solution. Then, we have to solve the numerical solutions. It is well known that the method of computer simulation is one of the best methods. There are some international mathematical softwares, 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 studying the complex dynamical systems.
References [1] Devaney RL. An introduction to chaotic dynamical systems. Menlo Park, CA: Benjamin/Cummings publishing company; 1986. [2] May RM. Necessity and chance: deterministic chaos in ecology and evolution. Am Math Soc 1995;32:291–308. [3] 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. [4] 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. [5] Gakkhar S, Naji MA. Order and chaos in predator to prey ratio-dependent food chain. Chaos, Solitons & Fractals 2003; 18:229–39. [6] Tang S, Chen L. Chaos in functional response host-parasitoid ecosystem models. Chaos, Solitons & Fractals 2002;13:875–84. [7] Zhang S, Chen L. Chaos in three species food chain system with impulsive perturbations. Chaos, Solitons & Fractals 2005;24: 73–83. [8] Varriale MC, Gomes AA. A study of a three species food chain. Ecol Moel 1998;110:119–33. [9] Kuznetsov YA, Inaldi SR. Remarks on food chains dynamics. Math Bios 1996;134:1–33. [10] McCann K, Yodis P. Bifurcation structure of a three-species food chains model. Theor Popu Biol 1995;48:93–125. [11] Gragnani A, DeFeo O, Rinaldi S. Food chain in the chemostat: relationships between mean yield and complex dynamics. Bull Math Biol 1998;60:703–19. [12] Upadhyay RK, Rai V. Crisis-limited chaotic dynamics in ecological system. Chaos, Solitons & Fractals 2001;12:205–18. [13] Rinaldi S, Casagrandi R, Gragnani A. Reduced order models for the prediction of the time of occurrence of extreme episodes. Chaos, Solitons & Fractals 2001;12:313–20. [14] Letellier C, Aziz-Alaoui MA. Analysis of the dynamics of a realistic ecological model. Chaos, Solitons & Fractals 2002;13:95–107. [15] Lenbury Y. Singular perturbation analysis of a model for a predator–prey system invaded by a parasite. Biosystems 1996;39: 251–62. [16] Ballinger G, Liu X. Permanence of population growth models with impulsive effects. Math Comput Model 1997;26:59–72. [17] Lenci S, Rega G. Periodic solutions and bifurcations in an impact inverted pendulum under impulsive excitation. Chaos, Solitons & Fractals 2000;11:2453–72. [18] Liu X, Chen L. Complex dynamics of Holling type II Lotka–Volterra predator–prey system with impulsive perturbations on the predator. Chaos, Solitons & Fractals 2003;16:311–20. [19] 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. [20] Zhang S, Dong L, Chen L. The study of predator–prey system with defensive ability of prey and impulsive perturbations on the predator. Chaos, Solitons & Fractals 2005;23:631–43. [21] Zhang S, Chen L. A Holling II functional response food chain model with impulsive perturbations. Chaos, Solitons & Fractals 2005;24:1269–78.
W. Wang et al. / Chaos, Solitons and Fractals 32 (2007) 1772–1785
1785
[22] Zhang S, Wang F, Chen L. A food chain model with impulsive perturbations and Holling IV functional response. Chaos, Solitons & Fractals 2005;26:855–66. [23] Zhang Y, Xiu Z, Chen L. Dynamic complexity of a two-prey one-predator system with impulsive effect. Chaos, Solitons & Fractals 2005;26:131–9. [24] Zhang Y, Xiu Z, Chen L. Chaos in a food chain chemostat with pulsed input and washout. Chaos, Solitons & Fractals 2005;26: 159–66. [25] 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. [26] Beddington JR. Mutual interference between parasites or predators and its effect on searching efficiency. J Animal Ecol 1975;44: 331–40. [27] Hwang TW. Global analysis of the predator–prey system with Beddington–DeAngelis functional response. J Math Anal Appl 2003;281:395–401. [28] Hwang TW. Uniqueness of limit cycles of the predator–prey system with Beddington–DeAngelis functional response. J Math Anal Appl 2004;290:113–22. [29] 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. [30] 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. [31] Dimitrov TD, Kojouharov VH. Complete mathematical analysis of predator–prey models with linear prey growth and Beddington–DeAngelis functional response. J Appl Math Comput 2005;162:523–38. [32] Li Z, Wang W, Wang H. The dynamics of a Beddington-type system with impulsive control strategy, Chaos, Solitons & Fractals, doi:10.1016/j.chaos.2005.08.195, in press. [33] Lakshmikantham V, Bainov DD, Simeonov PC. Theory of impulsive differential equations. Singapore: World Scientific; 1989. [34] Bainov DD, Simeonov PS. Impulsive differential equations: asymptotic properties of the solutions. Singapore: World Scientific; 1993. [35] Wang W. Computer algebra system and symbolic computation. Lanzhou, China: Gansu Sci-Tech press; 2004 (in Chinese). [36] Wang W, Lian X. A new algorithm for symbolic integral with application. Appl Math Comput 2005;162:949–68. [37] Wang W, Lin C. A new algorithm for integral of trigonometric functions with mechanization. Appl Math Comput 2005;164: 71–82. [38] Li Z, Wang W. Mechanization for solving SPP by reducing order method. Appl Math Comput 2005;169:1028–37.