- Email: [email protected]
Global stability and travelling waves of a predator-prey model with diffusion and nonlocal maturation delay
Commun Nonlinear Sci Numer Simulat 15 (2010) 3390–3401
Contents lists available at ScienceDirect
Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns
Global stability and travelling waves of a predator-prey model with diffusion and nonlocal maturation delay Xiao Zhang, Rui Xu * Institute of Applied Mathematics, Shijiazhuang Mechanical Engineering College, Shijiazhuang 050003, PR China
a r t i c l e
i n f o
Article history: Received 10 June 2009 Accepted 29 December 2009 Available online 4 January 2010 AMS subject classification: 35K57 35K55 92D25
a b s t r a c t In this paper, a diffusive predator-prey system with nonlocal maturation delay is investigated. By analyzing the corresponding characteristic equations, the local stability of each of uniform steady states of the system is discussed. Sufficient conditions are derived for the global stability of the positive steady state and the semi-trivial steady state of the system by using the method of upper–lower solutions and its associated monotone iteration scheme, respectively. The existence of travelling wave solution of the system is established by using the geometric singular perturbation theory. Numerical simulations are carried out to illustrate the theoretical results. Ó 2010 Elsevier B.V. All rights reserved.
Keywords: Diffusion Nonlocal delay Global stability Travelling wave Geometric singular perturbation
1. Introduction Predator-prey model is one of the basic models between different species in nature. In traditional predator-prey models the spatial content of the environment has been ignored. These models have been traditionally formulated in relation to the time evolution of uniform population distributions in the habitat. However, in many ecological systems, the species under consideration may disperse spatially as well as evolving in time ([12]). This spatial dispersal mainly due to resource limitation: in regions of high population density, food will become scarce, and individuals will tend to migrate to regions of lower population density. The effect of dispersion of a population in a bounded habitat has been researched by many authors, and in this situation the governing equations for the population densities are described by a system of reaction-diffusion equations. In realistic ecological models, any delays should be spatially inhomogeneous, that is, the delays affect both the temporal and spatial variables. This is due to the fact that any given individual may not necessarily have been at the same spatial location at precious times. Such delays are called spatio-temporal delays or nonlocal delays. In recent years, the effect of nonlocal delays on the dynamics of ecological models has been taken into account (see, for example, [2,5,6,9–16]). In [1], Al-Omari and Gourley studied the following single species population with stage structure and nonlocal maturation delay
* Corresponding author. E-mail addresses: [email protected] (X. Zhang), [email protected] (R. Xu). 1007-5704/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.cnsns.2009.12.031
X. Zhang, R. Xu / Commun Nonlinear Sci Numer Simulat 15 (2010) 3390–3401
Z þ1 Z þ1 @ui @ 2 ui Gðx; y; sÞkðsÞecs um ðy; t sÞdyds; ¼ d1 2 þ aum cui a @t @x 0 1 Z þ1 Z þ1 @um @ 2 um ¼ d2 þ a Gðx; y; sÞkðsÞecs um ðy; t sÞdyds bu2m @t @x2 0 1
3391
ð1:1Þ
for x 2 ð1; þ1Þ; t > 0, where Gðx; y; tÞ is defined by ðxyÞ2 1 Gðx; y; tÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi e 4d1 t : 4pd1 t
In system (1.1), ui ðx; tÞ and um ðx; tÞ represent the densities of immature and mature members, respectively. Here, a is the birth rate, assumed proportional to the number of adults and c is juvenile mortality; b represents adult mortality. The term ecs um ðy; t sÞ represents the number born at time t s and location y and is still alive at time t. G is a weighting function describing the distribution at past times of the individual of the species um who is at position x and time t. The kernel k satR1 isfies 0 kðsÞds ¼ 1. In [1], a detailed investigation of travelling front solutions connecting the extinction state and the positive steady state was carried out; much attention was paid on the minimum speed and the qualitative form of the profile, which appeared always to be monotone; a rigorous proof of existence was provided for a special, but realistic, choice of the probability distribution function representing the maturation delay. Motivated by the work of Al-Omari and Gourley in [1], in this paper, we are concerned with the global stability of each of feasible steady states and the existence of travelling wave solutions of a diffusive predator-prey model with nonlocal maturation delay. To this end, we consider the following reaction-diffusion predator-prey model
Z þ1 Z @u1 @ 2 u1 þ r1 Gðx; y; sÞkðsÞecs u1 ðy; t sÞdyds ¼ d1 2 @t @x 0 X a11 u21 ðx; tÞ a12 u1 ðx; tÞu2 ðx; tÞ;
ð1:2Þ
2
@u2 @ u2 ¼ d2 þ u2 ðx; tÞðr 2 þ a21 u1 ðx; tÞ a22 u2 ðx; tÞÞ; @t @x2 for t > 0; x 2 X, with homogeneous Neumann boundary conditions
@u1 @u2 ¼ ¼ 0; @n @n
x 2 @ X;
ð1:3Þ
and initial conditions
ui ðx; hÞ ¼ /i ðx; hÞ P 0 ði ¼ 1; 2Þ;
ðx; hÞ 2 X ð1; 0;
ð1:4Þ
where
Gðx; y; tÞ ¼
1
þ
1 2X
p p
2
ed1 n t cos nx cos ny;
n¼1
kðtÞ ¼
1
s
et=s ;
X ¼ ð0; pÞ;
G is the weighting function describing the distribution at past times of the individual of the prey u1 who is at position x and time t, and satisfies
Z
Gðx; y; tÞdy ¼ 1;
X
@G @2G ¼ d1 2 ; @t @x
Gðx; y; 0Þ ¼ dðx yÞ:
In system (1.2), u1 ðx; tÞ and u2 ðx; tÞ represent the densities of the prey and predator populations at location x and time t, respectively. The parameters d1 and d2 are the diffusion rates of the prey and predator populations, respectively; r1 is the intrinsic growth rate of the prey; r2 is the death rate of the predator; a11 is the intra-specific competition rate of the prey; a12 is the capture rate of the predator; a21 is the conversion rate of the predator by consuming prey; a22 is the intra-specific competition rate of the predator; c is defined as in (1.1). This paper is organized as follows. In Section 2, we discuss the local stability of steady states for system (1.2). In Section 3, we study the global stability of each of steady states of system (1.2). In Section 4, we are concerned with the existence of travelling waves of system (1.2) in the case X ¼ ð1; þ1Þ. Finally, numerical simulations are given to illustrate the theoretical results. 2. Local stability In this section, by analyzing the corresponding characteristic equations, we discuss the local stability of each of steady states of system (1.2). It is easy to show that system (1.2) has two steady states E00 ð0; 0Þ and E11 ða11 ðrc1sþ1Þ ; 0Þ. If the following holds
ðH1Þa21 r1 > a11 r 2 ðcs þ 1Þ;
3392
X. Zhang, R. Xu / Commun Nonlinear Sci Numer Simulat 15 (2010) 3390–3401
system (1.2) has a positive steady state E ðu1 ; u2 Þ, where
a22 r 1 þ a12 r2 ðcs þ 1Þ ; ða11 a22 þ a12 a21 Þðcs þ 1Þ
u1 ¼
u2 ¼
a21 r 1 a11 r 2 ðcs þ 1Þ : ða11 a22 þ a12 a21 Þðcs þ 1Þ
Letting t
Z p
1
0
Z
u3 ¼
Gðx; y; t sÞkðt sÞecðtsÞ u1 ðy; sÞdyds;
system (1.2) can be rewritten as
@u1 @ 2 u1 þ r 1 u3 a11 u21 a12 u1 u2 ; ¼ d1 @t @x2 @u2 @ 2 u2 þ u2 ðr 2 þ a21 u1 a22 u2 Þ; ¼ d2 @t @x2 2 @u3 @ u3 1 1 þ u1 u3 cu3 : ¼ d1 @t @x2 s s
ð2:1Þ
r1 The corresponding steady states of system (2.1) are E0 ð0; 0; 0Þ; E1 ða11 ð1þ csÞ ; 0; a
r1
u
sÞ2
11 ð1þc
Þ and E ðu1 ; u2 ; 1þ1scÞ.
Let 0 ¼ l1 < l2 < be the eigenvalues of the operator D on X with the homogeneous Neumann boundary conditions, and Eðli Þ be the eigenspace corresponding to li in C 1 ðXÞ. Let X ¼ ½C 1 ðXÞ3 ; f/ij ; j ¼ 1; . . . ; dimEðli Þg be an orthonormal basis of Eðli Þ, and Xij ¼ fc/ij jc 2 R3 g. Then 1
X ¼ Xi i¼0
and Xi ¼
dimEðli Þ
j¼1
Xij :
^1 ; u ^2 ; u ^ 3 Þ be any arbitrary steady state of system (2.1). Then Let D ¼ diagðd1 ; d2 ; d1 Þ; Lu ¼ DDu þ Gð b EÞu; b E ¼ ðu
0
^ 1 a12 u ^2 2a11 u B ^2 a21 u Gð b EÞ ¼ @
^1 a12 u ^1 2a22 u ^2 r 2 þ a21 u
r1
0
1s c
1
s
0
1 C A:
ð2:2Þ
The linearization of system (2.1) at b E is of the form ut ¼ Lu, and since Xi is invariant under the operator L for each i P 1, k EÞ. The characteristic equation of li D þ GðE0 Þ is an eigenvalue of L if and only if it is an eigenvalue of the matrix li D þ Gð b takes the form
1 1 r1 ðk þ li d2 þ r 2 Þ k2 þ 2li d1 þ þ c k þ li d1 li d1 þ þ c ¼ 0:
s
s
s
ð2:3Þ
Clearly, for i ¼ 1, Eq. (2.3) always has a positive real root. Therefore, there is a characteristic root k with positive real part in the spectrum of L. Accordingly, the trivial steady state E0 ð0; 0; 0Þ is always unstable. The characteristic equation of li D þ GðE1 Þ is
k þ li d2 þ r2
a21 r 1 1 2r1 1 2r 1 r li d1 þ 1 ¼ 0: k2 þ 2li d1 þ þ c þ k þ li d1 li d1 þ þ c þ a11 ð1 þ csÞ s 1 þ cs s 1 þ cs s
ð2:4Þ
21 r 1 If a21 r1 > a11 r2 ðcs þ 1Þ; for i ¼ 1, equation (2.4) always has a positive real root k ¼ a11að1þ csÞ r 2 , then the steady state E1 of system (2.1) is unstable. If a21 r1 < a11 r 2 ðcs þ 1Þ, we know all the roots of (2.4) have negative real parts. Therefore, the steady state E1 is stable. The characteristic equation of li D þ GðE Þ is of the form
k3 þ A2 k2 þ A1 k þ A0 ¼ 0;
ð2:5Þ
where
r1 1 þ c ðli d2 þ a22 u2 Þ > 0; þ ðli d1 þ a11 u1 Þ 1 þ cs s r 1 1 A1 ¼ li d1 li d1 þ a11 u1 þ þc þ ðli d1 þ a11 u1 Þ 1 þ cs s 1 r 1 þ ðli d2 þ a22 u2 Þ 2li d1 þ a11 u1 þ þ c þ þ a12 a21 u1 u2 > 0; s 1 þ cs 1 r1 A2 ¼ 2li d1 þ li d2 þ a11 u1 þ a12 u2 þ þ c þ > 0: s 1 þ cs A0 ¼
li d1 li d1 þ a11 u1 þ
ð2:6Þ
It is easy to see that A1 A2 > A0 . By Routh-Hurwitz criterion we know that the positive steady state E of system (2.1) is stable. In summary, we obtain the following results.
3393
X. Zhang, R. Xu / Commun Nonlinear Sci Numer Simulat 15 (2010) 3390–3401
Theorem 2.1. For system (2.1), we have
(i) The trivial steady state E0 ð0; 0Þ is always unstable. (ii) If a21 r 1 < a11 r2 ðcs þ 1Þ, the steady state E1 is stable; if a21 r1 > a11 r2 ðcs þ 1Þ, E1 is unstable. (iii) If a21 r 1 > a11 r2 ðcs þ 1Þ, the positive steady state E is always stable.
3. Global stability In this section, by using the method of upper–lower solutions and its associated monotone iteration scheme, we study the global stability of the positive steady state E and the boundary steady state E1 of system (2.1), respectively. We first investigate the asymptotic behavior of the solution of system (2.1) in relation to its corresponding constant steady state solution c1 ; c2 ; c3 of system (2.1) given by
r 1 c3 a11 c21 a12 c1 c2 ¼ 0; c2 ðr2 þ a21 c1 a22 c2 Þ ¼ 0; 1 1 c1 c3 cc3 ¼ 0:
s
ð3:1Þ
s
We deal with equations (3.1) as the following elliptic system
Dc1 ¼ r 1 c3 a11 c21 a12 c1 c2 ; Dc2 ¼ c2 ðr2 þ a21 c1 a22 c2 Þ; 1 1 Dc3 ¼ c1 c3 cc3 :
s
ð3:2Þ
s
Let ~c ¼ ð~c1 ; ~c2 ; ~c3 Þ and ^c ¼ ð^c1 ; ^c2 ; ^c3 Þ denote a pair of upper and lower solutions of system (3.2) satisfying ~c P ^c P 0 and
r 1 ~c3 a11 ~c21 a12 ~c1 ^c2 6 0 6 r 1 ^c3 a11 ^c21 a12 ^c1 ~c2 ; ~c2 ðr2 þ a21 ~c1 a22 ~c2 Þ 6 0 6 ^c2 ðr 2 þ a21 ^c1 a22 ^c2 Þ; 1 1 1 1 ~c1 ~c3 c~c3 6 0 6 ^c1 ^c3 c^c3 :
s
s
s
ð3:3Þ
s
It is easy to see that there is a positive constant K 0 ¼ K 0 ða11 ; a12 ; a21 ; a22 ; r 1 ; r 2 ; ~ci ; ^ci Þ such that the following Lipschitz condition holds
jr1 u3 a11 u21 a12 u1 u2 ðr 1 v 3 a11 v 21 a12 v 1 v 2 Þj 6 K 0 ðju1 v 1 j þ ju2 v 2 j þ ju3 v 3 jÞ; ju2 ðr 2 þ a21 u1 a22 u2 Þ v 2 ðr 2 þ a21 v 1 a22 v 2 Þj 6 K 0 ðju1 v 1 j þ ju2 v 2 jÞ; 1 u1 1 u3 cu3 ð1 v 1 1 v 3 cv 3 Þ 6 K 0 ðju1 v 1 j þ ju3 v 3 jÞ; s s s s
ð3:4Þ
for ^ci 6 ui ; v i 6 ~ci . ðmÞ ðmÞ ðmÞ ðmÞ ðmÞ ðmÞ We now construct two sequences cðmÞ ¼ ðc1 ; c2 ; c3 Þ and cðmÞ ¼ ðc1 ; c2 ; c3 Þ from the following iteration process ðm1Þ ðm1Þ ðm1Þ 2 ðm1Þ ðm1Þ cðmÞ ¼ c1 þ K10 ðr 1 c3 a11 ðc1 Þ a12 c1 c2 Þ; 1 ðm1Þ ðm1Þ ðm1Þ cðmÞ ¼ c2 þ K10 ðr2 þ a21 c1 a22 c2 Þ; 2 ðm1Þ ðm1Þ ðm1Þ ðm1Þ cðmÞ ¼ c3 þ K10 1s c1 1s c3 cc3 ; 3 ðmÞ
¼ c1
ðmÞ
¼ c2
c1 c2
ðmÞ
c3
ðm1Þ
ðm1Þ
þ K10 ðr 1 c3
ðm1Þ 2
a11 ðc1
ðm1Þ
ð3:5Þ
ðm1Þ ðm1Þ c2 Þ;
Þ a12 c1
ðm1Þ
ðm1Þ
þ K10 ðr2 þ a21 c1 a22 c2 Þ; ðm1Þ ðm1Þ 1 1 ðm1Þ 1 ðm1Þ ; ¼ c3 þ K 0 s c1 s c3 cc3
with initial data cð0Þ ¼ ~c; cð0Þ ¼ ^c, where K 0 is the Lipschitz constant above. Clearly, the sequences cðmÞ and cðmÞ possess the following property
^c 6 cðmÞ 6 cðmþ1Þ 6 cðmþ1Þ 6 cðmÞ 6 ~c;
m ¼ 1; 2; 3; . . . ;
where the ordering relation is in the componentwise sense. Therefore, the limits lim cðmÞ and lim cðmÞ exist. Denoting m!1
c ¼ lim cðmÞ ; m!1
c ¼ lim cðmÞ ; m!1
then c and c satisfy the following equations
m!1
ð3:6Þ
3394
X. Zhang, R. Xu / Commun Nonlinear Sci Numer Simulat 15 (2010) 3390–3401
r 1 c3 a11 c21 a12 c1 c2 ¼ 0 ¼ r 1 c3 a11 c21 a12 c1 c2 ; c2 ðr 2 þ a21 c1 a22 c2 Þ ¼ 0 ¼ c2 ðr2 þ a21 c1 a22 c2 Þ; 1
s
1 1 1 c1 c3 cc3 ¼ 0 ¼ c1 c3 cc3 :
s
s
ð3:7Þ
s
The constant vectors c and c are called quasi-solutions of system (3.1) in h^c; ~ci. Generally, the constant vectors c and c are not true solutions of system (3.1). However, if c ¼ c, then c (or c) is the unique solution of Eqs. (3.1) in h^c; ~ci. Next, we present a result on the relation between the solution of (2.1) and (3.1). Lemma 3.1. ([8], Theorems 2.1 and 2.2) Let ~c; ^c be a pair of upper and lower solutions of system (3.1) and let c and c be limits in (3.6) which are the quasi-solutions of (3.1) and satisfy (3.7). Then for any initial function satisfying ^ci 6 ui ðx; 0Þ 6 ~ci in ð0; pÞ for i ¼ 1; 2; 3, the corresponding solution u ¼ ðu1 ; u2 ; u3 Þ of system (2.1) possess the property
c 6 uðx; tÞ 6 c as t ! 1ðx 2 ½0; pÞ: Moreover, if c ¼ c then c (or c) is the unique solution of Eqs.(3.1) in h^c; ~ci and
lim uðx; tÞ ¼ c ðx 2 ½0; pÞ:
m!1
We now introduce two Lemmas. Lemma 3.2. (Hirsch[7]) Let F be a C 1 vector field in Rn , whose flow / preserves Rnþ for t P 0 and is strongly monotone in Rnþ . Assume that the origin is an steady state and that all trajectories in Rnþ are bounded. Suppose the matrix-valued map DF : Rn ! Rnn is strictly anti-monotone, in the sense that
if
x > y;
then DFðxÞ < DFðyÞ;
here, ‘‘<” and ‘‘>” represent the standard partial ordering in Rnþ tend to the origin, or else there is a unique steady state p 2 IntRnþ and all trajectories in Rnþ n f0g tend to p. We consider the following differential equations
u_ 1 ¼ r 1 u3 a11 u21 a12 Au1 ; 1 1 u_ 3 ¼ u1 u3 cu3 :
s
ð3:8Þ
s
System (3.8) always has a trivial equilibrium E0 ð0; 0Þ. If r 1 > a12 Að1 þ csÞ, then system (3.8) has a positive equilibrium 12 Að1þcsÞ r 1 a12 Að1þcsÞ ; a ð1þcsÞ2 Þ: E1 ðr1 a a11 ð1þcsÞ 11
By Lemma 3.2, it is easy to prove the following result. Lemma 3.3. If r 1 < a12 Að1 þ csÞ, the equilibrium E0 ð0; 0Þ of system (3.8) is globally asymptotically stable; if r1 > a12 Að1 þ csÞ, 12 Að1þcsÞ r 1 a12 Að1þcsÞ of system (3.8) is globally asymptotically stable. the equilibrium E1 r1 a 2 a11 ð1þcsÞ ; a11 ð1þcsÞ
From Lemmas 3.1–3.3, we have a result on the global stability of the positive steady state E of system (2.1). Theorem 3.1. Let ðu1 ðx; tÞ; u2 ðx; tÞ; u3 ðx; tÞÞ be a solution of system (2.1). If a11 a22 > a12 a21 and ðH1Þ hold, then
limðu1 ðx; tÞ; u2 ðx; tÞ; u3 ðx; tÞÞ ¼ u1 ; u2 ;
t!1
u1 1 þ cs
uniformly for x 2 ½0; p;
i.e., the positive steady state E of system (2.1) is globally asymptotically stable. Proof. Let ðu1 ðx; tÞ; u2 ðx; tÞ; u3 ðx; tÞÞ be a solution of system (2.1). First, we need to seek for a pair of upper and lower solutions ~c and ^c such that ~c > ^c > 0 and condition (3.3) holds. For e > 0 sufficiently small, we define
r1 a21 r1 a11 r 2 ð1 þ csÞ a21 e; þ e; ~c2 ¼ þ 1þ a22 a11 ð1 þ csÞ a11 ð1 þ csÞ a22 r1 ða11 a22 a12 a21 Þ þ a11 a12 r 2 ð1 þ csÞ a12 a12 a21 ¼ þ e; 1þ 2 a11 a11 a22 a11 a22 ð1 þ csÞ ða11 a22 a12 a21 Þðr 1 a21 a11 r 2 ð1 þ csÞÞ a21 a12 a21 a12 a221 ¼ þ þ e; 1 þ a22 a11 a22 a11 a222 a211 a222 ð1 þ csÞ ~c1 ^c1 ¼ ; ^c3 ¼ : 1 þ cs 1 þ cs
~c1 ¼ ^c1 ^c2 ~c3 Choose
ð3:9Þ
e > 0 sufficiently small satisfying
e < min
r 1 ða11 a22 a12 a21 Þ þ a11 a12 r 2 ð1 þ csÞ ða11 a22 a12 a21 Þðr 1 a21 a11 r 2 ð1 þ csÞÞ ; : a11 ð1 þ csÞða11 a22 þ a12 a22 þ a12 a21 Þ a11 ð1 þ csÞða11 a22 þ a12 a21 Þða22 þ a21 Þ
ð3:10Þ
X. Zhang, R. Xu / Commun Nonlinear Sci Numer Simulat 15 (2010) 3390–3401
3395
It follows from (3.9) and (3.10) that the pair ~c ¼ ð~c1 ; ~c2 ; ~c3 Þ and ^c ¼ ð^c1 ; ^c2 ; ^c3 Þ satisfy
r 1 ~c3 a11 ~c21 a12 ~c1 ^c2 < a11 ~c1 e < 0; r 2 þ a21 ~c1 a22 ~c2 ¼ a22 e < 0; 1
s
1 ~c1 ~c3 c~c3 ¼ 0;
s
ð3:11Þ
r 1 ^c3 a11 ^c21 a12 ^c1 ~c2 ¼ a11 ^c1 e > 0; r 2 þ a21 ^c1 a22 ^c2 ¼ a22 e > 0; 1
s
1 ^c1 ^c3 c^c3 ¼ 0:
s
This shows that the condition (3.3) holds. Therefore, ~c ¼ ð~c1 ; ~c2 ; ~c3 Þ and ^c ¼ ð^c1 ; ^c2 ; ^c3 Þ are the coupled upper and lower solutions of (3.2). Second, we claim that there exists a t such that
^ci 6 ui ðx; tÞ 6 ~ci ;
ðx; tÞ 2 ½0; p ½t ; þ1Þ;
i ¼ 1; 2; 3;
ð3:12Þ
where ð^c1 ; ^c2 ; ^c3 Þ and ð~c1 ; ~c2 ; ~c3 Þ are defined in (3.9). It follows from the first and third equations of system (2.1) that
@u1 @ 2 u1 þ r1 u3 a11 u21 ; 6 d1 @t @x2 @u3 @ 2 u3 1 1 þ u1 u3 cu3 : ¼ d1 @t @x2 s s
ð3:13Þ
Let ðu1 ðtÞ; u3 ðtÞÞ be the unique positive solution of the following system
u_ 1 ¼ r 1 u3 a11 u21 ; 1 1 u_ 3 ¼ u1 u3 cu3 :
s
ð3:14Þ
s
By Lemma 3.3, it follows that
lim u1 ðtÞ ¼
t!1
For
r1 ; a11 ð1 þ csÞ
lim u3 ðtÞ ¼
t!1
r1 a11 ð1 þ csÞ2
:
ð3:15Þ
e > 0 sufficiently small satisfying (3.10), by comparison we derive from (3.13)–(3.15) that there is a t1 > 0 such that u1 ðx; tÞ 6 u3 ðx; tÞ 6
r1 þ ¼ ~c1 ; a11 ð1þcsÞ r1 þ 1þecs ¼ a11 ð1þcsÞ2
e
~c3 ;
ðx; tÞ 2 ½0; p ½t 1 ; þ1Þ:
ð3:16Þ
It follows from the second equation of (2.1) that
@u2 @ 2 u2 þ u2 ðr 2 þ a21 ~c1 a22 u2 Þ: 6 d2 @t @x2 By comparison we derive that there is a t2 > t1 such that
u2 ðx; tÞ 6
a21 ~c1 r2 þ e ¼ ~c2 a22
ðx; tÞ 2 ½0; p ½t 2 ; þ1Þ:
ð3:17Þ
It follows from (3.17) and (2.1) that
@u1 @ 2 u1 þ r1 u3 a11 u21 a12 u1 ~c2 ; P d1 @t @x2 @u3 @ 2 u3 1 1 þ u1 u3 cu3 : ¼ d1 @t @x2 s s
ð3:18Þ
By comparison we derive that there is a t3 > t2 such that
u1 ðx; tÞ P
r 1 a12 ~c2 ð1þcsÞ a11 ð1þcsÞ
e ¼ ^c1 ;
u3 ðx; tÞ P
r 1 a12 ~c2 ð1þcsÞ a11 ð1þcsÞ2
1þecs ¼ ^c3 ;
ðx; tÞ 2 ½0; p ½t 3 ; þ1Þ:
It follows from the second equation of (2.1) that
@u2 @ 2 u2 P d2 þ u2 ðr 2 þ a21 ^c1 a22 u2 Þ: @t @x2 By comparison we derive that there is a t > t 3 such that
ð3:19Þ
3396
X. Zhang, R. Xu / Commun Nonlinear Sci Numer Simulat 15 (2010) 3390–3401
u2 ðx; tÞ P
a21 ^c1 r2 e ¼ ^c2 a22
ðx; tÞ 2 ½0; p ½t ; þ1Þ:
ð3:20Þ
Hence, (3.12) holds. u Last, we prove that the nonnegative solutions ðu1 ðx; tÞ; u2 ðx; tÞ; u3 ðx; tÞÞ of (2.1) converges uniformly to ðu1 ; u2 ; 1þ1csÞ. Since the asymptotic behavior of the solution of system (2.1) is the same as if we set the initial time at t ¼ t , without loss of generality, we assume c^i 6 /1 ðx; tÞ 6 ~ci ði ¼ 1; 2; 3Þ. From Lemma 3.1, there exist quasi-solution c ¼ ðc1 ; c2 ; c3 Þ and c ¼ ðc1 ; c2 ; c3 Þ such that the unique solution ðu1 ; u2 ; u3 Þ of (2.1) satisfies
0 < ^ci 6 ci 6 ui ðx; tÞ 6 ci 6 ~ci ; t ! 1ðx 2 ½0; pÞ;
i ¼ 1; 2; 3;
ð3:21Þ
and c and c satisfy (3.7). Therefore (3.7) is reduced to
r1 r1 a11 c1 a12 c2 ¼ 0 ¼ a11 c1 a12 c2 ; 1 þ cs 1 þ cs r 2 þ a21 c1 a22 c2 ¼ 0 ¼ r 2 þ a21 c1 a22 c2 ; 1 1 1 1 c1 c3 cc3 ¼ 0 ¼ c1 c3 cc3 :
s
s
s
ð3:22Þ
s
It follows from the first two equations of (3.22) that
a11 ðc1 c1 Þ a12 ðc2 c2 Þ ¼ 0; a21 ðc1 c1 Þ a22 ðc2 c2 Þ ¼ 0:
ð3:23Þ
Noting that a11 a22 > a12 a21 , we derive from (3.23) that c1 ¼ c1 ; c2 ¼ c2 . It therefore follows from (3.22) that c3 ¼ c3 . Applying Lemma 3.1 we have that
limðu1 ðx; tÞ; u2 ðx; tÞ; u3 ðx; tÞÞ ¼ u1 ; u2 ;
t!1
The proof is completed.
u1 1 þ cs
uniformly for x 2 ½0; p:
h
Theorem 3.2. Let ðu1 ; u2 ; u3 Þ be a solution of system (2.1). If a21 r1 < a11 r2 ðcs þ 1Þ, then
limðu1 ; u2 ; u3 Þ ¼
t!1
r1 r1 ; 0; a11 ð1 þ csÞ a11 ð1 þ csÞ2
! uniformly for
x 2 ½0; p;
i.e., the boundary steady state E1 of system (2.1) is globally asymptotically stable. Proof. Let ðu1 ; u2 ; u3 Þ be a solution of system (2.1). For
e < min
e > 0 sufficiently small satisfying
r1 r2 a11 ð1 þ csÞ a21 r1 ; ; ða12 þ a11 Þð1 þ csÞ a21 a11 ð1 þ csÞ
ð3:24Þ
we define r1 ~c1 ¼ a ð1þ csÞ þ e; 11
~c2 ¼ e;
~c3 ¼ 1þ~c1cs ;
r1 a12 e ^c1 ¼ a ð1þ csÞ a11 e; 11
^c2 ¼ 0;
^c3 ¼ 1þ^c1cs :
ð3:25Þ
It follows from (3.24) and (3.25) that the pair ~c ¼ ð~c1 ; ~c2 ; ~c3 Þ and ^c ¼ ð^c1 ; ^c2 ; ^c3 Þ satisfy
r1 ~c3 a11 ~c21 a12 ~c1 ^c2 ¼ a11 ~c1 e < 0; 21 r1 r 2 þ a21 ~c1 a22 ~c2 ¼ a að1þ csÞ r 2 a22 e þ a21 e < 0; 11
1~ s c1
1s ~c3 c~c3 ¼ 0; r1 ^c3 a11 ^c21 a12 ^c1 ~c2 ¼ a11 ^c1 e > 0; ^c2 ðr 2 þ a21 ^c1 a22 ^c2 Þ ¼ 0; 1^ 1^ ^ s c1 s c3 cc3 ¼ 0:
ð3:26Þ
This shows that the condition (3.3) holds. Therefore ~c ¼ ð~c1 ; ~c2 ; ~c3 Þ and ^c ¼ ð^c1 ; ^c2 ; ^c3 Þ are the coupled upper and lower solutions of (3.2). From Lemma 3.1, there exist quasi-solution c ¼ ðc1 ; c2 ; c3 Þ and c ¼ ðc1 ; c2 ; c3 Þ such that the unique solution ðu1 ; u2 ; u3 Þ of (2.1) satisfies
0 < ^ci 6 ci 6 ui 6 ci 6 ~ci ;
t ! 1;
i ¼ 1; 2; 3;
and c and c satisfy (3.7). Therefore (3.7) is reduced to
ð3:27Þ
3397
X. Zhang, R. Xu / Commun Nonlinear Sci Numer Simulat 15 (2010) 3390–3401
r1 r1 a11 c1 a12 c2 ¼ 0 ¼ a11 c1 a12 c2 ; 1 þ cs 1 þ cs c2 ðr 2 þ a21 c1 a22 c2 Þ ¼ 0 ¼ c2 ðr 2 þ a21 c1 a22 c2 Þ; 1 1 1 1 c1 c3 cc3 ¼ 0 ¼ c1 c3 cc3 :
s
s
s
ð3:28Þ
s
If c2 > 0, it then follows from the first two equations of (3.28) that r1 1þcs
a11 c1 a12 c2 ¼ 0;
ð3:29Þ
r 2 þ a21 c1 a22 c2 ¼ 0; which yields
a21 r 1 a11 r 2 ¼ a12 a21 c2 þ a11 a22 c2 : 1 þ cs Noting that a21 r 1 < a11 r2 ðcs þ 1Þ, a contradiction occurs. Hence, r1 r1 c1 ¼ c1 ¼ a ð1þ csÞ ; c3 ¼ c 3 ¼ a ð1þcsÞ2 . Applying Lemma 3.1 we have that 11 11
limðu1 ; u2 ; u3 Þ ¼
t!1
r1 r1 ; 0; a11 ð1 þ csÞ a11 ð1 þ csÞ2
c2 ¼ c2 ¼ 0.
We
derive
from
(3.28)
that
!
uniformly for x 2 ½0; p:
This completes the proof. h
4. The existence of travelling waves In this section, we discuss the existence of travelling waves of system (1.2) for the case of an infinite one-dimensional domain X ¼ ð1; þ1Þ by using the geometric singular perturbation theory. The infinite domain case is from a modelling point of view since there are no boundaries for individuals to interact with as they drift from their past to their present positions. The model assumes the form
@u1 @ 2 u1 þ r1 ¼ d1 @t @x2
Z
t
1
Z
þ1
Gðx y; t sÞkðt sÞecðtsÞ u1 ðy; sÞdyds
1
a11 u21 a12 u1 u2 ;
ð4:1Þ
2
@u2 @ u2 þ u2 ðr2 þ a21 u1 a22 u2 Þ: ¼ d2 @t @x2 where 2 1 x Gðx; tÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi e 4d1 t ; 4pd1 t
and G satisfies
@G @2G ¼ d1 2 ; @t @x
Gðx; 0Þ ¼ dðxÞ:
System (4.1) still has steady states E0 ð0; 0Þ; E1
r1 ;0 a11 ð1þcsÞ
and E ðu1 ; u2 Þ. Our interest in this section is the possibility of a tran-
sition between the steady states E1 and E in the form of a travelling wave solution. Theorem 4.1. For s sufficiently small, if a21 r 1 > a11 r2 ðcs þ 1Þ, then system (4.1) possesses a travelling wave solution connecting the steady states E1 and E . Proof. Denoting
u3 ¼
Z
t
1
Z
þ1
1
ðxyÞ2 1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e 4d1 ðtsÞ kðt sÞecðtsÞ u1 ðy; sÞdyds; 4pd1 ðt sÞ
system (4.1) can be rewritten as
@u1 @ 2 u1 þ r 1 u3 a11 u21 a12 u1 u2 ; ¼ d1 @t @x2 @u2 @ 2 u2 þ u2 ðr2 þ a21 u1 a22 u2 Þ; ¼ d2 @t @x2 2 @u3 @ u3 1 1 ¼ d1 þ u1 u3 cu3 : @t @x2 s s
ð4:2Þ
3398
X. Zhang, R. Xu / Commun Nonlinear Sci Numer Simulat 15 (2010) 3390–3401
Converting to travelling wave form, by writing
u1 ðx; tÞ ¼ u1 ðzÞ;
u2 ðx; tÞ ¼ u2 ðzÞ;
u3 ðx; tÞ ¼ u3 ðzÞ;
z ¼ x þ ct;
we derive that
€ 1 þ r 1 u3 a11 u21 a12 u1 u2 ; cu_ 1 ¼ d1 u € 2 þ u2 ðr 2 þ a21 u1 a22 u2 Þ; cu_ 2 ¼ d2 u 1 1 € 3 þ u1 u3 cu3 ; cu_ 3 ¼ d1 u
s
ð4:3Þ
s
where prime denotes differentiation with respect to z. We introduce
v 1 ¼ d1 u_ 1 ; v 2 ¼ d2 u_ 2 ; v 3 ¼ d1 u_ 3 : And, we replace
s with e2 s. System (4.3) becomes
1 v 1; d1 c v_ 1 ¼ v 1 ðr1 u3 a11 u21 a12 u1 u2 Þ; d1 1 u_ 2 ¼ v 2 ; d2 c v_ 2 ¼ v 2 u2 ðr2 þ a21 u1 a22 u2 Þ; d2 1 u_ 3 ¼ v 3 ; d1 e2 c 1 1 e2 v_ 3 ¼ v 3 u1 u3 þ e2 cu3 : d1 s s u_ 1 ¼
ð4:4Þ
By introducing new state variables
~ 1 ¼ u1 ; u
v~ 1 ¼ v 1 ;
~ 2 ¼ u2 ; u
v~ 2 ¼ v 2 ;
~ 3 ¼ u3 ; u
v~ 3 ¼ ev 3
and dropping the tildes, we have
1 v 1; d1 c v_ 1 ¼ v 1 ðr1 u3 a11 u21 a12 u1 u2 Þ; d1 1 u_ 2 ¼ v 2 ; d2 c v_ 2 ¼ v 2 u2 ðr2 þ a21 u1 a22 u2 Þ; d2 1 eu_ 3 ¼ v 3 ; d1 ec 1 1 ev_ 3 ¼ v 3 u1 u3 þ e2 cu3 : d1 s s u_ 1 ¼
ð4:5Þ
This system is called slow system. Introducing a new independent variable g defined by z ¼ eg, system (4.5) transforms into
1 v1; d1 c v_ 1 ¼ e v 1 ðr1 u3 a11 u21 a12 u1 u2 Þ ; d1 1 u_ 2 ¼ e v 2 ; d 2 c v_ 2 ¼ e v 2 u2 ðr2 þ a21 u1 a22 u2 Þ ; d2 1 u_ 3 ¼ v 3 ; d1 ec 1 1 v_ 3 ¼ v 3 u1 u3 þ e2 cu3 : d1 s s
u_ 1 ¼ e
ð4:6Þ
X. Zhang, R. Xu / Commun Nonlinear Sci Numer Simulat 15 (2010) 3390–3401
3399
This system is called fast system. In the slow system (4.5), if e ¼ 0, then the flow of that system is confined to the set
M0 ¼ fðu1 ; v 1 ; u2 ; v 2 ; u3 ; v 3 Þ 2 R6 : v 3 ¼ 0; u3 ¼ u1 g which is a four-dimensional manifold for system (4.5). By Theorem 9.1 [3], we know that if M 0 is normally hyperbolic, then for sufficiently small e > 0 we can obtain an invariant manifold of system (4.5). To verify normal hyperbolicity, we need to verify that the linearization of the fast system (4.6), restricted to M 0 , has precisely dim M0 eigenvalues on the imaginary axis with the remainder of the spectrum being hyperbolic. The linearization of the fast system (4.6) restricted to M0 is given by the matrix
0 B B B B B B B B B @
0 0 0 0 0 1s
1 0 0C C C 0 0 0 0 0C C C 0 0 0 0 0C C 1 C 0 0 0 0 d1 A 0 0 0 1s 0 0 0 0 0 0 0 0 0
pffiffiffiffiffiffiffiffi which has six eigenvalues f0; 0; 0; 0; 1= d1 sg. Thus, the manifold M 0 is normally hyperbolic. By the geometric singular perturbation theory, we know that the slow system (4.5) has a invariant manifold Me , which can be written as
Me ¼ fðu1 ; v 1 ; u2 ; v 2 ; u3 ; v 3 Þ 2 R6 : v 3 ¼ hðu1 ; v 1 ; u2 ; v 2 ; eÞ; u3 ¼ u1 þ gðu1 ; v 1 ; u2 ; v 2 ; eÞg where the functions h and g satisfy
hðu1 ; v 1 ; u2 ; v 2 ; 0Þ ¼ gðu1 ; v 1 ; u2 ; v 2 ; 0Þ ¼ 0: Thus the functions h and g can be expanded into the form of a Taylor series about e:
hðu1 ; v 1 ; u2 ; v 2 ; eÞ ¼ eh1 þ e2 h2 þ ;
ð4:7Þ
gðu1 ; v 1 ; u2 ; v 2 ; eÞ ¼ eg 1 þ e2 g 2 þ :
We determine the functions h and g. Note that M e is the invariant manifold for the flow of system (4.5). From w2 ¼ h; w1 ¼ u1 þ g and (4.5), we have
@g @g @g @g v_ 1 þ v_ 2 ; u_ 1 þ u_ 2 þ @u2 @v 1 @v 2 d1 @u1 e @h @h @h @h c v_ 1 þ v_ 2 h ecu1 g¼ u_ 1 þ u_ 2 þ 2 @u2 @v 1 @v 2 d1 1=s þ e c @u1 c 3 @h @h @h @h c ¼ ðe 6 4 e þ Þ v_ 1 þ v_ 2 h ecu1 : u_ 1 þ u_ 2 þ @u1 @u2 @v 1 @v 2 d1 s h ¼ ed1
v1
þ
Substituting (4.7) into (4.8) and comparing coefficients of
h1 ¼ v 1 ;
h2 ¼ 0;
g 1 ¼ 0;
ð4:8Þ
e and e2 , we get
g 2 ¼ r 1 u1 þ a11 u21 þ a12 u1 u2 cu1 :
ð4:9Þ
Then the slow system (2.4) restricted to Me is
1 v 1; d1 c v_ 1 ¼ v 1 u1 Mðu1 ; u2 Þ; d1 1 u_ 2 ¼ v 2 ; d2 c v_ 2 ¼ v 2 u2 Nðu1 ; u2 Þ; d2 u_ 1 ¼
ð4:10Þ
where
Mðu1 ; u2 Þ ¼ r 1 ð1 r 1 Þ a11 ð1 r 1 Þu1 a12 ð1 r 1 Þu2 r 1 c; Nðu1 ; u2 Þ ¼ r 2 þ a21 u1 a22 u2 : System (4.10) still possesses steady state points E1
r1 ; 0; 0; 0 a11 ð1þcsÞ
and E ðu1 ; 0; u2 ; 0Þ. Therefore, the results in Gardner [4]
are applicable, yielding a heteroclinic connection between the steady states E1 and E of (4.10). Thus, system (4.2) possesses travelling wave solution connecting the steady states E1 and E . This completes the proof. h
3400
X. Zhang, R. Xu / Commun Nonlinear Sci Numer Simulat 15 (2010) 3390–3401
5. Numerical simulations In this section, we give numerical simulations to illustrate the results of Section 4. In system (4.2), we let d1 ¼ 1; d 2 ¼ 2; a11 ¼ a12 ¼ a 22 ¼ 1; a21 ¼ 2; r 1 ¼ r 2 ¼ 1; c ¼ 1; s ¼ 0:1: Then system (4.2) admits 1 1 0:7 0:3 0:7 ; 0; 1:21 ; 1:1 ; 1:21 . By Theorem 4.1, we know that system (4.2) possesses travelling and E 1:1 steady states E1 ð0; 0; 0Þ; E2 1:1 wave solution connecting the steady states E1 and E (see Fig. 1). Next, we consider the effect of delay s on the wave speed. In system (4.2), we let d1 ¼ 1; d2 ¼ 2; a11 ¼ a12 ¼ a22 ¼ 1; a21 ¼ 2; r1 ¼ r 2 ¼ 1; c ¼ 1; s ¼ 0:5: Then the steady states of system (4.2) are
1 1 E1 ð0; 0; 0Þ; E2 23 ; 0; 2:25 and E 59 ; 19 ; 2:7 . Fig. 2 shows the temporal solution of system (4.2). In system (4.2), we let d1 ¼ 1; d2 ¼ 2; a11 ¼ a12 ¼ a22 ¼ 1; a21 ¼ 2; r1 ¼ r2 ¼ 1; c ¼ 1; s ¼ 0:8: Then the steady states of sys 5
1 ; 1 ; 0:7 . Fig. 3 shows the temporal solution of system (4.2). tem (4.2) are E1 ð0; 0; 0Þ; E2 9 ; 0; 3:24 and E 1:4 2:7 27 2:43 From Figs. 1–3, we can get the section drawing of the solutions u1 and u2 of system (4.2) (see Fig. 4). From Fig. 4, we see that as the delay s increases, the wave speed becomes more and more slower.
Fig. 1. The temporal solution found by numerical integration with
s ¼ 0:1.
Fig. 2. The temporal solution found by numerical integration with
s ¼ 0:5.
Fig. 3. The temporal solution found by numerical integration with
s ¼ 0:8.
3401
X. Zhang, R. Xu / Commun Nonlinear Sci Numer Simulat 15 (2010) 3390–3401 0.35
0.95 tau=0.1 tau=0.5 tau=0.8
0.9
0.25
0.8
solutions − u2
solutions − u1
0.85
tau=0.1 tau=0.5 tau=0.8
0.3
0.75 0.7 0.65
0.2 0.15 0.1
0.6 0.05
0.55 0.5
0
20
40
t − axis
60
80
100
0 0
20
40
t − axis
60
80
100
Fig. 4. The temporal solution found by numerical integration.
Acknowledgements This work was supported by the National Natural Science Foundation of China (No. 10671209). References [1] Al-Omari JFM, Gourley SA. Dynamics of a stage-structured population model incorporating a state-dependent maturation delay. Nonlinear Anal RWA 2005;6:13–33. [2] Britton NF. Spatial structures and periodic travelling waves in an integrodifferential reaction-diffusion population model. SIAM J Appl Math 1990;50:1663–88. [3] Fenichel N. Geometric singular perturbation theory for ordinary differential equations. J Differ Equat 1979;31:53–98. [4] Gardner R. Existence and stability of travelling wave solutions of competition models: a degree theoretic approch. J Differ Equat 1982;44:343–64. [5] Gourley SA, Britton NF. A predator-prey reaction-diffusion system with nonlocal effects. J Math Biol 1996;34:297–333. [6] Gourley SA, Ruan S. Convergence and travelling fronts in functional equations with nonlocal terms: a competition model. SIAM J Math Anal 2003;35:806–22. [7] Hirsch MW. The dynamical systems approach to differential equations. Bull Amer Math Soc 1984;11:1–64. [8] Pao CV. Convergence of solutions of reaction-diffusion systems with time delays. Nonlinear Anal TMA 2002;48:349–62. [9] Thieme HR, Zhao XQ. A non-local delayed and diffusive predator-prey model. Nonlinear Anal RWA 2001;2:145–60. [10] Wang ZC, Li WT, Ruan S. Travelling wave fronts of reaction-diffusion systems with spatio-temporal delays. J Differ Equat 2006;222:185–232. [11] Wang ZC, Li WT. Monotone travelling fronts of a food-limited population model with nonlocal delay. Nonlinear Anal RWA 2007;8:699–712. [12] Wu J. Theory and applications of partial functional differential equations. New York: Springer Verlag; 1996. [13] Xu R. A reaction-diffusion predator-prey model with stage structure and nonlocal delay. Appl Math Comput 2006;175:984–1006. [14] Xu R, Ma Z. Global stability of a reaction-diffusion predator-prey model with a nonlocal delay. Math Comput Model 2009;50:94–206. [15] Zhang GB, Li WT, Lin G. Traveling waves in delayed predator-prey systems with nonlocal diffusion and stage structure. Math Comput Model 2009;49:1021–9. [16] Zhang JM, Peng YH. Travelling waves of the diffusive Nicholson’s blowflioes equation with strong genetric delay kernel and non-local effect. Nonlinear Anal 2008;68:1263–70.
Contents lists available at ScienceDirect
Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns
Global stability and travelling waves of a predator-prey model with diffusion and nonlocal maturation delay Xiao Zhang, Rui Xu * Institute of Applied Mathematics, Shijiazhuang Mechanical Engineering College, Shijiazhuang 050003, PR China
a r t i c l e
i n f o
Article history: Received 10 June 2009 Accepted 29 December 2009 Available online 4 January 2010 AMS subject classification: 35K57 35K55 92D25
a b s t r a c t In this paper, a diffusive predator-prey system with nonlocal maturation delay is investigated. By analyzing the corresponding characteristic equations, the local stability of each of uniform steady states of the system is discussed. Sufficient conditions are derived for the global stability of the positive steady state and the semi-trivial steady state of the system by using the method of upper–lower solutions and its associated monotone iteration scheme, respectively. The existence of travelling wave solution of the system is established by using the geometric singular perturbation theory. Numerical simulations are carried out to illustrate the theoretical results. Ó 2010 Elsevier B.V. All rights reserved.
Keywords: Diffusion Nonlocal delay Global stability Travelling wave Geometric singular perturbation
1. Introduction Predator-prey model is one of the basic models between different species in nature. In traditional predator-prey models the spatial content of the environment has been ignored. These models have been traditionally formulated in relation to the time evolution of uniform population distributions in the habitat. However, in many ecological systems, the species under consideration may disperse spatially as well as evolving in time ([12]). This spatial dispersal mainly due to resource limitation: in regions of high population density, food will become scarce, and individuals will tend to migrate to regions of lower population density. The effect of dispersion of a population in a bounded habitat has been researched by many authors, and in this situation the governing equations for the population densities are described by a system of reaction-diffusion equations. In realistic ecological models, any delays should be spatially inhomogeneous, that is, the delays affect both the temporal and spatial variables. This is due to the fact that any given individual may not necessarily have been at the same spatial location at precious times. Such delays are called spatio-temporal delays or nonlocal delays. In recent years, the effect of nonlocal delays on the dynamics of ecological models has been taken into account (see, for example, [2,5,6,9–16]). In [1], Al-Omari and Gourley studied the following single species population with stage structure and nonlocal maturation delay
* Corresponding author. E-mail addresses: [email protected] (X. Zhang), [email protected] (R. Xu). 1007-5704/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.cnsns.2009.12.031
X. Zhang, R. Xu / Commun Nonlinear Sci Numer Simulat 15 (2010) 3390–3401
Z þ1 Z þ1 @ui @ 2 ui Gðx; y; sÞkðsÞecs um ðy; t sÞdyds; ¼ d1 2 þ aum cui a @t @x 0 1 Z þ1 Z þ1 @um @ 2 um ¼ d2 þ a Gðx; y; sÞkðsÞecs um ðy; t sÞdyds bu2m @t @x2 0 1
3391
ð1:1Þ
for x 2 ð1; þ1Þ; t > 0, where Gðx; y; tÞ is defined by ðxyÞ2 1 Gðx; y; tÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi e 4d1 t : 4pd1 t
In system (1.1), ui ðx; tÞ and um ðx; tÞ represent the densities of immature and mature members, respectively. Here, a is the birth rate, assumed proportional to the number of adults and c is juvenile mortality; b represents adult mortality. The term ecs um ðy; t sÞ represents the number born at time t s and location y and is still alive at time t. G is a weighting function describing the distribution at past times of the individual of the species um who is at position x and time t. The kernel k satR1 isfies 0 kðsÞds ¼ 1. In [1], a detailed investigation of travelling front solutions connecting the extinction state and the positive steady state was carried out; much attention was paid on the minimum speed and the qualitative form of the profile, which appeared always to be monotone; a rigorous proof of existence was provided for a special, but realistic, choice of the probability distribution function representing the maturation delay. Motivated by the work of Al-Omari and Gourley in [1], in this paper, we are concerned with the global stability of each of feasible steady states and the existence of travelling wave solutions of a diffusive predator-prey model with nonlocal maturation delay. To this end, we consider the following reaction-diffusion predator-prey model
Z þ1 Z @u1 @ 2 u1 þ r1 Gðx; y; sÞkðsÞecs u1 ðy; t sÞdyds ¼ d1 2 @t @x 0 X a11 u21 ðx; tÞ a12 u1 ðx; tÞu2 ðx; tÞ;
ð1:2Þ
2
@u2 @ u2 ¼ d2 þ u2 ðx; tÞðr 2 þ a21 u1 ðx; tÞ a22 u2 ðx; tÞÞ; @t @x2 for t > 0; x 2 X, with homogeneous Neumann boundary conditions
@u1 @u2 ¼ ¼ 0; @n @n
x 2 @ X;
ð1:3Þ
and initial conditions
ui ðx; hÞ ¼ /i ðx; hÞ P 0 ði ¼ 1; 2Þ;
ðx; hÞ 2 X ð1; 0;
ð1:4Þ
where
Gðx; y; tÞ ¼
1
þ
1 2X
p p
2
ed1 n t cos nx cos ny;
n¼1
kðtÞ ¼
1
s
et=s ;
X ¼ ð0; pÞ;
G is the weighting function describing the distribution at past times of the individual of the prey u1 who is at position x and time t, and satisfies
Z
Gðx; y; tÞdy ¼ 1;
X
@G @2G ¼ d1 2 ; @t @x
Gðx; y; 0Þ ¼ dðx yÞ:
In system (1.2), u1 ðx; tÞ and u2 ðx; tÞ represent the densities of the prey and predator populations at location x and time t, respectively. The parameters d1 and d2 are the diffusion rates of the prey and predator populations, respectively; r1 is the intrinsic growth rate of the prey; r2 is the death rate of the predator; a11 is the intra-specific competition rate of the prey; a12 is the capture rate of the predator; a21 is the conversion rate of the predator by consuming prey; a22 is the intra-specific competition rate of the predator; c is defined as in (1.1). This paper is organized as follows. In Section 2, we discuss the local stability of steady states for system (1.2). In Section 3, we study the global stability of each of steady states of system (1.2). In Section 4, we are concerned with the existence of travelling waves of system (1.2) in the case X ¼ ð1; þ1Þ. Finally, numerical simulations are given to illustrate the theoretical results. 2. Local stability In this section, by analyzing the corresponding characteristic equations, we discuss the local stability of each of steady states of system (1.2). It is easy to show that system (1.2) has two steady states E00 ð0; 0Þ and E11 ða11 ðrc1sþ1Þ ; 0Þ. If the following holds
ðH1Þa21 r1 > a11 r 2 ðcs þ 1Þ;
3392
X. Zhang, R. Xu / Commun Nonlinear Sci Numer Simulat 15 (2010) 3390–3401
system (1.2) has a positive steady state E ðu1 ; u2 Þ, where
a22 r 1 þ a12 r2 ðcs þ 1Þ ; ða11 a22 þ a12 a21 Þðcs þ 1Þ
u1 ¼
u2 ¼
a21 r 1 a11 r 2 ðcs þ 1Þ : ða11 a22 þ a12 a21 Þðcs þ 1Þ
Letting t
Z p
1
0
Z
u3 ¼
Gðx; y; t sÞkðt sÞecðtsÞ u1 ðy; sÞdyds;
system (1.2) can be rewritten as
@u1 @ 2 u1 þ r 1 u3 a11 u21 a12 u1 u2 ; ¼ d1 @t @x2 @u2 @ 2 u2 þ u2 ðr 2 þ a21 u1 a22 u2 Þ; ¼ d2 @t @x2 2 @u3 @ u3 1 1 þ u1 u3 cu3 : ¼ d1 @t @x2 s s
ð2:1Þ
r1 The corresponding steady states of system (2.1) are E0 ð0; 0; 0Þ; E1 ða11 ð1þ csÞ ; 0; a
r1
u
sÞ2
11 ð1þc
Þ and E ðu1 ; u2 ; 1þ1scÞ.
Let 0 ¼ l1 < l2 < be the eigenvalues of the operator D on X with the homogeneous Neumann boundary conditions, and Eðli Þ be the eigenspace corresponding to li in C 1 ðXÞ. Let X ¼ ½C 1 ðXÞ3 ; f/ij ; j ¼ 1; . . . ; dimEðli Þg be an orthonormal basis of Eðli Þ, and Xij ¼ fc/ij jc 2 R3 g. Then 1
X ¼ Xi i¼0
and Xi ¼
dimEðli Þ
j¼1
Xij :
^1 ; u ^2 ; u ^ 3 Þ be any arbitrary steady state of system (2.1). Then Let D ¼ diagðd1 ; d2 ; d1 Þ; Lu ¼ DDu þ Gð b EÞu; b E ¼ ðu
0
^ 1 a12 u ^2 2a11 u B ^2 a21 u Gð b EÞ ¼ @
^1 a12 u ^1 2a22 u ^2 r 2 þ a21 u
r1
0
1s c
1
s
0
1 C A:
ð2:2Þ
The linearization of system (2.1) at b E is of the form ut ¼ Lu, and since Xi is invariant under the operator L for each i P 1, k EÞ. The characteristic equation of li D þ GðE0 Þ is an eigenvalue of L if and only if it is an eigenvalue of the matrix li D þ Gð b takes the form
1 1 r1 ðk þ li d2 þ r 2 Þ k2 þ 2li d1 þ þ c k þ li d1 li d1 þ þ c ¼ 0:
s
s
s
ð2:3Þ
Clearly, for i ¼ 1, Eq. (2.3) always has a positive real root. Therefore, there is a characteristic root k with positive real part in the spectrum of L. Accordingly, the trivial steady state E0 ð0; 0; 0Þ is always unstable. The characteristic equation of li D þ GðE1 Þ is
k þ li d2 þ r2
a21 r 1 1 2r1 1 2r 1 r li d1 þ 1 ¼ 0: k2 þ 2li d1 þ þ c þ k þ li d1 li d1 þ þ c þ a11 ð1 þ csÞ s 1 þ cs s 1 þ cs s
ð2:4Þ
21 r 1 If a21 r1 > a11 r2 ðcs þ 1Þ; for i ¼ 1, equation (2.4) always has a positive real root k ¼ a11að1þ csÞ r 2 , then the steady state E1 of system (2.1) is unstable. If a21 r1 < a11 r 2 ðcs þ 1Þ, we know all the roots of (2.4) have negative real parts. Therefore, the steady state E1 is stable. The characteristic equation of li D þ GðE Þ is of the form
k3 þ A2 k2 þ A1 k þ A0 ¼ 0;
ð2:5Þ
where
r1 1 þ c ðli d2 þ a22 u2 Þ > 0; þ ðli d1 þ a11 u1 Þ 1 þ cs s r 1 1 A1 ¼ li d1 li d1 þ a11 u1 þ þc þ ðli d1 þ a11 u1 Þ 1 þ cs s 1 r 1 þ ðli d2 þ a22 u2 Þ 2li d1 þ a11 u1 þ þ c þ þ a12 a21 u1 u2 > 0; s 1 þ cs 1 r1 A2 ¼ 2li d1 þ li d2 þ a11 u1 þ a12 u2 þ þ c þ > 0: s 1 þ cs A0 ¼
li d1 li d1 þ a11 u1 þ
ð2:6Þ
It is easy to see that A1 A2 > A0 . By Routh-Hurwitz criterion we know that the positive steady state E of system (2.1) is stable. In summary, we obtain the following results.
3393
X. Zhang, R. Xu / Commun Nonlinear Sci Numer Simulat 15 (2010) 3390–3401
Theorem 2.1. For system (2.1), we have
(i) The trivial steady state E0 ð0; 0Þ is always unstable. (ii) If a21 r 1 < a11 r2 ðcs þ 1Þ, the steady state E1 is stable; if a21 r1 > a11 r2 ðcs þ 1Þ, E1 is unstable. (iii) If a21 r 1 > a11 r2 ðcs þ 1Þ, the positive steady state E is always stable.
3. Global stability In this section, by using the method of upper–lower solutions and its associated monotone iteration scheme, we study the global stability of the positive steady state E and the boundary steady state E1 of system (2.1), respectively. We first investigate the asymptotic behavior of the solution of system (2.1) in relation to its corresponding constant steady state solution c1 ; c2 ; c3 of system (2.1) given by
r 1 c3 a11 c21 a12 c1 c2 ¼ 0; c2 ðr2 þ a21 c1 a22 c2 Þ ¼ 0; 1 1 c1 c3 cc3 ¼ 0:
s
ð3:1Þ
s
We deal with equations (3.1) as the following elliptic system
Dc1 ¼ r 1 c3 a11 c21 a12 c1 c2 ; Dc2 ¼ c2 ðr2 þ a21 c1 a22 c2 Þ; 1 1 Dc3 ¼ c1 c3 cc3 :
s
ð3:2Þ
s
Let ~c ¼ ð~c1 ; ~c2 ; ~c3 Þ and ^c ¼ ð^c1 ; ^c2 ; ^c3 Þ denote a pair of upper and lower solutions of system (3.2) satisfying ~c P ^c P 0 and
r 1 ~c3 a11 ~c21 a12 ~c1 ^c2 6 0 6 r 1 ^c3 a11 ^c21 a12 ^c1 ~c2 ; ~c2 ðr2 þ a21 ~c1 a22 ~c2 Þ 6 0 6 ^c2 ðr 2 þ a21 ^c1 a22 ^c2 Þ; 1 1 1 1 ~c1 ~c3 c~c3 6 0 6 ^c1 ^c3 c^c3 :
s
s
s
ð3:3Þ
s
It is easy to see that there is a positive constant K 0 ¼ K 0 ða11 ; a12 ; a21 ; a22 ; r 1 ; r 2 ; ~ci ; ^ci Þ such that the following Lipschitz condition holds
jr1 u3 a11 u21 a12 u1 u2 ðr 1 v 3 a11 v 21 a12 v 1 v 2 Þj 6 K 0 ðju1 v 1 j þ ju2 v 2 j þ ju3 v 3 jÞ; ju2 ðr 2 þ a21 u1 a22 u2 Þ v 2 ðr 2 þ a21 v 1 a22 v 2 Þj 6 K 0 ðju1 v 1 j þ ju2 v 2 jÞ; 1 u1 1 u3 cu3 ð1 v 1 1 v 3 cv 3 Þ 6 K 0 ðju1 v 1 j þ ju3 v 3 jÞ; s s s s
ð3:4Þ
for ^ci 6 ui ; v i 6 ~ci . ðmÞ ðmÞ ðmÞ ðmÞ ðmÞ ðmÞ We now construct two sequences cðmÞ ¼ ðc1 ; c2 ; c3 Þ and cðmÞ ¼ ðc1 ; c2 ; c3 Þ from the following iteration process ðm1Þ ðm1Þ ðm1Þ 2 ðm1Þ ðm1Þ cðmÞ ¼ c1 þ K10 ðr 1 c3 a11 ðc1 Þ a12 c1 c2 Þ; 1 ðm1Þ ðm1Þ ðm1Þ cðmÞ ¼ c2 þ K10 ðr2 þ a21 c1 a22 c2 Þ; 2 ðm1Þ ðm1Þ ðm1Þ ðm1Þ cðmÞ ¼ c3 þ K10 1s c1 1s c3 cc3 ; 3 ðmÞ
¼ c1
ðmÞ
¼ c2
c1 c2
ðmÞ
c3
ðm1Þ
ðm1Þ
þ K10 ðr 1 c3
ðm1Þ 2
a11 ðc1
ðm1Þ
ð3:5Þ
ðm1Þ ðm1Þ c2 Þ;
Þ a12 c1
ðm1Þ
ðm1Þ
þ K10 ðr2 þ a21 c1 a22 c2 Þ; ðm1Þ ðm1Þ 1 1 ðm1Þ 1 ðm1Þ ; ¼ c3 þ K 0 s c1 s c3 cc3
with initial data cð0Þ ¼ ~c; cð0Þ ¼ ^c, where K 0 is the Lipschitz constant above. Clearly, the sequences cðmÞ and cðmÞ possess the following property
^c 6 cðmÞ 6 cðmþ1Þ 6 cðmþ1Þ 6 cðmÞ 6 ~c;
m ¼ 1; 2; 3; . . . ;
where the ordering relation is in the componentwise sense. Therefore, the limits lim cðmÞ and lim cðmÞ exist. Denoting m!1
c ¼ lim cðmÞ ; m!1
c ¼ lim cðmÞ ; m!1
then c and c satisfy the following equations
m!1
ð3:6Þ
3394
X. Zhang, R. Xu / Commun Nonlinear Sci Numer Simulat 15 (2010) 3390–3401
r 1 c3 a11 c21 a12 c1 c2 ¼ 0 ¼ r 1 c3 a11 c21 a12 c1 c2 ; c2 ðr 2 þ a21 c1 a22 c2 Þ ¼ 0 ¼ c2 ðr2 þ a21 c1 a22 c2 Þ; 1
s
1 1 1 c1 c3 cc3 ¼ 0 ¼ c1 c3 cc3 :
s
s
ð3:7Þ
s
The constant vectors c and c are called quasi-solutions of system (3.1) in h^c; ~ci. Generally, the constant vectors c and c are not true solutions of system (3.1). However, if c ¼ c, then c (or c) is the unique solution of Eqs. (3.1) in h^c; ~ci. Next, we present a result on the relation between the solution of (2.1) and (3.1). Lemma 3.1. ([8], Theorems 2.1 and 2.2) Let ~c; ^c be a pair of upper and lower solutions of system (3.1) and let c and c be limits in (3.6) which are the quasi-solutions of (3.1) and satisfy (3.7). Then for any initial function satisfying ^ci 6 ui ðx; 0Þ 6 ~ci in ð0; pÞ for i ¼ 1; 2; 3, the corresponding solution u ¼ ðu1 ; u2 ; u3 Þ of system (2.1) possess the property
c 6 uðx; tÞ 6 c as t ! 1ðx 2 ½0; pÞ: Moreover, if c ¼ c then c (or c) is the unique solution of Eqs.(3.1) in h^c; ~ci and
lim uðx; tÞ ¼ c ðx 2 ½0; pÞ:
m!1
We now introduce two Lemmas. Lemma 3.2. (Hirsch[7]) Let F be a C 1 vector field in Rn , whose flow / preserves Rnþ for t P 0 and is strongly monotone in Rnþ . Assume that the origin is an steady state and that all trajectories in Rnþ are bounded. Suppose the matrix-valued map DF : Rn ! Rnn is strictly anti-monotone, in the sense that
if
x > y;
then DFðxÞ < DFðyÞ;
here, ‘‘<” and ‘‘>” represent the standard partial ordering in Rnþ tend to the origin, or else there is a unique steady state p 2 IntRnþ and all trajectories in Rnþ n f0g tend to p. We consider the following differential equations
u_ 1 ¼ r 1 u3 a11 u21 a12 Au1 ; 1 1 u_ 3 ¼ u1 u3 cu3 :
s
ð3:8Þ
s
System (3.8) always has a trivial equilibrium E0 ð0; 0Þ. If r 1 > a12 Að1 þ csÞ, then system (3.8) has a positive equilibrium 12 Að1þcsÞ r 1 a12 Að1þcsÞ ; a ð1þcsÞ2 Þ: E1 ðr1 a a11 ð1þcsÞ 11
By Lemma 3.2, it is easy to prove the following result. Lemma 3.3. If r 1 < a12 Að1 þ csÞ, the equilibrium E0 ð0; 0Þ of system (3.8) is globally asymptotically stable; if r1 > a12 Að1 þ csÞ, 12 Að1þcsÞ r 1 a12 Að1þcsÞ of system (3.8) is globally asymptotically stable. the equilibrium E1 r1 a 2 a11 ð1þcsÞ ; a11 ð1þcsÞ
From Lemmas 3.1–3.3, we have a result on the global stability of the positive steady state E of system (2.1). Theorem 3.1. Let ðu1 ðx; tÞ; u2 ðx; tÞ; u3 ðx; tÞÞ be a solution of system (2.1). If a11 a22 > a12 a21 and ðH1Þ hold, then
limðu1 ðx; tÞ; u2 ðx; tÞ; u3 ðx; tÞÞ ¼ u1 ; u2 ;
t!1
u1 1 þ cs
uniformly for x 2 ½0; p;
i.e., the positive steady state E of system (2.1) is globally asymptotically stable. Proof. Let ðu1 ðx; tÞ; u2 ðx; tÞ; u3 ðx; tÞÞ be a solution of system (2.1). First, we need to seek for a pair of upper and lower solutions ~c and ^c such that ~c > ^c > 0 and condition (3.3) holds. For e > 0 sufficiently small, we define
r1 a21 r1 a11 r 2 ð1 þ csÞ a21 e; þ e; ~c2 ¼ þ 1þ a22 a11 ð1 þ csÞ a11 ð1 þ csÞ a22 r1 ða11 a22 a12 a21 Þ þ a11 a12 r 2 ð1 þ csÞ a12 a12 a21 ¼ þ e; 1þ 2 a11 a11 a22 a11 a22 ð1 þ csÞ ða11 a22 a12 a21 Þðr 1 a21 a11 r 2 ð1 þ csÞÞ a21 a12 a21 a12 a221 ¼ þ þ e; 1 þ a22 a11 a22 a11 a222 a211 a222 ð1 þ csÞ ~c1 ^c1 ¼ ; ^c3 ¼ : 1 þ cs 1 þ cs
~c1 ¼ ^c1 ^c2 ~c3 Choose
ð3:9Þ
e > 0 sufficiently small satisfying
e < min
r 1 ða11 a22 a12 a21 Þ þ a11 a12 r 2 ð1 þ csÞ ða11 a22 a12 a21 Þðr 1 a21 a11 r 2 ð1 þ csÞÞ ; : a11 ð1 þ csÞða11 a22 þ a12 a22 þ a12 a21 Þ a11 ð1 þ csÞða11 a22 þ a12 a21 Þða22 þ a21 Þ
ð3:10Þ
X. Zhang, R. Xu / Commun Nonlinear Sci Numer Simulat 15 (2010) 3390–3401
3395
It follows from (3.9) and (3.10) that the pair ~c ¼ ð~c1 ; ~c2 ; ~c3 Þ and ^c ¼ ð^c1 ; ^c2 ; ^c3 Þ satisfy
r 1 ~c3 a11 ~c21 a12 ~c1 ^c2 < a11 ~c1 e < 0; r 2 þ a21 ~c1 a22 ~c2 ¼ a22 e < 0; 1
s
1 ~c1 ~c3 c~c3 ¼ 0;
s
ð3:11Þ
r 1 ^c3 a11 ^c21 a12 ^c1 ~c2 ¼ a11 ^c1 e > 0; r 2 þ a21 ^c1 a22 ^c2 ¼ a22 e > 0; 1
s
1 ^c1 ^c3 c^c3 ¼ 0:
s
This shows that the condition (3.3) holds. Therefore, ~c ¼ ð~c1 ; ~c2 ; ~c3 Þ and ^c ¼ ð^c1 ; ^c2 ; ^c3 Þ are the coupled upper and lower solutions of (3.2). Second, we claim that there exists a t such that
^ci 6 ui ðx; tÞ 6 ~ci ;
ðx; tÞ 2 ½0; p ½t ; þ1Þ;
i ¼ 1; 2; 3;
ð3:12Þ
where ð^c1 ; ^c2 ; ^c3 Þ and ð~c1 ; ~c2 ; ~c3 Þ are defined in (3.9). It follows from the first and third equations of system (2.1) that
@u1 @ 2 u1 þ r1 u3 a11 u21 ; 6 d1 @t @x2 @u3 @ 2 u3 1 1 þ u1 u3 cu3 : ¼ d1 @t @x2 s s
ð3:13Þ
Let ðu1 ðtÞ; u3 ðtÞÞ be the unique positive solution of the following system
u_ 1 ¼ r 1 u3 a11 u21 ; 1 1 u_ 3 ¼ u1 u3 cu3 :
s
ð3:14Þ
s
By Lemma 3.3, it follows that
lim u1 ðtÞ ¼
t!1
For
r1 ; a11 ð1 þ csÞ
lim u3 ðtÞ ¼
t!1
r1 a11 ð1 þ csÞ2
:
ð3:15Þ
e > 0 sufficiently small satisfying (3.10), by comparison we derive from (3.13)–(3.15) that there is a t1 > 0 such that u1 ðx; tÞ 6 u3 ðx; tÞ 6
r1 þ ¼ ~c1 ; a11 ð1þcsÞ r1 þ 1þecs ¼ a11 ð1þcsÞ2
e
~c3 ;
ðx; tÞ 2 ½0; p ½t 1 ; þ1Þ:
ð3:16Þ
It follows from the second equation of (2.1) that
@u2 @ 2 u2 þ u2 ðr 2 þ a21 ~c1 a22 u2 Þ: 6 d2 @t @x2 By comparison we derive that there is a t2 > t1 such that
u2 ðx; tÞ 6
a21 ~c1 r2 þ e ¼ ~c2 a22
ðx; tÞ 2 ½0; p ½t 2 ; þ1Þ:
ð3:17Þ
It follows from (3.17) and (2.1) that
@u1 @ 2 u1 þ r1 u3 a11 u21 a12 u1 ~c2 ; P d1 @t @x2 @u3 @ 2 u3 1 1 þ u1 u3 cu3 : ¼ d1 @t @x2 s s
ð3:18Þ
By comparison we derive that there is a t3 > t2 such that
u1 ðx; tÞ P
r 1 a12 ~c2 ð1þcsÞ a11 ð1þcsÞ
e ¼ ^c1 ;
u3 ðx; tÞ P
r 1 a12 ~c2 ð1þcsÞ a11 ð1þcsÞ2
1þecs ¼ ^c3 ;
ðx; tÞ 2 ½0; p ½t 3 ; þ1Þ:
It follows from the second equation of (2.1) that
@u2 @ 2 u2 P d2 þ u2 ðr 2 þ a21 ^c1 a22 u2 Þ: @t @x2 By comparison we derive that there is a t > t 3 such that
ð3:19Þ
3396
X. Zhang, R. Xu / Commun Nonlinear Sci Numer Simulat 15 (2010) 3390–3401
u2 ðx; tÞ P
a21 ^c1 r2 e ¼ ^c2 a22
ðx; tÞ 2 ½0; p ½t ; þ1Þ:
ð3:20Þ
Hence, (3.12) holds. u Last, we prove that the nonnegative solutions ðu1 ðx; tÞ; u2 ðx; tÞ; u3 ðx; tÞÞ of (2.1) converges uniformly to ðu1 ; u2 ; 1þ1csÞ. Since the asymptotic behavior of the solution of system (2.1) is the same as if we set the initial time at t ¼ t , without loss of generality, we assume c^i 6 /1 ðx; tÞ 6 ~ci ði ¼ 1; 2; 3Þ. From Lemma 3.1, there exist quasi-solution c ¼ ðc1 ; c2 ; c3 Þ and c ¼ ðc1 ; c2 ; c3 Þ such that the unique solution ðu1 ; u2 ; u3 Þ of (2.1) satisfies
0 < ^ci 6 ci 6 ui ðx; tÞ 6 ci 6 ~ci ; t ! 1ðx 2 ½0; pÞ;
i ¼ 1; 2; 3;
ð3:21Þ
and c and c satisfy (3.7). Therefore (3.7) is reduced to
r1 r1 a11 c1 a12 c2 ¼ 0 ¼ a11 c1 a12 c2 ; 1 þ cs 1 þ cs r 2 þ a21 c1 a22 c2 ¼ 0 ¼ r 2 þ a21 c1 a22 c2 ; 1 1 1 1 c1 c3 cc3 ¼ 0 ¼ c1 c3 cc3 :
s
s
s
ð3:22Þ
s
It follows from the first two equations of (3.22) that
a11 ðc1 c1 Þ a12 ðc2 c2 Þ ¼ 0; a21 ðc1 c1 Þ a22 ðc2 c2 Þ ¼ 0:
ð3:23Þ
Noting that a11 a22 > a12 a21 , we derive from (3.23) that c1 ¼ c1 ; c2 ¼ c2 . It therefore follows from (3.22) that c3 ¼ c3 . Applying Lemma 3.1 we have that
limðu1 ðx; tÞ; u2 ðx; tÞ; u3 ðx; tÞÞ ¼ u1 ; u2 ;
t!1
The proof is completed.
u1 1 þ cs
uniformly for x 2 ½0; p:
h
Theorem 3.2. Let ðu1 ; u2 ; u3 Þ be a solution of system (2.1). If a21 r1 < a11 r2 ðcs þ 1Þ, then
limðu1 ; u2 ; u3 Þ ¼
t!1
r1 r1 ; 0; a11 ð1 þ csÞ a11 ð1 þ csÞ2
! uniformly for
x 2 ½0; p;
i.e., the boundary steady state E1 of system (2.1) is globally asymptotically stable. Proof. Let ðu1 ; u2 ; u3 Þ be a solution of system (2.1). For
e < min
e > 0 sufficiently small satisfying
r1 r2 a11 ð1 þ csÞ a21 r1 ; ; ða12 þ a11 Þð1 þ csÞ a21 a11 ð1 þ csÞ
ð3:24Þ
we define r1 ~c1 ¼ a ð1þ csÞ þ e; 11
~c2 ¼ e;
~c3 ¼ 1þ~c1cs ;
r1 a12 e ^c1 ¼ a ð1þ csÞ a11 e; 11
^c2 ¼ 0;
^c3 ¼ 1þ^c1cs :
ð3:25Þ
It follows from (3.24) and (3.25) that the pair ~c ¼ ð~c1 ; ~c2 ; ~c3 Þ and ^c ¼ ð^c1 ; ^c2 ; ^c3 Þ satisfy
r1 ~c3 a11 ~c21 a12 ~c1 ^c2 ¼ a11 ~c1 e < 0; 21 r1 r 2 þ a21 ~c1 a22 ~c2 ¼ a að1þ csÞ r 2 a22 e þ a21 e < 0; 11
1~ s c1
1s ~c3 c~c3 ¼ 0; r1 ^c3 a11 ^c21 a12 ^c1 ~c2 ¼ a11 ^c1 e > 0; ^c2 ðr 2 þ a21 ^c1 a22 ^c2 Þ ¼ 0; 1^ 1^ ^ s c1 s c3 cc3 ¼ 0:
ð3:26Þ
This shows that the condition (3.3) holds. Therefore ~c ¼ ð~c1 ; ~c2 ; ~c3 Þ and ^c ¼ ð^c1 ; ^c2 ; ^c3 Þ are the coupled upper and lower solutions of (3.2). From Lemma 3.1, there exist quasi-solution c ¼ ðc1 ; c2 ; c3 Þ and c ¼ ðc1 ; c2 ; c3 Þ such that the unique solution ðu1 ; u2 ; u3 Þ of (2.1) satisfies
0 < ^ci 6 ci 6 ui 6 ci 6 ~ci ;
t ! 1;
i ¼ 1; 2; 3;
and c and c satisfy (3.7). Therefore (3.7) is reduced to
ð3:27Þ
3397
X. Zhang, R. Xu / Commun Nonlinear Sci Numer Simulat 15 (2010) 3390–3401
r1 r1 a11 c1 a12 c2 ¼ 0 ¼ a11 c1 a12 c2 ; 1 þ cs 1 þ cs c2 ðr 2 þ a21 c1 a22 c2 Þ ¼ 0 ¼ c2 ðr 2 þ a21 c1 a22 c2 Þ; 1 1 1 1 c1 c3 cc3 ¼ 0 ¼ c1 c3 cc3 :
s
s
s
ð3:28Þ
s
If c2 > 0, it then follows from the first two equations of (3.28) that r1 1þcs
a11 c1 a12 c2 ¼ 0;
ð3:29Þ
r 2 þ a21 c1 a22 c2 ¼ 0; which yields
a21 r 1 a11 r 2 ¼ a12 a21 c2 þ a11 a22 c2 : 1 þ cs Noting that a21 r 1 < a11 r2 ðcs þ 1Þ, a contradiction occurs. Hence, r1 r1 c1 ¼ c1 ¼ a ð1þ csÞ ; c3 ¼ c 3 ¼ a ð1þcsÞ2 . Applying Lemma 3.1 we have that 11 11
limðu1 ; u2 ; u3 Þ ¼
t!1
r1 r1 ; 0; a11 ð1 þ csÞ a11 ð1 þ csÞ2
c2 ¼ c2 ¼ 0.
We
derive
from
(3.28)
that
!
uniformly for x 2 ½0; p:
This completes the proof. h
4. The existence of travelling waves In this section, we discuss the existence of travelling waves of system (1.2) for the case of an infinite one-dimensional domain X ¼ ð1; þ1Þ by using the geometric singular perturbation theory. The infinite domain case is from a modelling point of view since there are no boundaries for individuals to interact with as they drift from their past to their present positions. The model assumes the form
@u1 @ 2 u1 þ r1 ¼ d1 @t @x2
Z
t
1
Z
þ1
Gðx y; t sÞkðt sÞecðtsÞ u1 ðy; sÞdyds
1
a11 u21 a12 u1 u2 ;
ð4:1Þ
2
@u2 @ u2 þ u2 ðr2 þ a21 u1 a22 u2 Þ: ¼ d2 @t @x2 where 2 1 x Gðx; tÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi e 4d1 t ; 4pd1 t
and G satisfies
@G @2G ¼ d1 2 ; @t @x
Gðx; 0Þ ¼ dðxÞ:
System (4.1) still has steady states E0 ð0; 0Þ; E1
r1 ;0 a11 ð1þcsÞ
and E ðu1 ; u2 Þ. Our interest in this section is the possibility of a tran-
sition between the steady states E1 and E in the form of a travelling wave solution. Theorem 4.1. For s sufficiently small, if a21 r 1 > a11 r2 ðcs þ 1Þ, then system (4.1) possesses a travelling wave solution connecting the steady states E1 and E . Proof. Denoting
u3 ¼
Z
t
1
Z
þ1
1
ðxyÞ2 1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e 4d1 ðtsÞ kðt sÞecðtsÞ u1 ðy; sÞdyds; 4pd1 ðt sÞ
system (4.1) can be rewritten as
@u1 @ 2 u1 þ r 1 u3 a11 u21 a12 u1 u2 ; ¼ d1 @t @x2 @u2 @ 2 u2 þ u2 ðr2 þ a21 u1 a22 u2 Þ; ¼ d2 @t @x2 2 @u3 @ u3 1 1 ¼ d1 þ u1 u3 cu3 : @t @x2 s s
ð4:2Þ
3398
X. Zhang, R. Xu / Commun Nonlinear Sci Numer Simulat 15 (2010) 3390–3401
Converting to travelling wave form, by writing
u1 ðx; tÞ ¼ u1 ðzÞ;
u2 ðx; tÞ ¼ u2 ðzÞ;
u3 ðx; tÞ ¼ u3 ðzÞ;
z ¼ x þ ct;
we derive that
€ 1 þ r 1 u3 a11 u21 a12 u1 u2 ; cu_ 1 ¼ d1 u € 2 þ u2 ðr 2 þ a21 u1 a22 u2 Þ; cu_ 2 ¼ d2 u 1 1 € 3 þ u1 u3 cu3 ; cu_ 3 ¼ d1 u
s
ð4:3Þ
s
where prime denotes differentiation with respect to z. We introduce
v 1 ¼ d1 u_ 1 ; v 2 ¼ d2 u_ 2 ; v 3 ¼ d1 u_ 3 : And, we replace
s with e2 s. System (4.3) becomes
1 v 1; d1 c v_ 1 ¼ v 1 ðr1 u3 a11 u21 a12 u1 u2 Þ; d1 1 u_ 2 ¼ v 2 ; d2 c v_ 2 ¼ v 2 u2 ðr2 þ a21 u1 a22 u2 Þ; d2 1 u_ 3 ¼ v 3 ; d1 e2 c 1 1 e2 v_ 3 ¼ v 3 u1 u3 þ e2 cu3 : d1 s s u_ 1 ¼
ð4:4Þ
By introducing new state variables
~ 1 ¼ u1 ; u
v~ 1 ¼ v 1 ;
~ 2 ¼ u2 ; u
v~ 2 ¼ v 2 ;
~ 3 ¼ u3 ; u
v~ 3 ¼ ev 3
and dropping the tildes, we have
1 v 1; d1 c v_ 1 ¼ v 1 ðr1 u3 a11 u21 a12 u1 u2 Þ; d1 1 u_ 2 ¼ v 2 ; d2 c v_ 2 ¼ v 2 u2 ðr2 þ a21 u1 a22 u2 Þ; d2 1 eu_ 3 ¼ v 3 ; d1 ec 1 1 ev_ 3 ¼ v 3 u1 u3 þ e2 cu3 : d1 s s u_ 1 ¼
ð4:5Þ
This system is called slow system. Introducing a new independent variable g defined by z ¼ eg, system (4.5) transforms into
1 v1; d1 c v_ 1 ¼ e v 1 ðr1 u3 a11 u21 a12 u1 u2 Þ ; d1 1 u_ 2 ¼ e v 2 ; d 2 c v_ 2 ¼ e v 2 u2 ðr2 þ a21 u1 a22 u2 Þ ; d2 1 u_ 3 ¼ v 3 ; d1 ec 1 1 v_ 3 ¼ v 3 u1 u3 þ e2 cu3 : d1 s s
u_ 1 ¼ e
ð4:6Þ
X. Zhang, R. Xu / Commun Nonlinear Sci Numer Simulat 15 (2010) 3390–3401
3399
This system is called fast system. In the slow system (4.5), if e ¼ 0, then the flow of that system is confined to the set
M0 ¼ fðu1 ; v 1 ; u2 ; v 2 ; u3 ; v 3 Þ 2 R6 : v 3 ¼ 0; u3 ¼ u1 g which is a four-dimensional manifold for system (4.5). By Theorem 9.1 [3], we know that if M 0 is normally hyperbolic, then for sufficiently small e > 0 we can obtain an invariant manifold of system (4.5). To verify normal hyperbolicity, we need to verify that the linearization of the fast system (4.6), restricted to M 0 , has precisely dim M0 eigenvalues on the imaginary axis with the remainder of the spectrum being hyperbolic. The linearization of the fast system (4.6) restricted to M0 is given by the matrix
0 B B B B B B B B B @
0 0 0 0 0 1s
1 0 0C C C 0 0 0 0 0C C C 0 0 0 0 0C C 1 C 0 0 0 0 d1 A 0 0 0 1s 0 0 0 0 0 0 0 0 0
pffiffiffiffiffiffiffiffi which has six eigenvalues f0; 0; 0; 0; 1= d1 sg. Thus, the manifold M 0 is normally hyperbolic. By the geometric singular perturbation theory, we know that the slow system (4.5) has a invariant manifold Me , which can be written as
Me ¼ fðu1 ; v 1 ; u2 ; v 2 ; u3 ; v 3 Þ 2 R6 : v 3 ¼ hðu1 ; v 1 ; u2 ; v 2 ; eÞ; u3 ¼ u1 þ gðu1 ; v 1 ; u2 ; v 2 ; eÞg where the functions h and g satisfy
hðu1 ; v 1 ; u2 ; v 2 ; 0Þ ¼ gðu1 ; v 1 ; u2 ; v 2 ; 0Þ ¼ 0: Thus the functions h and g can be expanded into the form of a Taylor series about e:
hðu1 ; v 1 ; u2 ; v 2 ; eÞ ¼ eh1 þ e2 h2 þ ;
ð4:7Þ
gðu1 ; v 1 ; u2 ; v 2 ; eÞ ¼ eg 1 þ e2 g 2 þ :
We determine the functions h and g. Note that M e is the invariant manifold for the flow of system (4.5). From w2 ¼ h; w1 ¼ u1 þ g and (4.5), we have
@g @g @g @g v_ 1 þ v_ 2 ; u_ 1 þ u_ 2 þ @u2 @v 1 @v 2 d1 @u1 e @h @h @h @h c v_ 1 þ v_ 2 h ecu1 g¼ u_ 1 þ u_ 2 þ 2 @u2 @v 1 @v 2 d1 1=s þ e c @u1 c 3 @h @h @h @h c ¼ ðe 6 4 e þ Þ v_ 1 þ v_ 2 h ecu1 : u_ 1 þ u_ 2 þ @u1 @u2 @v 1 @v 2 d1 s h ¼ ed1
v1
þ
Substituting (4.7) into (4.8) and comparing coefficients of
h1 ¼ v 1 ;
h2 ¼ 0;
g 1 ¼ 0;
ð4:8Þ
e and e2 , we get
g 2 ¼ r 1 u1 þ a11 u21 þ a12 u1 u2 cu1 :
ð4:9Þ
Then the slow system (2.4) restricted to Me is
1 v 1; d1 c v_ 1 ¼ v 1 u1 Mðu1 ; u2 Þ; d1 1 u_ 2 ¼ v 2 ; d2 c v_ 2 ¼ v 2 u2 Nðu1 ; u2 Þ; d2 u_ 1 ¼
ð4:10Þ
where
Mðu1 ; u2 Þ ¼ r 1 ð1 r 1 Þ a11 ð1 r 1 Þu1 a12 ð1 r 1 Þu2 r 1 c; Nðu1 ; u2 Þ ¼ r 2 þ a21 u1 a22 u2 : System (4.10) still possesses steady state points E1
r1 ; 0; 0; 0 a11 ð1þcsÞ
and E ðu1 ; 0; u2 ; 0Þ. Therefore, the results in Gardner [4]
are applicable, yielding a heteroclinic connection between the steady states E1 and E of (4.10). Thus, system (4.2) possesses travelling wave solution connecting the steady states E1 and E . This completes the proof. h
3400
X. Zhang, R. Xu / Commun Nonlinear Sci Numer Simulat 15 (2010) 3390–3401
5. Numerical simulations In this section, we give numerical simulations to illustrate the results of Section 4. In system (4.2), we let d1 ¼ 1; d 2 ¼ 2; a11 ¼ a12 ¼ a 22 ¼ 1; a21 ¼ 2; r 1 ¼ r 2 ¼ 1; c ¼ 1; s ¼ 0:1: Then system (4.2) admits 1 1 0:7 0:3 0:7 ; 0; 1:21 ; 1:1 ; 1:21 . By Theorem 4.1, we know that system (4.2) possesses travelling and E 1:1 steady states E1 ð0; 0; 0Þ; E2 1:1 wave solution connecting the steady states E1 and E (see Fig. 1). Next, we consider the effect of delay s on the wave speed. In system (4.2), we let d1 ¼ 1; d2 ¼ 2; a11 ¼ a12 ¼ a22 ¼ 1; a21 ¼ 2; r1 ¼ r 2 ¼ 1; c ¼ 1; s ¼ 0:5: Then the steady states of system (4.2) are
1 1 E1 ð0; 0; 0Þ; E2 23 ; 0; 2:25 and E 59 ; 19 ; 2:7 . Fig. 2 shows the temporal solution of system (4.2). In system (4.2), we let d1 ¼ 1; d2 ¼ 2; a11 ¼ a12 ¼ a22 ¼ 1; a21 ¼ 2; r1 ¼ r2 ¼ 1; c ¼ 1; s ¼ 0:8: Then the steady states of sys 5
1 ; 1 ; 0:7 . Fig. 3 shows the temporal solution of system (4.2). tem (4.2) are E1 ð0; 0; 0Þ; E2 9 ; 0; 3:24 and E 1:4 2:7 27 2:43 From Figs. 1–3, we can get the section drawing of the solutions u1 and u2 of system (4.2) (see Fig. 4). From Fig. 4, we see that as the delay s increases, the wave speed becomes more and more slower.
Fig. 1. The temporal solution found by numerical integration with
s ¼ 0:1.
Fig. 2. The temporal solution found by numerical integration with
s ¼ 0:5.
Fig. 3. The temporal solution found by numerical integration with
s ¼ 0:8.
3401
X. Zhang, R. Xu / Commun Nonlinear Sci Numer Simulat 15 (2010) 3390–3401 0.35
0.95 tau=0.1 tau=0.5 tau=0.8
0.9
0.25
0.8
solutions − u2
solutions − u1
0.85
tau=0.1 tau=0.5 tau=0.8
0.3
0.75 0.7 0.65
0.2 0.15 0.1
0.6 0.05
0.55 0.5
0
20
40
t − axis
60
80
100
0 0
20
40
t − axis
60
80
100
Fig. 4. The temporal solution found by numerical integration.
Acknowledgements This work was supported by the National Natural Science Foundation of China (No. 10671209). References [1] Al-Omari JFM, Gourley SA. Dynamics of a stage-structured population model incorporating a state-dependent maturation delay. Nonlinear Anal RWA 2005;6:13–33. [2] Britton NF. Spatial structures and periodic travelling waves in an integrodifferential reaction-diffusion population model. SIAM J Appl Math 1990;50:1663–88. [3] Fenichel N. Geometric singular perturbation theory for ordinary differential equations. J Differ Equat 1979;31:53–98. [4] Gardner R. Existence and stability of travelling wave solutions of competition models: a degree theoretic approch. J Differ Equat 1982;44:343–64. [5] Gourley SA, Britton NF. A predator-prey reaction-diffusion system with nonlocal effects. J Math Biol 1996;34:297–333. [6] Gourley SA, Ruan S. Convergence and travelling fronts in functional equations with nonlocal terms: a competition model. SIAM J Math Anal 2003;35:806–22. [7] Hirsch MW. The dynamical systems approach to differential equations. Bull Amer Math Soc 1984;11:1–64. [8] Pao CV. Convergence of solutions of reaction-diffusion systems with time delays. Nonlinear Anal TMA 2002;48:349–62. [9] Thieme HR, Zhao XQ. A non-local delayed and diffusive predator-prey model. Nonlinear Anal RWA 2001;2:145–60. [10] Wang ZC, Li WT, Ruan S. Travelling wave fronts of reaction-diffusion systems with spatio-temporal delays. J Differ Equat 2006;222:185–232. [11] Wang ZC, Li WT. Monotone travelling fronts of a food-limited population model with nonlocal delay. Nonlinear Anal RWA 2007;8:699–712. [12] Wu J. Theory and applications of partial functional differential equations. New York: Springer Verlag; 1996. [13] Xu R. A reaction-diffusion predator-prey model with stage structure and nonlocal delay. Appl Math Comput 2006;175:984–1006. [14] Xu R, Ma Z. Global stability of a reaction-diffusion predator-prey model with a nonlocal delay. Math Comput Model 2009;50:94–206. [15] Zhang GB, Li WT, Lin G. Traveling waves in delayed predator-prey systems with nonlocal diffusion and stage structure. Math Comput Model 2009;49:1021–9. [16] Zhang JM, Peng YH. Travelling waves of the diffusive Nicholson’s blowflioes equation with strong genetric delay kernel and non-local effect. Nonlinear Anal 2008;68:1263–70.