Existence of traveling wave solutions for diffusive predator-prey type systems

(1)

Contents lists available atSciVerse ScienceDirect

Journal of Differential Equations

www.elsevier.com/locate/jde

Existence of traveling wave solutions for diffusive

predator–prey type systems

Cheng-Hsiung Hsu

a

,

1

, Chi-Ru Yang

b

, Ting-Hui Yang

c

,

2

, Tzi-Sheng Yang

d

,

∗,

2 a_{Department of Mathematics, National Central University, Chung-Li 32001, Taiwan}

b_{Department of Applied Mathematics, National Chiao Tung University, Hsin-Chu 30010, Taiwan} c_{Department of Mathematics, Tamkang University, Tamsui, Taipei County 25137, Taiwan} d_{Department of Mathematics, Tunghai University, Taichung 40704, Taiwan}

a r t i c l e

i n f o

a b s t r a c t

Article history:

Received 21 September 2010 Revised 18 October 2011 Available online 30 November 2011

Keywords:

Traveling wave Predator–prey Wazewski Theorem LaSalle’s Invariance Principle Lyapunov function Hopf bifurcation theory

In this work we investigate the existence of traveling wave solu-tions for a class of diffusive predator–prey type systems whose each nonlinear term can be separated as a product of suitable smooth functions satisfying some monotonic conditions. The pro-ﬁle equations for the above system can be reduced as a four-dimensional ODE system, and the traveling wave solutions which connect two different equilibria or the small amplitude traveling wave train solutions are equivalent to the heteroclinic orbits or small amplitude periodic solutions of the reduced system. Applying the methods of Wazewski Theorem, LaSalle’s Invariance Principle and Hopf bifurcation theory, we obtain the existence results. Our results can apply to various kinds of ecological models.

1. Introduction

This work concerns with the existence of traveling wave solutions for the following diffusive predator–prey type system:

ut

=

d1uxx

−

h

(

u

)

g

(

w

)

−

p

(

u

)

,

wt

=

d2wxx

− (

w

)

q

(

u

),

(1.1)

*

Corresponding author.

E-mail addresses:[email protected](C.-H. Hsu),[email protected](C.-R. Yang),[email protected] (T.-H. Yang),[email protected](T.-S. Yang).

1 _{Research supported in part by NSC, NCTS of Taiwan.} 2 _{Research supported in part by NSC of Taiwan.}

(2)

where d1

>

0, d2

>

0, p

(

u

)

, g

(

w

)

, h

(

u

)

,

(

w

)

and q

(

u

)

are smooth functions satisfying some

mono-tonic conditions which will be mentioned later. System (1.1) is a general form of the diffusive predator–prey system which contains many known models. Indeed, system (1.1) describes not only the interspecies relations for ecological and social models, but also the base block of more compli-cated food web, food chain and biochemical network structure. In ecology, the functions u

(

x

,

t

)

and

w

(

x

,

t

)

represent the species densities of the prey and predator, respectively; the constants d1 and d2

are the spatial diffusion rates of the two species; the function h

(

u

)

p

(

u

)

is the net growth rate of the prey in the absence of predator; the function h

(

u

)

is the predator functional response which describes consumption rate of prey by a unit number of predators; the graphs g

(

w

)

−

p

(

u

)

=

0 and q

(

u

)

=

0 are the prey nullcline and predator nullcline on the phase portrait, respectively. In the sequel, we will illustrate some models where the existence of traveling wave solutions has been studied in the past decades.

In 1983, Dunbar [4,5] considers the existence of traveling wave solutions for the following reaction–diffusion system based on the Lotka–Volterra differential equation model of a predator–prey interaction:

⎧

⎨

⎩

ut

=

d1uxx

+

Au

1

−

u K

−

Bu w

,

wt

=

d2wxx

−

C w

+

Du w

,

(1.2)

where d1, d2, A, B, C , D, K are positive constants. A is the intrinsic rate of increasing for the prey

species; C is the death rate for the predator in the absence of the prey; K is the carrying capacity of the environment; the predator functional response here is the identity function of u. By using the Wazewski Theorem (an extension of shooting argument in higher dimension) together with a Lya-punov function and LaSalle’s Invariance Principle, he proves the existence of traveling wave solutions. Dunbar [6] further considered the existence of traveling wave solutions for system (1.2) but with Holling type II functional response H2

(

u

)

=

1+uEu, i.e.,

⎧

⎨

⎩

ut

=

d1uxx

+

Au

1

−

u K

−

B H2

(

u

)

w

,

wt

=

d2wxx

−

C w

+

D H2

(

u

)

w

,

(1.3)

where E

>

0. System (1.3) includes the effects of predation satiation: the consumption rate of prey by a unit number of predators cannot continue to grow linearly with the number of prey available but must “saturate” at some value (see [8,9]). The parameter 1

/

E here is the satiation rate of predation.

Assume d1

=

0, Dunbar uses the method similar to that in [4,5] and the invariant manifold theory to

prove the existence of traveling wave train and traveling front solutions for system (1.3). The case for

d1

=

0 is then considered by Huang, Lu and Ruan [12]. Using the same shooting argument and the

Hopf bifurcation theory, they establish the existence of the traveling wave solutions connecting two rest states as well as the existence of small amplitude traveling wave train solutions.

Later, Li and Wu [16] also consider the system (1.3) but with Holling type III functional response

H3

(

u

)

=

u 2 1+Eu2, i.e.,

⎧

⎨

⎩

ut

=

d1uxx

