Nonlinear Analysis: Real World Applications 9 (2008) 2196 – 2213 www.elsevier.com/locate/na
Traveling wave solutions for a competitive reaction–diffusion system and their asymptotics Xiaojie Hou, Anthony W. Leung∗ Department of Mathematical Sciences, University of Cincinnati, Cincinnati, OH, USA Received 9 April 2007; accepted 20 July 2007
Abstract For a competitive PDE system, we study the existence, uniqueness, and asymptotic behaviors of the traveling wave solutions connecting a monoculture equilibrium to a co-existence equilibrium. We use the method of upper–lower solutions to prove the existence of traveling wave solutions, and investigate the asymptotic behavior of the traveling waves in relation to various interacting parameters of the system. By comparing with the upper solution, we obtain asymptotic description of the solution for large x or t in relation to the interacting parameters, and show the uniqueness of traveling wave solutions connecting the two equilibria with such asymptotic rates. Numerical results are also presented to illustrate the theoretical results. 䉷 2007 Elsevier Ltd. All rights reserved. MSC: primary 35B35 secondary 35K57; 35B40; 35P15 Keywords: Traveling wave; Existence; Asymptotics; Uniqueness; Numerical solutions
1. Introduction This article is concerned with the traveling wave solutions of system: ut = duxx + u(a1 − b1 u − c1 v), (x, t) ∈ R × R + , vt = dv xx + v(a2 − b2 u − c2 v),
(1.1)
where u = u(x, t), v = v(x, t); and d, ai , bi , ci are positive constants. The system describes two interacting species with diffusive effects. The quantities u and v are the concentrations of two competing populations. The parameters ai , i = 1, 2 are the intrinsic growth rates of the species; b1 , c2 are the crowding-effect coefficients. The parameters c1 and b2 are the coefficients describing competitions between the species; and d is the diffusion constant. It is a general belief among biologists that if the two competing species have different preferences in resource usage, then the intra-specific competition between the two species is stronger than the inter-specific competition, therefore, the co-existence of the two species is possible. ∗ Corresponding author. Tel./fax: +1 513 556 3417.
E-mail address:
[email protected] (A.W. Leung). 1468-1218/$ - see front matter 䉷 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.nonrwa.2007.07.007
X. Hou, A.W. Leung / Nonlinear Analysis: Real World Applications 9 (2008) 2196 – 2213
2197
We will show under some mild conditions, that the evolution from a monoculture state to the co-existence of two competing species can be proved rigorously. Moreover, this evolution is in the form of traveling wave solutions, it has certain speeds and a particular shape. We also investigate the uniqueness and the traveling wave the asymptotic rate of solutions a1 a1 for each given speed. The traveling wave solutions have the form u , c > 0 and d x + ca 1 t , v d x + ca 1 t a1 2 c1 a2 b1 −a1 b2 join the equilibria 0, ac22 and ab11 cc22 −a −b2 c1 , b1 c2 −b2 c1 as d x + ca 1 t moves from −∞ to +∞. This means that the first species moves from extinction to the positive co-existence state, while the second species moves from carrying capacity state to co-existence state. Throughout this paper, we will always assume the following hypotheses: a1 a2 [H1] < b 1 b2
and [H2]
a2 a1 < . c2 c1
In [21], the existence of traveling wave solutions moving from (0, 0) (rather than 0, ac22 ) to the positive coexistence state was proved. On the other hand, if both the inequalities in [H1] and [H2] are reversed, [4,7] proved the existence of traveling wave solutions moving from one (monoculture) equilibrium on one positive axis to another (monoculture) equilibrium on the other positive axis. In each of [21,4,7], dynamical system and ordinary differential equation methods are used to prove the existence of traveling wave solutions. More recently, [15] considers a special example of (1.1) satisfying [H1],[H2], and the existence of traveling wave connecting a monoculture state to co-existence state is proved by taking limit for recursion monotone maps as in [16]. The results and methods in this article are different. We use an alternative method of upper–lower solutions to prove the existence of traveling wave solutions. More importantly, we further investigate the asymptotic behavior of the traveling waves in relation to various interacting parameters of the system. We also consider the uniqueness of such traveling wave solutions. Related papers involving the study of traveling wave solutions for other interacting population systems can be found in e.g. [1,6,8,12,16,22]. Other boundary value problems involving system (1.1), not related traveling wave solutions, can be found in many references. Examples of related books are [9,10,14,17,19,20]. Other related papers concerning traveling waves can be found in [3,11,13,25]. In Section 2, the existence of traveling wave solution is proved by first reversing order in the second equation so that the resulting system becomes monotone. Then an appropriate upper solution is constructed, and the theory in [23] is used to obtain existence of traveling solutions with the desired limits at infinities. In Section 3, we make a more detailed examination of the asymptotic behavior of the traveling wave, by using comparison and upper–lower solution method in [24]. The construction of upper solution for the system is made by using traveling wave solutions of the well known KPP (Kolmogorov, Petrovsky and Piskunov) equation. The use of upper–lower solution method avoids the difficulty of analyzing the stable and unstable manifolds used in dynamical system method in [21,6]. The method also provides asymptotic description of the solution for large x or t in relation to the interacting parameters of the components, which is not considered in [15]. For the analysis of the asymptotics, appropriate linearized systems at infinities are used for approximation. The asymptotic analysis of the traveling wave solutions not only provides the necessary information to the uniqueness but also sheds light on the stability of the traveling wave solutions. However, due to a different nature of the stability analysis, we will address this issue in a separate paper. In Section 3, the uniqueness of traveling wave solutions with the same asymptotic properties is also proved by sliding method. The application of sliding method in the context of a system of PDE is new. Lastly, in Section 4, we numerically compute the traveling wave solutions which have the properties as described in Section 2. 2. Existence of traveling wave We first consider the ordinary differential equations corresponding to (1.1), i.e. d =0. Under conditions [H1] system of a2 a1 c2 −a2 c1 a2 b1 −a1 b2 and [H2], the points 0, c2 and b1 c2 −b2 c1 , b1 c2 −b2 c1 are equilibria, with the first equilibrium being unstable and the second stable. (SeeFig. 1.) We will prove in this section the of system (1.1), with existence of traveling wave solution d > 0, of the form u
a1 d x
+ ca 1 t , v
a1 d x
+ ca 1 t
joining the equilibria 0, ac22
and
a1 c2 −a2 c1 a2 b1 −a1 b2 b1 c2 −b2 c1 , b1 c2 −b2 c1
2198
X. Hou, A.W. Leung / Nonlinear Analysis: Real World Applications 9 (2008) 2196 – 2213
v a1 c1 a2 c2 (
a1c2 − a2c1 a2b1 − a1b2 , ) b1c2 − b2c1 b1c2 − b2c1
a1 b1
O
a2 b2
u
Fig. 1. Nullclines for ODE system corresponding to (1.1); the symbols · and ◦ represent stable and unstable equilibrium, respectively. √
a1 d x
as
+ ca 1 t moves from −∞ to +∞. Letting: = a1 t and x = a1−1 d x, ˜
Eq. (1.1) can be written as u = ux˜ x˜ + u(1 − a1−1 b1 u − a1−1 c1 v), v = vx˜ x˜ + v(a1−1 a2 − a1−1 b2 u − a1−1 c2 v),
(2.1)
(x, ˜ ) ∈ R × R.
(2.2)
Let u = kw,
v = qz,
(2.3)
where k = a1 b1−1 and q is a constant satisfying a2 c2−1 < q < a1 c1−1 .
(2.4)
System (2.2) becomes w = wx˜ x˜ + w(1 − w − rz), z = zx˜ x˜ + z(1 − bw − 1 (1 + 2 )z),
(2.5)
where r = a1−1 c1 q, b = b2 b1−1 ,
1 = a1−1 a2 , 2 = a2−1 c2 q − 1.
(2.6)
Note that from [H1] and (2.4), we have 0 < b < 1 ,
0 < r < 1,
2 > 0.
(2.7)
Observe that by choosing q close to a2 c2−1 in (2.4), 2 can be made arbitrarily small. Lemma 2.1. Consider system (1.1) under hypotheses [H1] and [H2]. The change of variables (2.1), (2.3) with q satisfying (2.4) transforms (1.1) into system (2.5). The parameters in (2.5) are related to those in (1.1) by (2.6). The parameters r, 1 , 2 and b satisfy the inequalities in (2.7).
X. Hou, A.W. Leung / Nonlinear Analysis: Real World Applications 9 (2008) 2196 – 2213
2199
One can readily verify that if (w(x, ˜ ), z(x, ˜ )) is a solution of (2.5), then (u(x, t), v(x, t)) = u a1−1 d x, ˜ a1−1 , v a1−1 d x, ˜ a1−1 = (kw(x, ˜ ), qz(x,)) ˜
(2.8)
is a solution of (1.1), where k, q are described in (2.3), (2.4). We will look for solution of (2.5) with the form (w(x, ˜ ), z(x, ˜ )) = (w(), z()), = x˜ + c, where c > 0 satisfying ⎧ 1 ⎪ , ⎪ lim (w(), z()) = 0, ⎨ →−∞ 1 + 2 (2.9) ⎪ 1 − b 1 (1 + 2 − r) ⎪ ⎩ lim→∞ (w(), z()) = , . 1 (1 + 2 ) − br 1 (1 + 2 ) − br Relating back to (2.5), we are thus looking for solutions of −w + cw = w(1 − w − rz), −z + cz = z(1 − bw − 1 (1 + 2 )z),
− ∞<<∞
(2.10)
with boundary conditions (2.9). Theorem 2.2. Under hypotheses [H1], [H2] and [H3] a2 < a1 , system (1.1) has a traveling wave solution of the form:
(u(x, t), v(x, t)) = kw
a1 a1 x + ca 1 t , qz x + ca 1 t d d
(2.11)
2 c1 for any c > 2 ba11 ab11 cc22 −a −b2 c1 . Here (w, z) is a function of one single variable satisfying (2.10) for −∞ < < ∞ and (2.9) as → ±∞. Moreover, w() and z() are positive monotonic functions for −∞ < < ∞. Remark 2.1. The function (w, z) in (2.11) now satisfies (2.10) with parameters satisfying (2.6), (2.7) as described above and 1 < 1 due to [H3] above. Moreover, we have lim (u(x, t), v(x, t)) = (0, a2 /c2 ),
t→−∞
lim (u(x, t), v(x, t)) =
t→+∞
a1 c2 − a2 c1 a2 b1 − a1 b2 , b1 c2 − b2 c1 b1 c2 − b2 c1
.
(2.12)
Proof. Note that by (2.6), we always have relation (1 + 2 − r)1 b 1 a1 c2 − a 2 c1 > 0. = 1 (1 + 2 ) − rb a1 b1 c2 − b2 c1
(2.13)
The last inequality above is due to [H1] and [H2]. After transforming (1.1) into (2.5) and (2.10) as described above, we next introduce the transformation u1 () = w(), u2 () =
1 − z() 1 + 2
(2.14)
2200
X. Hou, A.W. Leung / Nonlinear Analysis: Real World Applications 9 (2008) 2196 – 2213
u2 A
( 1 ) 1 + 2
B
II
O
u1 I
b(1+2 −r) (1+ −r)1 . I and II represent the nullclines 1− 1+r −u1 +ru2 =0 , (1+ )( (1+ Fig. 2. Nullclines for system (2.15); A= 1, 1+1 , B = (1+2 )−rb 2 )−rb) 2 2 1 2 2 1 and bu1 − 1 (1 + 2 )u2 = 0, respectively.
for −∞ < < ∞. This transforms (2.10) into 1 + 2 − r −(u1 ) + c(u1 ) = u1 − u1 + ru2 , 1 + 2 1 −(u2 ) + c(u2 ) = − u2 (bu1 − 1 (1 + 2 )u2 ) 1 + 2
(2.15)
for −∞ < < ∞, which is monotone for 0 u1 , 0 u2 1+12 . We illustrate the transformed system (2.15) by Fig. 2.
2 −r)1 , Note that the two expressions on the right-hand side of (2.15) now vanish simultaneously at (u1 , u2 ) = ( (1+ 1 (1+2 )−rb
b(1+2 −r) (1+2 )(1 (1+2 )−rb) )
in the first quadrant. The second component satisfies
b(1 + 2 − r) 1 < (1 + 2 )(1 (1 + 2 ) − rb) 1 + 2
(2.16)
because 0 < b < 1 . We now construct a pair of coupled solutions for system (2.15). Let u (), −∞ < < ∞, be the monotone solution of the following KPP equation: (1 + 2 − r)1 u − cu + u (2.17) − u = 0 1 (1 + 2 ) − rb 2 −r)1 with lim→−∞ u () = 0, lim→∞ u () = (1+ . From [14] and (2.13), such solution u () exists for c > 2 1 (1+2 )−rb b1 a1 c2 −a2 c1 a1 b1 c2 −b2 c1 . Define
u¯ 1 () = u (),
u¯ 2 () =
b u (). 1 (1 + 2 )
Observe that the system (2.15) is monotone in the region: 0 < u1 <
(1 + 2 − r)1 := u∗1 = lim u¯ 1 () 1 (1 + 2 ) − rb →∞
0 < u2 <
b(1 + 2 − r) := u∗2 = lim u¯ 2 (). (1 + 2 )(1 (1 + 2 ) − rb) →∞
(2.18)
X. Hou, A.W. Leung / Nonlinear Analysis: Real World Applications 9 (2008) 2196 – 2213
2201
For 0 u2 u¯ 2 (), one readily verifies that 1 + 2 − r − u¯ 1 + ru2 − (u¯ 1 ) + c(u¯ 1 ) − u¯ 1 1 + 2 (1 + 2 − r)1 1 + 2 − r u¯ 1 − u¯ 1 − + u¯ 1 − r u¯ 2 1 + 2 1 (1 + 2 ) − rb rb (1 + 2 − r)(1 + 2 )1 − 1 (1 + 2 − r)(1 + 2 ) + rb(1 + 2 − r) u¯ 1 − u∗1 (1 (1 + 2 ) − rb)(1 + 2 ) 1 (1 + 2 ) = u¯ 1
(1 + 2 − r)(1 + 2 )(1 − 1 ) = 0. (1 (1 + 2 ) − rb)(1 + 2 )
For 0 u1 u¯ 1 (), one verifies 1 − u¯ 2 (bu1 − 1 (1 + 2 )u¯ 2 ) − (u¯ 2 ) + c(u¯ 2 ) − 1 + 2 b (1 + 2 − r)1 1 b − u − u − u (bu − 1 (1 + 2 )u¯ 2 ) 1 (1 + 2 ) 1 (1 + 2 ) − rb 1 + 2 1 (1 + 2 ) b (1 + 2 − r)1 = u − u > 0. 1 (1 + 2 ) 1 (1 + 2 ) − rb
(2.19)
(2.20)
The pair of functions defined by ¯ 1 () = u¯ 1 (−),
¯ 2 () = u¯ 2 (−)
form a pair of upper solutions for the monotone system 1 + 2 − r − 1 + r2 = 0, (1 ) + c(1 ) + 1 1 + 2 1 − 2 (b1 − 1 (1 + 2 )2 ) = 0 (2 ) + c(2 ) + 1 + 2 for −∞ < < ∞, in the sense that 1 + 2 − r − ¯ 1 + r2 0, (¯ 1 ) + c(¯ 1 ) + ¯ 1 1 + 2 1 − ¯ 2 (b1 − 1 (1 + 2 )¯ 2 ) 0 (¯ 2 ) + c(¯ 2 ) + 1 + 2
(2.21)
(2.22)
(2.23)
for −∞ < < ∞, all 0 2 ¯ 2 (), 0 1 ¯ 1 (). (Note that the system (2.22) is monotone in the region: 0 1 u∗1 , 0 2 u∗2 , where u∗1 , u∗2 are defined in (2.18)). In particular, (2.23) is true when 1 = ¯ 1 () in the first equation, 2 = ¯ 2 () in the second equation for all −∞ < < ∞. Let 1 + 2 − r − 1 + r2 , F1 (1 , 2 ) = 1 1 + 2 1 − 2 (b1 − 1 (1 + 2 )2 ). F2 (1 , 2 ) = 1 + 2 We clearly have for i = 1, 2 that sb >0 Fi s, 21 (1 + 2 )
2202
X. Hou, A.W. Leung / Nonlinear Analysis: Real World Applications 9 (2008) 2196 – 2213
− → − → − → for s > 0 sufficiently small. We now apply Theorem 4.2 in [23, p. 176], with F = (F1 , F2 ), w+ = (0, 0), w− = − → (1+2 −r)1 b(1+2 −r) 2 1 (1+2 )−rb , (1+2 )(1 (1+2 )−rb) and K the class of vector-valued functions () = (1 (), 2 ()) ∈ C (−∞, ∞) − → − → decreasing monotonically and satisfying lim () = − w→ ¯ (), ¯ ()) ± . The existence of the function () = ( →±∞
1
2
satisfying (2.23) implies that
− → k + Fk ( ()) ∗ w := inf sup − → −k (} ∈K ,k is finite, and we have c w ∗ . We can thus apply Theorem 4.2 in [23] described above to assert that the system (2.22) has a solution (ˆ 1 (), ˆ 2 ()) with each component monotonically decreasing for −∞ < < ∞ and lim→±∞ (ˆ 1 (), ˆ 2 ())= − w→ ± . The function (uˆ 1 (), uˆ 2 ()) := (ˆ 1 (−), ˆ 2 (−))
(2.24)
is then a solution of the system (2.15). Finally, setting w() = uˆ 1 (), z() = 1+12 − uˆ 2 () for −∞ < < ∞ as in (2.14), then (u(x, t), v(x, t)) as defined in (2.11) is a traveling wave solution of system (1.1) for −∞ < x < ∞, −∞ < t < ∞, satisfying (2.12) as described in Remark 2.1 following the statement of Theorem 2.2. Remark 2.2.√ The wave speed of the solution (u(x, t), v(x, t)) for system (1.1) described in (2.11) in the original space variable is c a1 d. Interchanging the roles of a1 , b1 and c1 , respectively, with a2 , c2 and b2 , we obtain readily the following corollary from Theorem 2.2. Corollary 2.3. Assume hypotheses [H1], [H2] and [H3∗ ] a1 < a2 . 1 b2 Then for any c > 2 ac22 ac22 bb11 −a −c1 b2 , system (1.1) has a traveling wave solution of the form: a2 a2 (u(x, t), v(x, t)) = u˜ x + ca 2 t , v˜ x + ca 2 t . d d
(2.25)
Here u, ˜ v˜ ∈ C 2 (−∞, ∞) are monotonic functions of one single variable and lim (u(x, t), v(x, t)) = (a1 /b1 , 0),
t→−∞
lim (u(x, t), v(x, t)) =
t→+∞
a1 c2 − a2 c1 a2 b1 − a1 b2 , b 1 c2 − b 2 c 1 b 1 c2 − b 2 c1
.
(2.26)
Remark 2.3.√ The wave speed of the solution (u(x, t), v(x, t)) for system (1.1) described in (2.25) in the original space variable is c a2 d. 3. Asymptotic rates and uniqueness In order to obtain better description of the asymptotic behavior of the traveling wave solution of system (1.1) or (2.15) at ±∞, we will next establish some comparisons of the solution with appropriate upper and lower solutions. These comparisons will also be used for analyzing the uniqueness and stability of the traveling waves. Theorem 3.1. Under hypotheses [H1]–[H3], system (1.1) has a traveling wave solution of the form (2.11) for any 2 c1 c > 2 ba11 ab11 cc22 −a −b2 c1 . The function (w, z) in (2.11) is a function of one variable satisfying (2.10) for −∞ < < ∞
X. Hou, A.W. Leung / Nonlinear Analysis: Real World Applications 9 (2008) 2196 – 2213
2203
and (2.9) as → ±∞. Furthermore, w() and z() are monotonic functions for −∞ < < ∞, satisfying 0 < w() ≤ u(), ¯ 0<
1 b u(). ¯ − z() 1 + 2 1 (1 + 2 )
(3.1)
2 −r)1 ¯ = (1+ . ¯ = 0, lim→∞ u() Here, u() ¯ is the solution of Eq. (2.17) with lim→−∞ u() 1 (1+2 )−rb
Proof. The traveling wave solution constructed in Theorem 2.2 does not necessarily satisfy the right-hand side of inequalities (3.1). In order to obtain a solution with these additional properties, we first follow the proof of Theorem 2.2 to the end of Eq. (2.20), obtaining upper solution u¯ 1 (), u¯ 2 () for (2.15). Then we construct a pair of coupled lower solutions for (2.15) as follows: if 0, u1 () = (3.2) e if < 0 and u2 () =
u (). 1 + 2 1
(3.3)
Here ∈ (0, 1 − 1+r 2 ) and ∈ (c, ∞) are chosen as sufficiently small and large constants, respectively, so that we have the property u1 () u¯ 1 () = u(), ¯ u2 () u¯ 2 () =
b u() ¯ 1 (1 + 2 )
(3.4)
for −∞ < < ∞. The parameter in (3.3) will be reduced to a sufficiently small positive constant, 0 < < 1, to be determined below. For 0 u2 u¯ 2 (), > 0, we verify r r −u1 + cu1 − u1 1 − (3.5) − u1 + ru2 = − 1 − − + ru2 < 0. 1 + 2 1 + 2 For 0 u2 u¯ 2 (), < 0, we have r − u1 + cu1 − u1 1 − − u1 + ru2 1 + 2 r 2 = e − + c − 1 + + e − ru2 1 + 2 r 2 e − + c − 1 + + 1 + 2 e (− + c) < 0.
(3.6)
For u1 ()u1 u¯ 1 (), > 0, we verify 1 − u2 + cu2 − − u2 (bu1 − 1 (1 + 2 )u2 ) 1 + 2 1 − (b − 1 ) − 1 + 2 1 + 2 =−
1 (1 − )(b − 1 ) < 0. 1 + 2
The last inequality above is valid by choosing 0 < < b1 .
(3.7)
2204
X. Hou, A.W. Leung / Nonlinear Analysis: Real World Applications 9 (2008) 2196 – 2213
For u1 ()u1 u¯ 1 (), < 0, we have 1 − u2 + cu2 − − u2 (bu1 − 1 (1 + 2 )u2 ) 1 + 2 1 e (−2 + c) − (1 − e )(bu1 − 1 e ) 1 + 2 1 + 2 1 2 − + c − (1 − e )(b − 1 ) = e 1 + 2
e (− + c) < 0. 1 + 2
(3.8)
From (2.19), (2.20), (3.5)–(3.8), we see that the pairs (u¯ 1 (), u¯ 2 ()) and (u1 (), u2 ()) form coupled upper–lower solutions for the monotone system (2.15). We can apply Theorem 3.6 in [24], with zero delay, to conclude that system (2.15) has a traveling wave solution (u01 (), u02 ()) with lim→−∞ (u01 (), u02 ()) = (0, 0), lim→∞ (u01 (), u02 ()) = b(1+2 −r) (1+2 −)1 1 (1+2 )−rb , (1+2 )(1 (1+2 )−rb)
and
u1 ()u01 () u¯ 1 () = u(), ¯ u2 ()u02 () u¯ 2 () =
b u() ¯ 1 (1 + 2 )
(3.9)
for −∞ < < ∞. (Note that the differential inequalities for the lower solutions do not have to hold at = 0.) Moreover each u0i (), i = 1, 2 are monotone functions. The rectangle considered is described in (2.18). We thus obtain (3.1) by letting w() = u01 (), 1+12 − z() = u02 () for −∞ < < ∞. To obtain more asymptotic information of the traveling wave solution at ±∞, we let U () := (u01 (), u02 ())T
for − ∞ < < ∞,
(3.10)
which is the solution of (2.15) as described in the proof of Theorem 3.1. We differentiate (2.15) with respect to , and note that (U ()) := (1 , 2 )T satisfies r 0 0 (1 ) − c(1 ) + 1 − − 2u1 + ru2 1 + ru01 2 = 0, (3.11) 1 + 2 b (2 ) − c(2 ) + − bu02 1 + (21 (1 + 2 )u02 − bu01 − 1 )2 = 0. (3.12) 1 + 2 The limiting equation for (3.11) and (3.12) as → −∞ are, respectively, r ( 1 ) − c( 1 ) + 1 − 1 = 0, 1 + 2 ( 2 ) − c( 2 ) +
b − 1 2 = 0. 1 + 2 1
(3.13) (3.14)
Eq. (3.13) has two independent solutions of the form: √ √ 2 2 11 = e((c− c −4(1−r/(1+2 )))/2) , 21 = e((c+ c −4(1−r/(1+2 )))/2) . Relating (3.11) with (3.13), we deduce 1 has the following property as → −∞: √ √ 2 2 1 = [1 + o(1)]e((c− c −4(1−r/(1+2 )))/2) + [1 + o(1)]e((c+ c −4(1−r/(1+2 )))/2) for some constants and (see e.g. [5, Theorem 2, Chapter IV]). We next show that = 0.
(3.15)
X. Hou, A.W. Leung / Nonlinear Analysis: Real World Applications 9 (2008) 2196 – 2213
2205
In fact, let u˜ be the solution of (1 + 2 − r)1 −u =0 u − cu + u 1 (1 + 2 ) − rb 2 −r)1 as → +∞. From Theorem 3.1, we have that u01 () ≤ u(). ˜ Moreover, with u˜ → 0 as → −∞, and u˜ → (1+ 1 (1+2 )−rb (u) ˜ is the solution of (1 + 2 − r)1 − 2u˜ v = 0. (3.16) v − cv + 1 (1 + 2 ) − rb
From [18], we know it has the following asymptotic behavior
r(1 − b) 2 (u) ˜ = exp c − c − 4(1 − ) [cˆ + o(1)], 2 1 (1 + 2 ) − rb
cˆ = 0 as → −∞.
Integrating (3.15) and (3.17) from −∞ to , with (−) sufficiently large positive, and noticing r r(1 − b) 2 c− c −4 1− < c − c2 − 4 1 − , 1 + 2 1 (1 + 2 ) − rb
(3.17)
(3.18)
we conclude that = 0. (Here, we use the fact that u01 () ≤ u() ˜ for all ∈ R). Note that (3.18) is valid because 1 −b 1 < . Next, we consider Eq. (3.14). Writing (3.14) as 1+2 1 (1+2 )−rb ( 2 ) − c( 2 ) − 1 2 = − √
we notice that
c+
c2 −4(1−r/(1+2 )) 2
b , 1 + 2 1
(3.19)
is not a characteristic of
( 2 ) − c( 2 ) − 1 2 = 0.
(3.20)
Eq. (3.20) has two independent solutions of the form: √ √ 12 = e
c+
c2 +41 2
,
22 = e
c−
c2 +41 2
.
Thus, as → −∞, 2 has the property: √ √ 2 ¯ + o(1)]e((c− c2 +41 )/2) 2 = ¯ [1 + o(1)]e((c+ c +41 )/2) + [1 √ 2 + [1 + o(1)]e((c+ c −4(1−r/(1+2 )))/2) ¯ and = 0. Since lim→−∞ 2 () = 0, we obtain ¯ = 0. Hence 2 = [ + o(1)] for some constants ¯ , , √ 2 e((c+ c −4(1−r/(1+2 )))/2) , as → −∞. Therefore, as → −∞, we have the formula: √ 2 ( + o(1))e((c+ c −4(1−r/(1+2 )))/2) 1 √ , (3.21) = 2 2 ( + o(1))e((c+ c −4(1−r/(1+2 )))/2) where = 0 and = 0. Similarly, we consider (3.11) and (3.12) as → +∞, and obtain the limiting equation: ( 1 ) − c( 1 ) + A 1 + B 2 = 0, ( 2 ) − c( 2 ) + C 1 + D 2 = 0,
(3.22)
2206
X. Hou, A.W. Leung / Nonlinear Analysis: Real World Applications 9 (2008) 2196 – 2213
where A, B, C, D are given by A=
−(1 + 2 − r)1 , 1 (1 + 2 ) − rb
B=
r1 (1 + 2 − r) , 1 (1 + 2 ) − rb
C=
b(1 − b) , 1 (1 + 2 ) − rb
D=
1 b(1 + 2 − r) − 1 . 1 (1 + 2 ) − rb
(3.23)
˜ , i = 1, 2, we can write (3.22) as a first order system of ordinary differential equation in the four By setting ( i ) = i ˜ , , ) ˜ T with a constant 4 × 4 coefficient matrix. One can readily obtain that the general solution components ( 1 , 1 2 of the first order system corresponding to (3.22) has the following form: ˜ , , ˜ )T = 4 ci hi ei , ( 1 , 1 2 2 i=1 where 1 = 3 =
c+ c+
c2 + 41 , 2 c2 + 42 , 2
with 1 =
−(A + D) +
(3.24)
c2 + 41 , 2 c − c2 + 42 4 = , 2
2 =
c−
(A − D)2 + 4BC , 2
2 =
(3.25)
−(A + D) −
(A − D)2 + 4BC , 2
(3.26)
and hi , i = 1, 2, 3, 4 are eigenvectors of the constant 4 × 4 coefficient matrix mentioned above with i as corresponding eigenvalues. (Here, ci are arbitrary constants.) To calculate the sign of 1 and 2 , we note that A < 0, and from (2.7) D < b − 1 < 0. Hence we have A + D < 0 and we have BC = AD =
r1 (1 + 2 − r)b(1 − b) [1 (1 + 2 ) − rb]2 −(1 + 2 − r)1 [1 (1 + 2 ) − rb]2
BC − AD =
√
(A+D)−
(A−D)2 +4BC 2
√
< 0. We next show
(A+D)+
(A−D)2 +4BC 2
< 0. From (3.23),
,
1 (1 + 2 )(−b + 1 ),
(1 + 2 − r)1 (1 − b) [1 (1 + 2 ) − rb]2
[rb − 1 (1 + 2 )] < 0.
The last inequality is true due to (2.7). Hence, (A − D)2 + 4BC < (A + D)2 . √ 2 +4BC Consequently, we have (A+D)+ (A−D) < 0. 2 Thus from (3.26), we have shown that BC < AD,
1 > 0, 2 > 0 and 1 > 2 . Consequently, we find from (3.25) that 1 > 0, 2 < 0, 3 > 0, and 4 < 0. ˜ , , ˜ )T → 0 as → ∞, we obtain from (3.24) that c1 = c3 = 0, and Since ( 1 , 1 2 2 ˜ , , ˜ )T = c2 h2 e2 + c4 h4 e4 . ( 1 , 1 2 2
(3.27)
X. Hou, A.W. Leung / Nonlinear Analysis: Real World Applications 9 (2008) 2196 – 2213
We deduce as above the following asymptotic behavior as → +∞: √2 √2 1 (P1 + o(1))e((c− c +41 )/2) + 2 (P2 + o(1))e((c− c +42 )/2) 1 √2 √2 = , 2 1 (Q1 + o(1))e((c− c +41 )/2) + 2 (Q2 + o(1))e((c− c +42 )/2)
2207
(3.28)
where Pi , Qi , i , i = 1, 2 are constants, and 1 , 2 cannot be zero simultaneously. (See e.g. [4]). From (2.7) and (3.23), we readily verify that B = 0, and C = 0. Consider the solution hi ei of the first order system corresponding to (3.22). If one of the first and third components of the eigenvector hi is zero, the system (3.22) implies that the other component is also zero. This leads to the fact that hi = 0, which is a contradiction. Consequently, we conclude that Pi = 0, Qi = 0, i = 1, 2 in (3.28). Note that we are unable to show that 2 = 0, i.e. eliminate the dominant solution, by comparing with the upper solution as in the case when → −∞. This should be possible, if further assumptions are made in addition to [H1]–[H3]. Theorem 3.2. Under the hypotheses of Theorem 3.1, there exist positive constants A1 , A2 , B1 and B2 , such that system (2.15) has a traveling wave solution U with the following asymptotic properties: √ 2 (A1 + o(1))e((c+ c −4(1−r/(1+2 )))/2) √ (3.29) U () = 2 (A2 + o(1))e((c+ c −4(1−r/(1+2 )))/2) as → −∞ and ⎛
√ ⎞ (1 + 2 − r)1 2 − (B1 + o(1))e((c− c +4)/2) 1 (1 + 2 ) − rb ⎟ ⎜ U () = ⎝ ⎠ √ b(1 + 2 − r) 2 c +4 )/2) ((c− − (B2 + o(1))e (1 + 2 )(1 (1 + 2 ) − rb)
(3.30)
as → +∞. Here is one of 1 or 2 given by (3.26). Proof. The conclusions follow readily from formulas (3.21) and (3.28) by integration.
Remark 3.1. Expressing U () in terms of the original variables (x, t) the traveling wave solution of (1.1) have the following asymptotic forms:
a1 a a 1 c 1 2 1 u(x, t) = [A1 + o(1)] exp , c + c2 + 4 1 − x + ca 1 t b1 2 a1 c2 d
a1 1 c 1 a2 a2 2 − [qA2 + o(1)] exp , v(x, t) = c+ c +4 1− x + ca 1 t c2 2 a1 c2 d as
a1 d x
+ ca 1 t → −∞.
a1 a1 1 a1 c2 − a2 c1 2 − [B1 + o(1)] exp , (c − c + 4) x + ca 1 t u(x, t) = b 1 c 2 − b 2 c1 b1 2 d a1 1 a2 b1 − a1 b2 + [qB 2 + o(1)] exp , v(x, t) = (c − c2 + 4) x + ca 1 t b1 c2 − b 2 c1 2 d as
a1 d x
(3.31)
+ ca 1 t → ∞, where is one of the positive numbers 1 or 2 , with
1 =
−(A + D) +
(A − D)2 + 4BC , 2
2 =
−(A + D) −
(A − D)2 + 4BC , 2
(3.32)
2208
X. Hou, A.W. Leung / Nonlinear Analysis: Real World Applications 9 (2008) 2196 – 2213
A=
b1 (c1 a2 − c2 a1 ) , a1 (c2 b1 − c1 b2 )
C=
b2 (a2 b1 − a1 b2 ) , qb1 (c2 b1 − c1 b2 )
B=
qc1 b1 (a1 c2 − a2 c1 ) , a12 (c2 b1 − c1 b2 )
D=
b2 (c2 a1 − a2 c1 ) a2 − . a1 (c2 b1 − c1 b2 ) a1
1 c2 −a2 c1 ) Theorem 3.3. Assume all the hypotheses of Theorem 3.1, with c > 2 ba11 (a (b1 c2 −b2 c1 ) , then the monotone traveling wave solutions of (2.15) with asymptotic property as → ±∞ described by (3.29) and (3.30) are unique except for a possible shift of origin. That is, if U1 () and U2 () are two traveling wave solutions with such properties, then there exists ¯ ∈ R, such that U1 (¯ + ) = U2 (), ∈ R. Proof. We are going to use sliding method to prove the conclusion. (See [2] for the application of the method to a single equation.) Let U1 () and U2 () be traveling wave solutions of system (2.15), with speed c, satisfying the properties as stated above. There exist positive constants Ai , Bi , i = 1, 2, 3, 4 such that for < − N , √ 2 (A1 + o(1))e((c+ c −4(1−r/(1+2 )))/2) √ U1 () = , (3.33) 2 (A2 + o(1))e((c+ c −4(1−r/(1+2 )))/2) √ 2 (A3 + o(1))e(c+ c −4(1−r/(1+2 ))/2) √ U2 () = (3.34) 2 (A4 + o(1))e(c+ c −4(1−r/(1+2 ))/2) and for > N,
√ ⎞ (1 + 2 − r)1 2 − (B1 + o(1))e((c− c +4)/2) 1 (1 + 2 ) − rb ⎟ ⎜ U1 () = ⎝ ⎠, √ b(1 + 2 − r) 2 − (B2 + o(1))e((c− c +4)/2) (1 + 2 )(1 (1 + 2 ) − rb) √ ⎞ ⎛ (1 + 2 − r)1 2 − (B3 + o(1))e((c− c +4)/2) 1 (1 + 2 ) − rb ⎟ ⎜ U2 () = ⎝ ⎠. √ b(1 + 2 − r) 2 − (B4 + o(1))e((c− c +4)/2) (1 + 2 )(1 (1 + 2 ) − rb) ⎛
(3.35)
(3.36)
The traveling wave solutions of system (2.15) are translation invariant; thus for any > 0, U1 () := U1 ( + ) is also a traveling wave solution of (2.15). By (3.29) and (3.30), the solution U1 ( + ) has the following asymptotic behaviors: √ √ ((c+ c2 −4(1−r/(1+2 )))/2) e((c+ c2 −4(1−r/(1+2 )))/2) (A + o(1))e 1 √ √ U1 () = (3.37) 2 2 (A2 + o(1))e((c+ c −4(1−r/(1+2 )))/2) e((c+ c −4(1−r/(1+2 )))/2) for − N,
√ √ ⎞ (1 + 2 − r)1 2 2 − (B1 + o(1))e((c− c +4)/2) e((c− c +4)/2) 1 (1 + 2 ) − rb ⎟ ⎜ U1 () = ⎝ ⎠ √ √ b(1 + 2 − r) 2 2 − (B2 + o(1))e((c− c +4)/2) e((c− c +4)/2) (1 + 2 )(1 (1 + 2 ) − rb) for N. It is clear that for large enough, we have √ 2 A1 e((c+ c −4(1−r/(1+2 )))/2) > A3 , √ 2 A2 e((c+ c −4(1−r/(1+2 )))/2) > A4 , √ 2 B1 e((c− c +4)/2) < B3 , √ 2 B2 e((c− c +4)/2) < B4 . ⎛
(3.38)
(3.39) (3.40) (3.41) (3.42)
X. Hou, A.W. Leung / Nonlinear Analysis: Real World Applications 9 (2008) 2196 – 2213
2209
Formulas (3.33)–(3.42) imply that for large enough, U1 () > U2 ()
(3.43)
for ∈ (−∞, −N ]∪[N + ∞). Here, the inequality “>” in (3.43) is interpreted to hold for each component. Similar interpretation will be used for “ ” between two vectors. We now consider system (2.15) on [−N, +N]. First, suppose U1 ()U2 () on [−N, +N ], then the function W () := U1 () − U2 () satisfies for some ∈ (0, 1): ⎡ *F
1
(U2 + 1 (U1 − U2 )),
⎢ *u1 W − cW + ⎢ ⎣ *F 2 (U2 + 2 (U1 − U2 )), *u1
⎤ *F1 (U2 + 1 (U1 − U2 )) ⎥ *u2 ⎥W = 0 ⎦ *F2 (U2 + 2 (U1 − U2 )) *u2
(3.44)
for ∈ (−N, N), and W (−N ) > 0, W (+N ) > 0. − → Here, F = (F1 , F2 ) is a function of two variables as described by formulas following (2.23). Since the system above is monotone, we can readily deduce by maximum principle that W > 0, on [−N, N ]. Consequently, we have U1 () > U2 () on R in this case. Second, suppose there are some points in (−N, N ) such that U1 () < U2 (). We then increase , that is shift U1 () further left so that U1 (−N) > U2 (−N ), U1 (N ) > U2 (N ). By the monotonicity of U1 and U2 , we can ¯ > U2 (). Shifting U ( + ) ¯ back unfind ¯ ∈ (0, 2N) such that in the internal (−N, N ), we have U1 ( + ) 1 ¯ first touches one component of U2 () at some point ¯¯ ∈ (−N, N ). Then by til one component of U ( + ) 1
maximum principle for that component again, we find that component of U1 and U2 are identically equal for all ∈ [−N, N ] for a larger than the original one such that (3.43) holds. This is a contradiction. Therefore, we must have U1 () > U2 ()
for all ∈ R, where is the one chosen by means of (3.39)–(3.42) as described above. Now, decrease until one of the following situations happens. ¯
¯ 1. There exists 0, such that U1 () ≡ U2 (). In this case we have finished our proof. ¯ ¯ 2. For 0, there exists 1 ∈ R, such that one of the component of U and U2 is equal at the point 1 ; and for all ¯ ∈ R, we have U1 ()U2 (). We then consider the system (3.44) for large interval (−N, N ) and = ¯ in the ¯
definition for W. Suppose that the first component of U1 and U2 is equal at the point 1 . The maximum principle for ¯
the first component implies that the first component of U1 () is identically equal to that of U2 (). Also, we readily obtain that for large +, the limiting equation for (3.44) is the same as (3.22) as deduced from (3.11) and (3.12). Since the first component of W is identically zero, and the off diagonal limit coefficient B in (3.22) is not equal to zero, we conclude from the first equation in (3.44) that the second component of W must vanish for all large . By the maximum principle for the second equation, we conclude that the second component of W is also identicaly zero ¯ for all ∈ R. We next consider the case that the second component of U1 and U2 is equal at the point 1 . We can deduce similarly using the off diagonal coefficient C in (3.22) is not zero that both components of W are identically zero for all ∈ R as before. Consequently, in either situation, there exists ¯ 0, such that ¯
U1 () ≡ U2 () for all ∈ R.
2210
X. Hou, A.W. Leung / Nonlinear Analysis: Real World Applications 9 (2008) 2196 – 2213
4. Numerical simulations In this section, we will give numerical examples to illustrate the results in Section 2. To this end we will use the monotone system: ⎧ 1 + 2 − r ⎪ (u1 ) = (u1 )x˜ x˜ + u1 , − u + ru ⎪ 1 2 ⎨ 1 + 2 (4.1) ⎪ 1 ⎪ ⎩ (u2 ) = (u2 ) + − u2 (bu1 − 1 (1 + 2 )u2 ) x˜ x˜ 1 + 2 and system (2.22). Recall that system (1.1) is related to (4.1) by means of (2.1), (2.11), (2.14) and (2.15). Let (u1 (), u2 ()), = x˜ + c, be a solution of (4.1). As described in the proof of Theorem 2.2, the profile of the function 1 := u1 (−), 2 := u2 (−) can be computed as a solution of the following ODE system: (1 ) + c(1 ) + 1 (2 ) + c(2 ) +
1 + 2 − r − 1 + r2 = 0, 1 + 2
1 − 2 (b1 − 1 (1 + 2 )2 ) = 0. 1 + 2
(4.2)
We will choose the following values for the parameters: b = 0.2, 1 = 0.8, 2 = 0.3, r = 0.5. The parameters satisfy conditions given in (2.7). They are consequences of hypotheses [H1] and [H2]. Note also that 1 = a1−1 a2 = 0.8. Thus hypothesis [H3] is also satisfied. From the proof of Theorem 2.2, system (4.1) has a traveling wave solu-
2 −r)1 1 c2 −a2 c1 ) tion for every c > (1+ = 2 ba11 (a (b1 c2 −b2 c1 ) , connecting the equilibria (0, 0) and (0.6809, 0.1157). As an 1 (1+2 )−rb example for (1.1) with such corresponding parameters, we may set d = a1 = b1 = c1 = 1, a2 = 0.8, b2 = 0.2, and c3 = 2.08. In Figs. 3 and 4 , the initial value for (4.1) is chosen to be u1 (x, ˜ 0) = exp(−20x), ˜ u2 (x, ˜ 0) = exp(−20x) ˜ for x˜ ∈ [0, 100] and vanishing outside. After some time later, the profiles are shown in the three diagrams in Fig. 4. We use examples to analyze the stability of the traveling wave with various initial values. The smaller the initial values, the faster the solutions of (4.1) converges to a traveling wave solution. (See Figs. 5 and 6). This means the traveling wave solution obtained in Theorem 2.2 is unlikely to be stable in the usual Banach spaces such as L2 (R × R) or L∞ (R × R). We will study this issue in a forthcoming paper.
Fig. 3. Numerical traveling wave solutions, u1 and u2 ; t-axis represents , x-axis represents x. ˜
X. Hou, A.W. Leung / Nonlinear Analysis: Real World Applications 9 (2008) 2196 – 2213
2211
1 0.9
0.9
0.8
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0
0
10
20
30
40
50
60
70
80
90
0
100
0
10
20
30
40
60
70
80
90
100
50
60
70
80
90
100
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0
10
20
30
40
50
˜ The wave propagates to the left. Solid Fig. 4. Profiles of the traveling wave solutions at three different time , where the horizontal axis represents x. curve represents u1 and dashed curve represents u2 .
Fig. 5. Numerical traveling waves; the t-axis represents and the x-axis represents x. ˜
In deriving Fig. 6, we choose larger initial values to be (u1 (x, ˜ 0), u2 (x, ˜ 0))=(exp(−2x), ˜ exp(−2x)) ˜ for x˜ ∈ [0, 100] and vanishing outside. The four diagrams in Fig. 6 show that it takes longer for the profile to approach to that of the traveling wave in the fourth diagram.
2212
X. Hou, A.W. Leung / Nonlinear Analysis: Real World Applications 9 (2008) 2196 – 2213
1
1
0.9
0.9
0.8
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0
0
10
20
30
40
50
60
70
80
90
100
0
1
1
0.9
0.9
0.8
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0
0
10
20
30
40
50
60
70
80
90
100
0
0
10
20
30
40
50
60
70
80
90
100
0
10
20
30
40
50
60
70
80
90
100
Fig. 6. Profiles of the solution at different time with larger initial condition, where the horizontal axis represents x. ˜ The wave eventually propagates to the left. Solid curve represents u1 and dashed curve represents u2 .
References [1] S. Ahmad, A. Lazer, An elementary approach to traveling front solutions to a system of N competition-diffusion equations, Nonlinear Anal. 16 (10) (1991) 893–901. [2] H. Berestycki, B. Nicolaenko, B. Scheurer, Travelling wave solutions to combustion models and their singular limits, SIAM J. Math. 16 (6) (1985) 1207–1242. [3] H. Berestycki, L. Nirenberg, Travelling fronts in cylinders, Ann. Inst. H. Poincare Anal. Non Lineaire 9 (5) (1992) 497–572. [4] C. Conley, R. Gardner, An application of the generalized Morse index to traveling wave solutions of a competitive reaction diffusion model, Indiana Univ. Math. J. 33 (1984) 319–345. [5] W. Coppel, Stability and Asymptotic Behavior of Differential Equations, Heath Mathematical Monographs, D.C. Heath and Company, Boston, 1965. [6] S. Dunbar, Traveling wave solutions of diffusive Lotka–Volterra equations: a heteroclinic connection in R 4 , Trans. Am. Math. Soc. 286 (2) (1984) 557–594. [7] R. Gardner, Existence of traveling wave solutions of competing models, A degree theoretic approach, J. Differential Equations 44 (1982) 343–364. [8] T. Hallam, K. Lika, Traveling wave solutions of a nonlinear reaction–advection equation, J. Math. Biol. 38 (4) (1999) 346–358. [9] D. Henry, Geometric Theory of Semilinear Parabolic Equations, Lecture Notes in Mathematics, vol. 840, Springer, NY, 1981. [10] P. Hess, Periodic-Parabolic Boundary Value Problems and Positivity, Pitman, R.N.M., Longman, Harlow, 1991. [11] X. Hou, Y. Li, Local stability of traveling-wave solutions of nonlinear reaction–diffusion equations, Discrete Continuous Dynamical Syst. 15 (2) (2006) 681–701. [12] J. Huang, X. Zou, traveling wavefronts in diffusive and cooperative Lotka–Voltera system with delays, J. Math. Anal. Appl. 271 (2002) 455–466. [13] K. Kirchgassner, On the nonlinear dynamics of travelling fronts, J. Differential Equations 96 (1992) 256–278. [14] A. Leung, Systems of Nonlinear Partial Differential Equations: Applications to Biology and Engineering, MIA, Kuwer, Boston, 1989. [15] B. Li, H. Weinberger, M. Lewis, Spreading speeds as slowest wave speeds for cooperative systems, Math. Biosci. 196 (2005) 82–98.
X. Hou, A.W. Leung / Nonlinear Analysis: Real World Applications 9 (2008) 2196 – 2213 [16] [17] [18] [19] [20] [21] [22]
2213
R. Lui, Biological growth and spread modeled by systems of recursions, I. Mathematical Theory, Math. Biosci. 93 (1989) 269–295. C. Pao, Nonlinear Parabolic and Elliptic Equations, Plenum Press, NY, 1992. D.H. Sattinger, On the stability of waves of nonlinear parabolic systems, Adv. Math. 22 (1976) 312–355. H. Smith, Monotone Dynamical Systems: An Introduction to the Theory of Competitive and Cooperative Systems, AMS Providence, RI, 1995. J. Smoller, Shock Waves and Reaction-Diffusion Equations, Springer, NY, 1982. M. Tang, P. Fife, Propagating fronts for competing species with diffusion, Arch. Ration. Mech. Anal. 73 (1978) 69–77. J. Van Vuuren, The existence of traveling plane waves in a general class of competition–diffusive systems, IMA J. Appl. Math. 55 (1995) 135–148. [23] I. Volpert, V. Volpert, V. Volpert, Traveling Wave Solutions of Parabolic Systems, Translations of Mathematics Monograph, vol. 140, American Mathematical Society, Providence, RI, 1994. [24] J. Wu, X. Zou, Traveling wave fronts of reaction-diffusion systems with delay, J. Dyn. Differential Equations 13 (3) (2001). [25] J. Xin, Front propagation in heteogeneous media, SIAM Rev. 42 (2) (2000) 161–230.