- Email: [email protected]
Chaos and Hopf bifurcation analysis for a two species predator–prey system with prey refuge and diffusion
Nonlinear Analysis: Real World Applications 12 (2011) 1047–1061
Contents lists available at ScienceDirect
Nonlinear Analysis: Real World Applications journal homepage: www.elsevier.com/locate/nonrwa
Chaos and Hopf bifurcation analysis for a two species predator–prey system with prey refuge and diffusion✩ Xia Liu a , Maoan Han b,∗ a
College of Mathematics and Information Science, Henan Normal University, 453007, PR China
b
Department of Mathematics, Shanghai Normal University, Shanghai 200234, PR China
article
abstract
info
Article history: Received 8 May 2010 Accepted 23 August 2010
In this paper, a delayed predator–prey model with Holling type II functional response incorporating a constant prey refuge and diffusion is considered. By analyzing the characteristic equation of linearized system corresponding to the model, we study the local asymptotic stability of the positive equilibrium of the system. By choosing the time delay due to gestation as a bifurcation parameter, the existence of Hopf bifurcations at the positive equilibrium is established. By applying the normal form and the center manifold theory, an explicit algorithm to determine the direction of the Hopf bifurcation and the stability of the bifurcating periodic solutions is derived. Further, an example is presented to illustrate our main results. Finally, recurring to the numerical method, the influences of impulsive perturbations on the dynamics of the system are also investigated. © 2010 Elsevier Ltd. All rights reserved.
Keywords: Diffusion Hopf bifurcation Chaos Prey refuge
1. Introduction The predator–prey systems with prey refuge have been studied by several authors; one can refer to [1–5] and the references cited therein. Especially, González-Olivares et al. [1] proposed a predator–prey system with Holling type II functional response incorporating a constant prey refuge as follows
x˙ = α x 1 − y˙ = −d2 y +
x
−
β(x − m)y
k 1 + a(x − m) c β(x − m)y 1 + a( x − m )
, (1.1)
.
They investigated the local stability of equilibria and the existence of limit cycle of system (1.1). Motivated by the work of González-Olivares et al. [1], Chen et al. [2] further discussed the instability and global stability properties of the equilibria and the existence and uniqueness of limit cycle of system (1.1). Many papers have researched the predator–prey systems with diffusion [6–8]. As the authors in [9] noted, the dispersal between patches often occurs in ecological environments, and more realistic models should include the dispersal process.
✩ This work was supported by Shanghai Leading Academic Discipline Project (No. S30405).
∗
Corresponding author. E-mail address: [email protected] (M. Han).
1468-1218/$ – see front matter © 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.nonrwa.2010.08.027
1048
X. Liu, M. Han / Nonlinear Analysis: Real World Applications 12 (2011) 1047–1061
Based on the above discussion, by incorporating dispersal of prey and time delay into system (1.1), we establish the following model x1
x˙ 1 = α x1 1 −
k
−
β(x1 − m)y 1 + a(x1 − m)
+ D1 (x2 − x1 ),
x˙ 2 = −d1 x2 + D2 (x1 − x2 ), y˙ = −d2 y +
(1.2)
c β[x1 (t − τ ) − m]y 1 + a[ x 1 ( t − τ ) − m ]
,
with initial conditions x1 (θ ) = ϕ1 (θ ) ⩾ m,
x2 (θ ) = ϕ2 (θ ) ⩾ 0 ϕ2 (0) > 0, ψ(0) > 0,
ϕ1 (0) > m,
y(θ ) = ψ(θ ) ⩾ 0,
θ ∈ [−τ , 0],
(1.3)
where xi is the prey density in the ith patch, i = 1, 2, and y is the predator density; Di > 0 represents the dispersal rate of prey species, i = 1, 2; d1 > 0 denotes the maximal death rate of x2 in second patch; τ is a constant time delay due to gestation; the parameters α , k, d2 and c are positive constants and their ecological meaning can be found in Refs. [1,2]. The rest of the paper is organized as follows. The existence and local property of the positive equilibrium and the existence of local Hopf bifurcation of the model are discussed in Section 2. In Section 3, by applying the normal form theory and the center manifold reduction for functional differential equations, the conditions for determining the bifurcation direction and the stability of the bifurcating periodic solution are derived. We also give a numerical example to illustrate the validity of our results. Finally, in Section 4, by using the numerical method, we investigate the influences of impulsive perturbations on system (1.2). 2. Stability of the positive equilibrium and the existence of local Hopf bifurcation For simplicity, we take x¯ 1 = x1 − m, x¯ 2 = x2 and y¯ = β y so that system (1.2) can be rewritten as (still denote x¯ 1 , x¯ 2 and y¯ as x1 , x2 and y, respectively) x˙ 1 = p(x1 + m)(k − x1 − m) −
x1 y ax1 + 1
+ D1 (x2 − x1 − m),
x˙ 2 = −d1 x2 + D2 (x1 − x2 + m), y˙ = −d2 y +
c β x1 (t − τ )y ax1 (t − τ ) + 1
(2.1)
,
where p = αk . Then system (2.1) has an interior equilibrium E ∗ (x∗1 , x∗2 , y∗ ), where x∗1 = ∗
y =
d2 c β − ad2
,
p(ax∗1 + 1) x∗1
x∗2 =
D2 d1 + D2
(x∗1 + m),
(x1 + m) k − x1 − ∗
d1 D1
∗
p(D2 + d1 )
−m .
Obviously, x∗1 > 0 and x∗2 > 0 if and only if
(H1 ) c β − ad2 > 0 holds. Consider the quadratic function of m
R(m) = (x∗1 + m) k − x∗1 −
d1 D 1 p(D2 + d1 )
−m .
It is easy to see that if
(H2 ) 0 < m < k − x∗1 − then R(m) > 0, i.e. y∗ =
D 1 d1 p(d1 + D2 )
p(ax∗ +1 ) 1 R x∗ 1
(m) > 0. Hence, we have the following lemma.
Lemma 2.1. If (H1 ) and (H2 ) hold, then system (2.1) has a unique positive equilibrium E ∗ .
(2.2)
X. Liu, M. Han / Nonlinear Analysis: Real World Applications 12 (2011) 1047–1061
1049
In order to investigate the local stability of E ∗ of system (2.1), we linearize system (2.1) at E ∗ , obtaining a linear system of the form
[
x˙ 1 = p k − 2p(x1 + m) −
y∗
∗
(ax∗1 + 1)2
x∗1
] − D1 x1 + D1 x2 −
ax∗1 + 1
y,
x˙ 2 = D2 x1 − (d1 + D2 )x2 , c β y∗
y˙ =
(ax∗1 + 1)2
(2.3)
x1 (t − τ ).
The Jacobian matrix of system (2.3) for E ∗ is given by
A − D1 J (E ∗ ) = D2 C e−λτ
−B 0 ,
D1
−d 1 − D 2 0
(2.4)
0 x∗ 1
y∗
where A = p k − 2p(x∗1 + m) − (ax∗ +1)2 , B = 1 equation corresponding to (2.4)
+1 ax∗ 1
c β y∗ (ax∗1 +1)2
> 0 and C =
> 0. By det(λ I − J ) = 0, we get the characteristic
λ3 + (D2 + d1 + D1 − A)λ2 + [D1 d1 − A(d1 + D2 )]λ + BC (λ + D2 + d1 )e−λτ = 0.
(2.5)
Lemma 2.2. For τ = 0, if (H1 ) and
(H3 )
d1 D1
2
p(d1 + D2 )
1
max
k − 2x∗1 −
, 0 < m < k − x∗1 −
d1 D1 p(d1 + D2 )
are satisfied, then E ∗ of system (2.1) is locally asymptotically stable. Proof. For τ = 0, Eq. (2.5) becomes
λ3 + Q1 λ2 + Q2 λ + Q3 = 0,
(2.6)
where Q1 = d1 + D2 + D1 − A,
Q2 = D1 d1 + BC − A(d1 + D2 ),
Q3 = BC (d1 + D2 ).
First, we will prove that Q1 > 0 holds. Clearly, if (H3 ) holds then (H2 ) holds. By (H2 ), we have p(x∗ + m) D1 d1 Q1 = d1 + D2 + D1 − p k + 2p(x∗1 + m) + ∗ 1 ∗ −m k − x∗1 − x1 (ax1 + 1) p(d1 + D2 )
[
]
> d1 + D2 + D1 − p k + 2p(x∗1 + m). To prove Q1 > 0, it suffices to prove D1 − p k + 2p(x∗1 + m) > 0 which is equivalent to m > by
D1 p
> 1 2
d1 D1 , p(d1 +D2 )
1 2
(k − x∗1 −
D1 p
)−
x∗ 1 2
. Note that
we have
k − x∗1 −
D1
p
−
x∗1 2
<
1 2
k − x∗1 −
Considering (H3 ) there exists m > 12 (k − x∗1 − Further, let g (m) = Q1 Q2 − Q3 . Then
d1 D1 p(d1 + D2 ) D1 p
)−
x∗ 1 2
−
x∗1 2
< k − x∗1 −
d1 D 1 p(d1 + D2 )
.
. Hence, under condition (H3 ) there exists Q1 > 0.
g (m) = (d1 + D2 )(D1 d1 − A(d1 + D2 )) + (D1 − A)(D1 d1 + BC − A(d1 + D2 ))
[ ] D1 d1 − A + (D1 − A) (d1 + D2 ) − A + BC d1 + D 2 d1 + D 2 ∗ D d y D1 y∗ 1 1 ∗ = p(d1 + D2 )2 − k + 2(x∗1 + m) + + p − k + 2 ( x + m ) + 1 p(d1 + D2 ) p(ax∗1 + 1)2 p p(ax∗1 + 1)2 [ ] ∗ y D1 d1 × p(d1 + D2 ) − k + 2(x∗1 + m) + + BC . ∗ p(d1 + D2 ) p(ax1 + 1)2 = (d1 + D2 )2
D1 d1
It is easy to verify that if (H3 ) holds, it has D1 d1 p(d1 + D2 )
− k + 2(x∗1 + m) > 0,
which implies g (m) > 0.
D1 p
− k + 2(x∗1 + m) > 0,
1050
X. Liu, M. Han / Nonlinear Analysis: Real World Applications 12 (2011) 1047–1061
By Q1 > 0 and g (m) > 0, together with Q3 > 0, applying Routh–Hurwitz criterion, all the roots of Eq. (2.6) have negative real parts. This ends the proof. For τ > 0, we investigate the existence of purely imaginary roots of Eq. (2.5). Obviously, iω (ω > 0) is a root of Eq. (2.5) if and only if iω satisfies the following equation
− iω3 − (D2 + d1 + D1 − A)ω2 + i[D1 d1 − A(d1 + D2 )]ω + BC (iω + D2 + d1 )(cos ωτ − i sin ωτ ) = 0.
(2.7)
Separating the real and imaginary parts, we get
−ω3 + [D1 d1 − A(d1 + D2 )]ω = BC [(d1 + D2 ) sin ωτ − ω cos ωτ ],
(2.8)
(d1 + D2 + D1 − A)ω2 = BC [(d1 + D2 ) cos ωτ + ω sin ωτ ], which yields
ω6 + G1 ω4 + G2 ω2 + G3 = 0,
(2.9)
where G1 = (d1 + D2 )2 + (A − D1 )2 + 2D1 D2 > 0, G2 = [D1 d1 − A(d1 + D2 )]2 − B2 C 2 , G3 = −B2 C 2 (d1 + D2 )2 < 0. G1 > 0 and G3 < 0 together indicate that Eq. (2.9) has a unique positive root ω02 . Substituting ω02 into (2.8) and solving for τ , we obtain
ω02 [ω02 + (d1 + D2 )2 + D2 D1 ] 2jπ τj = arccos , + ω0 BC [ω0 2 + (d1 + D2 )2 ] ω0
1
j = 0, 1, 2, 3, . . . .
(2.10)
Theorem 2.1. Let (H1 ) and (H3 ) hold. Then (i) if 0 ⩽ τ < τ0 , the positive equilibrium E ∗ of system (2.1) is locally asymptotically stable; (ii) if τ > τ0 , the positive equilibrium E ∗ is unstable. Furthermore, system (2.1) undergoes Hopf bifurcation at E ∗ when τ = τj , j = 0, 1, 2, 3, . . .. Proof. Define by λ(τ ) = υ(τ ) + iω(τ ) the root of Eq. (2.5) such that
υ(τj ) = 0,
ω(τj ) = ω0 ,
j = 0, 1, 2, . . . .
Substituting λ(τ ) into (2.5) and taking the derivative with respect to τ , we obtain
[
dλ
] −1
[3λ2 + 2(D2 + d1 + D1 − A)λ + D1 d1 − A(d1 + D2 ) + BC e−λτ ] τ − BC λ(λ + D2 + d1 )e−λτ λ 2λ3 + (D2 + d1 + D1 − A)λ2 − BC (d1 + D2 )e−λτ τ = − BC λ2 (λ + D2 + d1 )e−λτ λ [ ] 1 2λ3 + (D2 + d1 + D1 − A)λ2 d1 + D2 τ =− 2 + − . 3 2 λ λ + (D2 + d1 + D1 − A)λ + [D1 d1 − A(d1 + D2 )]λ λ + d1 + D2 λ
=
dτ
When λ = iω0 , then
[ Re
dλ
2ω02 + [(d1 + D2 )2 + (A − D1 )2 + 2D1 D2 ]
]−1
dτ τ =τj
=
]
[(D1 d1 − A(d1 + D2 )) − ω ] + (d1 + D2 + D1 − A) ω 2 2 0
2
2 0
+
(d1 + D2 )2 > 0. ω [ω02 + (d1 + D2 )2 ] 2 0
Hence, sign Re
[
dλ
dτ τ =τj
= sign Re
[
dλ
] −1
dτ τ =τj
= 1.
(2.11)
Therefore, according to Lemmas 2.1 and 2.2, the transversality condition (2.11) and the Hopf bifurcation theorem for functional differential equations [10], we obtain the conclusions. 3. Direction and stability of Hopf bifurcation In Section 2, we have shown that system (2.1) undergoes Hopf bifurcations at τj , j = 0, 1, 2, . . .. Note that Theorem 2.1 only indicates the existence of the periodic solutions. In this section, by referring the works in [11–13] and applying the ideas of Hassard et al. [14], we derive some explicit formulae to determine the direction of the Hopf bifurcation and stability of the bifurcating periodic solutions of system (2.1) at τ = τj .
X. Liu, M. Han / Nonlinear Analysis: Real World Applications 12 (2011) 1047–1061
1051
Let x˜ 1 = x1 − x∗1 , x˜ 2 = x2 − x∗2 , x˜ 3 = y − y∗ and t˜ = τt , τ = τj + µ. Then system (2.1) can be written as a functional differential equation in C = C ([−1, 0], R3 ) x˙ (t ) = Lµ (xt ) + F (µ, xt ),
(3.1)
where x(t ) = (˜x1 (t ), x˜ 2 (t ), x˜ 3 (t )) ∈ R , Lµ : C → R and F : R × C → R are defined, respectively, by T
3
[p k − 2p(x∗1 + m) − y∗ q2 − D1 ]φ1 (0) + D1 φ2 (0) − x∗1 qφ3 (0) , D2 φ1 (0) − (d1 + D2 )φ2 (0) Lµ (φ) = (τj + µ) c β y∗ q2 φ1 (−1)
and
F (µ, φ) = (τj + µ)
− (−1)j aj−2 qj (ay∗ qφ1j (0) − φ1j−1 (0)φ3 (0)) j⩾2 0 , − j j −1 j −1 j −2 j ∗ cβ (−1) a q (ay qφ1 (−1) − φ1 (−1)φ3 (0))
−pφ12 (0) +
j⩾2
where q =
1 . ax∗ +1 1
Note that µ = τ − τj , thus, µ = 0 is a Hopf bifurcation value of system (3.1).
By the Riesz representation theorem, there exists a 3 × 3 matrix function ρ(θ , µ) whose components are bounded variation for θ ∈ [−1, 0], such that Lµ (φ) =
∫
0
dρ(θ , µ)φ(θ ), −1
here, we choose
p k − 2p(x∗1 + m) − y∗ q2 − D1 ρ(θ , µ) = (τj + µ) D2 0
D1
−d1 − D2 0
−x∗1 q 0 δ(θ ) − 0
0 0 c β y∗ q2
0 0 0
0 0 δ(θ + 1) , 0
where δ(θ ) is delta function. For φ ∈ C = C 1 ([−1, 0], R3 ), define
dφ(θ ) , dθ ∫ M (µ)φ = 0 dρ(v, µ)φ(v),
θ ∈ [−1, 0), R(µ)φ =
θ = 0,
0, F (µ, φ),
θ ∈ [−1, 0), θ = 0.
(3.2)
−1
Then system (3.1) can be written as x˙ t = M (µ)xt + R(µ)xt ,
(3.3)
where xt = x(t + θ ) for θ ∈ [−1, 0]. For φ ∈ C = C 1 ([−1, 0], R3 ), ϕ ∈ C ∗ = C 1 ([0, 1], R3∗ ), we define the adjoint operator M∗ of M (0) by
dϕ(s) , − ds M∗ ϕ(s) = ∫ 0 dρ T (t , 0)ϕ(−t ),
s ∈ (0, 1], s = 0,
−1
and the adjoint bilinear form by
⟨ϕ(s), φ(s)⟩ = ϕ( ¯ 0)φ(0) −
0
∫
∫
−1
θ
υ=0
ϕ(υ ¯ − θ )dρ(θ, 0)φ(υ)dυ.
(3.4)
In terms of the discussion of Section 2, we see that ± iω0 τj are eigenvalues of M (0) and M∗ . Next we calculate the eigenvector η(θ ) of M (0) corresponding to iω0 τj and the eigenvector η∗ (s) of M∗ belonging to −iω0 τj , respectively. By the definition of M (0) and M∗ , it is not difficult to obtain that T iω0 τj θ
η(θ ) = (1, b1 , b2 ) e
= 1,
D2 iω0 + d1 + D2
,
c β y∗ q2 e−iω0 τj
T
iω0
and
η∗ (s) = D(1, b∗1 , b∗2 )eiω0 τj s = D 1,
D2
,
x∗1 q
−iω0 + d1 + D2 iω0
eiω0 τj s .
eiω0 τj θ
1052
X. Liu, M. Han / Nonlinear Analysis: Real World Applications 12 (2011) 1047–1061
It follows from (3.4) that
∫
0
θ
∫
η¯ ∗ (ξ − θ )dρ(θ , 0)η(ξ )dξ ⟨η , η⟩ = η¯ (0)η(0) − −1 ξ =0 0 0 0 = D¯ 1 + b1 b¯ ∗1 + b2 b¯ ∗2 + τj (1, b¯ ∗1 , b¯ ∗2 ) 0 c β y∗ q2 0 ∗
∗
0 1 0 b1 b2 0
e−iω0 τj
= D¯ (1 + b1 b¯ ∗1 + b2 b¯ ∗2 + τj c β y∗ q2 b¯ ∗2 e−iω0 τj ). ¯ = 0. We choose D = (1 + b¯ 1 b∗1 + b¯ 2 b∗2 + τj c β y∗ q2 b∗2 eiω0 τj )−1 , such that ⟨η∗ , η⟩ = 1 and ⟨η∗ , η⟩ Following the algorithms in [14], we compute the coordinates describing center manifold C0 at µ = 0. Define z (t ) = ⟨η∗ , xt ⟩,
W (t , θ ) = xt (θ ) − z (t )η(θ ) − z (t )η(θ ),
(3.5)
therefore, xt (θ ) = (˜x1t (θ ), x˜ 2t (θ ), x˜ 3t (θ ))T
= W (t , θ ) + z (1, b1 , b2 )T eiθω0 τj + z¯ (1, b¯ 1 , b¯ 2 )T e−iθ ω0 τj .
(3.6)
On the center manifold C0 , we know that W (t , θ ) = W (z (t ), z¯ (t ), θ) = W20
z2 2
+ W11 z z¯ + W02
z¯ 2 2
+ ···,
(3.7)
where z and z¯ are local coordinates for the center manifold C0 in the direction of η∗ and η¯ ∗ . Note that W is real if xt is real. For the solution xt ∈ C0 of (3.3), when µ = 0, it follows from (3.2)–(3.5) that z˙ (t ) = ⟨η∗ , x˙ t ⟩ = iω0 τj z (t ) + ⟨η∗ (θ ), R(0)(W (z , z¯ , θ ) + 2Re{z (t )η(θ )})⟩
= iω0 τj z (t ) + η¯ ∗ (0)F (0)(W (z , z¯ , 0) + 2Re{z (t )η(0)}) , iω0 τj z (t ) + g (z , z¯ ),
(3.8)
where z2
z¯ 2
z 2 z¯
+ ···. (3.9) 2 2 2 From [14], we know that the properties of bifurcating periodic solutions at the critical value τj are determined by the following values g (z , z¯ ) = g20
c1 (0) =
+ g11 z z¯ + g02
i 2ω0 τj
+ g21
2
g20 g11 − 2|g11 | −
β2 = 2Re{c1 (0)},
T2 = −
|g02 |2 3
+
g21 2
,
µ2 = −
Im{c1 (0)} + µ2 Im{λ′ (τj )}
ω0 τ j
Re{c1 (0)} Re{λ′ (τj )}
,
.
Therefore, we need to compute the values of g20 , g11 , g02 and g21 . Firstly, from (3.6) and (3.7), we obtain (1)
z2
(3)
2 z2
x˜ 1t (0) = W20 (0) x˜ 3t (0) = W20 (0) (1)
2
(1) (1) + W11 (0)z z¯ + W02 (0) (3) (3) + W11 (0)z z¯ + W02 (0)
x˜ 1t (−1) = W20 (−1)
z2 2
z¯ 2 2 z¯ 2 2
+ z + z¯ + O(|(z , z¯ )|3 ); + b2 z + b¯ 2 z¯ + O(|(z , z¯ )|3 );
(1) (1) + W11 (−1)z z¯ + W02 (−1)
z¯ 2 2
+ ze−iω0 τj + z¯ eiω0 τj + O(|(z , z¯ )|3 ).
On the other hand, in terms of the definition of F (µ, xt ), we have g (z , z¯ ) = η¯ ∗ (0)F (0)(W (z , z¯ , 0) + 2Re{z (t )η(0)}) = η¯ ∗ (0)F0 (z , z¯ )
− j−1 −px˜ 21t (0) + (−1)j aj−2 qj (ay∗ qx˜ j1t (0) − x˜ 1t (0)˜x3t (0)) j⩾2 0 = τj D¯ (1, b¯ ∗1 , b¯ ∗2 )T − j −1 j−1 j−2 j ∗ j cβ (−1) a q (ay qx˜ 1t (−1) − x˜ 1t (−1)˜x3t (0)) j⩾2
= τj D¯ {z 2 (−p + aq3 y∗ − q2 b2 − b¯ ∗2 c β aq3 y∗ e−iω0 τj + b¯ ∗2 b2 c β q2 e−iω0 τj ) + 2z z¯ (−p + aq3 y∗ − Re{b2 }q2 − b¯ ∗2 c β aq3 y∗ + b¯ ∗2 c β q2 Re{b2 eiω0 τj })
(3.10)
X. Liu, M. Han / Nonlinear Analysis: Real World Applications 12 (2011) 1047–1061
1053
+ z¯ 2 (−p + aq3 y∗ − q2 b¯ 2 − b¯ ∗2 c β aq3 y∗ eiω0 τj + b¯ ∗2 b¯ 2 c β q2 eiω0 τj ) z 2 z¯ (3) (1) (1) (−1) [2(−p + ay∗ q3 )(W20 (0) + 2W11 (0)) − 2b¯ ∗2 c β ay∗ q3 eiω0 τj W20 + 2
(1) (1) (3) (3) − q2 (W20 (0)b¯ 2 + 2W11 (0)b2 + 2W11 (0) + W20 (0)) ( 3 ) ( 3 ) (1) (1) ∗ 2 i ω τ − i ω τ + b¯ 2 c β q (W20 (0)e 0 j + 2W11 (0)e 0 j + 2b2 W11 (−1) + b¯ 2 W20 (−1))
+ 2aq3 (b¯ 2 + 2b2 ) − 2b¯ ∗2 c β aq3 (2b2 + b¯ 2 e−2iω0 τj ) + 6a2 q4 y∗ (b¯ ∗2 c β e−iω0 τj − 1)] + · · ·}. By comparing coefficients with (3.9), we get
¯ (−p + aq3 y∗ − q2 b2 − b¯ ∗2 c β aq3 y∗ e−iω0 τj + b¯ ∗2 b2 c β q2 e−iω0 τj ); g20 = 2τj D ¯ (−p + aq3 y∗ − Re {b2 } q2 − b¯ ∗2 c β aq3 y∗ + b¯ ∗2 c β q2 Re{b2 eiω0 τj }); g11 = 2τj D ¯ (−p + aq3 y∗ − q2 b¯ 2 − b¯ ∗2 c β aq3 y∗ eiω0 τj + b¯ ∗2 b¯ 2 c β q2 eiω0 τj ); g02 = 2τj D (1)
(3)
(1)
¯ [2(−p + ay∗ q3 )(W20 (0) + 2W11 (0)) − 2b¯ ∗2 c β ay∗ q3 eiω0 τj W20 (−1) g21 = τj D (1)
(1)
(3)
(3)
− q (W20 (0)b¯ 2 + 2W11 (0)b2 + 2W11 (0) + W20 (0)) (3) (3) (1) (1) + b¯ ∗2 c β q2 (W20 (0)eiω0 τj + 2W11 (0)e−iω0 τj + 2b2 W11 (−1) + b¯ 2 W20 (−1)) + 2aq3 (b¯ 2 + 2b2 ) − 2b¯ ∗2 c β aq3 (2b2 + b¯ 2 e−2iω0 τj ) + 6a2 q4 y∗ (b¯ ∗2 c β e−iω0 τj − 1)]. 2
Obviously, the value of g21 is determined by W20 (θ ) and W11 (θ ). From (3.3) and (3.5), we have
˙ (t , θ ) = x˙ t (θ ) − 2Re{˙z (t )η(θ )} W = M (0)W − 2Re{g (z , z¯ )η(θ )} + R(0)(W + 2Re{z (t )η(θ )}) , M (0)W + H (z , z¯ , θ ),
(3.11)
where H (z , z¯ , θ ) = −2Re{g (z , z¯ )η(θ )} + R(0)(W + 2Re{z (t )η(θ )}) , H20 (θ )
z2
+ H11 (θ )z z¯ + H02 (θ )
z¯ 2
2 2 It follows from (3.12) and the definition of R(0) that H (z , z¯ , θ ) =
−2Re{g (z , z¯ )η(θ )}, −2Re{g (z , z¯ )η(0)} + F0 (z , z¯ ),
= H20 (θ )
z2
+ H11 (θ )z z¯ + H02 (θ )
z¯ 2
2 2 On the center manifold C0 near the origin, we have
+ ···.
(3.12)
θ ∈ [−1, 0) θ =0 + ···,
(3.13)
˙ = Wz z˙ + Wz¯ z˙¯ = W20 z z˙ + W11 z¯ z˙ + W11 z z˙¯ + W02 z¯ z˙¯ + · · · . W
(3.14)
Substituting z˙ = iω0 τj z (t ) + g (z , z¯ ) into (3.14) and comparing the coefficients with (3.11), we obtain
(M (0) − 2iω0 τj )W20 (θ ) = −H20 (θ ),
M (0)W11 (θ ) = −H11 (θ ).
(3.15)
It follows from (3.9) and (3.13) that for θ ∈ [−1, 0), H20 (θ ) = −g20 η(θ ) − g¯20 η(θ ¯ ),
H11 (θ ) = −g11 η(θ ) − g¯11 η(θ ¯ ).
Substituting H20 (θ ) and H11 (θ ) into (3.15) and using the definition of M (0) for θ ∈ [−1, 0), we obtain
˙ 20 (θ ) = 2iω0 τj W20 (θ ) + g20 η(θ ) + g¯20 η(θ W ¯ ),
˙ 11 (θ ) = g11 η(θ ) + g¯11 η(θ W ¯ ).
A straightforward computation yields W20 (θ ) =
ig20 η(0) iω0 τj θ ig¯20 η( ¯ 0) −iω0 τj θ e + e + S1 e2iω0 τj θ , ω0 τj 3ω0 τj
W11 (θ ) = −
ig11 η(0) iω0 τj θ ig¯11 η( ¯ 0) −iω0 τj θ e + e + S2 ,
ω0 τj
(3.16)
ω0 τj
3
where S1 , S2 ∈ R are undetermined constant vectors. In the following, we will determine the values of S1 and S2 . 0 Note that M (0)W20 = −1 dρ(v, 0)W20 (v) for θ = 0, then by (3.15), we have
∫
0
dη(v, 0)W20 (v) = 2iω0 τj W20 (0) − H20 (0), −1
∫
0
dη(v, 0)W20 (v) = −H11 (0). −1
(3.17)
1054
X. Liu, M. Han / Nonlinear Analysis: Real World Applications 12 (2011) 1047–1061
It follows from (3.13) that for θ = 0,
−p + aq3 y∗ − q2 b2
, 0 −c β aq3 y∗ e−iω0 τj + b2 c β q2 e−iω0 τj −p + aq3 y∗ − Re {b2 } q2 0 H11 (0) = −g11 η(0) − g¯11 η( ¯ 0) + 2τj iω τ . 3 ∗ 2 0 j −c β aq y + c β q Re b2 e
H20 (0) = −g20 η(0) − g¯20 η( ¯ 0) + 2τj
(3.18)
Substituting H20 (0) defined by (3.18) and W20 (θ ) defined by (3.16) in (3.17) respectively, and noting that
iω0 τj I −
∫
0
e
iω0 τj v
dρ(v, 0) η(0) = 0,
0
∫ −iω0 τj I −
e
−iω0 τj v
¯ 0) = 0, dρ(v, 0) η(
−1
−1
we derive
2iω0 τj I −
∫
0
e2iω0 τj v dρ(v, 0) S1 = 2τj
−1
−p + aq3 y∗ − q2 b2 0
−c β aq3 y∗ e−iω0 τj + b2 c β q2 e−iω0 τj
.
Solving the above linear equations for S1 , we may obtain 2iω0 τj − p k − 2p(x∗1 + m) + y∗ q2 + D1 −D 2 S1 = 2 −c β y∗ q2 e−2iω0 τj
×
3 ∗
0
−1
x∗1 q 0 2iω0 τj
2
−p + aq y − q b2 0
−D1 2iω0 τj + d1 + D2
−c β aq3 y∗ e−iω0 τj + b2 c β q2 e−iω0 τj
.
Similarly, from (3.16)–(3.18), we may get S2 which is defined by D1 − p k − 2p(x∗1 + m) + y∗ q2 −D2 S2 = 2 −c β y∗ q2
−D1 d1 + D 2 0
−1
x∗1 q 0 0
aq3 y∗ − p − Re {b2 } q2 . 0 c β q2 Re b2 eiω0 τj − c β aq3 y∗
Substituting S1 and S2 into (3.16), we may calculate W20 (θ ) and W11 (θ ). Furthermore, g21 can be determined. Thus, by (3.10), we can compute the values of c1 (0), µ2 , β2 and T2 . It is clear that µ2 determines the direction of the Hopf bifurcation: if µ2 > 0 (µ2 < 0), then the Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solutions exist for τ > τj (τ < τj ); β2 determines the stability of bifurcating periodic solutions: if β2 < 0 (β2 > 0) then the bifurcating periodic solutions are orbitally asymptotically stable (unstable); and T2 determines the period of the bifurcating periodic solutions: the period increases (decreases) if T2 > 0 (T2 < 0). To illustrate the feasibility of our main results, we give an example and carry out numerical simulations by using Matlab. Example. In (2.1), let α = 1.8, β = 0.8, k = 2.6, a = 1.2, c = 0.9, d1 = 0.1, d2 = 0.2, D1 = 0.8, D2 = 0.3 and D d m = 0.3. Clearly, c β − ad2 = 0.48 > 0 and 0 < m < k − x∗1 − p(d 1+1D ) = 1.8944. Then conditions (H1 ) and (H3 ) in 1 2 Theorem 2.1 are satisfied. System (2.1) has a unique positive equilibrium E ∗ (0.4167, 0.5375, 2.8479). Solving (2.9) for ω, we ∗ obtain ω0 ≈ 0.2565334393, then τ0 ≈ 3.916807176. By Theorem 2.1-(i) E is locally asymptotically stable for 0 ⩽ τ < τ0 . One can see Fig. 1(a)–(d). Furthermore, when the delay τ > τ0 , the positive equilibrium E ∗ loses its stability and the system goes into oscillations. The system undergoes a Hopf bifurcation at τ = τ0 . By computation we obtain g20 = 3.470469910 + 2.439350744i; g11 = 9.520322552 + 8.634501512i; g02 = 11.83851897 + 14.74884396i; g21 = −36.6992317 + 326.174705i; c1 (0) = −44.817374 − 54.6846853i;
µ2 = 806.3833986 > 0; β2 = −89.6347475 < 0; T2 = 61.45127048 > 0. Then we can conclude that the Hopf bifurcation is supercritical. Moreover, the bifurcating periodic solution for τ > τ0 is asymptotically stable (see Fig. 2(a)–(d)).
X. Liu, M. Han / Nonlinear Analysis: Real World Applications 12 (2011) 1047–1061
1055
Fig. 1. The time course of a solution tending to the equilibrium (0.4167, 0.5375, 2.8479) of system (2.1) with τ = 3.4 < τ0 , where α = 1.8, β = 0.8, k = 2.6, a = 1.2, c = 0.9, d1 = 0.1, d2 = 0.2, D1 = 0.8, D2 = 0.3, m = 0.3 and initial data x1 (0) = 0.6, x2 (0) = 0.6, y(0) = 2.76. (a) Time-state t − x1 (t ) relation graph; (b) Phase portrait of t − x2 (t ) relation graph; (c) Time-state t − y(t ) relation graph; (d) Phase portrait of x1 (t ), x2 (t ) and y(t ).
4. Numerical analysis of system (1.2) with impulsive In this section, by using numerical method, we will discuss the effect of impulsive perturbations on system (1.2). As we know, impulsive differential equations can be seen in almost every domain of applied science and studied extensively (see, for example [15–18]). In general, the impulsive perturbations will bring to the system complexity because of their propensity for chaos. Motivated by Meng et al. [17], Tang and Chen [18], by introducing periodic impulsive proportional harvesting (releasing) on the prey x2 and predator y, respectively, we develop model (1.2) as follows
β(x1 − m)y + D1 (x2 − x1 ), k 1 + a(x1 − m) x˙ 2 = −d1 x2 + D2 (x1 − x2 ), c β[x1 (t − τ ) − m]y , y˙ = −d2 y + 1 + a[x1 (t− τ ) − m] x1 (t + ) = x1 (t ), x2 (t + ) = (1 − p)x2 (t ), t = nT , n ∈ N , y(t + ) = (1 − δ)y(t ),
x˙ 1 = α x1 1 −
x1
−
t ̸= nT , n ∈ N , (4.1)
where 0 < p < 1 and 0 < δ < 1 (p < 0 and δ < 0) represent partial impulsive harvest (release) to preys x2 and predators, respectively, at t = nT , n ∈ N. T is the period of impulsive effect.
1056
X. Liu, M. Han / Nonlinear Analysis: Real World Applications 12 (2011) 1047–1061
Fig. 2. The trajectory graph of system (2.1) with the same parameters as in Fig. 1, when τ = 3.98 > τ0 , it shows the bifurcating periodic solution from (0.4167, 0.5375, 2.8479) which is asymptotically stable. (a) Time-state t − x1 (t ) relation graph; (b) Time-state t − x2 (t ) relation graph; (c) Phase portrait of t − y(t ); (d) Phase portrait of x1 (t ), x2 (t ) and y(t ).
Firstly, we investigate the effect of period of impulsive control T on dynamical behavior of system (1.2). For convenience, we keep the parameters α = 1.8, β = 0.8, k = 2.6, a = 1.2, c = 0.9, d1 = 0.1, d2 = 0.2, D1 = 0.8, D2 = 0.3, m = 0.3 fixed, and let p = 0.2, δ = 0.32 and τ = 6. Then the influence of T may be documented by stroboscopically sampling some of the variables over a range of T . Here, we carry out numerical simulations with respect to T in the range 0.5 ⩽ T ⩽ 44. And the resulting bifurcation diagrams provide a description of essential dynamical behavior of the system. With the increasing of T , the system experiences process of cycles, period-doubling bifurcation, period-halving bifurcation, intermittency (intermittency occurs whenever the behavior of a system seems to switch back and forth between two qualitatively different behaviors as a control parameter is changed), chaotic band, narrow periodic window, wide periodic window and non-unique dynamics, i.e., several attractors may coexist with the same T , and every attractor is corresponding to the initial values. In order to investigate the existence of the predator-free solution, we give Fig. 3, firstly. And it shows that the predator population goes into extinction if T < Tmin ≈ 1.8, where Tmin can be regarded as the minimum period of impulsive control to prevent the predator from extinction. From the bifurcation diagrams Fig. 4, we observe that a period-doubling bifurcation occurs when T ≈ 10.66. It means that T -periodic solution disappears suddenly at this point and 2T -periodic solution appears. Then the behavior of 2T -periodicity is kept for 10.66 < T < 16.6256. When T closes to 16.626, the 2T -periodic solution becomes unstable and shrinks to a chaotic attractor, rapidly (see Fig. 6). Following, system (4.1) comes into a chaotic region with some narrow periodic
X. Liu, M. Han / Nonlinear Analysis: Real World Applications 12 (2011) 1047–1061
1057
Fig. 3. (a) The bifurcation diagrams of the predator y for T over [0.5, 5]; (b) Time series of the prey x1 when T = 1.8; (c) Time series of the prey x2 when T = 1.8; (d) Time series of the predator y when T = 1.8.
10 9.5 9 8.5 8 7.5 7 6.5 6 5.5 5 4.5 4 3.5 3
b
Max(y)
Max(y)
a
5
6
7
8
9 10 11 12 13 14 15 16 17 18 19 20
10 9.5 9 8.5 8 7.5 7 6.5 6 5.5 5 4.5 4 3.5 3 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
T
c
T
d
9 8.5 8 7
Max(y)
Max(y)
7.5 6.5 6 5.5 5 4.5 4 35
35.5
36
36.5
37
37.5 T
38
38.5
39
39.5
40
10 9.5 9 8.5 8 7.5 7 6.5 6 5.5 5 4.5 4 3.5 3 40
40.4 40.8 41.2 41.6
42
42.4 42.8 43.2 43.6
44
T
Fig. 4. Bifurcation diagrams of system (5.1) showing the influence of T with α = 1.8, β = 0.8, k = 2.6, a = 1.2, c = 0.9, d1 = 0.1, d2 = 0.2, D1 = 0.8, D2 = 0.3, m = 0.3, p = 0.2 and δ = 0.32 when varying T from 5 to 44. (a) The predator y is plotted for T over interval [5, 20]. (b) The predator y is plotted for T over interval [20, 35]. (c) The predator y is plotted for T over interval [35, 40]. (d) The predator y is plotted for T over interval [40, 44].
windows. When T ⩾ 19.25 and nearby 19.25, the chaos bands disappear and a 3T -periodic solution occurs. And when T = 20.84, a non-unique attractor appears, i.e., a 3T -periodic solution coexists with a T -periodic solution (see Fig. 7). when T ≈ 31.49, the system displays another non-unique attractor (two different types of T -periodic solutions coexist). When varying T in the range [31.5, 36.91], there is a transition from periodicity to quasi-periodicity to chaos to periodicity again (see Fig. 5(b)). As T increases further, the evidence for 9T -periodic solutions leading to chaos can be observed for 36.91 < T < 37.25. And the chaotic behavior is kept until T closes to the value 37.38, then the system displays a 4T -periodicity behavior. When 37.38 ⩽ T ⩽ 38.798, there is 4T -periodic solution with twice period-doubling bifurcations and period-halving bifurcations. When T = 38.798, the non-unique attractor appears again (a 4T -periodic solution and a
1058
X. Liu, M. Han / Nonlinear Analysis: Real World Applications 12 (2011) 1047–1061
a
b
8.5
8.5
8
8
7.5
7.5
7
7 6.5
Max(y)
Max(y)
6.5 6 5.5
6 5.5
5 4.5
5
4
4.5
3.5
16.5
4 36.2
16.8 17 17.2 17.4 17.6 17.8 18 18.2 18.4 18.6 18.8 19
36.3
36.4
36.5
36.6
T
36.7
36.8
36.9
37
T
Fig. 5. The magnified parts of Fig. 4. (a) The predator y is plotted for T over interval [16.5, 19]. (b) The predator y is plotted for T over interval [36.2, 37].
a
7
b
T=16.625
6
6
T=16.627
5 4
4
y
y
5
3
3
2
2
1 0
1 0.5
0.5 x1
1 1.5
0.4
2 0.2
0.8
0.6
0.75
1.2
1
x1 1
x2
1.25 1.5 0.2
0.4
0.6
0.8
1
1.2
x2
Fig. 6. The dynamical behaviors of system (5.1) corresponding to T = 16.625 and T = 16.627, respectively. (a) A 2T -periodic solution corresponding to T = 16.625. (b) A chaotic attractor corresponding to T = 16.627. 8
b
T=20.84
6
6
5
4
4
y
y
a
T=20.84
3
2
2 0
0 0 0.5 x1
1 1.5 0.2
0.4
0.6
0.8
1
1.2
0.5 x1
x2
1 1.5 0.4
0.5
0.6
0.7 x2
0.8
0.9
Fig. 7. Coexistence of 3T -periodic solution with a T -periodic solution when the period T = 20.84. (a) With initial value x1 (0) = 0.6, x2 (0) = 0.6, y(0) = 2.76, the solution will finally tend to a 3T -periodic solution. (b) With initial value x1 (0) = 0.4, x2 (0) = 0.1, y(0) = 0.1, the solution will finally tend to a T -periodic solution.
3T -periodic solution may coexists; see Fig. 8). When T nearby 38.9, the behavior of period-doubling bifurcation and periodhalving bifurcation is given second. When T ⩾ 39.78, the 3T -periodic solution occurs again, and we note that until T = 42.3 no chaos occurs. Following the chaotic area, the behavior of cascade of period-halving bifurcations from chaos to cycles may also be seen for 42.85 < T < 43.6. In addition, the phenomenon of intermittency also can be illustrated by Fig. 5(a) which shows a series of switch from chaotic band to periodic window to chaotic window back and forth. On the other hand, in order to further investigate the dynamical behaviors of system (4.1), the bifurcation diagrams with respect to the parameters δ , p and m are also carried out, respectively. Now, consider the following set of parameters
α = 1.8, D1 = 0.8,
β = 0.8, D2 = 0.3,
k = 2.6, m = 0.3,
a = 1.2,
τ = 6,
c = 0.9,
d1 = 0.1,
d2 = 0.2,
T = 18.
Let p = 0.1, then the bifurcation diagrams with respect to parameter δ over the range [−0.26, 0.4] are shown in Fig. 9. From these bifurcation diagrams, we can see that system (4.1) displays complex dynamics with increasing δ from −0.26 to 0.4, gradually. For example, some strange attractors are captured when δ = −0.1, −0.075, 0.315, and 0.35, respectively (see Figs. 10 and 11).
X. Liu, M. Han / Nonlinear Analysis: Real World Applications 12 (2011) 1047–1061
a
10
b
T=38.798
8
1059
T=38.798
7
8
6 6 y
y
5 4
4
3 2
2
0 0
1 0 0.5 1
x1
0.4
1.5 0.2
0.8 x2
0.6
1
0.5
1.2
1
x1
0.6
0.4
1.5 0.2
0.8
1
1.2
x2
Fig. 8. Coexistence of 4T -period solution with a 3T -periodic solution when the period T = 38.798. (a) With initial value x1 (0) = 0.4, x2 (0) = 0.1, y(0) = 0.1, the solution will finally tend to a 4T -periodic solution. (b) With initial value x1 (0) = 0.6, x2 (0) = 0.6, y(0) = 2.76, the solution will finally tend to a 3T -periodic solution.
a
b
6.5
7 6.5
5.5
6
5
5.5
Max(y)
Max(y)
6
4.5
5
4
4.5
3.5
4
3
3.5
-0.26 -0.235 -0.21 -0.185 -0.16 -0.135 -0.11 -0.085 -0.06 -0.035 -0.010
0
0.05
0.1
0.15
0.2
δ
0.25
0.3
0.35
0.4
δ
Fig. 9. The bifurcation diagrams of system (5.1) with respect to the parameter δ in the range −0.26 ⩽ δ ⩽ 0.4. (a) −0.26 ⩽ δ ⩽ 0; (b) 0 ⩽ δ ⩽ 0.4.
a
6
b
δ=-0.11
5.5
6
δ=-0.075
5.5
y
5 4.5
y
5 4.5 4
4
3.5
3.5
3
3
2.5 0
2.5 0 0.5 x1
1 1.5
0.4
0.5
0.6 x2
0.7
0.8
0.5 1
x1
1.5
0.5
0.4
0.6 x2
0.7
0.8
Fig. 10. The dynamical behaviors of system (5.1) when α = 1.8, β = 0.8, k = 2.6, a = 1.2, c = 0.9, d1 = 0.1, d2 = 0.2, D1 = 0.8, D2 = 0.3, m = 0.3, p = 0.1 and T = 18. (a) A strange attractor corresponding to δ = −0.11. (b) A strange attractor corresponding to δ = −0.075.
a
b
6 δ=0.315
5.5 δ=0.35
5
5
y
y
4.5 4
4 3.5 3 2.5
3 2 0.5 1 x1
1.5 0.2
0.4
0.6 x2
0.8
1
2 0.5 1 x1
1.5 0.4
0.5
0.6
0.7
0.8
0.9
x2
Fig. 11. The dynamical behaviors of system (5.1) when α = 1.8, β = 0.8, k = 2.6, a = 1.2, c = 0.9, d1 = 0.1, d2 = 0.2, D1 = 0.8, D2 = 0.3, m = 0.3, p = 0.1 and T = 18. (a) A strange attractor corresponding to δ = 0.315. (b) A strange attractor corresponding to δ = 0.35.
1060
X. Liu, M. Han / Nonlinear Analysis: Real World Applications 12 (2011) 1047–1061
a
b
6 5.8
6.5 6
5.6 5.5 Max(y)
Max(y)
5.4 5.2 5 4.8
5 4.5 4
4.6 3.5
4.4
3
4.2 p
p
Fig. 12. The bifurcation diagrams of system (5.1) with respect to the parameter p in the range −0.8 ⩽ p ⩽ 0.85. (a) −0.8 ⩽ p ⩽ 0; (b) 0 ⩽ p ⩽ 0.85.
16 15 14 13 12 11 10 9 8 7 6 5 4 3 2
b
Max(y)
Max(y)
a
0
0.05
0.1
0.15
0.2 m
0.25
0.3
0.35
0.4
12 11.5 11 10.5 10 9.5 9 8.5 8 7.5 7 6.5 6 5.5 5 4.5 4 3.5 3
0.15
0.16
0.17
0.18
0.19
0.2
0.21
0.225
m
Fig. 13. The bifurcation diagrams of system (5.1) with respect to the parameter p in the range 0 ⩽ m ⩽ 0.4. (a) 0 ⩽ m ⩽ 0.4; (b) 0.15 ⩽ m ⩽ 0.225.
Similarly, if we keep the values of α , β , k, a, c, d1 , d2 , D1 , D2 , m, T and τ fixed, and let δ = 0.3. For −0.8 ⩽ p ⩽ 0.85 the effect of p on the dynamics of system (4.1) can be shown in Fig. 12. From the discussion of Section 2, it is clear that the behavior of prey refuge also has an important influence on the system. For this reason, we plot the bifurcation diagrams of system (4.1) with respect to m in the range 0 ⩽ m ⩽ 0.4 (see Fig. 13). By choosing T , δ , p and m as bifurcation parameters, respectively, the rich and complex dynamical behaviors of system (4.1) have been observed. All these results show that the dynamical behaviors of system (4.1) will change evidently with varying the values of some parameters. References [1] E. Gonlez-Olivares, R. Ramos-Jiliberto, Dynamic consequences of prey refuges in a simple model system: more prey, fewer predators and enhanced stability, Ecological Modelling 166 (2003) 135–146. [2] L. Chen, F. Chen, L. Chen, Qualitative analysis of a predator–prey model with Holling type II functional response incorporating a constant prey refuge, Nonlinear Analysis: Real World Applications 11 (2010) 246–252. [3] T.K. Kar, Modelling and analysis of a harvested prey–predator system incorporating a prey refuge, Journal of Computational and Applied Mathematics 185 (2006) 19–33. [4] L. Ji, C. Wu, Qualitative analysis of a predator–prey model with constant-rate prey harvesting incorparating a constant prey refuge, Nonlinear Analysis: Real World Applications 11 (2010) 2285–2295. [5] Z. Ma, et al., Effects of prey refuges on a predator–prey model with a class of functional responses: the role of refuges, Mathematical Biosciences 218 (2009) 73–79. [6] Y. Kuang, Y. Takeuchi, Predator–prey dynamics in models of prey dispersal in two-patch environments, Mathematical Biosciences 120 (1994) 77–98. [7] J. Cui, L. Chen, The effect of diffusion on the time varying logistic population growth, Computers & Mathematics with Applications 36 (1998) 1–9. [8] J. Cui, L. Chen, Permanence and extinction in logistic and Lotka–Volterra systems with diffusion, Journal of Mathematical Analysis and Applications 258 (2001) 512–535. [9] R. Xu, L. Chen, Persistence and stability for a two-species ratio-dependent predator–prey system with time delay in a two-patch environment, Computers and Mathematics with Applications 40 (2000) 577–588. [10] J. Hale, S.V. Lunel, Introduction to Functional Differential Equations, Springer-Verlag, New York, 1993. [11] F. Lian, Y. Xu, Hopf bifurcation analysis of a predator prey system with Holling type IV functional response and time delay, Applied Mathematics and Computation 215 (2009) 1484–1495. [12] Y. Song, S. Yuan, Bifurcations analysis in a predator–prey system with time delay, Nonlinear Analysis: Real World Applications 7 (2006) 265–284. [13] R. Xu, Q. Gan, Z. Ma, Stability and bifurcation analysis on a ratio-dependent predator–prey model with time delay, Journal of Computational and Applied Mathematics 230 (2009) 187–203. [14] B.D. Hassard, N.D. Kazarinoff, Y.H. Wan, Theory and Applications of Hopf Bifurcation, Cambridge University Press, Cambridge, 1981. [15] S. Zhang, L. Chen, The study of predator–prey system with defensive ability of prey and impulsive perturbations on the predator, Chaos, Solitons & Fractals 23 (2005) 311–320.
X. Liu, M. Han / Nonlinear Analysis: Real World Applications 12 (2011) 1047–1061
1061
[16] Y. Xia, Positive periodic solutions for a neutral impulsive delayed Lotka–Volterra competition system with the effect of toxic substance, Nonlinear Analysis: Real World Applications 8 (2007) 204–221. [17] X. Meng, L. Chen, H. Cheng, Two profitless delays for the SEIRS epidemic disease model with nonlinear incidence and pulse vaccination, Applied Mathematics and Computation 186 (2007) 516–529. [18] S. Tang, L. Chen, Density-dependent birth rate, birth pulse and their population dynamic consequences, Journal of Mathematical Biology 44 (2002) 185–199.
Contents lists available at ScienceDirect
Nonlinear Analysis: Real World Applications journal homepage: www.elsevier.com/locate/nonrwa
Chaos and Hopf bifurcation analysis for a two species predator–prey system with prey refuge and diffusion✩ Xia Liu a , Maoan Han b,∗ a
College of Mathematics and Information Science, Henan Normal University, 453007, PR China
b
Department of Mathematics, Shanghai Normal University, Shanghai 200234, PR China
article
abstract
info
Article history: Received 8 May 2010 Accepted 23 August 2010
In this paper, a delayed predator–prey model with Holling type II functional response incorporating a constant prey refuge and diffusion is considered. By analyzing the characteristic equation of linearized system corresponding to the model, we study the local asymptotic stability of the positive equilibrium of the system. By choosing the time delay due to gestation as a bifurcation parameter, the existence of Hopf bifurcations at the positive equilibrium is established. By applying the normal form and the center manifold theory, an explicit algorithm to determine the direction of the Hopf bifurcation and the stability of the bifurcating periodic solutions is derived. Further, an example is presented to illustrate our main results. Finally, recurring to the numerical method, the influences of impulsive perturbations on the dynamics of the system are also investigated. © 2010 Elsevier Ltd. All rights reserved.
Keywords: Diffusion Hopf bifurcation Chaos Prey refuge
1. Introduction The predator–prey systems with prey refuge have been studied by several authors; one can refer to [1–5] and the references cited therein. Especially, González-Olivares et al. [1] proposed a predator–prey system with Holling type II functional response incorporating a constant prey refuge as follows
x˙ = α x 1 − y˙ = −d2 y +
x
−
β(x − m)y
k 1 + a(x − m) c β(x − m)y 1 + a( x − m )
, (1.1)
.
They investigated the local stability of equilibria and the existence of limit cycle of system (1.1). Motivated by the work of González-Olivares et al. [1], Chen et al. [2] further discussed the instability and global stability properties of the equilibria and the existence and uniqueness of limit cycle of system (1.1). Many papers have researched the predator–prey systems with diffusion [6–8]. As the authors in [9] noted, the dispersal between patches often occurs in ecological environments, and more realistic models should include the dispersal process.
✩ This work was supported by Shanghai Leading Academic Discipline Project (No. S30405).
∗
Corresponding author. E-mail address: [email protected] (M. Han).
1468-1218/$ – see front matter © 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.nonrwa.2010.08.027
1048
X. Liu, M. Han / Nonlinear Analysis: Real World Applications 12 (2011) 1047–1061
Based on the above discussion, by incorporating dispersal of prey and time delay into system (1.1), we establish the following model x1
x˙ 1 = α x1 1 −
k
−
β(x1 − m)y 1 + a(x1 − m)
+ D1 (x2 − x1 ),
x˙ 2 = −d1 x2 + D2 (x1 − x2 ), y˙ = −d2 y +
(1.2)
c β[x1 (t − τ ) − m]y 1 + a[ x 1 ( t − τ ) − m ]
,
with initial conditions x1 (θ ) = ϕ1 (θ ) ⩾ m,
x2 (θ ) = ϕ2 (θ ) ⩾ 0 ϕ2 (0) > 0, ψ(0) > 0,
ϕ1 (0) > m,
y(θ ) = ψ(θ ) ⩾ 0,
θ ∈ [−τ , 0],
(1.3)
where xi is the prey density in the ith patch, i = 1, 2, and y is the predator density; Di > 0 represents the dispersal rate of prey species, i = 1, 2; d1 > 0 denotes the maximal death rate of x2 in second patch; τ is a constant time delay due to gestation; the parameters α , k, d2 and c are positive constants and their ecological meaning can be found in Refs. [1,2]. The rest of the paper is organized as follows. The existence and local property of the positive equilibrium and the existence of local Hopf bifurcation of the model are discussed in Section 2. In Section 3, by applying the normal form theory and the center manifold reduction for functional differential equations, the conditions for determining the bifurcation direction and the stability of the bifurcating periodic solution are derived. We also give a numerical example to illustrate the validity of our results. Finally, in Section 4, by using the numerical method, we investigate the influences of impulsive perturbations on system (1.2). 2. Stability of the positive equilibrium and the existence of local Hopf bifurcation For simplicity, we take x¯ 1 = x1 − m, x¯ 2 = x2 and y¯ = β y so that system (1.2) can be rewritten as (still denote x¯ 1 , x¯ 2 and y¯ as x1 , x2 and y, respectively) x˙ 1 = p(x1 + m)(k − x1 − m) −
x1 y ax1 + 1
+ D1 (x2 − x1 − m),
x˙ 2 = −d1 x2 + D2 (x1 − x2 + m), y˙ = −d2 y +
c β x1 (t − τ )y ax1 (t − τ ) + 1
(2.1)
,
where p = αk . Then system (2.1) has an interior equilibrium E ∗ (x∗1 , x∗2 , y∗ ), where x∗1 = ∗
y =
d2 c β − ad2
,
p(ax∗1 + 1) x∗1
x∗2 =
D2 d1 + D2
(x∗1 + m),
(x1 + m) k − x1 − ∗
d1 D1
∗
p(D2 + d1 )
−m .
Obviously, x∗1 > 0 and x∗2 > 0 if and only if
(H1 ) c β − ad2 > 0 holds. Consider the quadratic function of m
R(m) = (x∗1 + m) k − x∗1 −
d1 D 1 p(D2 + d1 )
−m .
It is easy to see that if
(H2 ) 0 < m < k − x∗1 − then R(m) > 0, i.e. y∗ =
D 1 d1 p(d1 + D2 )
p(ax∗ +1 ) 1 R x∗ 1
(m) > 0. Hence, we have the following lemma.
Lemma 2.1. If (H1 ) and (H2 ) hold, then system (2.1) has a unique positive equilibrium E ∗ .
(2.2)
X. Liu, M. Han / Nonlinear Analysis: Real World Applications 12 (2011) 1047–1061
1049
In order to investigate the local stability of E ∗ of system (2.1), we linearize system (2.1) at E ∗ , obtaining a linear system of the form
[
x˙ 1 = p k − 2p(x1 + m) −
y∗
∗
(ax∗1 + 1)2
x∗1
] − D1 x1 + D1 x2 −
ax∗1 + 1
y,
x˙ 2 = D2 x1 − (d1 + D2 )x2 , c β y∗
y˙ =
(ax∗1 + 1)2
(2.3)
x1 (t − τ ).
The Jacobian matrix of system (2.3) for E ∗ is given by
A − D1 J (E ∗ ) = D2 C e−λτ
−B 0 ,
D1
−d 1 − D 2 0
(2.4)
0 x∗ 1
y∗
where A = p k − 2p(x∗1 + m) − (ax∗ +1)2 , B = 1 equation corresponding to (2.4)
+1 ax∗ 1
c β y∗ (ax∗1 +1)2
> 0 and C =
> 0. By det(λ I − J ) = 0, we get the characteristic
λ3 + (D2 + d1 + D1 − A)λ2 + [D1 d1 − A(d1 + D2 )]λ + BC (λ + D2 + d1 )e−λτ = 0.
(2.5)
Lemma 2.2. For τ = 0, if (H1 ) and
(H3 )
d1 D1
2
p(d1 + D2 )
1
max
k − 2x∗1 −
, 0 < m < k − x∗1 −
d1 D1 p(d1 + D2 )
are satisfied, then E ∗ of system (2.1) is locally asymptotically stable. Proof. For τ = 0, Eq. (2.5) becomes
λ3 + Q1 λ2 + Q2 λ + Q3 = 0,
(2.6)
where Q1 = d1 + D2 + D1 − A,
Q2 = D1 d1 + BC − A(d1 + D2 ),
Q3 = BC (d1 + D2 ).
First, we will prove that Q1 > 0 holds. Clearly, if (H3 ) holds then (H2 ) holds. By (H2 ), we have p(x∗ + m) D1 d1 Q1 = d1 + D2 + D1 − p k + 2p(x∗1 + m) + ∗ 1 ∗ −m k − x∗1 − x1 (ax1 + 1) p(d1 + D2 )
[
]
> d1 + D2 + D1 − p k + 2p(x∗1 + m). To prove Q1 > 0, it suffices to prove D1 − p k + 2p(x∗1 + m) > 0 which is equivalent to m > by
D1 p
> 1 2
d1 D1 , p(d1 +D2 )
1 2
(k − x∗1 −
D1 p
)−
x∗ 1 2
. Note that
we have
k − x∗1 −
D1
p
−
x∗1 2
<
1 2
k − x∗1 −
Considering (H3 ) there exists m > 12 (k − x∗1 − Further, let g (m) = Q1 Q2 − Q3 . Then
d1 D1 p(d1 + D2 ) D1 p
)−
x∗ 1 2
−
x∗1 2
< k − x∗1 −
d1 D 1 p(d1 + D2 )
.
. Hence, under condition (H3 ) there exists Q1 > 0.
g (m) = (d1 + D2 )(D1 d1 − A(d1 + D2 )) + (D1 − A)(D1 d1 + BC − A(d1 + D2 ))
[ ] D1 d1 − A + (D1 − A) (d1 + D2 ) − A + BC d1 + D 2 d1 + D 2 ∗ D d y D1 y∗ 1 1 ∗ = p(d1 + D2 )2 − k + 2(x∗1 + m) + + p − k + 2 ( x + m ) + 1 p(d1 + D2 ) p(ax∗1 + 1)2 p p(ax∗1 + 1)2 [ ] ∗ y D1 d1 × p(d1 + D2 ) − k + 2(x∗1 + m) + + BC . ∗ p(d1 + D2 ) p(ax1 + 1)2 = (d1 + D2 )2
D1 d1
It is easy to verify that if (H3 ) holds, it has D1 d1 p(d1 + D2 )
− k + 2(x∗1 + m) > 0,
which implies g (m) > 0.
D1 p
− k + 2(x∗1 + m) > 0,
1050
X. Liu, M. Han / Nonlinear Analysis: Real World Applications 12 (2011) 1047–1061
By Q1 > 0 and g (m) > 0, together with Q3 > 0, applying Routh–Hurwitz criterion, all the roots of Eq. (2.6) have negative real parts. This ends the proof. For τ > 0, we investigate the existence of purely imaginary roots of Eq. (2.5). Obviously, iω (ω > 0) is a root of Eq. (2.5) if and only if iω satisfies the following equation
− iω3 − (D2 + d1 + D1 − A)ω2 + i[D1 d1 − A(d1 + D2 )]ω + BC (iω + D2 + d1 )(cos ωτ − i sin ωτ ) = 0.
(2.7)
Separating the real and imaginary parts, we get
−ω3 + [D1 d1 − A(d1 + D2 )]ω = BC [(d1 + D2 ) sin ωτ − ω cos ωτ ],
(2.8)
(d1 + D2 + D1 − A)ω2 = BC [(d1 + D2 ) cos ωτ + ω sin ωτ ], which yields
ω6 + G1 ω4 + G2 ω2 + G3 = 0,
(2.9)
where G1 = (d1 + D2 )2 + (A − D1 )2 + 2D1 D2 > 0, G2 = [D1 d1 − A(d1 + D2 )]2 − B2 C 2 , G3 = −B2 C 2 (d1 + D2 )2 < 0. G1 > 0 and G3 < 0 together indicate that Eq. (2.9) has a unique positive root ω02 . Substituting ω02 into (2.8) and solving for τ , we obtain
ω02 [ω02 + (d1 + D2 )2 + D2 D1 ] 2jπ τj = arccos , + ω0 BC [ω0 2 + (d1 + D2 )2 ] ω0
1
j = 0, 1, 2, 3, . . . .
(2.10)
Theorem 2.1. Let (H1 ) and (H3 ) hold. Then (i) if 0 ⩽ τ < τ0 , the positive equilibrium E ∗ of system (2.1) is locally asymptotically stable; (ii) if τ > τ0 , the positive equilibrium E ∗ is unstable. Furthermore, system (2.1) undergoes Hopf bifurcation at E ∗ when τ = τj , j = 0, 1, 2, 3, . . .. Proof. Define by λ(τ ) = υ(τ ) + iω(τ ) the root of Eq. (2.5) such that
υ(τj ) = 0,
ω(τj ) = ω0 ,
j = 0, 1, 2, . . . .
Substituting λ(τ ) into (2.5) and taking the derivative with respect to τ , we obtain
[
dλ
] −1
[3λ2 + 2(D2 + d1 + D1 − A)λ + D1 d1 − A(d1 + D2 ) + BC e−λτ ] τ − BC λ(λ + D2 + d1 )e−λτ λ 2λ3 + (D2 + d1 + D1 − A)λ2 − BC (d1 + D2 )e−λτ τ = − BC λ2 (λ + D2 + d1 )e−λτ λ [ ] 1 2λ3 + (D2 + d1 + D1 − A)λ2 d1 + D2 τ =− 2 + − . 3 2 λ λ + (D2 + d1 + D1 − A)λ + [D1 d1 − A(d1 + D2 )]λ λ + d1 + D2 λ
=
dτ
When λ = iω0 , then
[ Re
dλ
2ω02 + [(d1 + D2 )2 + (A − D1 )2 + 2D1 D2 ]
]−1
dτ τ =τj
=
]
[(D1 d1 − A(d1 + D2 )) − ω ] + (d1 + D2 + D1 − A) ω 2 2 0
2
2 0
+
(d1 + D2 )2 > 0. ω [ω02 + (d1 + D2 )2 ] 2 0
Hence, sign Re
[
dλ
dτ τ =τj
= sign Re
[
dλ
] −1
dτ τ =τj
= 1.
(2.11)
Therefore, according to Lemmas 2.1 and 2.2, the transversality condition (2.11) and the Hopf bifurcation theorem for functional differential equations [10], we obtain the conclusions. 3. Direction and stability of Hopf bifurcation In Section 2, we have shown that system (2.1) undergoes Hopf bifurcations at τj , j = 0, 1, 2, . . .. Note that Theorem 2.1 only indicates the existence of the periodic solutions. In this section, by referring the works in [11–13] and applying the ideas of Hassard et al. [14], we derive some explicit formulae to determine the direction of the Hopf bifurcation and stability of the bifurcating periodic solutions of system (2.1) at τ = τj .
X. Liu, M. Han / Nonlinear Analysis: Real World Applications 12 (2011) 1047–1061
1051
Let x˜ 1 = x1 − x∗1 , x˜ 2 = x2 − x∗2 , x˜ 3 = y − y∗ and t˜ = τt , τ = τj + µ. Then system (2.1) can be written as a functional differential equation in C = C ([−1, 0], R3 ) x˙ (t ) = Lµ (xt ) + F (µ, xt ),
(3.1)
where x(t ) = (˜x1 (t ), x˜ 2 (t ), x˜ 3 (t )) ∈ R , Lµ : C → R and F : R × C → R are defined, respectively, by T
3
[p k − 2p(x∗1 + m) − y∗ q2 − D1 ]φ1 (0) + D1 φ2 (0) − x∗1 qφ3 (0) , D2 φ1 (0) − (d1 + D2 )φ2 (0) Lµ (φ) = (τj + µ) c β y∗ q2 φ1 (−1)
and
F (µ, φ) = (τj + µ)
− (−1)j aj−2 qj (ay∗ qφ1j (0) − φ1j−1 (0)φ3 (0)) j⩾2 0 , − j j −1 j −1 j −2 j ∗ cβ (−1) a q (ay qφ1 (−1) − φ1 (−1)φ3 (0))
−pφ12 (0) +
j⩾2
where q =
1 . ax∗ +1 1
Note that µ = τ − τj , thus, µ = 0 is a Hopf bifurcation value of system (3.1).
By the Riesz representation theorem, there exists a 3 × 3 matrix function ρ(θ , µ) whose components are bounded variation for θ ∈ [−1, 0], such that Lµ (φ) =
∫
0
dρ(θ , µ)φ(θ ), −1
here, we choose
p k − 2p(x∗1 + m) − y∗ q2 − D1 ρ(θ , µ) = (τj + µ) D2 0
D1
−d1 − D2 0
−x∗1 q 0 δ(θ ) − 0
0 0 c β y∗ q2
0 0 0
0 0 δ(θ + 1) , 0
where δ(θ ) is delta function. For φ ∈ C = C 1 ([−1, 0], R3 ), define
dφ(θ ) , dθ ∫ M (µ)φ = 0 dρ(v, µ)φ(v),
θ ∈ [−1, 0), R(µ)φ =
θ = 0,
0, F (µ, φ),
θ ∈ [−1, 0), θ = 0.
(3.2)
−1
Then system (3.1) can be written as x˙ t = M (µ)xt + R(µ)xt ,
(3.3)
where xt = x(t + θ ) for θ ∈ [−1, 0]. For φ ∈ C = C 1 ([−1, 0], R3 ), ϕ ∈ C ∗ = C 1 ([0, 1], R3∗ ), we define the adjoint operator M∗ of M (0) by
dϕ(s) , − ds M∗ ϕ(s) = ∫ 0 dρ T (t , 0)ϕ(−t ),
s ∈ (0, 1], s = 0,
−1
and the adjoint bilinear form by
⟨ϕ(s), φ(s)⟩ = ϕ( ¯ 0)φ(0) −
0
∫
∫
−1
θ
υ=0
ϕ(υ ¯ − θ )dρ(θ, 0)φ(υ)dυ.
(3.4)
In terms of the discussion of Section 2, we see that ± iω0 τj are eigenvalues of M (0) and M∗ . Next we calculate the eigenvector η(θ ) of M (0) corresponding to iω0 τj and the eigenvector η∗ (s) of M∗ belonging to −iω0 τj , respectively. By the definition of M (0) and M∗ , it is not difficult to obtain that T iω0 τj θ
η(θ ) = (1, b1 , b2 ) e
= 1,
D2 iω0 + d1 + D2
,
c β y∗ q2 e−iω0 τj
T
iω0
and
η∗ (s) = D(1, b∗1 , b∗2 )eiω0 τj s = D 1,
D2
,
x∗1 q
−iω0 + d1 + D2 iω0
eiω0 τj s .
eiω0 τj θ
1052
X. Liu, M. Han / Nonlinear Analysis: Real World Applications 12 (2011) 1047–1061
It follows from (3.4) that
∫
0
θ
∫
η¯ ∗ (ξ − θ )dρ(θ , 0)η(ξ )dξ ⟨η , η⟩ = η¯ (0)η(0) − −1 ξ =0 0 0 0 = D¯ 1 + b1 b¯ ∗1 + b2 b¯ ∗2 + τj (1, b¯ ∗1 , b¯ ∗2 ) 0 c β y∗ q2 0 ∗
∗
0 1 0 b1 b2 0
e−iω0 τj
= D¯ (1 + b1 b¯ ∗1 + b2 b¯ ∗2 + τj c β y∗ q2 b¯ ∗2 e−iω0 τj ). ¯ = 0. We choose D = (1 + b¯ 1 b∗1 + b¯ 2 b∗2 + τj c β y∗ q2 b∗2 eiω0 τj )−1 , such that ⟨η∗ , η⟩ = 1 and ⟨η∗ , η⟩ Following the algorithms in [14], we compute the coordinates describing center manifold C0 at µ = 0. Define z (t ) = ⟨η∗ , xt ⟩,
W (t , θ ) = xt (θ ) − z (t )η(θ ) − z (t )η(θ ),
(3.5)
therefore, xt (θ ) = (˜x1t (θ ), x˜ 2t (θ ), x˜ 3t (θ ))T
= W (t , θ ) + z (1, b1 , b2 )T eiθω0 τj + z¯ (1, b¯ 1 , b¯ 2 )T e−iθ ω0 τj .
(3.6)
On the center manifold C0 , we know that W (t , θ ) = W (z (t ), z¯ (t ), θ) = W20
z2 2
+ W11 z z¯ + W02
z¯ 2 2
+ ···,
(3.7)
where z and z¯ are local coordinates for the center manifold C0 in the direction of η∗ and η¯ ∗ . Note that W is real if xt is real. For the solution xt ∈ C0 of (3.3), when µ = 0, it follows from (3.2)–(3.5) that z˙ (t ) = ⟨η∗ , x˙ t ⟩ = iω0 τj z (t ) + ⟨η∗ (θ ), R(0)(W (z , z¯ , θ ) + 2Re{z (t )η(θ )})⟩
= iω0 τj z (t ) + η¯ ∗ (0)F (0)(W (z , z¯ , 0) + 2Re{z (t )η(0)}) , iω0 τj z (t ) + g (z , z¯ ),
(3.8)
where z2
z¯ 2
z 2 z¯
+ ···. (3.9) 2 2 2 From [14], we know that the properties of bifurcating periodic solutions at the critical value τj are determined by the following values g (z , z¯ ) = g20
c1 (0) =
+ g11 z z¯ + g02
i 2ω0 τj
+ g21
2
g20 g11 − 2|g11 | −
β2 = 2Re{c1 (0)},
T2 = −
|g02 |2 3
+
g21 2
,
µ2 = −
Im{c1 (0)} + µ2 Im{λ′ (τj )}
ω0 τ j
Re{c1 (0)} Re{λ′ (τj )}
,
.
Therefore, we need to compute the values of g20 , g11 , g02 and g21 . Firstly, from (3.6) and (3.7), we obtain (1)
z2
(3)
2 z2
x˜ 1t (0) = W20 (0) x˜ 3t (0) = W20 (0) (1)
2
(1) (1) + W11 (0)z z¯ + W02 (0) (3) (3) + W11 (0)z z¯ + W02 (0)
x˜ 1t (−1) = W20 (−1)
z2 2
z¯ 2 2 z¯ 2 2
+ z + z¯ + O(|(z , z¯ )|3 ); + b2 z + b¯ 2 z¯ + O(|(z , z¯ )|3 );
(1) (1) + W11 (−1)z z¯ + W02 (−1)
z¯ 2 2
+ ze−iω0 τj + z¯ eiω0 τj + O(|(z , z¯ )|3 ).
On the other hand, in terms of the definition of F (µ, xt ), we have g (z , z¯ ) = η¯ ∗ (0)F (0)(W (z , z¯ , 0) + 2Re{z (t )η(0)}) = η¯ ∗ (0)F0 (z , z¯ )
− j−1 −px˜ 21t (0) + (−1)j aj−2 qj (ay∗ qx˜ j1t (0) − x˜ 1t (0)˜x3t (0)) j⩾2 0 = τj D¯ (1, b¯ ∗1 , b¯ ∗2 )T − j −1 j−1 j−2 j ∗ j cβ (−1) a q (ay qx˜ 1t (−1) − x˜ 1t (−1)˜x3t (0)) j⩾2
= τj D¯ {z 2 (−p + aq3 y∗ − q2 b2 − b¯ ∗2 c β aq3 y∗ e−iω0 τj + b¯ ∗2 b2 c β q2 e−iω0 τj ) + 2z z¯ (−p + aq3 y∗ − Re{b2 }q2 − b¯ ∗2 c β aq3 y∗ + b¯ ∗2 c β q2 Re{b2 eiω0 τj })
(3.10)
X. Liu, M. Han / Nonlinear Analysis: Real World Applications 12 (2011) 1047–1061
1053
+ z¯ 2 (−p + aq3 y∗ − q2 b¯ 2 − b¯ ∗2 c β aq3 y∗ eiω0 τj + b¯ ∗2 b¯ 2 c β q2 eiω0 τj ) z 2 z¯ (3) (1) (1) (−1) [2(−p + ay∗ q3 )(W20 (0) + 2W11 (0)) − 2b¯ ∗2 c β ay∗ q3 eiω0 τj W20 + 2
(1) (1) (3) (3) − q2 (W20 (0)b¯ 2 + 2W11 (0)b2 + 2W11 (0) + W20 (0)) ( 3 ) ( 3 ) (1) (1) ∗ 2 i ω τ − i ω τ + b¯ 2 c β q (W20 (0)e 0 j + 2W11 (0)e 0 j + 2b2 W11 (−1) + b¯ 2 W20 (−1))
+ 2aq3 (b¯ 2 + 2b2 ) − 2b¯ ∗2 c β aq3 (2b2 + b¯ 2 e−2iω0 τj ) + 6a2 q4 y∗ (b¯ ∗2 c β e−iω0 τj − 1)] + · · ·}. By comparing coefficients with (3.9), we get
¯ (−p + aq3 y∗ − q2 b2 − b¯ ∗2 c β aq3 y∗ e−iω0 τj + b¯ ∗2 b2 c β q2 e−iω0 τj ); g20 = 2τj D ¯ (−p + aq3 y∗ − Re {b2 } q2 − b¯ ∗2 c β aq3 y∗ + b¯ ∗2 c β q2 Re{b2 eiω0 τj }); g11 = 2τj D ¯ (−p + aq3 y∗ − q2 b¯ 2 − b¯ ∗2 c β aq3 y∗ eiω0 τj + b¯ ∗2 b¯ 2 c β q2 eiω0 τj ); g02 = 2τj D (1)
(3)
(1)
¯ [2(−p + ay∗ q3 )(W20 (0) + 2W11 (0)) − 2b¯ ∗2 c β ay∗ q3 eiω0 τj W20 (−1) g21 = τj D (1)
(1)
(3)
(3)
− q (W20 (0)b¯ 2 + 2W11 (0)b2 + 2W11 (0) + W20 (0)) (3) (3) (1) (1) + b¯ ∗2 c β q2 (W20 (0)eiω0 τj + 2W11 (0)e−iω0 τj + 2b2 W11 (−1) + b¯ 2 W20 (−1)) + 2aq3 (b¯ 2 + 2b2 ) − 2b¯ ∗2 c β aq3 (2b2 + b¯ 2 e−2iω0 τj ) + 6a2 q4 y∗ (b¯ ∗2 c β e−iω0 τj − 1)]. 2
Obviously, the value of g21 is determined by W20 (θ ) and W11 (θ ). From (3.3) and (3.5), we have
˙ (t , θ ) = x˙ t (θ ) − 2Re{˙z (t )η(θ )} W = M (0)W − 2Re{g (z , z¯ )η(θ )} + R(0)(W + 2Re{z (t )η(θ )}) , M (0)W + H (z , z¯ , θ ),
(3.11)
where H (z , z¯ , θ ) = −2Re{g (z , z¯ )η(θ )} + R(0)(W + 2Re{z (t )η(θ )}) , H20 (θ )
z2
+ H11 (θ )z z¯ + H02 (θ )
z¯ 2
2 2 It follows from (3.12) and the definition of R(0) that H (z , z¯ , θ ) =
−2Re{g (z , z¯ )η(θ )}, −2Re{g (z , z¯ )η(0)} + F0 (z , z¯ ),
= H20 (θ )
z2
+ H11 (θ )z z¯ + H02 (θ )
z¯ 2
2 2 On the center manifold C0 near the origin, we have
+ ···.
(3.12)
θ ∈ [−1, 0) θ =0 + ···,
(3.13)
˙ = Wz z˙ + Wz¯ z˙¯ = W20 z z˙ + W11 z¯ z˙ + W11 z z˙¯ + W02 z¯ z˙¯ + · · · . W
(3.14)
Substituting z˙ = iω0 τj z (t ) + g (z , z¯ ) into (3.14) and comparing the coefficients with (3.11), we obtain
(M (0) − 2iω0 τj )W20 (θ ) = −H20 (θ ),
M (0)W11 (θ ) = −H11 (θ ).
(3.15)
It follows from (3.9) and (3.13) that for θ ∈ [−1, 0), H20 (θ ) = −g20 η(θ ) − g¯20 η(θ ¯ ),
H11 (θ ) = −g11 η(θ ) − g¯11 η(θ ¯ ).
Substituting H20 (θ ) and H11 (θ ) into (3.15) and using the definition of M (0) for θ ∈ [−1, 0), we obtain
˙ 20 (θ ) = 2iω0 τj W20 (θ ) + g20 η(θ ) + g¯20 η(θ W ¯ ),
˙ 11 (θ ) = g11 η(θ ) + g¯11 η(θ W ¯ ).
A straightforward computation yields W20 (θ ) =
ig20 η(0) iω0 τj θ ig¯20 η( ¯ 0) −iω0 τj θ e + e + S1 e2iω0 τj θ , ω0 τj 3ω0 τj
W11 (θ ) = −
ig11 η(0) iω0 τj θ ig¯11 η( ¯ 0) −iω0 τj θ e + e + S2 ,
ω0 τj
(3.16)
ω0 τj
3
where S1 , S2 ∈ R are undetermined constant vectors. In the following, we will determine the values of S1 and S2 . 0 Note that M (0)W20 = −1 dρ(v, 0)W20 (v) for θ = 0, then by (3.15), we have
∫
0
dη(v, 0)W20 (v) = 2iω0 τj W20 (0) − H20 (0), −1
∫
0
dη(v, 0)W20 (v) = −H11 (0). −1
(3.17)
1054
X. Liu, M. Han / Nonlinear Analysis: Real World Applications 12 (2011) 1047–1061
It follows from (3.13) that for θ = 0,
−p + aq3 y∗ − q2 b2
, 0 −c β aq3 y∗ e−iω0 τj + b2 c β q2 e−iω0 τj −p + aq3 y∗ − Re {b2 } q2 0 H11 (0) = −g11 η(0) − g¯11 η( ¯ 0) + 2τj iω τ . 3 ∗ 2 0 j −c β aq y + c β q Re b2 e
H20 (0) = −g20 η(0) − g¯20 η( ¯ 0) + 2τj
(3.18)
Substituting H20 (0) defined by (3.18) and W20 (θ ) defined by (3.16) in (3.17) respectively, and noting that
iω0 τj I −
∫
0
e
iω0 τj v
dρ(v, 0) η(0) = 0,
0
∫ −iω0 τj I −
e
−iω0 τj v
¯ 0) = 0, dρ(v, 0) η(
−1
−1
we derive
2iω0 τj I −
∫
0
e2iω0 τj v dρ(v, 0) S1 = 2τj
−1
−p + aq3 y∗ − q2 b2 0
−c β aq3 y∗ e−iω0 τj + b2 c β q2 e−iω0 τj
.
Solving the above linear equations for S1 , we may obtain 2iω0 τj − p k − 2p(x∗1 + m) + y∗ q2 + D1 −D 2 S1 = 2 −c β y∗ q2 e−2iω0 τj
×
3 ∗
0
−1
x∗1 q 0 2iω0 τj
2
−p + aq y − q b2 0
−D1 2iω0 τj + d1 + D2
−c β aq3 y∗ e−iω0 τj + b2 c β q2 e−iω0 τj
.
Similarly, from (3.16)–(3.18), we may get S2 which is defined by D1 − p k − 2p(x∗1 + m) + y∗ q2 −D2 S2 = 2 −c β y∗ q2
−D1 d1 + D 2 0
−1
x∗1 q 0 0
aq3 y∗ − p − Re {b2 } q2 . 0 c β q2 Re b2 eiω0 τj − c β aq3 y∗
Substituting S1 and S2 into (3.16), we may calculate W20 (θ ) and W11 (θ ). Furthermore, g21 can be determined. Thus, by (3.10), we can compute the values of c1 (0), µ2 , β2 and T2 . It is clear that µ2 determines the direction of the Hopf bifurcation: if µ2 > 0 (µ2 < 0), then the Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solutions exist for τ > τj (τ < τj ); β2 determines the stability of bifurcating periodic solutions: if β2 < 0 (β2 > 0) then the bifurcating periodic solutions are orbitally asymptotically stable (unstable); and T2 determines the period of the bifurcating periodic solutions: the period increases (decreases) if T2 > 0 (T2 < 0). To illustrate the feasibility of our main results, we give an example and carry out numerical simulations by using Matlab. Example. In (2.1), let α = 1.8, β = 0.8, k = 2.6, a = 1.2, c = 0.9, d1 = 0.1, d2 = 0.2, D1 = 0.8, D2 = 0.3 and D d m = 0.3. Clearly, c β − ad2 = 0.48 > 0 and 0 < m < k − x∗1 − p(d 1+1D ) = 1.8944. Then conditions (H1 ) and (H3 ) in 1 2 Theorem 2.1 are satisfied. System (2.1) has a unique positive equilibrium E ∗ (0.4167, 0.5375, 2.8479). Solving (2.9) for ω, we ∗ obtain ω0 ≈ 0.2565334393, then τ0 ≈ 3.916807176. By Theorem 2.1-(i) E is locally asymptotically stable for 0 ⩽ τ < τ0 . One can see Fig. 1(a)–(d). Furthermore, when the delay τ > τ0 , the positive equilibrium E ∗ loses its stability and the system goes into oscillations. The system undergoes a Hopf bifurcation at τ = τ0 . By computation we obtain g20 = 3.470469910 + 2.439350744i; g11 = 9.520322552 + 8.634501512i; g02 = 11.83851897 + 14.74884396i; g21 = −36.6992317 + 326.174705i; c1 (0) = −44.817374 − 54.6846853i;
µ2 = 806.3833986 > 0; β2 = −89.6347475 < 0; T2 = 61.45127048 > 0. Then we can conclude that the Hopf bifurcation is supercritical. Moreover, the bifurcating periodic solution for τ > τ0 is asymptotically stable (see Fig. 2(a)–(d)).
X. Liu, M. Han / Nonlinear Analysis: Real World Applications 12 (2011) 1047–1061
1055
Fig. 1. The time course of a solution tending to the equilibrium (0.4167, 0.5375, 2.8479) of system (2.1) with τ = 3.4 < τ0 , where α = 1.8, β = 0.8, k = 2.6, a = 1.2, c = 0.9, d1 = 0.1, d2 = 0.2, D1 = 0.8, D2 = 0.3, m = 0.3 and initial data x1 (0) = 0.6, x2 (0) = 0.6, y(0) = 2.76. (a) Time-state t − x1 (t ) relation graph; (b) Phase portrait of t − x2 (t ) relation graph; (c) Time-state t − y(t ) relation graph; (d) Phase portrait of x1 (t ), x2 (t ) and y(t ).
4. Numerical analysis of system (1.2) with impulsive In this section, by using numerical method, we will discuss the effect of impulsive perturbations on system (1.2). As we know, impulsive differential equations can be seen in almost every domain of applied science and studied extensively (see, for example [15–18]). In general, the impulsive perturbations will bring to the system complexity because of their propensity for chaos. Motivated by Meng et al. [17], Tang and Chen [18], by introducing periodic impulsive proportional harvesting (releasing) on the prey x2 and predator y, respectively, we develop model (1.2) as follows
β(x1 − m)y + D1 (x2 − x1 ), k 1 + a(x1 − m) x˙ 2 = −d1 x2 + D2 (x1 − x2 ), c β[x1 (t − τ ) − m]y , y˙ = −d2 y + 1 + a[x1 (t− τ ) − m] x1 (t + ) = x1 (t ), x2 (t + ) = (1 − p)x2 (t ), t = nT , n ∈ N , y(t + ) = (1 − δ)y(t ),
x˙ 1 = α x1 1 −
x1
−
t ̸= nT , n ∈ N , (4.1)
where 0 < p < 1 and 0 < δ < 1 (p < 0 and δ < 0) represent partial impulsive harvest (release) to preys x2 and predators, respectively, at t = nT , n ∈ N. T is the period of impulsive effect.
1056
X. Liu, M. Han / Nonlinear Analysis: Real World Applications 12 (2011) 1047–1061
Fig. 2. The trajectory graph of system (2.1) with the same parameters as in Fig. 1, when τ = 3.98 > τ0 , it shows the bifurcating periodic solution from (0.4167, 0.5375, 2.8479) which is asymptotically stable. (a) Time-state t − x1 (t ) relation graph; (b) Time-state t − x2 (t ) relation graph; (c) Phase portrait of t − y(t ); (d) Phase portrait of x1 (t ), x2 (t ) and y(t ).
Firstly, we investigate the effect of period of impulsive control T on dynamical behavior of system (1.2). For convenience, we keep the parameters α = 1.8, β = 0.8, k = 2.6, a = 1.2, c = 0.9, d1 = 0.1, d2 = 0.2, D1 = 0.8, D2 = 0.3, m = 0.3 fixed, and let p = 0.2, δ = 0.32 and τ = 6. Then the influence of T may be documented by stroboscopically sampling some of the variables over a range of T . Here, we carry out numerical simulations with respect to T in the range 0.5 ⩽ T ⩽ 44. And the resulting bifurcation diagrams provide a description of essential dynamical behavior of the system. With the increasing of T , the system experiences process of cycles, period-doubling bifurcation, period-halving bifurcation, intermittency (intermittency occurs whenever the behavior of a system seems to switch back and forth between two qualitatively different behaviors as a control parameter is changed), chaotic band, narrow periodic window, wide periodic window and non-unique dynamics, i.e., several attractors may coexist with the same T , and every attractor is corresponding to the initial values. In order to investigate the existence of the predator-free solution, we give Fig. 3, firstly. And it shows that the predator population goes into extinction if T < Tmin ≈ 1.8, where Tmin can be regarded as the minimum period of impulsive control to prevent the predator from extinction. From the bifurcation diagrams Fig. 4, we observe that a period-doubling bifurcation occurs when T ≈ 10.66. It means that T -periodic solution disappears suddenly at this point and 2T -periodic solution appears. Then the behavior of 2T -periodicity is kept for 10.66 < T < 16.6256. When T closes to 16.626, the 2T -periodic solution becomes unstable and shrinks to a chaotic attractor, rapidly (see Fig. 6). Following, system (4.1) comes into a chaotic region with some narrow periodic
X. Liu, M. Han / Nonlinear Analysis: Real World Applications 12 (2011) 1047–1061
1057
Fig. 3. (a) The bifurcation diagrams of the predator y for T over [0.5, 5]; (b) Time series of the prey x1 when T = 1.8; (c) Time series of the prey x2 when T = 1.8; (d) Time series of the predator y when T = 1.8.
10 9.5 9 8.5 8 7.5 7 6.5 6 5.5 5 4.5 4 3.5 3
b
Max(y)
Max(y)
a
5
6
7
8
9 10 11 12 13 14 15 16 17 18 19 20
10 9.5 9 8.5 8 7.5 7 6.5 6 5.5 5 4.5 4 3.5 3 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
T
c
T
d
9 8.5 8 7
Max(y)
Max(y)
7.5 6.5 6 5.5 5 4.5 4 35
35.5
36
36.5
37
37.5 T
38
38.5
39
39.5
40
10 9.5 9 8.5 8 7.5 7 6.5 6 5.5 5 4.5 4 3.5 3 40
40.4 40.8 41.2 41.6
42
42.4 42.8 43.2 43.6
44
T
Fig. 4. Bifurcation diagrams of system (5.1) showing the influence of T with α = 1.8, β = 0.8, k = 2.6, a = 1.2, c = 0.9, d1 = 0.1, d2 = 0.2, D1 = 0.8, D2 = 0.3, m = 0.3, p = 0.2 and δ = 0.32 when varying T from 5 to 44. (a) The predator y is plotted for T over interval [5, 20]. (b) The predator y is plotted for T over interval [20, 35]. (c) The predator y is plotted for T over interval [35, 40]. (d) The predator y is plotted for T over interval [40, 44].
windows. When T ⩾ 19.25 and nearby 19.25, the chaos bands disappear and a 3T -periodic solution occurs. And when T = 20.84, a non-unique attractor appears, i.e., a 3T -periodic solution coexists with a T -periodic solution (see Fig. 7). when T ≈ 31.49, the system displays another non-unique attractor (two different types of T -periodic solutions coexist). When varying T in the range [31.5, 36.91], there is a transition from periodicity to quasi-periodicity to chaos to periodicity again (see Fig. 5(b)). As T increases further, the evidence for 9T -periodic solutions leading to chaos can be observed for 36.91 < T < 37.25. And the chaotic behavior is kept until T closes to the value 37.38, then the system displays a 4T -periodicity behavior. When 37.38 ⩽ T ⩽ 38.798, there is 4T -periodic solution with twice period-doubling bifurcations and period-halving bifurcations. When T = 38.798, the non-unique attractor appears again (a 4T -periodic solution and a
1058
X. Liu, M. Han / Nonlinear Analysis: Real World Applications 12 (2011) 1047–1061
a
b
8.5
8.5
8
8
7.5
7.5
7
7 6.5
Max(y)
Max(y)
6.5 6 5.5
6 5.5
5 4.5
5
4
4.5
3.5
16.5
4 36.2
16.8 17 17.2 17.4 17.6 17.8 18 18.2 18.4 18.6 18.8 19
36.3
36.4
36.5
36.6
T
36.7
36.8
36.9
37
T
Fig. 5. The magnified parts of Fig. 4. (a) The predator y is plotted for T over interval [16.5, 19]. (b) The predator y is plotted for T over interval [36.2, 37].
a
7
b
T=16.625
6
6
T=16.627
5 4
4
y
y
5
3
3
2
2
1 0
1 0.5
0.5 x1
1 1.5
0.4
2 0.2
0.8
0.6
0.75
1.2
1
x1 1
x2
1.25 1.5 0.2
0.4
0.6
0.8
1
1.2
x2
Fig. 6. The dynamical behaviors of system (5.1) corresponding to T = 16.625 and T = 16.627, respectively. (a) A 2T -periodic solution corresponding to T = 16.625. (b) A chaotic attractor corresponding to T = 16.627. 8
b
T=20.84
6
6
5
4
4
y
y
a
T=20.84
3
2
2 0
0 0 0.5 x1
1 1.5 0.2
0.4
0.6
0.8
1
1.2
0.5 x1
x2
1 1.5 0.4
0.5
0.6
0.7 x2
0.8
0.9
Fig. 7. Coexistence of 3T -periodic solution with a T -periodic solution when the period T = 20.84. (a) With initial value x1 (0) = 0.6, x2 (0) = 0.6, y(0) = 2.76, the solution will finally tend to a 3T -periodic solution. (b) With initial value x1 (0) = 0.4, x2 (0) = 0.1, y(0) = 0.1, the solution will finally tend to a T -periodic solution.
3T -periodic solution may coexists; see Fig. 8). When T nearby 38.9, the behavior of period-doubling bifurcation and periodhalving bifurcation is given second. When T ⩾ 39.78, the 3T -periodic solution occurs again, and we note that until T = 42.3 no chaos occurs. Following the chaotic area, the behavior of cascade of period-halving bifurcations from chaos to cycles may also be seen for 42.85 < T < 43.6. In addition, the phenomenon of intermittency also can be illustrated by Fig. 5(a) which shows a series of switch from chaotic band to periodic window to chaotic window back and forth. On the other hand, in order to further investigate the dynamical behaviors of system (4.1), the bifurcation diagrams with respect to the parameters δ , p and m are also carried out, respectively. Now, consider the following set of parameters
α = 1.8, D1 = 0.8,
β = 0.8, D2 = 0.3,
k = 2.6, m = 0.3,
a = 1.2,
τ = 6,
c = 0.9,
d1 = 0.1,
d2 = 0.2,
T = 18.
Let p = 0.1, then the bifurcation diagrams with respect to parameter δ over the range [−0.26, 0.4] are shown in Fig. 9. From these bifurcation diagrams, we can see that system (4.1) displays complex dynamics with increasing δ from −0.26 to 0.4, gradually. For example, some strange attractors are captured when δ = −0.1, −0.075, 0.315, and 0.35, respectively (see Figs. 10 and 11).
X. Liu, M. Han / Nonlinear Analysis: Real World Applications 12 (2011) 1047–1061
a
10
b
T=38.798
8
1059
T=38.798
7
8
6 6 y
y
5 4
4
3 2
2
0 0
1 0 0.5 1
x1
0.4
1.5 0.2
0.8 x2
0.6
1
0.5
1.2
1
x1
0.6
0.4
1.5 0.2
0.8
1
1.2
x2
Fig. 8. Coexistence of 4T -period solution with a 3T -periodic solution when the period T = 38.798. (a) With initial value x1 (0) = 0.4, x2 (0) = 0.1, y(0) = 0.1, the solution will finally tend to a 4T -periodic solution. (b) With initial value x1 (0) = 0.6, x2 (0) = 0.6, y(0) = 2.76, the solution will finally tend to a 3T -periodic solution.
a
b
6.5
7 6.5
5.5
6
5
5.5
Max(y)
Max(y)
6
4.5
5
4
4.5
3.5
4
3
3.5
-0.26 -0.235 -0.21 -0.185 -0.16 -0.135 -0.11 -0.085 -0.06 -0.035 -0.010
0
0.05
0.1
0.15
0.2
δ
0.25
0.3
0.35
0.4
δ
Fig. 9. The bifurcation diagrams of system (5.1) with respect to the parameter δ in the range −0.26 ⩽ δ ⩽ 0.4. (a) −0.26 ⩽ δ ⩽ 0; (b) 0 ⩽ δ ⩽ 0.4.
a
6
b
δ=-0.11
5.5
6
δ=-0.075
5.5
y
5 4.5
y
5 4.5 4
4
3.5
3.5
3
3
2.5 0
2.5 0 0.5 x1
1 1.5
0.4
0.5
0.6 x2
0.7
0.8
0.5 1
x1
1.5
0.5
0.4
0.6 x2
0.7
0.8
Fig. 10. The dynamical behaviors of system (5.1) when α = 1.8, β = 0.8, k = 2.6, a = 1.2, c = 0.9, d1 = 0.1, d2 = 0.2, D1 = 0.8, D2 = 0.3, m = 0.3, p = 0.1 and T = 18. (a) A strange attractor corresponding to δ = −0.11. (b) A strange attractor corresponding to δ = −0.075.
a
b
6 δ=0.315
5.5 δ=0.35
5
5
y
y
4.5 4
4 3.5 3 2.5
3 2 0.5 1 x1
1.5 0.2
0.4
0.6 x2
0.8
1
2 0.5 1 x1
1.5 0.4
0.5
0.6
0.7
0.8
0.9
x2
Fig. 11. The dynamical behaviors of system (5.1) when α = 1.8, β = 0.8, k = 2.6, a = 1.2, c = 0.9, d1 = 0.1, d2 = 0.2, D1 = 0.8, D2 = 0.3, m = 0.3, p = 0.1 and T = 18. (a) A strange attractor corresponding to δ = 0.315. (b) A strange attractor corresponding to δ = 0.35.
1060
X. Liu, M. Han / Nonlinear Analysis: Real World Applications 12 (2011) 1047–1061
a
b
6 5.8
6.5 6
5.6 5.5 Max(y)
Max(y)
5.4 5.2 5 4.8
5 4.5 4
4.6 3.5
4.4
3
4.2 p
p
Fig. 12. The bifurcation diagrams of system (5.1) with respect to the parameter p in the range −0.8 ⩽ p ⩽ 0.85. (a) −0.8 ⩽ p ⩽ 0; (b) 0 ⩽ p ⩽ 0.85.
16 15 14 13 12 11 10 9 8 7 6 5 4 3 2
b
Max(y)
Max(y)
a
0
0.05
0.1
0.15
0.2 m
0.25
0.3
0.35
0.4
12 11.5 11 10.5 10 9.5 9 8.5 8 7.5 7 6.5 6 5.5 5 4.5 4 3.5 3
0.15
0.16
0.17
0.18
0.19
0.2
0.21
0.225
m
Fig. 13. The bifurcation diagrams of system (5.1) with respect to the parameter p in the range 0 ⩽ m ⩽ 0.4. (a) 0 ⩽ m ⩽ 0.4; (b) 0.15 ⩽ m ⩽ 0.225.
Similarly, if we keep the values of α , β , k, a, c, d1 , d2 , D1 , D2 , m, T and τ fixed, and let δ = 0.3. For −0.8 ⩽ p ⩽ 0.85 the effect of p on the dynamics of system (4.1) can be shown in Fig. 12. From the discussion of Section 2, it is clear that the behavior of prey refuge also has an important influence on the system. For this reason, we plot the bifurcation diagrams of system (4.1) with respect to m in the range 0 ⩽ m ⩽ 0.4 (see Fig. 13). By choosing T , δ , p and m as bifurcation parameters, respectively, the rich and complex dynamical behaviors of system (4.1) have been observed. All these results show that the dynamical behaviors of system (4.1) will change evidently with varying the values of some parameters. References [1] E. Gonlez-Olivares, R. Ramos-Jiliberto, Dynamic consequences of prey refuges in a simple model system: more prey, fewer predators and enhanced stability, Ecological Modelling 166 (2003) 135–146. [2] L. Chen, F. Chen, L. Chen, Qualitative analysis of a predator–prey model with Holling type II functional response incorporating a constant prey refuge, Nonlinear Analysis: Real World Applications 11 (2010) 246–252. [3] T.K. Kar, Modelling and analysis of a harvested prey–predator system incorporating a prey refuge, Journal of Computational and Applied Mathematics 185 (2006) 19–33. [4] L. Ji, C. Wu, Qualitative analysis of a predator–prey model with constant-rate prey harvesting incorparating a constant prey refuge, Nonlinear Analysis: Real World Applications 11 (2010) 2285–2295. [5] Z. Ma, et al., Effects of prey refuges on a predator–prey model with a class of functional responses: the role of refuges, Mathematical Biosciences 218 (2009) 73–79. [6] Y. Kuang, Y. Takeuchi, Predator–prey dynamics in models of prey dispersal in two-patch environments, Mathematical Biosciences 120 (1994) 77–98. [7] J. Cui, L. Chen, The effect of diffusion on the time varying logistic population growth, Computers & Mathematics with Applications 36 (1998) 1–9. [8] J. Cui, L. Chen, Permanence and extinction in logistic and Lotka–Volterra systems with diffusion, Journal of Mathematical Analysis and Applications 258 (2001) 512–535. [9] R. Xu, L. Chen, Persistence and stability for a two-species ratio-dependent predator–prey system with time delay in a two-patch environment, Computers and Mathematics with Applications 40 (2000) 577–588. [10] J. Hale, S.V. Lunel, Introduction to Functional Differential Equations, Springer-Verlag, New York, 1993. [11] F. Lian, Y. Xu, Hopf bifurcation analysis of a predator prey system with Holling type IV functional response and time delay, Applied Mathematics and Computation 215 (2009) 1484–1495. [12] Y. Song, S. Yuan, Bifurcations analysis in a predator–prey system with time delay, Nonlinear Analysis: Real World Applications 7 (2006) 265–284. [13] R. Xu, Q. Gan, Z. Ma, Stability and bifurcation analysis on a ratio-dependent predator–prey model with time delay, Journal of Computational and Applied Mathematics 230 (2009) 187–203. [14] B.D. Hassard, N.D. Kazarinoff, Y.H. Wan, Theory and Applications of Hopf Bifurcation, Cambridge University Press, Cambridge, 1981. [15] S. Zhang, L. Chen, The study of predator–prey system with defensive ability of prey and impulsive perturbations on the predator, Chaos, Solitons & Fractals 23 (2005) 311–320.
X. Liu, M. Han / Nonlinear Analysis: Real World Applications 12 (2011) 1047–1061
1061
[16] Y. Xia, Positive periodic solutions for a neutral impulsive delayed Lotka–Volterra competition system with the effect of toxic substance, Nonlinear Analysis: Real World Applications 8 (2007) 204–221. [17] X. Meng, L. Chen, H. Cheng, Two profitless delays for the SEIRS epidemic disease model with nonlinear incidence and pulse vaccination, Applied Mathematics and Computation 186 (2007) 516–529. [18] S. Tang, L. Chen, Density-dependent birth rate, birth pulse and their population dynamic consequences, Journal of Mathematical Biology 44 (2002) 185–199.