+

Au

1

−

u K

−

B H3

(

u

)

w

,

wt

=

d2wxx

−

C w

+

D H3

(

u

)

w

.

(1.4)

By using the similar methods of [4,5], they establish the existence of traveling wave solutions of (1.4) for the case d1

=

0. In this work, we generalize the results of [16] to the case d1

=

0.

In addition to the previous Holling types of functional responses, Ivlev in 1961 [14] introduces another functional response H4

(

u

)

=

E

(

1

−

e−Mu

)

where E represents the maximum rate of predation

(3)

and M is a constant representing the decrease in motivation to hunt. The diffusive predator–prey model with logistic growth rate of prey and Ivlev type functional response is described by

⎧

⎨

⎩

ut

=

d1uxx

+

Au

1

−

u K

−

B H4

(

u

)

w

,

wt

=

d2wxx

−

C w

+

D H4

(

u

)

w

.

(1.5)

If d1

=

d2

=

0, system (1.5) is studied by many authors, see [1,2,15,17,18,20,22,24,25]. Most of these

papers concentrate on the existence and stability of limit cycle. Recently, in [23], Wang, Shi and Wei also study the global bifurcation of a class of more general predator–prey models with a strong Allee effect in prey population. On the other hand, if d1

=

0 and d2

=

0, there seems no results for the

existence of traveling wave solution of system (1.5). In Section 5.4 of this work, we will apply our main theorem to obtain the new existence results for traveling wave solutions of system (1.5).

For other examples, Owen and Lewis [19] consider the following general system

ut

=

εα

0uxx

+

α

1u f1

(

u

)

−

α

2w f2

(

u

),

wt

=

α

0wxx

+

α

3w f2

(

u

)

−

α

4w

,

(1.6)

where

ε

≈

0 and

α

i’s are positive constants. They study the mechanism for which predation pressure can slow, stall or reverse a spatial invasion of prey. Some numerical results of traveling wave solutions are demonstrated in [19] for speciﬁc fi’s described below. The function f1is given by f1

(

u

)

= (

1

−

u

)

or f1

(

u

)

=

k

(

1

−

u

)(

u

−

a

)

for some constants k and a; while f2is given by Holling type I ( f2

(

u

)

=

u),

type II, or type III functional response. However there is no theoretical proof for their numerical results.

Motivated by the above models, throughout this article, we consider p

(

u

)

, g

(

w

)

, h

(

u

)

,

(

w

)

and

q

(

u

)

to be C1 _{functions satisfying the following assumptions:}

(A1) p

(

u

) <

0 for u

>

0, and p

(

K

)

=

0 for some u

=

K

>

0. (A2) q

(

u

) <

0 for u

>

0, and q

(

u_∗

)

=

0 for some u_∗

∈ (

0

,

K

)

.

(A3) g

(

w

) >

0,

(

w

) >

0,

(

w

)

0 for w

∈ R

, g

(

0

)

= (

0

)

=

0 and g

(

∞) = (∞) = ∞

. (A4) h

(

0

)

=

0 and h

(

u

) >

0 for u

∈ R

.

Note that (A1)–(A4) hold for the systems (1.2)–(1.6) provided the corresponding parameters lying in suitable regions. For example, let p

(

u

)

=

A

(

1

−

u

/

K

)

, g

(

w

)

=

B w, h

(

u

)

=

u,

(

w

)

=

w and q

(

u

)

=

C

−

Du for (1.2), then (A1)–(A4) hold if C

/

D

<

K .

For further simpliﬁcation, we introduce the parameter d

=

d1

/

d2and rescale the spatial variable x

by

˜

x

=

x

/

√

d2

.

Then system (1.1) is recast as (still using x instead ofx)

˜

ut

=

duxx

−

h

(

u

)

g

(

w

)

−

p

(

u

)

,

wt

=

wxx

− (

w

)

q

(

u

).

(1.7) According to assumptions (A1)–(A4), it is easy to see that system (1.7) has three spatially uniform equilibria: E0

= (

0

,

0

)

, E1

= (

K

,

0

)

, and E2

= (

u∗

,

w_∗

)

where

w_∗

=

g−1

◦

p

(

u_∗

) >

0

.

Note that E0 corresponds to the absence of both species; E1 corresponds to the prey being at the

environment carrying capacity in the absence of the predator; and E2 corresponds to the coexistence

of the two species. The purpose of this work is to establish the traveling wave solutions of system (1.7) connecting the equilibria E1 and E2, which is called the “wave of invasion”, cf. [3].

A traveling wave solution of (1.7) is a solution of the form

(4)

where the constant c

>

0 is the wave speed; s

=

x

+

ct is called the moving coordinate.

Substitut-ing (1.8) into (1.7), we have the followSubstitut-ing proﬁle equations:

cu

=

du

−

h

(

u

)

g

(

w

)

−

p

(

u

)

,

c w

=

w

− (

w

)

q

(

u

),

(1.9)

where denotes the differentiation with respect to s. It is required that u and w of system (1.7) are nonnegative for natural ecological restriction. Then we look for the nonnegative solutions of (1.9) connecting the equilibria E1 and E2, i.e., satisfying the following boundary conditions:

u

(

−∞) =

K

,

w

(

−∞) =

0

,

u

(

∞) =

u_∗

,

and w

(

∞) =

w_∗

.

(1.10) Our main results are stated as follows.

Theorem 1.1. Assume (A1)–(A4) hold, and let d

<

1, c∗

:=

√

−

4

