- Email: [email protected]
Stability and Hopf bifurcation analysis on a delayed Leslie-Gower predator–prey system incorporating a prey refuge
Applied Mathematics and Computation 219 (2013) 4576–4589
Contents lists available at SciVerse ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Stability and Hopf bifurcation analysis on a delayed Leslie-Gower predator–prey system incorporating a prey refuge q Yongkun Li a,⇑, Changzhao Li b a b
Department of Mathematics, Yunnan University, Kunming, Yunnan 650091, People’s Republic of China Department of Applied Mathematics, Kunming University of Science and Technology, Kunming 650093, Yunnan, People’s Republic of China
a r t i c l e
i n f o
a b s t r a c t In this paper, a modified Leslie-Gower predator–prey system with time delays is investigated, where the time delays are regarded as bifurcation parameters. By analyzing the corresponding characteristic equation, the local stability of a positive equilibrium is considered. Moreover, we show that Hopf bifurcations occur when time delay crosses some critical values. By deriving the equation describing the flow on the center manifold, we give the formula for determining the direction of the Hopf bifurcation and the stability of bifurcating periodic solutions. In addition, we also try on the global existence of periodic solutions by using the global Hopf bifurcation result of Wu [J. Wu, Symmetric functional differential equations and neural networks with memory, Trans. Amer. Math. Soc. 350 (1998) 4799–4838.] for functional differential equations. Numerical simulations are carried out to illustrate the theoretical results and they show that the time delays in the system under consideration can destroy the stability of the system. Ó 2012 Elsevier Inc. All rights reserved.
Keywords: Predator–prey system Time delay Hopf bifurcation Leslie-Gower Stability
1. Introduction The dynamic relationship between predators and their preys has long been and will continue to be one of the dominant themes in both ecology and mathematical ecology due to its universal existence and importance [1,2]. Considering the reduction in a predator population has a reciprocal relationship with per capita availability of its preferred food, Leslie [3] introduced the following predator–prey model where the ’’carrying capacity’’ of the predator’s environment is proportional to the number of prey:
dH ¼ ðr 1 a1 P qðHÞÞH; dt
dP ¼ dt
r 2 a2
P P; H
ð1:1Þ
where H is the density of prey species at time t, and P is the density of the predator species at time t, qðHÞ is the so-called predator functional response to prey and satisfies the usual properties that make system (1.1) a predator–prey one. In the last decade, based on system (1.1) with different predator functional response to prey, many excellent works have been done, for example, see [1,2,4–7] and the references cited therein. For more results on predator–prey systems, one could refer to [8–10] for details. As was pointed out by Tapan Kumar Kar [11], mite predator–prey interactions often exhibit spatial refugia which afford the prey some degree of protection from predation and reduce the chance of extinction due to predation. By virtue of this, in Ref. [12], the authors considered the following Leslie-Gower predator–prey model incorporating a prey refuge: q
This work is supported by the National Natural Sciences Foundation of People’s Republic of China under Grant 10971183.
⇑ Corresponding author.
E-mail address: [email protected] (Y. Li). 0096-3003/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2012.10.069
Y. Li, C. Li / Applied Mathematics and Computation 219 (2013) 4576–4589
dH ¼ ðr 1 b1 HÞH a1 ð1 mÞPH; dt
dP ¼ dt
r 2 a2
4577
P P: ð1 mÞH
They found that the prey refuge could influence the densities of both prey and predator species greatly. And they got the global stability and persistence of the system under proper conditions. Time delays are often incorporated into population models for resource regeneration times, maturing times or gestation periods [13,14]. It is generally recognized that some kinds of time delays are inevitable in population interactions and tend to be destabilizing in the sense that longer delays may destroy the stability of positive equilibria. Time delay due to the gestation is a common example, because generally the consumption of prey by the predator throughout its past history governs the present birth rate of the predator. Recently, great attention has been received and a large body of work has been carried out on the existence of Hopf bifurcations in delayed population models (see, for example, [15–17] and references cited therein). The stability of positive equilibria and the existence and the direction of Hopf bifurcations were discussed respectively in the references mentioned above. In pioneering work [18], Yuan and Song considered the following delayed Leslie-Gower predator–prey system:
8 sÞ > < x0 ðtÞ ¼ r1 xðtÞ 1 xðt mxðtÞyðtÞ; K > : y0 ðtÞ ¼ r 2 yðtÞ 1 cyðtÞ : xðtÞ
ð1:2Þ
They investigated the stability and Hopf bifurcation of the above system without considering the effects of time delays on predator. As is known to all, periodic solutions can arise through the Hopf bifurcation in delay differential equations. However, these bifurcating periodic solutions are generally local, i.e., they remain valid only in a small neighborhood of the critical value. Therefore, we mean to known if these non-constant periodic solutions obtained through local Hopf bifurcations exist globally. Motivated by the above discussion, in this paper, by choosing the time delay as a bifurcation parameter, we investigate a modified Leslie-Gower predator–prey system with time delays described by the following system:
8 sÞ > < x0 ðtÞ ¼ r1 xðtÞ 1 xðt að1 mÞxðtÞyðtÞ; K yðt s Þ > : y0 ðtÞ ¼ r 2 yðtÞ 1 cð1mÞxðt sÞ ;
ð1:3Þ
where xðtÞ represents the density of the prey at time t; yðtÞ represents the density of the predator at time t; the parameters a; K; m; r 1 ; r 2 ; c are positive constants in which r1 and r 2 are the intrinsic growth rates of the prey and the predator, respectively. The value K is the carrying capacity of the prey and cx takes on the role of a prey-dependent carrying capacity for the predator. The parameter c is a measure of the quality of the prey as food for the predator; mH is a refuge protecting of the prey, where m 2 ½0; 1Þ is a constant. This leaves ð1 mÞH of the prey available to the predator. The initial conditions for system (1.3) take the form
xðhÞ ¼ /ðhÞ;
yðhÞ ¼ wðhÞ;
/ðhÞ P 0;
wðhÞ P 0;
/ð0Þ > 0;
wð0Þ > 0;
h 2 ½s; 0;
ð1:4Þ
where ð/ðhÞ; wðhÞÞ 2 Cð½s; 0; R2þ0 Þ; R2þ0 ¼ fðx1 ; x2 Þ : xi P 0; i ¼ 1; 2g. It is well known by the fundamental theory of functional differential equations [19] that system (1.3) has a unique solution ðxðtÞ; yðtÞÞ satisfying initial conditions (1.4). It is easy to show that all solutions of system (1.3) corresponding to initial conditions (1.4) are defined on ½0; þ1Þ and remain positive for all t P 0. Remark 1.1. If m ¼ 0, then model (1.3) with s ¼ 0 in the second equation turns to be the model (1.2). The organization of this paper is as follows. In the next section, by choosing the time delay s as a parameter and analyzing the associated characteristic equation of a linearized system, we investigate the linear stability of the positive equilibrium for system (1.3). In addition, we get sufficient conditions for the existence of Hopf bifurcations. In Section 3, we derive formulae to determine the direction of bifurcations and the stability of bifurcating periodic solutions by using the normal form theory and center manifold theorem. In Section 4, numerical simulations are carried out to illustrate the theoretical results. In Section 5, based on the global Hopf bifurcation theorem for general functional differential equations, we investigate the global existence of periodic solutions by using degree theory methods. 2. Local stability and Hopf bifurcation In this section, we discuss the stability of a positive equilibrium and the existence of Hopf bifurcation for system (1.3) with time delay s as a parameter. 1 Obviously, system (1.3) always has a unique equilibrium E ðx ; y Þ with x ¼ r þaKKr ; y ¼ cð1 mÞx . cð1mÞ2 1
4578
Y. Li, C. Li / Applied Mathematics and Computation 219 (2013) 4576–4589
¼ y y . Dropping the bars, system (1.3) becomes Let x ¼ x x ; y
8 _ ¼ a1 xðt sÞ a2 yðtÞ a3 xðtÞxðt sÞ a4 xðtÞyðtÞ; < xðtÞ X ð1Þ 1 _ f xi ðt sÞyj ðt sÞyk ðtÞ; ¼ cð1 mÞr 2 xðt sÞ r 2 yðt sÞ þ i!j!k! ijk : yðtÞ
ð2:1Þ
iþjþkP2
where
r 1 x r1 ; a2 ¼ að1 mÞx ; a3 ¼ ; a4 ¼ að1 mÞ; K K @ iþjþk f ð1Þ y ð1Þ : ¼ i j k kðx ;y ;y Þ ; f ðx; y; y1 Þ ¼ r 2 y1 1 cð1 mÞx @x @y @y1
a1 ¼ ð1Þ
fijk
We then obtain the linearized system
_ xðtÞ ¼ a1 xðt sÞ a2 yðtÞ; _ yðtÞ ¼ cð1 mÞr 2 xðt sÞ r2 yðt sÞ:
From the above linearized system, we obtain the characteristic equation
det
k þ a1 eks
a2
cr2 ð1 mÞeks
k þ r2 eks
¼ 0;
i.e.,
k2 þ ½ðr2 þ a1 Þk þ ca2 r 2 ð1 mÞeks þ a1 r2 e2ks ¼ 0: Clearly, k ¼ 0 is not a root of (2.2). When
ð2:2Þ
s ¼ 0, (2.2) becomes
2
k þ ðr 2 þ a1 Þk þ ca2 r 2 ð1 mÞ þ a1 r2 ¼ 0:
ð2:3Þ
Noting that (2.3) must have negative real parts. Now we rewrite (2.2) as
k2 eks þ ½ðr2 þ a1 Þk þ ca2 r 2 ð1 mÞ þ a1 r 2 eks ¼ 0:
ð2:4Þ
Then to study the stability of the equilibrium ðx ; y Þ of (1.3) and Hopf bifurcation, it is equivalent to study the distribution of the roots of (2.4). When s > 0, noting that ixðx > 0Þ is a root of (2.4) if and only if x satisfies
(
ðx2 a1 r 2 Þ cos xs ¼ ca2 r 2 ð1 mÞ; ðx2 þ a1 r 2 Þ sin xs ¼ ðr 2 þ a1 Þx:
That is,
sin xs ¼
ðr 2 þ a1 Þx ca2 r2 ð1 mÞ ; cos xs ¼ : x2 þ a1 r2 x2 a1 r2
ð2:5Þ
2
Since sin xs þ cos2 xs ¼ 1, we obtain
x8 þ p3 x6 þ p2 x4 þ p1 x2 þ p0 ¼ 0;
ð2:6Þ
where
p0 ¼ a41 r42 c2 ð1 mÞ2 a21 a22 r 42 ;
p1 ¼ ½2a1 c2 ð1 mÞ2 a22 r 32 þ ðr 2 þ a1 Þ2 a21 r22 ;
p2 ¼ ½2a21 r 22 þ c2 ð1 mÞ2 a22 r 22 2a1 r 2 ðr 2 þ a1 Þ2 ;
p3 ¼ ðr 2 þ a1 Þ2 :
Denote u ¼ x2 , (1.1) turns to be
u4 þ p3 u3 þ p2 u2 þ p1 u þ p0 ¼ 0:
ð2:7Þ
Let
/ðuÞ ¼ u4 þ p3 u3 þ p2 u2 þ p1 u þ p0 : If all the coefficients of system (1.3) are given, we can easily work out the roots of (2.7) by computer. Since limu!þ1 /ðuÞ ¼ þ1, we claim that if p0 < 0, then (2.7) has at least one positive real root, i.e., ðH1 Þ If r1 < aK cð1 mÞ2 , (2.7) has at least one positive real root. Without loss of generality, assuming that it has four positive real roots, denoted by u1 ; u2 ; u3 ; u4 , respectively. Then (2.6) has four positive roots
Y. Li, C. Li / Applied Mathematics and Computation 219 (2013) 4576–4589
pffiffiffiffiffi
pffiffiffiffiffi
sin xk s ¼
ðr 2 þ a1 Þxk ; x2k þ a1 r2
pffiffiffiffiffi
4579
pffiffiffiffiffi
x1 ¼ u1 ; x2 ¼ u2 ; x3 ¼ u3 ; x4 ¼ u4 : So we have
k ¼ 1; 2; 3; 4:
Thus, denoting
sjk ¼
ðr 2 þ a1 Þxk þ 2j ; arcsin p xk x2k þ a1 r2 1
ð2:8Þ
s ¼ skj , the characteristic (2.2) has a pair of purely imaginary roots ixk .
where k ¼ 1; 2; 3; 4; j ¼ 0; 1; 2; . . .. Thus, when Define
s0 ¼ s0k0 ¼
min
k2fk¼1;2;3;4g
fs0k g;
x0 ¼ xk0 :
ð2:9Þ
Taking the derivative of k with respect to s in (2.4), it is easy to have
2keks
dk dk ks dk dk e þ ðr 2 þ a1 Þ ¼ 0: þ k2 k þ s a1 r 2 eks k þ s ds ds ds ds
Hence
dkðsÞ k3 eks þ a1 r 2 keks ¼ : k s ds 2ke þ k2 seks þ r2 þ a1 a1 r2 seks To simplify our expression, we denote x0 ; sk by x; s, respectively, then
dk ds
1 ¼
2keks þ k2 seks þ r 2 þ a1 a1 r 2 seks kðk2 eks
a1 r2
eks Þ
¼
2keks þ r2 þ a1 kðk2 eks
ks
¼
2ke þ r2 þ a1 ðr 2 þ a1 Þk2 þ cð1 mÞa2 r2 k þ 2a1 r 2 eks
s k
¼
a1 r 2
eks Þ
s k
A þ iB is þ ; C þ iD x
where
A ¼ r 2 þ a1 2x sin xs; B ¼ 2x cos xs; C ¼ 2a1 r2 cos xs x2 ðr 2 þ a1 Þ; D ¼ cð1 mÞa2 r2 x 2a1 r2 sin xs: Let Q ¼ C 2 þ D2 > 0, then
1 dk Q Re ¼ AC þ BD: ds Noting that
2 # 1 dk dk 4 ¼ sign Re sign Re ds s¼s0 ds
3
"
5:
s¼s0
To obtain the transversal condition, we also need the condition as follows. ðH2 Þ ReðddksÞjs¼s0 – 0, i.e., AC þ BD – 0, where A; B; C; D are defined above. Till now, we will apply the following lemma cited from [15,20] to analyze (2.2). Lemma 2.1 ([15,20]). Consider the exponential polynomial ð0Þ
ð0Þ
ð1Þ
ð1Þ
n1 ks1 Pðk; eks1 ; . . . ; eksm Þ ¼ kn þ p1 kn1 þ þ pn1 k þ pð0Þ þ þ pn1 k þ pð1Þ þ n þ ½p1 k n e ðmÞ
ðmÞ
þ ½p1 kn1 þ þ pn1 k þ pnðmÞ eksm ; ðiÞ
where s P 0 ði ¼ 1; 2; . . . ; mÞ and pj ði ¼ 0; 1; 2; . . . ; m; j ¼ 1; 2; . . . ; nÞ are constants. As ðs1 ; s2 ; . . . ; sm Þ vary, the sum of the order of the zeros of Pðk; eks1 ; . . . ; eksm Þ on the open right half plane can change only if a zero appears on or crosses the imaginary axis. From Lemma 2.1, it is easy to obtain the following theorem:
4580
Y. Li, C. Li / Applied Mathematics and Computation 219 (2013) 4576–4589
Theorem 2.1. Assume that ðH1 Þ ðH2 Þ hold, then the following results hold: (i) For system (1.3), its zero solution is asymptotically stable for s 2 ½0; s0 Þ. (ii) System (1.3) undergoes a Hopf bifurcation at the origin when s ¼ s0 . That is, system (1.3) has a branch of periodic solutions bifurcating from the zero solution near s ¼ s0 . 3. Direction and stability of Hopf bifurcations In this section, we study the direction of bifurcations and the stability of bifurcating periodic solutions of system (1.3) at
s0 by using the method based on the normal form theory and center manifold theory introduced by Hassard et al. in [21]. ^ðsÞ ¼ yðssÞ; s ¼ s0 þ l; l 2 R, and still denote by Now, we re-scale the time by t ¼ ss; ^ xðsÞ ¼ xðssÞ; y ^ðsÞ, then system (2.1) can be written as xðtÞ ¼ ^ xðsÞ; yðtÞ ¼ y
8_ xðtÞ ¼ ðs0 þ lÞ½a1 xðt 1Þ a2 yðtÞ a3 xðtÞxðt 1Þ a4 xðtÞyðtÞ; > < " # X ð1Þ i 1 j k _ > f x ðt 1Þy ðt 1Þy ðtÞ : : yðtÞ ¼ ðs0 þ lÞ cð1 mÞr 2 xðt 1Þ r2 yðt 1Þ þ i!j!k! ijk
ð3:1Þ
iþjþkP2
T
2
For u ¼ ðu1 ; u2 Þ 2 C½1; 0 ¼ Cð½1; 0; R Þ, define a family of operators
Ll u ¼ B1 uð0Þ þ B2 uð1Þ;
ð3:2Þ
where
B1 ¼ ðs0 þ lÞ
0 a2 0
;
0
B2 ¼ ðs0 þ lÞ
a1
0
cð1 mÞr2 r2
and define
! a3 u1 ð0Þu1 ð1Þ a4 u1 ð0Þu2 ð0Þ f ðl; uÞ ¼ ðs0 þ lÞ P : ð1Þ i j 1 k iþjþkP2 i!j!k! fijk u1 ð1Þu2 ð1Þu2 ð0Þ By the Riesz representation theorem, there exists a matrix whose components are bounded variation functions gðh; lÞ : ½1; 0 ! R2 such that
Ll u ¼
Z
0
dgðh; lÞuðhÞ:
1
In fact, choosing
gðh; lÞ ¼ B1 dðhÞ þ B2 dðh þ 1Þ;
ð3:3Þ
1; 0;
h¼0 where dðhÞ ¼ is a Dirac delta function, then (3.2) is satisfied. h–0 For u ¼ ðu1 ; u2 ÞT 2 C½1; 0 ¼ Cð½1; 0; R2 Þ, define
(
AðlÞu ¼
u_ ðhÞ; h 2 ½1; 0Þ 0; h 2 ½1; 0Þ; ; RðlÞu ¼ R0 l ; u Þ; h ¼ 0: f ð d g ðs; l Þ u ðsÞ; h ¼ 0 1
ð3:4Þ
Hence, (3.1) can be rewritten as
U_ t ¼ AðlÞU t þ RðlÞU t ;
ð3:5Þ
T
where U ¼ ðx; yÞ and U t ðhÞ ¼ Uðt þ hÞ; h 2 ð1; 0. For w 2 C 1 ½0; 1, define the adjoint operator A of A as
(
A wðsÞ ¼
_ wðsÞ; s 2 ð0; 1; R0 T dg ðt; 0ÞwðtÞ; s ¼ 0; 1
ð3:6Þ
where gT is the transpose of the matrix g. For u 2 Cð½1; 0; R2 Þ and w 2 Cð½0; 1; ðR2 Þ Þ, in order to normalize the eigenvectors of operator A and adjoint operator A , we define a bilinear inner product
hw; ui ¼ wT ð0Þuð0Þ
Z
0
h¼1
Z
n¼0
where gðhÞ ¼ gðh; 0Þ. Then as usual,
hw; Aui ¼ hA w; ui
h
wT ðn hÞ dgðhÞuðnÞ dn;
Y. Li, C. Li / Applied Mathematics and Computation 219 (2013) 4576–4589
4581
for ðu; wÞ 2 DðAÞ DðA Þ. Normalizing q and q by the condition
hq ; qi ¼ 1;
i ¼ 0: hq ; q
By the discussion in Section 2 and transformation t ¼ ss, we know that is0 x0 are eigenvalues of A and other eigenvalues have strictly negative real parts. Thus, they are also eigenvalues of A . Next we calculate the eigenvector q of A associated to the eigenvalue is0 x0 and the eigenvector q of A associated to the eigenvalue is0 x0 . Let
qðhÞ ¼
1 ix0 s0 h e ;
1 6 h < 0:
q
ð3:7Þ
From the above discussion, it is easy to know that
Aqð0Þ ¼ is0 x0 qð0Þ; i.e.,
s0
ix0 þ a1 eix0 s0 cr 2 ð1 mÞeix0 s0
! 1
a2 ix0 þ r2 eix0 s0
q
¼
0 : 0
ð3:8Þ
So we have
q¼
a1 eis0 x0 þ ix0 : a2
Similarly, assume that the eigenvector q of A is
q ðsÞ ¼
1 D
1
r
eix0 s0 s ;
0 6 s < 1:
ð3:9Þ
Then we obtain
A q ð0Þ ¼ is0 x0 q ð0Þ; i.e.,
s0
ix0 þ a1 eix0 s0
cr 2 ð1 mÞeix0 s0
a2
ix0 þ r 2 eix0 s0
!
1
r
¼
0 : 0
ð3:10Þ
Hence, we obtain
r¼
ix0 þ a1 eix0 s0 : cr 2 ð1 mÞeix0 s0
By hq ; qi ¼ 1, we can get q, that is,
ð0ÞÞT qð0Þ hq ; qi ¼ ðq ¼
1 D
Z
0 h¼1
Z
Z Þ 1 s0 ð 1r
q
h
ÞT ðn hÞ dgðhÞqðnÞ dn ðq
n¼0 0
Z
h¼1
h
ÞeiðnhÞx0 dgðhÞ ð 1r
n¼0
1 inx0 e dn
q
Z 0 1 1 s0 þ s0 ða1 r Þ ¼ 1: Þheix0 s0 h dgðhÞ 1 cð1 mÞr 2 r 2 qr ¼ ½1 þ qr 1 þ qr ¼ ð 1r q D D 1 Hence, we have
þ s0 ða1 r Þ; cð1 mÞr 2 r 2 qr D ¼ 1 þ qr such that
hq ðhÞ; qðhÞi ¼ 1;
ðhÞi ¼ 0: hq ðhÞ; q
So q and q are obtained. In the remainder of this section, we use the same notations as those in Ref. [21]. We first compute the coordinates to describe the center manifold C0 at l ¼ 0. Let U t be the solution of (3.5) when l ¼ 0, and define
zðtÞ ¼ hq ; U t i;
Wðt; hÞ ¼ U t ðhÞ 2RefzðtÞqðhÞg:
ð3:11Þ
On the center manifold C0 , we have Wðt; hÞ ¼ WðzðtÞ; zðtÞ; hÞ, where
Wðz; z; hÞ ¼ W 20 ðhÞ
z2 z2 þ W 11 ðhÞzz þ W 02 ðhÞ þ ; 2 2
ð3:12Þ
4582
Y. Li, C. Li / Applied Mathematics and Computation 219 (2013) 4576–4589
. Note that W is real if U t is real, we consider z and z are local coordinates for center manifold C0 in the direction of q and q only real solutions. For the solution U t 2 C0 , since l ¼ 0, then
z_ ðtÞ ¼ hq ; U_ t i ¼ hq ; AU t þ RU t i ¼ hq ; AU t i þ hq ; RU t i ¼ hA q ; U t i þ hq ; RU t i ð0Þf ð0; Wðt; 0Þ þ 2Re½zðtÞqð0ÞÞ: ¼ is0 x0 zðtÞ þ q
ð3:13Þ
We rewrite this equation as
z_ ðtÞ ¼ is0 x0 zðtÞ þ gðz; zÞ
ð3:14Þ
with
gðz; zÞ ¼ g 20
z2 z2 z2 z þ þ g 11 zz þ g 02 þ g 21 2 2 2
In the following, the values of
ð3:15Þ
l2 and b2 presented by Hassard at [21] are obtained. Substituting (3.5) and (3.13) into
_ ¼ U_ t z_ q z_ q ; W we have
_ ¼ U_ t z_ q z_ q ¼ AU t þ RU t ½is0 x0 z þ q ð0Þf ðz; zÞq ½is0 x0 z þ q ð0Þf ðz; zÞq W ð0Þf ðz; zÞqðhÞ 2Re½is0 x0 zqðhÞ ¼ AW 2Re½q ð0Þf ðz; zÞqðhÞ þ RU t ¼ AW þ 2AReðzqÞ þ RU t 2Re½q 8 ð0Þf ðz; zÞqðhÞ; 1 6 h < 0; < AW 2Re½q ð0Þf ðz; zÞqðhÞ þ f ðl; uÞ; h ¼ 0 ¼ AW 2Re½q : def ¼ AW þ Hðz; z; hÞ;
ð3:16Þ
where
Hðz; z; hÞ ¼ h20 ðhÞ
z2 z2 þ h11 ðhÞzz þ h02 ðhÞ þ 2 2
ð3:17Þ
Taking the derivative of W with respect to t in (3.12), we have
_ ¼ W z z_ þ W z z_ : W
ð3:18Þ
From (3.12), (3.14) and (3.18), we obtain
_ ¼ ðW 20 z þ W 11 z þ Þðis0 x0 z þ gÞ þ ðW 11 z þ W 02 z þ Þðis0 x0 z þ gÞ: W
ð3:19Þ
Then substituting (3.12) and (3.17) into (3.16), we have 2 2 _ ¼ ðAW 20 þ h20 Þ z þ ðAW 11 þ h11 Þzz þ ðAW 02 þ h02 Þ z þ W 2 2
ð3:20Þ
By comparing the coefficients of (3.19) with (3.20), the following hold:
ðA 2is0 x0 ÞW 20 ðhÞ ¼ h20 ðhÞ;
ð3:21Þ
AW 11 ðhÞ ¼ h11 ðhÞ:
ð3:22Þ
According to (3.13) and (3.14), we obtain
ð0Þf ðz; zÞ ¼ gðz; zÞ ¼ q
s0 D
Þ P ð1 r
a3 U 1t ð0ÞU 1t ð1Þa4 U 1t ð0ÞU 2t ð0Þ ð1Þ i j k 1 iþjþkP2 i!j!k! fijk U 1t ð1ÞU 2t ð1ÞU 2t ð0Þ
ðhÞ and qðhÞ ¼ ð 1 where U t ðhÞ ¼ ðU 1t ðhÞ; U 2t ðhÞ; U 3t ðhÞÞT ¼ Wðt; hÞ þ zqðhÞ þ zq
! ;
q ÞT eix0 h . Noting that
z2 z2 ð1Þ ð1Þ þ W 11 ð0Þzz þ W 02 ð0Þ þ ; 2 2 2 z2 z ð2Þ ð2Þ z þ W ð2Þ þ W 11 ð0Þzz þ W 02 ð0Þ þ ; U 2t ð0Þ ¼ qz þ q 20 ð0Þ 2 2 2 z2 z ð1Þ ð1Þ ð1Þ U 1t ð1Þ ¼ zeis0 x0 þ zeis0 x0 þ W 20 ð1Þ þ W 11 ð1Þzz þ W 02 ð1Þ þ ; 2 2 2 z2 z ð2Þ ð2Þ zeis0 x0 þ W ð2Þ þ W 11 ð1Þzz þ W 02 ð1Þ þ : U 2t ð1Þ ¼ qzeis0 x0 þ q 20 ð1Þ 2 2 ð1Þ
U 1t ð0Þ ¼ z þ z þ W 20 ð0Þ
ð3:23Þ
4583
Y. Li, C. Li / Applied Mathematics and Computation 219 (2013) 4576–4589
Therefore, we have
(" !# ð1Þ f 1 ð1Þ ð1Þ ð1Þ 2 ix0 s0 z2 a4 q a3 eix0 s0 þ r 200 e2ix0 s0 þ f110 qe2ix0 s0 þ f101 qeix0 s0 þ f011 qe 2 D h i ð1Þ ð1Þ ð1Þ ð1Þ f200 þ qÞ þ r Þ þ f101 eix0 s0 þ qeix0 s0 Þ þ f011 qeix0 s0 Þ zz eix0 s0 þ q þ a3 ðeix0 s0 þ eix0 s0 Þa4 ðq r þ f110 ðq þ q ðq ðqq " ! ! ð1Þ ð1Þ ð2Þ ð1Þ W 20 ð1Þ W 20 ð0Þ ix0 s0 W 20 ð0Þ W 20 ð0Þ ð1Þ ð1Þ ð2Þ ð1Þ ix0 s0 þ W 11 ð0Þe q þ qW 11 ð0Þ þ e þ þ a3 W 11 ð1Þ þ a4 W 11 ð0Þ þ 2 2 2 2
gðz; zÞ ¼
f200 ð1Þ ð1Þ 2W 11 ð1Þeix0 s0 þ W 20 ð1Þeix0 s0 2 ð1Þ
þr
ð2Þ
ð1Þ
ð1Þ
W 20 ð1Þ ix0 s0 W 20 ð1Þ ix0 s0 ð1Þ e þ q e þ W 11 ð1Þqeix0 s0 2 2 ! ð2Þ ð1Þ W 20 ð0Þ ix0 s0 W 20 ð1Þ ð2Þ ð1Þ ix0 s0 W 11 ð0Þe þ þ q þ W 11 ð1Þq e 2 2 ! ð2Þ ð2Þ W ð1Þ W ð0Þ ð2Þ ix0 s0 W 11 ð1Þq þ 20 q þ 20 q eix0 s0 þ W ð2Þ ð0Þ q e 11 2 2
!
ð2Þ
þ f110 W 11 ð1Þeix0 s0 þ ð1Þ
þ f101 ð1Þ
þ f011
!# ð1Þ ð1Þ ð1Þ f300 ix0 s0 f210 f ð1Þ þ 2qÞeix0 s0 þ 201 ðq e2ix0 s0 þ 2qÞ þ f111 þ qðq þ q þ q e2ix0 s0 Þ z2z e ðq 2 2 2 " !# ) ð1Þ f200 2ix0 s0 ð1Þ ð1Þ ð1Þ 2 ix0 s0 ix0 s0 2i x s i x s 0 0 0 0 z2 þ : þr e e e a4 q þ f110 q þ f101 q þ f011 q e þ a3 e 2 þ
Comparing the coefficients of (3.15) with the above equality, it follows that
" !# ð1Þ f200 2ix0 s0 2s0 ð1Þ ð1Þ ð1Þ 2 ix0 s0 ix0 s0 2i x s i x s 0 0 0 0 ; a4 q a3 e g 20 ¼ þr þ f110 qe þ f101 qe þ f011 q e e 2 D h i s0 ð1Þ ð1Þ ð1Þ ð1Þ f200 eix0 s0 þ q þ qÞ þ r Þ þ f101 eix0 s0 þ qeix0 s0 Þ þ f011 qeix0 s0 Þ ; g 11 ¼ r þ f110 ðq þ q ðq ðqq a3 ðeix0 s0 þ eix0 s0 Þ a4 ðq D " !# ð1Þ f200 2ix0 s0 2s0 ð1Þ ð1Þ ð1Þ 2 ix0 s0 ix0 s0 2i x s i x s 0 0 0 0 þr e e e ; a4 q þ f110 q þ f101 q þ f011 q a3 e e g 02 ¼ 2 D ! " ! ð1Þ ð1Þ ð2Þ ð1Þ W ð1Þ W 20 ð0Þ ix0 s0 W ð0Þ W 20 ð0Þ 2s0 ð1Þ ð1Þ ð2Þ g 21 ¼ a3 W 11 ð1Þ þ 20 þ W 11 ð0Þeix0 s0 a4 W 11 ð0Þ þ 20 q þ qW ð1Þ ð0Þ þ e þ 11 2 2 2 2 D f200 ð1Þ ð1Þ 2W 11 ð1Þeix0 s0 þ W 20 ð1Þeix0 s0 2 ð1Þ
þr ð1Þ
! ð2Þ ð1Þ W 20 ð1Þ ix0 s0 W 20 ð1Þ ix0 s0 ð1Þ þ q e þ W 11 ð1Þqeix0 s0 e 2 2 ! ð2Þ ð1Þ W 20 ð0Þ ix0 s0 W 20 ð1Þ ð2Þ ð1Þ ix0 s0 e W 11 ð0Þe þ þ q þ W 11 ð1Þq 2 2 ! ð2Þ ð2Þ ð1Þ ð1Þ W 20 ð1Þ W 20 ð0Þ ix0 s0 f f ð2Þ ð2Þ ix0 s0 þ 2qÞeix0 s0 þ 300 eix0 s0 þ 210 ðq W 11 ð1Þq þ qþ qe þ W 11 ð0Þqe 2 2 2 2 !# ð2Þ
þf110 W 11 ð1Þeix0 s0 þ ð1Þ
þf101 ð1Þ
þf011 ð1Þ
þ
f201 ð1Þ e2ix0 s0 þ 2qÞ þ f111 qðq þ q þ q e2ix0 s0 Þ ðq 2
:
Now we compute W 20 ðhÞ and W 11 ðhÞ. From (3.16), we have
ð0Þf ðz; zÞqðhÞÞ þ RU t ¼ gqðhÞ gq ðhÞ þ RU t Hðz; z; hÞ ¼ 2Reðq 1 1 1 1 ðhÞ þ RU t : g20 z2 þ g11 zz þ g02 z2 þ q ¼ g 20 z2 þ g 11 zz þ g 02 z2 þ qðhÞ 2 2 2 2
ð3:25Þ
Comparing the coefficients with (3.17), we can obtain, for h 2 ½1; 0Þ,
ðhÞ; h11 ðhÞ ¼ g 11 qðhÞ g11 q ðhÞ: h20 ðhÞ ¼ g 20 qðhÞ g02 q
ð3:26Þ
4584
Y. Li, C. Li / Applied Mathematics and Computation 219 (2013) 4576–4589
Substituting the above equalities into (3.21) and (3.22), respectively, it follows that
(
_ 20 ðhÞ ¼ 2is0 x0 W 20 ðhÞ þ g qðhÞ þ g02 q ðhÞ; W 20 _ 11 ðhÞ ¼ g qðhÞ þ g11 q ðhÞ: W
ð3:27Þ
11
Solving (3.27), we have
W 20 ðhÞ ¼
ig 20
s0 x0
qð0Þeis0 x0 h þ
ig02 ð0Þeis0 x0 h þ Ee2is0 x0 h q 3s 0 x0
ð3:28Þ
and
ig 11
W 11 ðhÞ ¼
s0 x0
qð0Þeis0 x0 h þ
g11 ð0Þeis0 x0 h þ F: q is0 x0
ð3:29Þ
In what follows, we seek appropriate E and F. From (3.21) and (3.22), we can get
AW 20 ð0Þ ¼ 2is0 x0 W 20 ð0Þ h20 ð0Þ
ð3:30Þ
AW 11 ð0Þ ¼ h11 ð0Þ:
ð3:31Þ
and
From the definition of A in (3.4), we obtain
Z
0
dgðhÞW 20 ðhÞ ¼ 2is0 x0 W 20 ð0Þ h20 ð0Þ
ð3:32Þ
dgðhÞW 11 ðhÞ ¼ h11 ð0Þ:
ð3:33Þ
1
and
Z
0
1
From (3.26) and the definition of f ðl; sÞ, we have
ð0Þ þ s0 h20 ð0Þ ¼ g 20 qð0Þ g02 q
d1
ð3:34Þ
d2
and
ð0Þ þ s0 h11 ð0Þ ¼ g 11 qð0Þ g11 q
d3 ; d4
ð3:35Þ
where
d1 ¼ 2a4 q 2a3 eix0 s0 ; ð1Þ
ð1Þ
ð1Þ
ð1Þ
d2 ¼ f200 e2ix0 s0 þ 2f110 qe2ix0 s0 þ 2f101 qeix0 s0 þ 2f011 q2 eix0 s0 ; þ qÞ; d3 ¼ a3 ðeix0 s0 þ eix0 s0 Þ a4 ðq ð1Þ
ð1Þ
ð1Þ
ð1Þ
eix0 s0 þ q Þ þ f101 ðq eix0 s0 þ qeix0 s0 Þ þ f011 ðqq qeix0 s0 Þ: d4 ¼ f200 þ f110 ðq þ q Substituting (3.28) and (3.34) into (3.32) and noticing that
Z is0 x0 I
0
eis0 x0 h dgðhÞ qð0Þ ¼ 0
1
and
Z is0 x0 I
ð0Þ ¼ 0; eis0 x0 h dgðhÞ q
0
1
we obtain
d1 ; e2is0 x0 h dgðhÞ E ¼ s0 d2 1
Z 2is0 x0 I
0
ð3:36Þ
which leads to
2ix0 þ a1 e2ix0 s0 2ix0 s0
cr 2 ð1 mÞe It follows that
!
a2 2ix0 s0
2ix0 þ r2 e
E¼
d1 d2
:
4585
Y. Li, C. Li / Applied Mathematics and Computation 219 (2013) 4576–4589
E¼
!1
2ix0 þ a1 e2ix0 s0
a2
cr2 ð1 mÞe2ix0 s0
2ix0 þ r 2 e2ix0 s0
d1
ð3:37Þ
:
d2
Similarly, substituting (3.29) and (3.35) into (3.31), we have
F¼
a1
a2
cr2 ð1 mÞ
r2
1
d3 d4
:
ð3:38Þ
As a result, we know that W 20 ðhÞ and W 11 ðhÞ, then g ij in (3.24) is determined by the parameters and the time delay in (2.2). Thus, we can compute the following quantities
C1 ð0Þ ¼
1 g g 20 g 11 2jg 11 j2 jg 02 j2 þ 21 ; 2s 0 x0 3 2 i
l2 ¼
ReC1 ð0Þ ; Rek0 ðs0 Þ
t2 ¼
ImC1 ð0Þ þ l2 Imk0 ðs0 Þ
s0 x0
b2 ¼ 2ReC1 ð0Þ:
;
ð3:39Þ From the expression of C1 ð0Þ in (3.39), it is easy to get the values of l2 ; b2 and t 2 . Thus, if system (1.3) is given, then we will analyze the direction of the Hopf bifurcation and stability of the bifurcating periodic orbits at s ¼ s0 according to the following results: (i) (ii)
l ¼ 0 is Hopf bifurcation value of system (3.1); l2 determines the direction of the Hopf bifurcation: if l2 > 0ð< 0Þ, then the Hopf bifurcation is supercritical ðsubcriticalÞ and the bifurcating periodic solutions exist for s > s0 ð< s0 Þ;
(iii) b2 determines the stability of the bifurcating periodic solutions: if b2 < 0ð> 0Þ the bifurcating periodic solutions are stable (unstable); (iv) t2 determines the period of the bifurcating periodic solutions: the period increase (decrease) if t 2 > 0ð< 0Þ.
1.02
1 0.98 0.96
0.96
0.94
0.94
0.92
y
0.98
0.92
0.9
0.9
0.88
0.88
0.86
0.86
0.84
0.84
0
20
40
60
80
0.82
100
t
0
20
40
t
60
80
1 0.98 0.96 0.94
y
x
1
0.92 0.9 0.88 0.86 0.84 0.82 0.84 0.86 0.88 0.9
Fig. 1. Behavior and phase portrait of system (1.3) with
0.92 0.94 0.96 0.98
x
1
1.02
s ¼ 0:5 < s0 , the positive equilibrium E ðx ; y Þ ¼ ð0:8696; 0:8696Þ is stable.
100
4586
Y. Li, C. Li / Applied Mathematics and Computation 219 (2013) 4576–4589
4. Numerical example To illustrate the theoretical results, let us give some numerical simulations. Example 4.1. Consider the following system.
(
sÞ x0 ðtÞ ¼ xðtÞð1 xðt Þ 1:3ð1 0:5ÞxðtÞyðtÞ; 2
ð4:1Þ
yðtsÞ y0 ðtÞ ¼ 1:5yðtÞð1 2ð10:5Þxðt sÞÞ;
where r 1 ¼ 1; r 2 ¼ 1:5; K ¼ 2; a ¼ 1:3; c ¼ 2; m ¼ 0:5. The positive equilibrium is E ¼ ð0:8696; 0:8696Þ. (2.6) has only one positive root. Furthermore, we can obtain that s0 ¼ 0:6729. Simple calculations show that conditions of Theorem 2.1 hold. We can see from Figs. 1 and 2 that E is asymptotically stable when s ¼ 0:5 < s0 while it loses its stability when s ¼ 0:7 > s0 . Remark 4.1. In Ref. [12], the authors considered system (1.3) with s ¼ 0, they obtained the system is globally stable. However, system (1.3) is not always stable, which is illustrated in the above example. That is, time delays even if very little ones could disturb the stability of dynamical systems.
5. Global continuation of local Hopf bifurcation In this section, we study the existence of global Hopf bifurcations. The method we used here is based on the global Hopf bifurcating theorem for general functional differential equations introduced by Wu in [22]. For convenience, we write system (1.3) as the following form:
z_ ðtÞ ¼ Fðzt ; s; pÞ;
ð5:1Þ
where z ¼ ðx; y1 ; y2 ÞT ; zt ðhÞ 2 Cð½s; 0; R3 Þ. Following the work of Wu [19], we define
1
1.05 1
0.95
0.95 0.9
y
x
0.9 0.85
0.85
0.8 0.8
0.75
0
0.75 20
40
60
t
80
0.7
100
0
20
40
t
60
80
100
1.1
y
1 0.9 0.8 0.7 1 0.9 0.8
x Fig. 2. Behavior and phase portrait of system (1.3) with occurs.
0.7
0
20
40
60
80
100
t
s ¼ 0:7 > s0 , the positive equilibrium E ð0:8696; 0:8696Þ loses its stability and Hopf bifurcation
Y. Li, C. Li / Applied Mathematics and Computation 219 (2013) 4576–4589
4587
X ¼ Cð½s; 0; R3 Þ;
R ¼ ClfðzðtÞ; s; pÞ 2 X R Rþ ; zðtÞis a T-periodic solution of Ref: ð5:1Þg; N ¼ fðz; s; pÞ; Fðz; s; pÞ ¼ 0g: 2p Let lðz ;sj ; 2p Þ be the connected component of ðz ; sj ; x Þ in R, where 0 k
x0
sj and x0 are defined in (2.9).
Lemma 5.1. If s is bounded, all solutions of system (1.3) are uniformly bounded. Proof. Let ðxðtÞ; yðtÞÞT is a solution of system (1.3) through ð/; wÞ at t ¼ 0 with /ð0Þ > 0; wð0Þ > 0. Then we know that ðxðtÞ; yðtÞÞT is positive on the maximum existence interval of solution. h If ðxðtÞ; yðtÞÞT is a solution of (1.3) with xðtÞ > 0; yðtÞ > 0, then it follows from the first equation of (1.3) that
xðt sÞ : x0 ðtÞ < r 1 xðtÞ 1 K
ð5:2Þ
Clearly, xðtÞ < xðt sÞer1 s for t > s, which implies that xðt sÞ > xðtÞer1 s . This, together with (5.2), leads to
xðtÞer1 s ; x0 ðtÞ < r 1 xðtÞ 1 K
for; t > s;
which implies that
lim sup xðtÞ 6 Ker1 s : t!þ1
Thus, for any
e > 0, there exists a T > 0ðT > sÞ such that when t > T, we have
xðtÞ 6 Ker1 s þ e :¼ d: It follows from the second equation of (1.3) that, for t > T þ s,
yðt sÞ y0 ðtÞ < r 2 yðtÞ 1 : cð1 mÞd
ð5:3Þ
Then by similar discussion as above, we can easily have
lim sup yðtÞ 6 cð1 mÞder2 s : t!þ1
Thus the possible positive solutions of system (1.3) must be uniformly bounded. This completes the proof of the lemma. Lemma 5.2. System (1.3) has no nontrivial periodic solutions of period s.
Proof. Assume system (1.3) has a nontrivial periodic solution of period s, then the following differential system, which has the same equilibrium as (1.3),
(
x0 ðtÞ ¼ r 1 xðtÞð1 xðtÞ Þ að1 mÞxðtÞyðtÞ; K
ð5:4Þ
yðtÞ y0 ðtÞ ¼ r 2 yðtÞð1 cð1mÞxðtÞ Þ;
has periodic solutions. Due to Lemma 5.1, we restrict our attention to 0 < x < K; 0 < yðtÞ < cKð1 mÞ, respectively. System (5.4) also has the equilibrium E ¼ ðx ; y Þ. Define
Vðx; yÞ ¼ ln
x x ax cð1 mÞ2 y y þ ln þ þ : x x y y r2
Obviously, Vðx; yÞ is well defined and continuous for all xðtÞ > 0; yðtÞ > 0. Calculating the derivative of V along the solution of system (1.3), we have
dV r1 að1 mÞ ¼ ðx x Þ2 ðy y Þ2 : dt y Kx Obviously, dV < 0 strictly for all x; y > 0 except the positive equilibrium ðx ; y Þ where dt Thus, Vðx; yÞ satisfies Lyapunov’s asymptotic stability theorem, we conclude that
dV dt
¼ 0.
limðxðtÞ; yðtÞÞ ¼ ðx ; y Þ;
t!1
which contradicts the fact that system (5.4) has periodic solutions. This ends the proof. h
4588
Y. Li, C. Li / Applied Mathematics and Computation 219 (2013) 4576–4589
Theorem 5.1. Suppose that ðH1 Þ ðH2 Þ hold, then for each s > sjk ðj ¼ 1; 2; . . .Þ, system (1.3) has at least j þ 1 periodic solutions, where sjk is defined in (2.8). Proof. The characteristic matrix of system (1.3) at the positive equilibrium z is of the form
Dðz ; s; pÞðkÞ ¼
k þ a1 eks
a2
cr 2 ð1 mÞeks
k þ r2 e
: ks
ð5:5Þ
2p From the discussion of Section 2, it can be verified that ðz ; sj ; x Þ; j ¼ 1; 2; . . . are isolated centers. 0 Let
2p X;x2p ¼ ðg; pÞ : 0 < g < ;
p
< :
x0
0
Clearly, if js sj j 6 d and ðg; pÞ 2 @ X , then the necessary and sufficient conditions for Dðz ; s; pÞðg þ i 2ppÞ ¼ 0 are
g ¼ 0; s ¼ sj and p ¼ x2p0 . 2p Þðg; pÞ ¼ Dðz ; sj d; pÞðg þ i 2ppÞ, then we have the transversal number Define H ðz ; sj ; x 0
c z ; sj ;
2p
x0
2p 2p ¼ degB H z ; sj ; ; X;x2p degB Hþ z ; sj ; ; X;x2p ¼ 1:
x0
x0
0
0
2p By Theorem 2.1 of Wu [22], we conclude that the connected component lðz ;sj ;x2p Þ through ðz ; sj ; x Þ in R is nonempty. Mean0 0 while, we have
X ðz;s;pÞ2l
cðz; s; pÞ < 0
2p Þ ðz ;sj ;x 0
and hence lðz ;sj ;x2p Þ is unbounded. 0
From (2.8), we see that, for j P 1,
2p < sjk : x0 ; 1Þ, where s < sjk . Clearly, it follows from Now, we are in a position to prove that the projection of lðz ;sj ; 2p Þ onto s-space is ½s k
x0
the proof of Lemma 5.2 that system (1.3) with s ¼ 0 has no nontrivial periodic solution. Hence, the projection of lðz ;sj ; 2p Þ onto k
x0
s-space is away from zero. For a contradiction, we suppose that the projection of lðz ;sj ; 2p Þ onto s-space is bounded. This means that the projection of k
x0
2p lðz ;sj ; 2p Þ onto s-space is included in an interval ð0; s Þ. Noting that x < sjk and applying Lemma 5.2, we have p < s for ðz; s; pÞ 0 k
x0
belonging to lðz ;sj ; 2p Þ . This implies that the projection of the connected component lðz ;sj ; 2p Þ onto p-space is bounded. k
x0
k
x0
Moreover, from Lemma 5.1, we obtain that the projection of lðz ;sj ; 2p Þ onto z-space is bounded if the projection of lðz ;sj ; 2p Þ onto k
x0
k
x0
s-space is bounded. Thus, the connected component lðz ;sj ; 2p Þ crossing through ðE ; sjk ; x2p0 ÞÞ is bounded, which is a k
x0
; 1Þ for each j P 1, where s
of the proof.
x0
h
References [1] C. Ji, D. Jiang, N. Shi, Analysis of a predator–prey model with modified Leslie-Gower and Holling-type II schemes with stochastic perturbation, J. Math. Anal. Appl. 359 (2009) 482–498. [2] Y. Song, S. Yuan, J. Zhang, Bifurcation analysis in the delayed Leslie-Gower predator–prey system, Appl. Math. Model. 33 (2009) 4049–4061. [3] P.H. Leslie, Some further notes on the use of matrices in population mathematics, Biometrika 35 (1948) 213–245. [4] Y. Chen, Z. Liu, M. Haque, Analysis of a Leslie-Gower-type prey-predator model with periodic impulsive perturbations, Commun. Nonlinear Sci. Numer. Simulat. 14 (2009) 3412–3423. [5] P. Aguirre et al, Two limit cycles in a Leslie-Gower predator–prey model with additive Allee effect, Nonlinear Anal. Real World Appl. 10 (2009) 1401– 1416. [6] A. Korobeinikov, A Lyapunov function for Leslie-Gower predator–prey models, Appl. Math. Lett. 14 (2001) 697–699. [7] A.F. Nindjin, M.A. Aziz-Alaoui, M. Cadivel, Analysis of a predator–prey model with modified Leslie-Gower and Holling-type II schemes with time delay, Nonlinear Anal. Real World Appl. 7 (5) (2007) 1104–1118. [8] Z. Jing, J. Yang, Bifurcation and chaos in discrete-time predator–prey system, Chaos Solitons Fract. 27 (2006) 259–277. [9] W. Wang, L. Chen, A predator–prey system with stage structure for predator, Comput. Math. Appl. 33 (1997) 83–91. [10] Y. Yang, Hopf bifurcation in a two-competitor, one-prey system with time delay, Appl. Math. Comput. 214 (2009) 228–235. [11] T. Kumar Kar, Stability analysis of a prey-predator model incorporating a prey refuge, Commun. Nonlinear Sci. Numer. Simulat. 10 (2005) 681–691. [12] F. Chen, L. Chen, X. Xie, On a Leslie-Gower predator–prey model incorporating a prey refuge, Nonlinear Anal. Real World Appl. 10 (2009) 2905–2908. [13] J.M. Cushing, Integrodifferential Equations and Delay Models in Population Dynamics, Springer Verlag, Heidelberg, 1977.
Y. Li, C. Li / Applied Mathematics and Computation 219 (2013) 4576–4589
4589
[14] S.A. Gourley, M.V. Bartuccelli, Parameter domains for instability of uniform states in system with many delays, J. Math. Biol. 35 (1997) 843–867. [15] W. Yu, J. Cao, Stability and Hopf bifurcation analysis on a four-neuron BAM neural networks with time delays, Phy. Lett. A 351 (2006) 64–78. [16] S. Toaha, M.A. Hassan, Stability analysis of predator–prey population model with time delay and constant rate of harvesting, Punjab Univ. J. Math. 40 (2008) 37–48. [17] R. Xu, Z. Ma, Stability and Hopf bifurcation in a predator–prey model with stage structure for the predator, Nonlinear Anal. Real World Appl. 9 (2008) 1444–1460. [18] S. Yuan, Y. Song, Stability and Hopf bifurcations in a delayed Leslie-Gower predator–prey system, J. Math. Anal. Appl. 355 (2009) 82–100. [19] J. Hale, Theory of Functional Differential Equations, Springer-Verlag, Heidelberg, 1977. [20] S. Ruan, J. Wei, On the zeros of transcendental functions with applications to stability of delay differential equations with two delays, Dynam. Contin. Discrete Impuls. Syst. Ser. A Math. Anal. 10 (2003) 863–874. [21] B. Hassard, N. Kazarinoff, Y.H. Wan, Theory and Applications of Hopf Bifurcation, London Mathematical Society Lecture Note Series, vol. 41, Cambridge University Press, Cambridge, 1981. [22] J. Wu, Symmetric functional differential equations and neural networks with memory, Trans. Amer. Math. Soc. 350 (1998) 4799–4838.
Contents lists available at SciVerse ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Stability and Hopf bifurcation analysis on a delayed Leslie-Gower predator–prey system incorporating a prey refuge q Yongkun Li a,⇑, Changzhao Li b a b
Department of Mathematics, Yunnan University, Kunming, Yunnan 650091, People’s Republic of China Department of Applied Mathematics, Kunming University of Science and Technology, Kunming 650093, Yunnan, People’s Republic of China
a r t i c l e
i n f o
a b s t r a c t In this paper, a modified Leslie-Gower predator–prey system with time delays is investigated, where the time delays are regarded as bifurcation parameters. By analyzing the corresponding characteristic equation, the local stability of a positive equilibrium is considered. Moreover, we show that Hopf bifurcations occur when time delay crosses some critical values. By deriving the equation describing the flow on the center manifold, we give the formula for determining the direction of the Hopf bifurcation and the stability of bifurcating periodic solutions. In addition, we also try on the global existence of periodic solutions by using the global Hopf bifurcation result of Wu [J. Wu, Symmetric functional differential equations and neural networks with memory, Trans. Amer. Math. Soc. 350 (1998) 4799–4838.] for functional differential equations. Numerical simulations are carried out to illustrate the theoretical results and they show that the time delays in the system under consideration can destroy the stability of the system. Ó 2012 Elsevier Inc. All rights reserved.
Keywords: Predator–prey system Time delay Hopf bifurcation Leslie-Gower Stability
1. Introduction The dynamic relationship between predators and their preys has long been and will continue to be one of the dominant themes in both ecology and mathematical ecology due to its universal existence and importance [1,2]. Considering the reduction in a predator population has a reciprocal relationship with per capita availability of its preferred food, Leslie [3] introduced the following predator–prey model where the ’’carrying capacity’’ of the predator’s environment is proportional to the number of prey:
dH ¼ ðr 1 a1 P qðHÞÞH; dt
dP ¼ dt
r 2 a2
P P; H
ð1:1Þ
where H is the density of prey species at time t, and P is the density of the predator species at time t, qðHÞ is the so-called predator functional response to prey and satisfies the usual properties that make system (1.1) a predator–prey one. In the last decade, based on system (1.1) with different predator functional response to prey, many excellent works have been done, for example, see [1,2,4–7] and the references cited therein. For more results on predator–prey systems, one could refer to [8–10] for details. As was pointed out by Tapan Kumar Kar [11], mite predator–prey interactions often exhibit spatial refugia which afford the prey some degree of protection from predation and reduce the chance of extinction due to predation. By virtue of this, in Ref. [12], the authors considered the following Leslie-Gower predator–prey model incorporating a prey refuge: q
This work is supported by the National Natural Sciences Foundation of People’s Republic of China under Grant 10971183.
⇑ Corresponding author.
E-mail address: [email protected] (Y. Li). 0096-3003/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2012.10.069
Y. Li, C. Li / Applied Mathematics and Computation 219 (2013) 4576–4589
dH ¼ ðr 1 b1 HÞH a1 ð1 mÞPH; dt
dP ¼ dt
r 2 a2
4577
P P: ð1 mÞH
They found that the prey refuge could influence the densities of both prey and predator species greatly. And they got the global stability and persistence of the system under proper conditions. Time delays are often incorporated into population models for resource regeneration times, maturing times or gestation periods [13,14]. It is generally recognized that some kinds of time delays are inevitable in population interactions and tend to be destabilizing in the sense that longer delays may destroy the stability of positive equilibria. Time delay due to the gestation is a common example, because generally the consumption of prey by the predator throughout its past history governs the present birth rate of the predator. Recently, great attention has been received and a large body of work has been carried out on the existence of Hopf bifurcations in delayed population models (see, for example, [15–17] and references cited therein). The stability of positive equilibria and the existence and the direction of Hopf bifurcations were discussed respectively in the references mentioned above. In pioneering work [18], Yuan and Song considered the following delayed Leslie-Gower predator–prey system:
8 sÞ > < x0 ðtÞ ¼ r1 xðtÞ 1 xðt mxðtÞyðtÞ; K > : y0 ðtÞ ¼ r 2 yðtÞ 1 cyðtÞ : xðtÞ
ð1:2Þ
They investigated the stability and Hopf bifurcation of the above system without considering the effects of time delays on predator. As is known to all, periodic solutions can arise through the Hopf bifurcation in delay differential equations. However, these bifurcating periodic solutions are generally local, i.e., they remain valid only in a small neighborhood of the critical value. Therefore, we mean to known if these non-constant periodic solutions obtained through local Hopf bifurcations exist globally. Motivated by the above discussion, in this paper, by choosing the time delay as a bifurcation parameter, we investigate a modified Leslie-Gower predator–prey system with time delays described by the following system:
8 sÞ > < x0 ðtÞ ¼ r1 xðtÞ 1 xðt að1 mÞxðtÞyðtÞ; K yðt s Þ > : y0 ðtÞ ¼ r 2 yðtÞ 1 cð1mÞxðt sÞ ;
ð1:3Þ
where xðtÞ represents the density of the prey at time t; yðtÞ represents the density of the predator at time t; the parameters a; K; m; r 1 ; r 2 ; c are positive constants in which r1 and r 2 are the intrinsic growth rates of the prey and the predator, respectively. The value K is the carrying capacity of the prey and cx takes on the role of a prey-dependent carrying capacity for the predator. The parameter c is a measure of the quality of the prey as food for the predator; mH is a refuge protecting of the prey, where m 2 ½0; 1Þ is a constant. This leaves ð1 mÞH of the prey available to the predator. The initial conditions for system (1.3) take the form
xðhÞ ¼ /ðhÞ;
yðhÞ ¼ wðhÞ;
/ðhÞ P 0;
wðhÞ P 0;
/ð0Þ > 0;
wð0Þ > 0;
h 2 ½s; 0;
ð1:4Þ
where ð/ðhÞ; wðhÞÞ 2 Cð½s; 0; R2þ0 Þ; R2þ0 ¼ fðx1 ; x2 Þ : xi P 0; i ¼ 1; 2g. It is well known by the fundamental theory of functional differential equations [19] that system (1.3) has a unique solution ðxðtÞ; yðtÞÞ satisfying initial conditions (1.4). It is easy to show that all solutions of system (1.3) corresponding to initial conditions (1.4) are defined on ½0; þ1Þ and remain positive for all t P 0. Remark 1.1. If m ¼ 0, then model (1.3) with s ¼ 0 in the second equation turns to be the model (1.2). The organization of this paper is as follows. In the next section, by choosing the time delay s as a parameter and analyzing the associated characteristic equation of a linearized system, we investigate the linear stability of the positive equilibrium for system (1.3). In addition, we get sufficient conditions for the existence of Hopf bifurcations. In Section 3, we derive formulae to determine the direction of bifurcations and the stability of bifurcating periodic solutions by using the normal form theory and center manifold theorem. In Section 4, numerical simulations are carried out to illustrate the theoretical results. In Section 5, based on the global Hopf bifurcation theorem for general functional differential equations, we investigate the global existence of periodic solutions by using degree theory methods. 2. Local stability and Hopf bifurcation In this section, we discuss the stability of a positive equilibrium and the existence of Hopf bifurcation for system (1.3) with time delay s as a parameter. 1 Obviously, system (1.3) always has a unique equilibrium E ðx ; y Þ with x ¼ r þaKKr ; y ¼ cð1 mÞx . cð1mÞ2 1
4578
Y. Li, C. Li / Applied Mathematics and Computation 219 (2013) 4576–4589
¼ y y . Dropping the bars, system (1.3) becomes Let x ¼ x x ; y
8 _ ¼ a1 xðt sÞ a2 yðtÞ a3 xðtÞxðt sÞ a4 xðtÞyðtÞ; < xðtÞ X ð1Þ 1 _ f xi ðt sÞyj ðt sÞyk ðtÞ; ¼ cð1 mÞr 2 xðt sÞ r 2 yðt sÞ þ i!j!k! ijk : yðtÞ
ð2:1Þ
iþjþkP2
where
r 1 x r1 ; a2 ¼ að1 mÞx ; a3 ¼ ; a4 ¼ að1 mÞ; K K @ iþjþk f ð1Þ y ð1Þ : ¼ i j k kðx ;y ;y Þ ; f ðx; y; y1 Þ ¼ r 2 y1 1 cð1 mÞx @x @y @y1
a1 ¼ ð1Þ
fijk
We then obtain the linearized system
_ xðtÞ ¼ a1 xðt sÞ a2 yðtÞ; _ yðtÞ ¼ cð1 mÞr 2 xðt sÞ r2 yðt sÞ:
From the above linearized system, we obtain the characteristic equation
det
k þ a1 eks
a2
cr2 ð1 mÞeks
k þ r2 eks
¼ 0;
i.e.,
k2 þ ½ðr2 þ a1 Þk þ ca2 r 2 ð1 mÞeks þ a1 r2 e2ks ¼ 0: Clearly, k ¼ 0 is not a root of (2.2). When
ð2:2Þ
s ¼ 0, (2.2) becomes
2
k þ ðr 2 þ a1 Þk þ ca2 r 2 ð1 mÞ þ a1 r2 ¼ 0:
ð2:3Þ
Noting that (2.3) must have negative real parts. Now we rewrite (2.2) as
k2 eks þ ½ðr2 þ a1 Þk þ ca2 r 2 ð1 mÞ þ a1 r 2 eks ¼ 0:
ð2:4Þ
Then to study the stability of the equilibrium ðx ; y Þ of (1.3) and Hopf bifurcation, it is equivalent to study the distribution of the roots of (2.4). When s > 0, noting that ixðx > 0Þ is a root of (2.4) if and only if x satisfies
(
ðx2 a1 r 2 Þ cos xs ¼ ca2 r 2 ð1 mÞ; ðx2 þ a1 r 2 Þ sin xs ¼ ðr 2 þ a1 Þx:
That is,
sin xs ¼
ðr 2 þ a1 Þx ca2 r2 ð1 mÞ ; cos xs ¼ : x2 þ a1 r2 x2 a1 r2
ð2:5Þ
2
Since sin xs þ cos2 xs ¼ 1, we obtain
x8 þ p3 x6 þ p2 x4 þ p1 x2 þ p0 ¼ 0;
ð2:6Þ
where
p0 ¼ a41 r42 c2 ð1 mÞ2 a21 a22 r 42 ;
p1 ¼ ½2a1 c2 ð1 mÞ2 a22 r 32 þ ðr 2 þ a1 Þ2 a21 r22 ;
p2 ¼ ½2a21 r 22 þ c2 ð1 mÞ2 a22 r 22 2a1 r 2 ðr 2 þ a1 Þ2 ;
p3 ¼ ðr 2 þ a1 Þ2 :
Denote u ¼ x2 , (1.1) turns to be
u4 þ p3 u3 þ p2 u2 þ p1 u þ p0 ¼ 0:
ð2:7Þ
Let
/ðuÞ ¼ u4 þ p3 u3 þ p2 u2 þ p1 u þ p0 : If all the coefficients of system (1.3) are given, we can easily work out the roots of (2.7) by computer. Since limu!þ1 /ðuÞ ¼ þ1, we claim that if p0 < 0, then (2.7) has at least one positive real root, i.e., ðH1 Þ If r1 < aK cð1 mÞ2 , (2.7) has at least one positive real root. Without loss of generality, assuming that it has four positive real roots, denoted by u1 ; u2 ; u3 ; u4 , respectively. Then (2.6) has four positive roots
Y. Li, C. Li / Applied Mathematics and Computation 219 (2013) 4576–4589
pffiffiffiffiffi
pffiffiffiffiffi
sin xk s ¼
ðr 2 þ a1 Þxk ; x2k þ a1 r2
pffiffiffiffiffi
4579
pffiffiffiffiffi
x1 ¼ u1 ; x2 ¼ u2 ; x3 ¼ u3 ; x4 ¼ u4 : So we have
k ¼ 1; 2; 3; 4:
Thus, denoting
sjk ¼
ðr 2 þ a1 Þxk þ 2j ; arcsin p xk x2k þ a1 r2 1
ð2:8Þ
s ¼ skj , the characteristic (2.2) has a pair of purely imaginary roots ixk .
where k ¼ 1; 2; 3; 4; j ¼ 0; 1; 2; . . .. Thus, when Define
s0 ¼ s0k0 ¼
min
k2fk¼1;2;3;4g
fs0k g;
x0 ¼ xk0 :
ð2:9Þ
Taking the derivative of k with respect to s in (2.4), it is easy to have
2keks
dk dk ks dk dk e þ ðr 2 þ a1 Þ ¼ 0: þ k2 k þ s a1 r 2 eks k þ s ds ds ds ds
Hence
dkðsÞ k3 eks þ a1 r 2 keks ¼ : k s ds 2ke þ k2 seks þ r2 þ a1 a1 r2 seks To simplify our expression, we denote x0 ; sk by x; s, respectively, then
dk ds
1 ¼
2keks þ k2 seks þ r 2 þ a1 a1 r 2 seks kðk2 eks
a1 r2
eks Þ
¼
2keks þ r2 þ a1 kðk2 eks
ks
¼
2ke þ r2 þ a1 ðr 2 þ a1 Þk2 þ cð1 mÞa2 r2 k þ 2a1 r 2 eks
s k
¼
a1 r 2
eks Þ
s k
A þ iB is þ ; C þ iD x
where
A ¼ r 2 þ a1 2x sin xs; B ¼ 2x cos xs; C ¼ 2a1 r2 cos xs x2 ðr 2 þ a1 Þ; D ¼ cð1 mÞa2 r2 x 2a1 r2 sin xs: Let Q ¼ C 2 þ D2 > 0, then
1 dk Q Re ¼ AC þ BD: ds Noting that
2 # 1 dk dk 4 ¼ sign Re sign Re ds s¼s0 ds
3
"
5:
s¼s0
To obtain the transversal condition, we also need the condition as follows. ðH2 Þ ReðddksÞjs¼s0 – 0, i.e., AC þ BD – 0, where A; B; C; D are defined above. Till now, we will apply the following lemma cited from [15,20] to analyze (2.2). Lemma 2.1 ([15,20]). Consider the exponential polynomial ð0Þ
ð0Þ
ð1Þ
ð1Þ
n1 ks1 Pðk; eks1 ; . . . ; eksm Þ ¼ kn þ p1 kn1 þ þ pn1 k þ pð0Þ þ þ pn1 k þ pð1Þ þ n þ ½p1 k n e ðmÞ
ðmÞ
þ ½p1 kn1 þ þ pn1 k þ pnðmÞ eksm ; ðiÞ
where s P 0 ði ¼ 1; 2; . . . ; mÞ and pj ði ¼ 0; 1; 2; . . . ; m; j ¼ 1; 2; . . . ; nÞ are constants. As ðs1 ; s2 ; . . . ; sm Þ vary, the sum of the order of the zeros of Pðk; eks1 ; . . . ; eksm Þ on the open right half plane can change only if a zero appears on or crosses the imaginary axis. From Lemma 2.1, it is easy to obtain the following theorem:
4580
Y. Li, C. Li / Applied Mathematics and Computation 219 (2013) 4576–4589
Theorem 2.1. Assume that ðH1 Þ ðH2 Þ hold, then the following results hold: (i) For system (1.3), its zero solution is asymptotically stable for s 2 ½0; s0 Þ. (ii) System (1.3) undergoes a Hopf bifurcation at the origin when s ¼ s0 . That is, system (1.3) has a branch of periodic solutions bifurcating from the zero solution near s ¼ s0 . 3. Direction and stability of Hopf bifurcations In this section, we study the direction of bifurcations and the stability of bifurcating periodic solutions of system (1.3) at
s0 by using the method based on the normal form theory and center manifold theory introduced by Hassard et al. in [21]. ^ðsÞ ¼ yðssÞ; s ¼ s0 þ l; l 2 R, and still denote by Now, we re-scale the time by t ¼ ss; ^ xðsÞ ¼ xðssÞ; y ^ðsÞ, then system (2.1) can be written as xðtÞ ¼ ^ xðsÞ; yðtÞ ¼ y
8_ xðtÞ ¼ ðs0 þ lÞ½a1 xðt 1Þ a2 yðtÞ a3 xðtÞxðt 1Þ a4 xðtÞyðtÞ; > < " # X ð1Þ i 1 j k _ > f x ðt 1Þy ðt 1Þy ðtÞ : : yðtÞ ¼ ðs0 þ lÞ cð1 mÞr 2 xðt 1Þ r2 yðt 1Þ þ i!j!k! ijk
ð3:1Þ
iþjþkP2
T
2
For u ¼ ðu1 ; u2 Þ 2 C½1; 0 ¼ Cð½1; 0; R Þ, define a family of operators
Ll u ¼ B1 uð0Þ þ B2 uð1Þ;
ð3:2Þ
where
B1 ¼ ðs0 þ lÞ
0 a2 0
;
0
B2 ¼ ðs0 þ lÞ
a1
0
cð1 mÞr2 r2
and define
! a3 u1 ð0Þu1 ð1Þ a4 u1 ð0Þu2 ð0Þ f ðl; uÞ ¼ ðs0 þ lÞ P : ð1Þ i j 1 k iþjþkP2 i!j!k! fijk u1 ð1Þu2 ð1Þu2 ð0Þ By the Riesz representation theorem, there exists a matrix whose components are bounded variation functions gðh; lÞ : ½1; 0 ! R2 such that
Ll u ¼
Z
0
dgðh; lÞuðhÞ:
1
In fact, choosing
gðh; lÞ ¼ B1 dðhÞ þ B2 dðh þ 1Þ;
ð3:3Þ
1; 0;
h¼0 where dðhÞ ¼ is a Dirac delta function, then (3.2) is satisfied. h–0 For u ¼ ðu1 ; u2 ÞT 2 C½1; 0 ¼ Cð½1; 0; R2 Þ, define
(
AðlÞu ¼
u_ ðhÞ; h 2 ½1; 0Þ 0; h 2 ½1; 0Þ; ; RðlÞu ¼ R0 l ; u Þ; h ¼ 0: f ð d g ðs; l Þ u ðsÞ; h ¼ 0 1
ð3:4Þ
Hence, (3.1) can be rewritten as
U_ t ¼ AðlÞU t þ RðlÞU t ;
ð3:5Þ
T
where U ¼ ðx; yÞ and U t ðhÞ ¼ Uðt þ hÞ; h 2 ð1; 0. For w 2 C 1 ½0; 1, define the adjoint operator A of A as
(
A wðsÞ ¼
_ wðsÞ; s 2 ð0; 1; R0 T dg ðt; 0ÞwðtÞ; s ¼ 0; 1
ð3:6Þ
where gT is the transpose of the matrix g. For u 2 Cð½1; 0; R2 Þ and w 2 Cð½0; 1; ðR2 Þ Þ, in order to normalize the eigenvectors of operator A and adjoint operator A , we define a bilinear inner product
hw; ui ¼ wT ð0Þuð0Þ
Z
0
h¼1
Z
n¼0
where gðhÞ ¼ gðh; 0Þ. Then as usual,
hw; Aui ¼ hA w; ui
h
wT ðn hÞ dgðhÞuðnÞ dn;
Y. Li, C. Li / Applied Mathematics and Computation 219 (2013) 4576–4589
4581
for ðu; wÞ 2 DðAÞ DðA Þ. Normalizing q and q by the condition
hq ; qi ¼ 1;
i ¼ 0: hq ; q
By the discussion in Section 2 and transformation t ¼ ss, we know that is0 x0 are eigenvalues of A and other eigenvalues have strictly negative real parts. Thus, they are also eigenvalues of A . Next we calculate the eigenvector q of A associated to the eigenvalue is0 x0 and the eigenvector q of A associated to the eigenvalue is0 x0 . Let
qðhÞ ¼
1 ix0 s0 h e ;
1 6 h < 0:
q
ð3:7Þ
From the above discussion, it is easy to know that
Aqð0Þ ¼ is0 x0 qð0Þ; i.e.,
s0
ix0 þ a1 eix0 s0 cr 2 ð1 mÞeix0 s0
! 1
a2 ix0 þ r2 eix0 s0
q
¼
0 : 0
ð3:8Þ
So we have
q¼
a1 eis0 x0 þ ix0 : a2
Similarly, assume that the eigenvector q of A is
q ðsÞ ¼
1 D
1
r
eix0 s0 s ;
0 6 s < 1:
ð3:9Þ
Then we obtain
A q ð0Þ ¼ is0 x0 q ð0Þ; i.e.,
s0
ix0 þ a1 eix0 s0
cr 2 ð1 mÞeix0 s0
a2
ix0 þ r 2 eix0 s0
!
1
r
¼
0 : 0
ð3:10Þ
Hence, we obtain
r¼
ix0 þ a1 eix0 s0 : cr 2 ð1 mÞeix0 s0
By hq ; qi ¼ 1, we can get q, that is,
ð0ÞÞT qð0Þ hq ; qi ¼ ðq ¼
1 D
Z
0 h¼1
Z
Z Þ 1 s0 ð 1r
q
h
ÞT ðn hÞ dgðhÞqðnÞ dn ðq
n¼0 0
Z
h¼1
h
ÞeiðnhÞx0 dgðhÞ ð 1r
n¼0
1 inx0 e dn
q
Z 0 1 1 s0 þ s0 ða1 r Þ ¼ 1: Þheix0 s0 h dgðhÞ 1 cð1 mÞr 2 r 2 qr ¼ ½1 þ qr 1 þ qr ¼ ð 1r q D D 1 Hence, we have
þ s0 ða1 r Þ; cð1 mÞr 2 r 2 qr D ¼ 1 þ qr such that
hq ðhÞ; qðhÞi ¼ 1;
ðhÞi ¼ 0: hq ðhÞ; q
So q and q are obtained. In the remainder of this section, we use the same notations as those in Ref. [21]. We first compute the coordinates to describe the center manifold C0 at l ¼ 0. Let U t be the solution of (3.5) when l ¼ 0, and define
zðtÞ ¼ hq ; U t i;
Wðt; hÞ ¼ U t ðhÞ 2RefzðtÞqðhÞg:
ð3:11Þ
On the center manifold C0 , we have Wðt; hÞ ¼ WðzðtÞ; zðtÞ; hÞ, where
Wðz; z; hÞ ¼ W 20 ðhÞ
z2 z2 þ W 11 ðhÞzz þ W 02 ðhÞ þ ; 2 2
ð3:12Þ
4582
Y. Li, C. Li / Applied Mathematics and Computation 219 (2013) 4576–4589
. Note that W is real if U t is real, we consider z and z are local coordinates for center manifold C0 in the direction of q and q only real solutions. For the solution U t 2 C0 , since l ¼ 0, then
z_ ðtÞ ¼ hq ; U_ t i ¼ hq ; AU t þ RU t i ¼ hq ; AU t i þ hq ; RU t i ¼ hA q ; U t i þ hq ; RU t i ð0Þf ð0; Wðt; 0Þ þ 2Re½zðtÞqð0ÞÞ: ¼ is0 x0 zðtÞ þ q
ð3:13Þ
We rewrite this equation as
z_ ðtÞ ¼ is0 x0 zðtÞ þ gðz; zÞ
ð3:14Þ
with
gðz; zÞ ¼ g 20
z2 z2 z2 z þ þ g 11 zz þ g 02 þ g 21 2 2 2
In the following, the values of
ð3:15Þ
l2 and b2 presented by Hassard at [21] are obtained. Substituting (3.5) and (3.13) into
_ ¼ U_ t z_ q z_ q ; W we have
_ ¼ U_ t z_ q z_ q ¼ AU t þ RU t ½is0 x0 z þ q ð0Þf ðz; zÞq ½is0 x0 z þ q ð0Þf ðz; zÞq W ð0Þf ðz; zÞqðhÞ 2Re½is0 x0 zqðhÞ ¼ AW 2Re½q ð0Þf ðz; zÞqðhÞ þ RU t ¼ AW þ 2AReðzqÞ þ RU t 2Re½q 8 ð0Þf ðz; zÞqðhÞ; 1 6 h < 0; < AW 2Re½q ð0Þf ðz; zÞqðhÞ þ f ðl; uÞ; h ¼ 0 ¼ AW 2Re½q : def ¼ AW þ Hðz; z; hÞ;
ð3:16Þ
where
Hðz; z; hÞ ¼ h20 ðhÞ
z2 z2 þ h11 ðhÞzz þ h02 ðhÞ þ 2 2
ð3:17Þ
Taking the derivative of W with respect to t in (3.12), we have
_ ¼ W z z_ þ W z z_ : W
ð3:18Þ
From (3.12), (3.14) and (3.18), we obtain
_ ¼ ðW 20 z þ W 11 z þ Þðis0 x0 z þ gÞ þ ðW 11 z þ W 02 z þ Þðis0 x0 z þ gÞ: W
ð3:19Þ
Then substituting (3.12) and (3.17) into (3.16), we have 2 2 _ ¼ ðAW 20 þ h20 Þ z þ ðAW 11 þ h11 Þzz þ ðAW 02 þ h02 Þ z þ W 2 2
ð3:20Þ
By comparing the coefficients of (3.19) with (3.20), the following hold:
ðA 2is0 x0 ÞW 20 ðhÞ ¼ h20 ðhÞ;
ð3:21Þ
AW 11 ðhÞ ¼ h11 ðhÞ:
ð3:22Þ
According to (3.13) and (3.14), we obtain
ð0Þf ðz; zÞ ¼ gðz; zÞ ¼ q
s0 D
Þ P ð1 r
a3 U 1t ð0ÞU 1t ð1Þa4 U 1t ð0ÞU 2t ð0Þ ð1Þ i j k 1 iþjþkP2 i!j!k! fijk U 1t ð1ÞU 2t ð1ÞU 2t ð0Þ
ðhÞ and qðhÞ ¼ ð 1 where U t ðhÞ ¼ ðU 1t ðhÞ; U 2t ðhÞ; U 3t ðhÞÞT ¼ Wðt; hÞ þ zqðhÞ þ zq
! ;
q ÞT eix0 h . Noting that
z2 z2 ð1Þ ð1Þ þ W 11 ð0Þzz þ W 02 ð0Þ þ ; 2 2 2 z2 z ð2Þ ð2Þ z þ W ð2Þ þ W 11 ð0Þzz þ W 02 ð0Þ þ ; U 2t ð0Þ ¼ qz þ q 20 ð0Þ 2 2 2 z2 z ð1Þ ð1Þ ð1Þ U 1t ð1Þ ¼ zeis0 x0 þ zeis0 x0 þ W 20 ð1Þ þ W 11 ð1Þzz þ W 02 ð1Þ þ ; 2 2 2 z2 z ð2Þ ð2Þ zeis0 x0 þ W ð2Þ þ W 11 ð1Þzz þ W 02 ð1Þ þ : U 2t ð1Þ ¼ qzeis0 x0 þ q 20 ð1Þ 2 2 ð1Þ
U 1t ð0Þ ¼ z þ z þ W 20 ð0Þ
ð3:23Þ
4583
Y. Li, C. Li / Applied Mathematics and Computation 219 (2013) 4576–4589
Therefore, we have
(" !# ð1Þ f 1 ð1Þ ð1Þ ð1Þ 2 ix0 s0 z2 a4 q a3 eix0 s0 þ r 200 e2ix0 s0 þ f110 qe2ix0 s0 þ f101 qeix0 s0 þ f011 qe 2 D h i ð1Þ ð1Þ ð1Þ ð1Þ f200 þ qÞ þ r Þ þ f101 eix0 s0 þ qeix0 s0 Þ þ f011 qeix0 s0 Þ zz eix0 s0 þ q þ a3 ðeix0 s0 þ eix0 s0 Þa4 ðq r þ f110 ðq þ q ðq ðqq " ! ! ð1Þ ð1Þ ð2Þ ð1Þ W 20 ð1Þ W 20 ð0Þ ix0 s0 W 20 ð0Þ W 20 ð0Þ ð1Þ ð1Þ ð2Þ ð1Þ ix0 s0 þ W 11 ð0Þe q þ qW 11 ð0Þ þ e þ þ a3 W 11 ð1Þ þ a4 W 11 ð0Þ þ 2 2 2 2
gðz; zÞ ¼
f200 ð1Þ ð1Þ 2W 11 ð1Þeix0 s0 þ W 20 ð1Þeix0 s0 2 ð1Þ
þr
ð2Þ
ð1Þ
ð1Þ
W 20 ð1Þ ix0 s0 W 20 ð1Þ ix0 s0 ð1Þ e þ q e þ W 11 ð1Þqeix0 s0 2 2 ! ð2Þ ð1Þ W 20 ð0Þ ix0 s0 W 20 ð1Þ ð2Þ ð1Þ ix0 s0 W 11 ð0Þe þ þ q þ W 11 ð1Þq e 2 2 ! ð2Þ ð2Þ W ð1Þ W ð0Þ ð2Þ ix0 s0 W 11 ð1Þq þ 20 q þ 20 q eix0 s0 þ W ð2Þ ð0Þ q e 11 2 2
!
ð2Þ
þ f110 W 11 ð1Þeix0 s0 þ ð1Þ
þ f101 ð1Þ
þ f011
!# ð1Þ ð1Þ ð1Þ f300 ix0 s0 f210 f ð1Þ þ 2qÞeix0 s0 þ 201 ðq e2ix0 s0 þ 2qÞ þ f111 þ qðq þ q þ q e2ix0 s0 Þ z2z e ðq 2 2 2 " !# ) ð1Þ f200 2ix0 s0 ð1Þ ð1Þ ð1Þ 2 ix0 s0 ix0 s0 2i x s i x s 0 0 0 0 z2 þ : þr e e e a4 q þ f110 q þ f101 q þ f011 q e þ a3 e 2 þ
Comparing the coefficients of (3.15) with the above equality, it follows that
" !# ð1Þ f200 2ix0 s0 2s0 ð1Þ ð1Þ ð1Þ 2 ix0 s0 ix0 s0 2i x s i x s 0 0 0 0 ; a4 q a3 e g 20 ¼ þr þ f110 qe þ f101 qe þ f011 q e e 2 D h i s0 ð1Þ ð1Þ ð1Þ ð1Þ f200 eix0 s0 þ q þ qÞ þ r Þ þ f101 eix0 s0 þ qeix0 s0 Þ þ f011 qeix0 s0 Þ ; g 11 ¼ r þ f110 ðq þ q ðq ðqq a3 ðeix0 s0 þ eix0 s0 Þ a4 ðq D " !# ð1Þ f200 2ix0 s0 2s0 ð1Þ ð1Þ ð1Þ 2 ix0 s0 ix0 s0 2i x s i x s 0 0 0 0 þr e e e ; a4 q þ f110 q þ f101 q þ f011 q a3 e e g 02 ¼ 2 D ! " ! ð1Þ ð1Þ ð2Þ ð1Þ W ð1Þ W 20 ð0Þ ix0 s0 W ð0Þ W 20 ð0Þ 2s0 ð1Þ ð1Þ ð2Þ g 21 ¼ a3 W 11 ð1Þ þ 20 þ W 11 ð0Þeix0 s0 a4 W 11 ð0Þ þ 20 q þ qW ð1Þ ð0Þ þ e þ 11 2 2 2 2 D f200 ð1Þ ð1Þ 2W 11 ð1Þeix0 s0 þ W 20 ð1Þeix0 s0 2 ð1Þ
þr ð1Þ
! ð2Þ ð1Þ W 20 ð1Þ ix0 s0 W 20 ð1Þ ix0 s0 ð1Þ þ q e þ W 11 ð1Þqeix0 s0 e 2 2 ! ð2Þ ð1Þ W 20 ð0Þ ix0 s0 W 20 ð1Þ ð2Þ ð1Þ ix0 s0 e W 11 ð0Þe þ þ q þ W 11 ð1Þq 2 2 ! ð2Þ ð2Þ ð1Þ ð1Þ W 20 ð1Þ W 20 ð0Þ ix0 s0 f f ð2Þ ð2Þ ix0 s0 þ 2qÞeix0 s0 þ 300 eix0 s0 þ 210 ðq W 11 ð1Þq þ qþ qe þ W 11 ð0Þqe 2 2 2 2 !# ð2Þ
þf110 W 11 ð1Þeix0 s0 þ ð1Þ
þf101 ð1Þ
þf011 ð1Þ
þ
f201 ð1Þ e2ix0 s0 þ 2qÞ þ f111 qðq þ q þ q e2ix0 s0 Þ ðq 2
:
Now we compute W 20 ðhÞ and W 11 ðhÞ. From (3.16), we have
ð0Þf ðz; zÞqðhÞÞ þ RU t ¼ gqðhÞ gq ðhÞ þ RU t Hðz; z; hÞ ¼ 2Reðq 1 1 1 1 ðhÞ þ RU t : g20 z2 þ g11 zz þ g02 z2 þ q ¼ g 20 z2 þ g 11 zz þ g 02 z2 þ qðhÞ 2 2 2 2
ð3:25Þ
Comparing the coefficients with (3.17), we can obtain, for h 2 ½1; 0Þ,
ðhÞ; h11 ðhÞ ¼ g 11 qðhÞ g11 q ðhÞ: h20 ðhÞ ¼ g 20 qðhÞ g02 q
ð3:26Þ
4584
Y. Li, C. Li / Applied Mathematics and Computation 219 (2013) 4576–4589
Substituting the above equalities into (3.21) and (3.22), respectively, it follows that
(
_ 20 ðhÞ ¼ 2is0 x0 W 20 ðhÞ þ g qðhÞ þ g02 q ðhÞ; W 20 _ 11 ðhÞ ¼ g qðhÞ þ g11 q ðhÞ: W
ð3:27Þ
11
Solving (3.27), we have
W 20 ðhÞ ¼
ig 20
s0 x0
qð0Þeis0 x0 h þ
ig02 ð0Þeis0 x0 h þ Ee2is0 x0 h q 3s 0 x0
ð3:28Þ
and
ig 11
W 11 ðhÞ ¼
s0 x0
qð0Þeis0 x0 h þ
g11 ð0Þeis0 x0 h þ F: q is0 x0
ð3:29Þ
In what follows, we seek appropriate E and F. From (3.21) and (3.22), we can get
AW 20 ð0Þ ¼ 2is0 x0 W 20 ð0Þ h20 ð0Þ
ð3:30Þ
AW 11 ð0Þ ¼ h11 ð0Þ:
ð3:31Þ
and
From the definition of A in (3.4), we obtain
Z
0
dgðhÞW 20 ðhÞ ¼ 2is0 x0 W 20 ð0Þ h20 ð0Þ
ð3:32Þ
dgðhÞW 11 ðhÞ ¼ h11 ð0Þ:
ð3:33Þ
1
and
Z
0
1
From (3.26) and the definition of f ðl; sÞ, we have
ð0Þ þ s0 h20 ð0Þ ¼ g 20 qð0Þ g02 q
d1
ð3:34Þ
d2
and
ð0Þ þ s0 h11 ð0Þ ¼ g 11 qð0Þ g11 q
d3 ; d4
ð3:35Þ
where
d1 ¼ 2a4 q 2a3 eix0 s0 ; ð1Þ
ð1Þ
ð1Þ
ð1Þ
d2 ¼ f200 e2ix0 s0 þ 2f110 qe2ix0 s0 þ 2f101 qeix0 s0 þ 2f011 q2 eix0 s0 ; þ qÞ; d3 ¼ a3 ðeix0 s0 þ eix0 s0 Þ a4 ðq ð1Þ
ð1Þ
ð1Þ
ð1Þ
eix0 s0 þ q Þ þ f101 ðq eix0 s0 þ qeix0 s0 Þ þ f011 ðqq qeix0 s0 Þ: d4 ¼ f200 þ f110 ðq þ q Substituting (3.28) and (3.34) into (3.32) and noticing that
Z is0 x0 I
0
eis0 x0 h dgðhÞ qð0Þ ¼ 0
1
and
Z is0 x0 I
ð0Þ ¼ 0; eis0 x0 h dgðhÞ q
0
1
we obtain
d1 ; e2is0 x0 h dgðhÞ E ¼ s0 d2 1
Z 2is0 x0 I
0
ð3:36Þ
which leads to
2ix0 þ a1 e2ix0 s0 2ix0 s0
cr 2 ð1 mÞe It follows that
!
a2 2ix0 s0
2ix0 þ r2 e
E¼
d1 d2
:
4585
Y. Li, C. Li / Applied Mathematics and Computation 219 (2013) 4576–4589
E¼
!1
2ix0 þ a1 e2ix0 s0
a2
cr2 ð1 mÞe2ix0 s0
2ix0 þ r 2 e2ix0 s0
d1
ð3:37Þ
:
d2
Similarly, substituting (3.29) and (3.35) into (3.31), we have
F¼
a1
a2
cr2 ð1 mÞ
r2
1
d3 d4
:
ð3:38Þ
As a result, we know that W 20 ðhÞ and W 11 ðhÞ, then g ij in (3.24) is determined by the parameters and the time delay in (2.2). Thus, we can compute the following quantities
C1 ð0Þ ¼
1 g g 20 g 11 2jg 11 j2 jg 02 j2 þ 21 ; 2s 0 x0 3 2 i
l2 ¼
ReC1 ð0Þ ; Rek0 ðs0 Þ
t2 ¼
ImC1 ð0Þ þ l2 Imk0 ðs0 Þ
s0 x0
b2 ¼ 2ReC1 ð0Þ:
;
ð3:39Þ From the expression of C1 ð0Þ in (3.39), it is easy to get the values of l2 ; b2 and t 2 . Thus, if system (1.3) is given, then we will analyze the direction of the Hopf bifurcation and stability of the bifurcating periodic orbits at s ¼ s0 according to the following results: (i) (ii)
l ¼ 0 is Hopf bifurcation value of system (3.1); l2 determines the direction of the Hopf bifurcation: if l2 > 0ð< 0Þ, then the Hopf bifurcation is supercritical ðsubcriticalÞ and the bifurcating periodic solutions exist for s > s0 ð< s0 Þ;
(iii) b2 determines the stability of the bifurcating periodic solutions: if b2 < 0ð> 0Þ the bifurcating periodic solutions are stable (unstable); (iv) t2 determines the period of the bifurcating periodic solutions: the period increase (decrease) if t 2 > 0ð< 0Þ.
1.02
1 0.98 0.96
0.96
0.94
0.94
0.92
y
0.98
0.92
0.9
0.9
0.88
0.88
0.86
0.86
0.84
0.84
0
20
40
60
80
0.82
100
t
0
20
40
t
60
80
1 0.98 0.96 0.94
y
x
1
0.92 0.9 0.88 0.86 0.84 0.82 0.84 0.86 0.88 0.9
Fig. 1. Behavior and phase portrait of system (1.3) with
0.92 0.94 0.96 0.98
x
1
1.02
s ¼ 0:5 < s0 , the positive equilibrium E ðx ; y Þ ¼ ð0:8696; 0:8696Þ is stable.
100
4586
Y. Li, C. Li / Applied Mathematics and Computation 219 (2013) 4576–4589
4. Numerical example To illustrate the theoretical results, let us give some numerical simulations. Example 4.1. Consider the following system.
(
sÞ x0 ðtÞ ¼ xðtÞð1 xðt Þ 1:3ð1 0:5ÞxðtÞyðtÞ; 2
ð4:1Þ
yðtsÞ y0 ðtÞ ¼ 1:5yðtÞð1 2ð10:5Þxðt sÞÞ;
where r 1 ¼ 1; r 2 ¼ 1:5; K ¼ 2; a ¼ 1:3; c ¼ 2; m ¼ 0:5. The positive equilibrium is E ¼ ð0:8696; 0:8696Þ. (2.6) has only one positive root. Furthermore, we can obtain that s0 ¼ 0:6729. Simple calculations show that conditions of Theorem 2.1 hold. We can see from Figs. 1 and 2 that E is asymptotically stable when s ¼ 0:5 < s0 while it loses its stability when s ¼ 0:7 > s0 . Remark 4.1. In Ref. [12], the authors considered system (1.3) with s ¼ 0, they obtained the system is globally stable. However, system (1.3) is not always stable, which is illustrated in the above example. That is, time delays even if very little ones could disturb the stability of dynamical systems.
5. Global continuation of local Hopf bifurcation In this section, we study the existence of global Hopf bifurcations. The method we used here is based on the global Hopf bifurcating theorem for general functional differential equations introduced by Wu in [22]. For convenience, we write system (1.3) as the following form:
z_ ðtÞ ¼ Fðzt ; s; pÞ;
ð5:1Þ
where z ¼ ðx; y1 ; y2 ÞT ; zt ðhÞ 2 Cð½s; 0; R3 Þ. Following the work of Wu [19], we define
1
1.05 1
0.95
0.95 0.9
y
x
0.9 0.85
0.85
0.8 0.8
0.75
0
0.75 20
40
60
t
80
0.7
100
0
20
40
t
60
80
100
1.1
y
1 0.9 0.8 0.7 1 0.9 0.8
x Fig. 2. Behavior and phase portrait of system (1.3) with occurs.
0.7
0
20
40
60
80
100
t
s ¼ 0:7 > s0 , the positive equilibrium E ð0:8696; 0:8696Þ loses its stability and Hopf bifurcation
Y. Li, C. Li / Applied Mathematics and Computation 219 (2013) 4576–4589
4587
X ¼ Cð½s; 0; R3 Þ;
R ¼ ClfðzðtÞ; s; pÞ 2 X R Rþ ; zðtÞis a T-periodic solution of Ref: ð5:1Þg; N ¼ fðz; s; pÞ; Fðz; s; pÞ ¼ 0g: 2p Let lðz ;sj ; 2p Þ be the connected component of ðz ; sj ; x Þ in R, where 0 k
x0
sj and x0 are defined in (2.9).
Lemma 5.1. If s is bounded, all solutions of system (1.3) are uniformly bounded. Proof. Let ðxðtÞ; yðtÞÞT is a solution of system (1.3) through ð/; wÞ at t ¼ 0 with /ð0Þ > 0; wð0Þ > 0. Then we know that ðxðtÞ; yðtÞÞT is positive on the maximum existence interval of solution. h If ðxðtÞ; yðtÞÞT is a solution of (1.3) with xðtÞ > 0; yðtÞ > 0, then it follows from the first equation of (1.3) that
xðt sÞ : x0 ðtÞ < r 1 xðtÞ 1 K
ð5:2Þ
Clearly, xðtÞ < xðt sÞer1 s for t > s, which implies that xðt sÞ > xðtÞer1 s . This, together with (5.2), leads to
xðtÞer1 s ; x0 ðtÞ < r 1 xðtÞ 1 K
for; t > s;
which implies that
lim sup xðtÞ 6 Ker1 s : t!þ1
Thus, for any
e > 0, there exists a T > 0ðT > sÞ such that when t > T, we have
xðtÞ 6 Ker1 s þ e :¼ d: It follows from the second equation of (1.3) that, for t > T þ s,
yðt sÞ y0 ðtÞ < r 2 yðtÞ 1 : cð1 mÞd
ð5:3Þ
Then by similar discussion as above, we can easily have
lim sup yðtÞ 6 cð1 mÞder2 s : t!þ1
Thus the possible positive solutions of system (1.3) must be uniformly bounded. This completes the proof of the lemma. Lemma 5.2. System (1.3) has no nontrivial periodic solutions of period s.
Proof. Assume system (1.3) has a nontrivial periodic solution of period s, then the following differential system, which has the same equilibrium as (1.3),
(
x0 ðtÞ ¼ r 1 xðtÞð1 xðtÞ Þ að1 mÞxðtÞyðtÞ; K
ð5:4Þ
yðtÞ y0 ðtÞ ¼ r 2 yðtÞð1 cð1mÞxðtÞ Þ;
has periodic solutions. Due to Lemma 5.1, we restrict our attention to 0 < x < K; 0 < yðtÞ < cKð1 mÞ, respectively. System (5.4) also has the equilibrium E ¼ ðx ; y Þ. Define
Vðx; yÞ ¼ ln
x x ax cð1 mÞ2 y y þ ln þ þ : x x y y r2
Obviously, Vðx; yÞ is well defined and continuous for all xðtÞ > 0; yðtÞ > 0. Calculating the derivative of V along the solution of system (1.3), we have
dV r1 að1 mÞ ¼ ðx x Þ2 ðy y Þ2 : dt y Kx Obviously, dV < 0 strictly for all x; y > 0 except the positive equilibrium ðx ; y Þ where dt Thus, Vðx; yÞ satisfies Lyapunov’s asymptotic stability theorem, we conclude that
dV dt
¼ 0.
limðxðtÞ; yðtÞÞ ¼ ðx ; y Þ;
t!1
which contradicts the fact that system (5.4) has periodic solutions. This ends the proof. h
4588
Y. Li, C. Li / Applied Mathematics and Computation 219 (2013) 4576–4589
Theorem 5.1. Suppose that ðH1 Þ ðH2 Þ hold, then for each s > sjk ðj ¼ 1; 2; . . .Þ, system (1.3) has at least j þ 1 periodic solutions, where sjk is defined in (2.8). Proof. The characteristic matrix of system (1.3) at the positive equilibrium z is of the form
Dðz ; s; pÞðkÞ ¼
k þ a1 eks
a2
cr 2 ð1 mÞeks
k þ r2 e
: ks
ð5:5Þ
2p From the discussion of Section 2, it can be verified that ðz ; sj ; x Þ; j ¼ 1; 2; . . . are isolated centers. 0 Let
2p X;x2p ¼ ðg; pÞ : 0 < g < ;
p
< :
x0
0
Clearly, if js sj j 6 d and ðg; pÞ 2 @ X , then the necessary and sufficient conditions for Dðz ; s; pÞðg þ i 2ppÞ ¼ 0 are
g ¼ 0; s ¼ sj and p ¼ x2p0 . 2p Þðg; pÞ ¼ Dðz ; sj d; pÞðg þ i 2ppÞ, then we have the transversal number Define H ðz ; sj ; x 0
c z ; sj ;
2p
x0
2p 2p ¼ degB H z ; sj ; ; X;x2p degB Hþ z ; sj ; ; X;x2p ¼ 1:
x0
x0
0
0
2p By Theorem 2.1 of Wu [22], we conclude that the connected component lðz ;sj ;x2p Þ through ðz ; sj ; x Þ in R is nonempty. Mean0 0 while, we have
X ðz;s;pÞ2l
cðz; s; pÞ < 0
2p Þ ðz ;sj ;x 0
and hence lðz ;sj ;x2p Þ is unbounded. 0
From (2.8), we see that, for j P 1,
2p < sjk : x0 ; 1Þ, where s < sjk . Clearly, it follows from Now, we are in a position to prove that the projection of lðz ;sj ; 2p Þ onto s-space is ½s k
x0
the proof of Lemma 5.2 that system (1.3) with s ¼ 0 has no nontrivial periodic solution. Hence, the projection of lðz ;sj ; 2p Þ onto k
x0
s-space is away from zero. For a contradiction, we suppose that the projection of lðz ;sj ; 2p Þ onto s-space is bounded. This means that the projection of k
x0
2p lðz ;sj ; 2p Þ onto s-space is included in an interval ð0; s Þ. Noting that x < sjk and applying Lemma 5.2, we have p < s for ðz; s; pÞ 0 k
x0
belonging to lðz ;sj ; 2p Þ . This implies that the projection of the connected component lðz ;sj ; 2p Þ onto p-space is bounded. k
x0
k
x0
Moreover, from Lemma 5.1, we obtain that the projection of lðz ;sj ; 2p Þ onto z-space is bounded if the projection of lðz ;sj ; 2p Þ onto k
x0
k
x0
s-space is bounded. Thus, the connected component lðz ;sj ; 2p Þ crossing through ðE ; sjk ; x2p0 ÞÞ is bounded, which is a k
x0
; 1Þ for each j P 1, where s
of the proof.
x0
h
References [1] C. Ji, D. Jiang, N. Shi, Analysis of a predator–prey model with modified Leslie-Gower and Holling-type II schemes with stochastic perturbation, J. Math. Anal. Appl. 359 (2009) 482–498. [2] Y. Song, S. Yuan, J. Zhang, Bifurcation analysis in the delayed Leslie-Gower predator–prey system, Appl. Math. Model. 33 (2009) 4049–4061. [3] P.H. Leslie, Some further notes on the use of matrices in population mathematics, Biometrika 35 (1948) 213–245. [4] Y. Chen, Z. Liu, M. Haque, Analysis of a Leslie-Gower-type prey-predator model with periodic impulsive perturbations, Commun. Nonlinear Sci. Numer. Simulat. 14 (2009) 3412–3423. [5] P. Aguirre et al, Two limit cycles in a Leslie-Gower predator–prey model with additive Allee effect, Nonlinear Anal. Real World Appl. 10 (2009) 1401– 1416. [6] A. Korobeinikov, A Lyapunov function for Leslie-Gower predator–prey models, Appl. Math. Lett. 14 (2001) 697–699. [7] A.F. Nindjin, M.A. Aziz-Alaoui, M. Cadivel, Analysis of a predator–prey model with modified Leslie-Gower and Holling-type II schemes with time delay, Nonlinear Anal. Real World Appl. 7 (5) (2007) 1104–1118. [8] Z. Jing, J. Yang, Bifurcation and chaos in discrete-time predator–prey system, Chaos Solitons Fract. 27 (2006) 259–277. [9] W. Wang, L. Chen, A predator–prey system with stage structure for predator, Comput. Math. Appl. 33 (1997) 83–91. [10] Y. Yang, Hopf bifurcation in a two-competitor, one-prey system with time delay, Appl. Math. Comput. 214 (2009) 228–235. [11] T. Kumar Kar, Stability analysis of a prey-predator model incorporating a prey refuge, Commun. Nonlinear Sci. Numer. Simulat. 10 (2005) 681–691. [12] F. Chen, L. Chen, X. Xie, On a Leslie-Gower predator–prey model incorporating a prey refuge, Nonlinear Anal. Real World Appl. 10 (2009) 2905–2908. [13] J.M. Cushing, Integrodifferential Equations and Delay Models in Population Dynamics, Springer Verlag, Heidelberg, 1977.
Y. Li, C. Li / Applied Mathematics and Computation 219 (2013) 4576–4589
4589
[14] S.A. Gourley, M.V. Bartuccelli, Parameter domains for instability of uniform states in system with many delays, J. Math. Biol. 35 (1997) 843–867. [15] W. Yu, J. Cao, Stability and Hopf bifurcation analysis on a four-neuron BAM neural networks with time delays, Phy. Lett. A 351 (2006) 64–78. [16] S. Toaha, M.A. Hassan, Stability analysis of predator–prey population model with time delay and constant rate of harvesting, Punjab Univ. J. Math. 40 (2008) 37–48. [17] R. Xu, Z. Ma, Stability and Hopf bifurcation in a predator–prey model with stage structure for the predator, Nonlinear Anal. Real World Appl. 9 (2008) 1444–1460. [18] S. Yuan, Y. Song, Stability and Hopf bifurcations in a delayed Leslie-Gower predator–prey system, J. Math. Anal. Appl. 355 (2009) 82–100. [19] J. Hale, Theory of Functional Differential Equations, Springer-Verlag, Heidelberg, 1977. [20] S. Ruan, J. Wei, On the zeros of transcendental functions with applications to stability of delay differential equations with two delays, Dynam. Contin. Discrete Impuls. Syst. Ser. A Math. Anal. 10 (2003) 863–874. [21] B. Hassard, N. Kazarinoff, Y.H. Wan, Theory and Applications of Hopf Bifurcation, London Mathematical Society Lecture Note Series, vol. 41, Cambridge University Press, Cambridge, 1981. [22] J. Wu, Symmetric functional differential equations and neural networks with memory, Trans. Amer. Math. Soc. 350 (1998) 4799–4838.