(

0

)

q

(

K

).

(i) If 0

<

c

<

c∗, then there is no nonnegative traveling wave solution of system (1.7) connecting the equilibria E1and E2.

(ii) If c

>

c∗,

(

w

)

=

αg

(

w

)

and q

(

u

)

= β(

h

(

u

)

−

h

(

u_∗

))

for some

α

>

0 and

β <

0, then there is a

nonneg-ative traveling wave solution of (1.7) connecting the equilibria E1and E2. Furthermore, there exists a

σ

∗

>

0 such that

(1) if

|

α

β

| <

σ

∗, then the traveling wave solutions approach E2monotonically for large s;

(2) if

|

α

β

| >

σ

∗, then the traveling wave solutions have exponentially damped oscillations about E2for large s.

Extending the ideas of [4,5], we apply the Wazewski Theorem (see Theorem 2.3) together with LaSalle’s Invariance Principle (see [11]) to prove Theorem 1.1. Note that although we apply the tech-niques similar to those of [4,5], there are some differences. First, the model that we consider is more general, and our results contain (or extend) all the results of [4,5,12,19] and some other models, e.g., the predator–prey model with Ivlev’s functional response (1.5) and some typical S.I.R. models, such as Kermack–McKendrick model (cf. [7]). Second, due to the general setting of system (1.1), the con-struction of Wazewski set is more complicated than those of [4,5], and it’s more diﬃcult to ﬁnd an invariant orbit of system (1.9) in the Wazewski set. Third, we construct the Lyapunov function for system (1.1) more generally to prove the existence results.

According to Theorem 1.1, we know that

c∗

=

2

√

D K

−

C

,

2

D H2

(

K

)

−

C

,

2

D H3

(

K

)

−

C

for systems (1.2), (1.3) and (1.4) respectively. Note that for speciﬁc form of system (1.2), Dunbar [4] pointed out that c∗ is a distinguished speed dividing the positive traveling wave solutions into two types: wave of speed c

<

c∗ being one type connecting E0 and E2, wave of speed c

c∗ being of

the other type connecting E1 and E2. In our case, the existence of positive traveling wave solutions

connecting E0 and E2 is still open, and will be in our further study.

This paper is organized as follows. In Section 2, we recall a variant of Wazewski Theorem and construct the Wazewski set. Then we use the standard Stable Manifold Theorem to investigate the behavior of solutions for system (1.9) in the 4-dimensional phase space and prove that there is an invariant solution orbit in the Wazewski set. In Section 3, we construct the Lyapunov function for the invariant orbit. In Section 4, we prove Theorem 1.1 by using LaSalle’s Invariance Principle. In Section 5, we apply our main theorem to systems (1.2)–(1.5). We further investigate the existence of traveling wave train solutions for these systems by using the Hopf bifurcation theory. The technical proofs for Proposition 2.4 and Lemma 2.18 are given in Appendices A and B respectively.

(5)

2. Construction of Wazewski set and invariant orbit

In this section, we will apply the Wazewski Theorem to prove that there is an orbit invariant in a bounded region containing E1 and E2. First, let’s rewrite system (1.9) as a system of ﬁrst order ODEs

in

R

4_,

⎧

⎪

⎨

⎪

⎩

u

=

v

,

dv

=

cv

+

h

(

u

)

g

(

w

)

−

p

(

u

)

,

w

=

z

,

z

=

cz

+ (

w

)

q

(

u

).

(2.1)

Then the boundary conditions (1.10) yield

u

(

−∞) =

K

,

v

(

−∞) =

0

,

w

(

−∞) =

0

,

z

(

−∞) =

0

,

u

(

∞) =

u_∗

,

v

(

∞) =

0

,

w

(

∞) =

w_∗

,

z

(

∞) =

0

.

(2.2)

It’s obvious that

H

:=

(

u

,

v

,

w

,

z

)

: u

=

v

=

0

and

V

:=

(

u

,

v

,

w

,

z

)

: w

=

z

=

0

are invariant manifolds of (2.1). The eigenvalues of the linearization of (2.1) at

(

K

,

0

,

0

,

0

)

are

λ

1

=

c

+

c2

₋

_4dh

₍

_K

₎

_p

₍

_K

₎

2d

>

0

,

λ

2

=

c

+

c2

₊

₄

₍

₀

₎

_q

₍

_K

₎

2

,

λ

3

=

c

−

c2

₊

₄

₍

₀

₎

_q

₍

_K

₎

2

,

λ

4

=

c

−

c2

₋

_4dh

₍

_K

₎

_p

₍

_K

₎

2d

<

0

.

The corresponding eigenvectors are given by

e1

= (−1

,

−λ

1

,

0

,

0

),

e2

=

−

1

,

−λ

2

,

−ψ(λ

2

),

−λ

2

ψ (λ

2

)

,

e3

=

−

1

,

−λ

3

,

−ψ(λ

3

),

−λ

3

ψ (λ

3

)

,

e4

= (−

1

,

−λ

4

,

0

,

0

),

(2.3) where

ψ (λ)

=

1 g

(

0

)

h

(

K

)

d

λ

2

−

c

λ

+

h

(

K

)

p

(

K

)

.

(2.4)

Note that

λ

1 and

λ

4 satisfy the equation

d

λ

2

−

c

λ

+

h

(

K

)

p

(

K

)

=

0

;

(2.5)

λ

2 and

λ

3 satisfy the equation

λ

2

−

c

λ

−

(

0

)

q

(

K

)

=

0

.

Let d

<

1. If c2

_<

₋

₄

₍

₀

₎

_q

₍

_K

₎

_{, then}

_λ

2 and

λ

3 are complex conjugate eigenvalues and

λ

1

>

Re

λ

2

=

Re

λ

3

>

0

> λ

4. Thus there is a 1-dimensional strongest unstable manifold, which is tangent to e1 at

(

K

,

0

,

0

,

0

)

. This manifold is actually contained in the invariant manifold

V

. Therefore a solution of

(6)

(2.1)–(2.2) cannot lie in the strongest unstable manifold. It follows that a solution of (2.1)–(2.2) must tend spirally to

(

K

,

0

,

0

,

0

)

. Hence w

(

s

) <

0 for some s. Therefore, there is no nonnegative solution of (2.1)–(2.2). The part (i) of Theorem 1.1 is then proved.

On the other hand, if c2

>

−

4

(

0

)

q

(

K

)

, then it’s obvious that

λ

1

> λ

2

> λ

3

>

0

> λ

4. Note that

ψ(λ

2

) <

0 and

ψ(λ

3

) <

0.

To investigate the structure of the eigenvalues at

(

u_∗

,

0

,

w_∗

,

0

)

, we recall the Routh–Hurwitz Sta-bility Criterion. Consider the polynomial equation

anxn

+

an−1xn−1

+ · · · +

a1x

+

a0

=

0

.

The Routh array for the above equation is deﬁned by

⎛

⎜

⎝

an an−2 an−4 an−6

. . .

an−1 an−3 an−5 an−7

. . .

b1 b2 b3 b4

. . .

c1 c2 c3 c4

. . .

..

.

..

.

..

.

..

.

. . .

⎞

⎟

⎠

where bk

= −

1 an−1

an an−2k an−1 an−2k−1

,

ck

= −

1 b1

an−1 an−2k−1 b1 bk+1

and so on. For example, the Routh array for a four-degree polynomial (n

=

4) is given by

⎛

⎜

⎝

a4 a2 a0 0 a3 a1 0 0 b1 b2 0 0 c1 c2 0 0 d1 0 0 0

⎞

⎟

⎠

where b1

= −

1 a3

a4 a2 a3 a1

,

b2

= −

1 a3

a4 a0 a3 0

,

c1

= −

1 b1

a3 a1 b1 b2

,

c2

= −

1 b1

a3 0 b1 0

=

0

,

d1

= −

1 c1

b1 b2 c1 c2

.

With the Routh array, the Routh–Hurwitz Stability Criterion [10,13,21] tells us how many roots having positive real parts.

Proposition 2.1. The number of sign changes in the ﬁrst column of the Routh array equals to the number of

roots with positive real parts.

Now we consider the characteristic equation of the linearization of (2.1) at

(

u_∗

,

0

,

w_∗

,

0

)

, i.e.,

λ

4

−

c

+

c d

λ

3

+

c 2

_{− ξ}

∗ d

λ

2

₊

c

ξ

∗ d

λ

+

ζ

_∗ d

=

0

,

(2.6)

where

ξ

_∗

= −

h

(

u_∗

)

p

(

u_∗

) >

0 and

ζ

_∗

= −(

w_∗

)

h

(

u_∗

)

q

(

u_∗

)

g

(

w_∗

) >

0. Applying Proposition 2.1, we have the following lemma.

(7)

Lemma 2.2. Eq. (2.6) has two eigenvalues with positive real parts and two eigenvalues with negative real

parts.

Proof. After simple computation, we have the following Routh array for Eq. (2.6)

⎛

⎜

⎝

1

(

c2

_{− ξ}

∗

)/

d

ζ

∗

/

d 0

−

c

−

c

/

d c

ξ

_∗

/

d 0 0

(

c2

− ξ

_∗

)/

d

− ξ

_∗

/(

d

+

1

)

ζ

_∗

/

d 0 0 c

ξ

∗

/

d

+

c

(

1

+

d

)ζ

∗

/(

b1d2

)

0 0 0

ζ

_∗

/

d 0 0 0

⎞

⎟

⎠

,

where b1

= (

c2

− ξ∗

)/

d

− ξ∗

/(

d

+

1

)

. It can be veriﬁed that the signs of ﬁrst column always change

twice. Hence Eq. (2.6) has two roots with positive real parts. On the other hand, if we replace

λ

by

iωin Eq. (2.6) then we have

ω

2

= −ξ

∗

/(

1

+

d

) <

0

,

which is a contradiction. Hence there is no pure imaginary roots. The proof is complete.

2

2.1. Wazewski Theorem

We now recall a variant of Wazewski Theorem which is a formalization and extension of the shooting method in higher dimension (see Proposition 2 of [5]).

Let us consider the differential equation:

y

(

s

)

=

f

y

(

s

)

,

(2.7)

where f :

R

n

→ R

nis a Lipschitz continuous function. Denote y

(

s

;

y0

)

as the unique solution of (2.7)

with initial value y

(

0

)

=

y0. For convenience, the notation y0

·

s stands for y

(

s

;

y0

)

and y0

·

S for the

set of points y

·

s with s

∈

S

⊂ R

. Now we deﬁne the following sets.

◦

Given W

⊆ R

n, we deﬁne the immediate exit set W− of W by

W−

=

y0

∈

W :

∀

s

>

0

,

y0

· [0

,

s

)

W

.

◦

Given

Σ

⊆

W , we set

Σ

0

_{= {}

_y

0

∈ Σ

:

∃

s0

>

0 such that y0

·

s0

∈

/

W

}.

◦

For y0

∈ Σ

0, we deﬁne the exit time T

(

y0

)

of y0by

T

(

y0

)

=

sup

s: y0

· [

0

,

s

] ⊂

W

.

Note that y0

·

T

(

y0

)

∈

W− and T

(

y0

)

=

0 if and only if y0

∈

W−. The Wazewski Theorem is stated as

the following.

Theorem 2.3. Consider Eq. (2.7). Suppose that

(i) if y0

∈ Σ

and y0

· [

0

,

s

] ⊆

c

(

W

)

, then y0

· [

0

,

s

] ⊆

W ;

(ii) if y0

∈ Σ

, y0

·

s

∈

W and y0

·

s

∈

/

W−, then there is an open set Vsabout y0

·

s disjoint from W−;

(iii)

Σ

= Σ

0,

Σ

is a compact set and intersects a trajectory of y

=

f

(

y

)

only once. Then the mapping F

(

y0

)

=

y0

·

T

(

y0

)

is a homeomorphism from

Σ

to its image on W−.

(8)

Fig. 1. The projection of P , Q , R, and S in the u w-plane.

2.2. The exit set W−

According to Theorem 2.3, the idea for choosing a Wazewski set for (2.1) is to exclude the region where the trajectories will go to inﬁnity. The vector ﬁeld of system (2.1) leads us to exclude the region where v and v(or z and z, resp.) are both positive or negative. Thus, we set W (see Fig. 1) by

W

= R

+

⊕ R

3

\ (

P

∪

Q

∪

R

∪

S

),

(2.8) where P

=

(

u

,

v

,

w

,

z

)

: 0

<

u

<

u_∗

,

w

>

w_∗

,

z

>

0

,

Q

=

(

u

,

v

,

w

,

z

)

: u

>

u_∗

,

w

<

w_∗

,

z

<

0

,

R

=

(

u

,

v

,

w

,

z

)

: 0

<

u

<

u_∗

,

g

(

w

)

−

p

(

u

) <

0

,

v

<

0

,

S

=

(

u

,

v

,

w

,

z

)

: u

>

u_∗

,

g

(

w

)

−

p

(

u

) >

0

,

v

>

0

.

Note that in the block P (or Q

∩ {

w

>

0

}

, resp.) z

→ ∞

(or z

→ −∞

, resp.); in the block S (or R, resp.) v

→ ∞

(or v

→ −∞

, resp.); the set W is the complement of the four blocks P , Q , R, S in

R

+

_{⊕ R}

3_{. It is easy to see that}

∂

W

= (∂

P

\

R

)

∪ (∂

Q

\

S

)

∪ (∂

S

\

Q

)

∪ (∂

R

\

P

),

since P

∩

R

= ∅

, and Q

∩

S

= ∅

. Using the phase space analysis, the structure of W− is described in the following proposition.

Proposition 2.4. The exit set W−is given by

W−

= ∂

W

\

(

u_∗

,

0

,

w_∗

,

0

)

∪ (

K

,

0

,

0

,

0

)

∪

J1

∪

J2

,

where J1

=

J10

∪

J11

∪

J12

∪

J13

,

J10

=

(

u

,

v

,

w

,

z

)

: u

u_∗

,

v

>

0

,

w

=

z

=

0

,

J11

=

(

u

,

v

,

w

,

z

)

: u

=

u_∗

,

v

>

0

,

w

<

0

,

z

=

0

,

(9)

J12

=

(

u

,

v

,

w

,

z

)

: u

>

u_∗

,

v

<

0

,

w

<

0

,

z

=

0

,

J13

=

(

u

,

v

,

w

,

z

)

: u

>

u_∗

,

v

0

,

w

<

0

,

z

=

0

,

g

(

w

)

−

p

(

u

) <

0

,

J2

=

(

u

,

v

,

w

,

z

)

: u

=

v

=

0

,

w

∈ R,

z

∈ R

.

Proof. The proof is tedious and illustrated in Appendix A.

2

2.3. Construction of

Σ

By the standard Stable Manifold Theorem, there is a 1-dimensional strongest unstable manifold

Ω

1

tangent to e1at

(

K

,

0

,

0

,

0

)

, and a parametric representation for this manifold in a small neighborhood

of

(

K

,

0

,

0

,

0

)

given by F1

(

ε

1

)

= (

K

,

0

,

0

,

0

)

+

ε

1e1

+

O

|

ε

1

|

2

.

There is also a 2-dimensional strongly unstable manifold

Ω

2 tangent to the linear subspace spanned

by e1 and e2at

(

K

,

0

,

0

,

0

)

, and a parametric representation for this manifold in a small neighborhood

of

(

K

,

0

,

0

,

0

)

given by F2

(

ε

1

,

ε

2

)

= (

K

,

0

,

0

,

0

)

+

ε

1e1

+

ε

2e2

+

O

|

ε

1

|

2

+ |

ε

2

|

2

.

Finally, the 3-dimensional unstable manifold

Ω

3 at

(

K

,

0

,

0

,

0

)

has a parametric representation in a

small neighborhood of

(

K

,

0

,

0

,

0

)

given by

F3

(

ε

1

,

ε

2

,

ε

3

)

= (

K

,

0

,

0

,

0

)

+

ε

1e1

+

ε

2e2

+

ε

3e3

+

O

|

ε

1

|

2

+ |

ε

2

|

2

+ |

ε

3

|

2

.

Throughout the rest of this article, y

(

s

;

y0

)

stands for the solution of (2.1) with initial value y0

=

(

u0

,

v0

,

w0

,

z0

)

; u

(

s

;

y0

)

stands for the u coordinate of y

(

s

;

y0

)

, and similarly for the other three

coordinates of y.

For y0

∈ Ω

1, we have the following properties.

Lemma 2.5. Let y

(

s

;

y0

)

be the solution of (2.1) with y0

∈ Ω

1and 0

<

u0

<

K . Then there is a ﬁnite s0

>

0 such that u

(

s0

;

y0

) <

u∗and v

(

s

;

y0

) <

0 for s

∈ [

0

,

s0

]

. That is, the solution enters region R.

Proof. Since e1

∈

V

is an invariant manifold, it follows that

Ω

1

⊂

V

. Thus, to investigate the dynamics

of solutions on

Ω

1, we may let w

=

z

=

0 in (2.1). Let us ﬁx a y0

∈ Ω

1 closed to

(

K

,

0

,

0

,

0

)

. The

parametrization F1 of

Ω

1 implies that there exists m

>

n

>

0 such that y0 lies between the two

curves: v

=

m

(

h

(

u

)

−

h

(

K

))

and v

=

n

(

h

(

u

)

−

h

(

K

))

. If m and n are large and small enough respectively, then we claim that y

(

s

;

y0

)

always lies between those two curves until u

=

u∗. We prove the claim

by contradiction. Suppose that there is an s1

>

0 such that v

=

m

(

h

(

u

)

−

h

(

K

))

and

(

v

−

m

(

h

(

u

)

−

h

(

K

)))

0 at s

=

s1(where u∗

<

u

(

s1

) <

K ), then we have

0

v

(

s1

)

−

mh

u

(

s1

)

v

(

s1

)

=

cv

(

s1

)

−

h

u

(

s1

)

p

u

(

s1

)

−

mh

u

(

s1

)

v

(

s1

)

= −

h

u

(

s1

)

h

u

(

s1

)

−

h

(

K

)

m2

+

c

h

u

(

s1

)

−

h

(

K

)

m

−

p

u

(

s1

)

h

u

(

s1

)

.

However, the above inequality cannot hold when m is large enough. Therefore, the trajectory y

(

s

;

y0

)

with s

>

0 cannot lie below the curve v

=

m

(

h

(

u

)

−

h

(

K

))

whenever u_∗

<

u

(

s

;

y0

) <

K .

Similarly, suppose that there is an s2

>

0 such that v

=

n

(

h

(

u

)

−

h

(

K

))

and

(

v

−

n

(

h

(

u

)

−

h

(

K

)))

0

(10)

0

v

(

s2

)

−

nh

u

(

s2

)

v

(

s2

)

=

cv

(

s2

)

−

h

u

(

s2

)

p

u

(

s2

)

−

nh

u

(

s2

)

v

(

s2

)

= −

h

u

(

s2

)

h

u

(

s2

)

−

h

(

K

)

n2

+

c

h

u

(

s2

)

−

h

(

K

)

n

−

h

u

(

s2

)

p

u

(

s2

)

.

The above inequality also cannot hold when n is small enough. Therefore, y

(

s

;

y0

)

with s

>

0 cannot

lie above the curve v

=

n

(

h

(

u

)

−

h

(

K

))

whenever u_∗

<

u

(

s

;

y0

) <

K .

Since y

(

s

;

y0

)

is bounded by the curves v

=

m

(

h

(

u

)

−

h

(

K

))

and v

=

n

(

h

(

u

)

−

h

(

K

))

, it follows that v

(

s

;

y0

) <

0 and u

(

s

;

y0

)

decreases until u

(

s

;

y0

) <

u∗. The proof is complete.

2

Since the invariant manifold

Ω

1has w

=

0 and z

=

0, we immediately have the following lemma.

Lemma 2.6. Any trajectory y

(

s

;

y0

)

with y0

∈ Ω

1, u0

>

K and v0

>

0 will stay in the region

{

u

>

K

,

v

>

0

}

for s

>

0.

Proof. Let y0

∈ Ω

1 be near

(

K

,

0

,

0

,

0

)

, then w

(

s

;

y0

)

=

0 for all s. Since u0

>

K and v0

>

0, we have v

(

s

;

y0

) >

0 for s

>

0. Hence the assertion follows.

2

Lemma 2.7. Any trajectory y

(

s

;

y0

)

with 0

<

u0

<

K , w0

>

0, and z0

>

c2w0 will stay in the region

{

w

>

0

,

z

>

c₂w

}

whenever 0

<

u

(

s

;

y0

) <

K .

Proof. Assume the assertion of this lemma is false. Let s1

>

0 be the ﬁrst time that y

(

s

;

y0

)

leaves

the region

{

w

>

0

,

z

>

c₂w

}

with 0

<

u

(

s1

,

y0

) <

K . Then for s

∈ [

0

,

s1

)

, we have

w

(

s

)

=

z

(

s

) >

c

2w

(

s

)

with w

(

0

) >

0

,

which implies w

(

s1

) >

0. Since

z

(

s1

)

=

c w

(

s1

)/

2 and z

(

s1

)

− (

c

/

2

)

w

(

s1

)

0

,

we have 0

cz

(

s1

)

+

w

(

s1

)

q

u

(

s1

)

−

c 2z

(

s1

)

c2 4 w

(

s1

)

+

w

(

s1

)

q

(

K

)

c2 4

+

₍

₀

₎

_q

₍

_K

₎

_w

₍

_s 1

).

This contradicts the assumption c

>

c∗. The proof is complete.

2

On

Ω

2, let’s parameterize a small circle centered at

(

K

,

0

,

0

,

0

)

by

G

(θ )

=

F2

ε

cos

(θ

+ ψ

0

),

ε

sin

(θ

+ ψ

0

)

,

(2.9)

where

θ

∈ [

0

,

2

π

]

and the constant phase

ψ

0is chosen such that G

(

0

)

lies in

Ω

1with u

<

K . Set

A

:=

θ

∈ [

0

,

2

π

)

:

∃

s0

>

0 satisfying u

s0

;

G

(θ )

=

u_∗and v

s

;

G

(θ )

<

0 on s

∈ (

0

,

s0

]

.

By Lemma 2.5, A is nonempty since

θ

=

0

∈

A. Denote

θ

1

:=

sup

(11)

Remark 2.8.

(i)

ψ

0is close to zero provided

ε

≈

0.

(ii) According to Lemma 2.5, there exists an s0

>

0 such that u

(

s0

;

G

(

0

)) <

u∗ and v

(

s

;

G

(

0

)) <

0 for

s

∈ [

0

,

s0

].

The continuous dependence of a solution on initial condition implies that

θ

1

>

0.

(iii) Since v

(

0

;

G

(θ ))

0 for

θ

∈

A, we have A

⊂ [

0

,

3

π

/

4

− ψ

0

)

. If

θ

∈ [

0

,

3

π

/

4

− ψ

0

)

, then the

components u and w of G

(θ )

satisfy 0

<

u

<

K and w

>

0. Thus, we have w

(

0

;

y1

) >

0.

Lemma 2.9. Let

ε

>

0 be small enough. If

θ

∈ [

0

,

3

π

/

4

− ψ

0

)

, then the trajectory y

(

s

;

G

(θ ))

with s

0 will stay in the region

{

w

>

0

,

z

>

c w

/

2

}

whenever 0

<

u

(

s

;

G

(θ )) <

K .

Proof. Let y0

=

G

(θ )

∈ Ω

2. From (2.9), the w and z coordinates of y0 satisfy w

>

0 and z

≈ λ

2w

>

c w

/

2. Then the assertion follows by Lemma 2.7.

2

Lemma 2.10. Suppose y0

=

G

(θ )

for some

θ

∈ (

0

, θ

1

)

. Then y

(

s

;

y0

)

will leave W and enter the region R or P .

Proof. Fix a

θ

∈ (

0

, θ

1

)

, then there exists s0such that

u

s0

;

G

(θ )

=

u_∗ and v

s

;

G

(θ )

<

0 for s

∈ (

0

,

s0

].

If

(

g

(

w

)

−

p

(

u

))

s=s0

<

0, we have dv

(

s0

)

=

cv

+

h

(

u

)

g

(

w

)

−

p

(

u

)

_s₌_s 0

<

0

,

which implies v

(

s+₀

) <

0 and u

(

s+₀

) <

u_∗. That is, the trajectory enters region R.

If

(

g

(

w

)

−

p

(

u

))

s=s0

0, then w

(

s0

)

w∗ by u

(

s0

)

=

u∗. Since v

(

s0

) <

0, we have u

(

s+0

) <

u∗.

By Lemma 2.9, we have w

(

s0

) >

0 and z

(

s0

) >

c₂w

(

s0

) >

0. Thus w

(

s+₀

) >

w∗. That is, the trajectory enters region P . The proof is complete.

2

The next lemma shows that there is a “last” trajectory on

Ω

2 such that u

(

s

)

decreases to the value u

=

u_∗.

Lemma 2.11. There exists an s0

>

0 such that u

(

s0

;

y1

)

=

u∗and v

(

s0

;

y1

)

=

0, see Fig. 2. Moreover, we have

g

w

(

s0

;

y1

)

−

p

u

(

s0

;

y1

)

>

0 and w

(

s0

;

y1

) >

w∗

.

Proof. Recall that u_∗

<

u

(

0

;

y1

) <

K , v

(

0

;

y1

)

0 and w

(

0

;

y1

) >

0. The proof consists of several steps

as follows.

(1) We claim that u

(

s

;

y1

)

u∗ or v

(

s

;

y1

)

0 for some s

>

0

.

Suppose the claim is false, i.e., u

(

s

;

y1

) >

u∗ and v

(

s

;

y1

) <

0 for all s

>

0

.

Then u

(

s

;

y1

)

decreases

monotonically to u

(

∞;

y1

)

u∗ and v

(

∞;

y1

)

=

0. By Lemma 2.9, we have

w

(

s

;

y1

)

=

z

(

s

;

y1

) >

c

/

2w

(

s

;

y1

),

which implies w

(

∞;

y1

)

= ∞

. Then it follows that

dv

(

s

;

y1

)

=

cv

(

s

;

y1

)

+

h

u

(

s

;

y1

)

g

w

(

s

;

y1

)

−

p

u

(

s

;

y1

)

→ ∞,

as s

→ ∞

. However, this fact contradicts v

(

∞;

y1

)

=

0. Hence the claim follows.

(2) Let s0 be the ﬁrst time that u

(

s

;

y1

)

=

u∗ or v

(

s

;

y1

)

=

0. We claim that v

(

s0

;

y1

)

=

0, v

(

s

;

y1

) <

0 for s

∈ (

0

,

s0

)

, and u

(

s0

;

y1

)

u∗

.

(12)

Fig. 2. Projection of the trajectory y(s;y1)in the u w-plane.

Suppose the claim is false, i.e.,

v

(

s

;

y1

) <

0 for s

∈ (

0

,

s0

]

and u

(

s0

;

y1

)

=

u∗

.

Then, by the Implicit Function Theorem, there exists a function s0

(θ )

with

θ

≈ θ

1such that

u

s0

(θ )

;

G

(θ )

=

u_∗

.

By the continuous dependence of the solution on

θ

, we have for

θ

≈ θ

1

v

s

;

G

(θ )

<

0 on s

∈

0

,

s0

(θ

1

)

+ δ

.

Also, by continuity of the function s0

(θ )

, we have s0

(θ )

∈ (

0

,

s0

(θ

1

)

+ δ]

for

θ

≈ θ

1. Therefore, there

are

θ

1satisfying v

s0

;

G

(θ )

<

0 on

0

,

s0

(θ )

and u

s0

(θ )

;

G

(θ )

=

u_∗

.

This fact contradicts the deﬁnition of

θ

1. Thus the claim follows.

(3) We claim that g

(

w

(

s0

;

y1

))

−

p

(

u

(

s0

;

y1

)) >

0 and v

(

s0

;

y1

) >

0.

Indeed, since v

(

s0

;

y1

)

=

0 and v

(

s

;

y1

) <

0 on s

∈ (

0

,

s0

)

, we have v

(

s0

;

y1

)

0 and

dv

(

s0

;

y1

)

=

h

u

(

s0

;

y1

)

g

w

(

s0

;

y1

)

−

p

u

(

s0

;

y1

)

0

.

Thus g

(

w

(

s0

;

y1

))

−

p

(

u

(

s0

;

y1

))

0. Suppose g

(

w

(

s0

;

y1

))

−

p

(

u

(

s0

;

y1

))

=

0, then

dv

(

s0

;

y1

)

h

u

(

s0

;

y1

)

g

w

(

s0

;

y1

)

z

(

s0

;

y1

),

which leads to dv

(

s0

;

y1

) >

h

u

(

s0

;

y1

)

g

w

(

s0

;

y1

)

c w

(

s0

;

y1

)/

2

>

0

by Lemma 2.9. This implies that v

(

s

;

y1

)

0 for s

≈

s0, which contradicts the deﬁnition of s0.

There-fore g

(

w

(

s0

;

y1

))

−

p

(

u

(

s0

;

y1

)) >

0 and v

(

s0

;

y1

) >

0.

(13)

Since v

(

s0

;

y1

)

=

0 and v

(

s0

;

y1

) >

0, by the Implicit Function Theorem, there exists a function s0

(θ )

for

θ

≈ θ

1 such that v

(

s0

(θ )

;

G

(θ ))

=

0. Suppose u

(

s0

;

y1

) >

u∗. Then, the continuous depen-dence of the solution on

θ

implies

v

s0

(θ )

;

G

(θ )

=

0

,

v

s0

(θ )

;

G

(θ )

>

0

;

v

s

;

G

(θ )

<

0 on s

∈

0

,

s0

(θ )

;

and u

(

s0

(θ ),

G

(θ )) >

u∗ for

θ

≈ θ

1. Thus

θ /

∈

A for

θ

≈ θ

1, a contradiction. Hence u

(

s0

;

y1

)

=

u∗. It follows from g

(

w

(

s0

;

y1

))

−

p

(

u

(

s0

;

y1

)) >

0 that w

(

s0

;

y1

) >

w∗. The proof is complete.

2

Lemma 2.12. There exists a

θ

2

∈ [θ

1

,

3

π

/

4

− ψ

0

)

such that the v coordinate of y2

:=

G

(θ

2

)

is equal to zero.

Proof. By (2.9), the v coordinate of G

(θ )

is given by

v

= −

ε

λ

2₁

+ λ

2₂sin

(θ

+ ψ

0

+ ψ

1

)

+

O

ε

2

,

where sin

ψ

1

= λ

1

/

λ

2₁

+ λ

2₂and

ψ

1

∈ (

π

/

4

,

π

/

2

)

. Obviously v

=

0 at

θ

2

:=

π

− ψ

0

− ψ

1

+

O

(

ε

)

∈ (

0

,

3

π

/

4

− ψ

0

).

Recall that the v coordinate of G

(θ

1

)

is non-positive. It follows that

θ

2

θ

1. The proof is complete.

2

On

Ω

3, we consider a small sphere centered at

(

K

,

0

,

0

,

0

)

with radius

ε

, which is parameterized

by

U

(θ, φ)

=

F3

ε

cos

(θ

+ ψ

0

)

sin

φ,

ε

sin

(θ

+ ψ

0

)

sin

φ,

ε

cos

φ

,

(2.10)

where

θ

∈ [

0

,

2

π

]

and

φ

∈ [

0

,

π

]

. The constant phase

ψ

0 is the one in (2.9). This sphere contains the

arc G

(θ )

=

U

(θ,

π

/

2

)

. According to Lemma 2.12 we know that the sphere intersects the hyperplane

v

=

0 at least one point U

(θ

2

,

π

/

2

)

. The next lemma shows that the intersection is a smooth closed

curve.

Lemma 2.13. The intersection of the sphere deﬁned by (2.10) and the hyperplane v

Existence of traveling wave solutions for diffusive predator-prey type systems

Journal of Differential Equations