Bounded traveling waves of the Burgers–Huxley equation

Nonlinear Analysis 74 (2011) 1047–1060

Contents lists available at ScienceDirect

Nonlinear Analysis journal homepage: www.elsevier.com/locate/na

Bounded traveling waves of the Burgers–Huxley equation✩ Yuqian Zhou a,b,∗ , Qian Liu c , Weinian Zhang d a

School of Mathematics and Statistics, Yunnan University, Kunming, Yunnan 650091, PR China

b

School of Mathematics, Chengdu University of Information Technology, Chengdu, Sichuan 610225, PR China

c

School of Computer Science Technology, Southwest University for Nationalities, Chengdu, Sichuan 610041, PR China

d

School of Mathematics, Sichuan University, Chengdu, Sichuan 610064, PR China

article

info

Article history: Received 12 March 2009 Accepted 5 September 2010 MSC: 35B32 37J20 35Q51 Keywords: Burgers–Huxley equation Solitary wave Kink wave Periodic wave Bifurcation

abstract In order to investigate bounded traveling waves of the Burgers–Huxley equation, bifurcations of codimension 1 and 2 are discussed for its traveling wave system. By reduction to center manifolds and normal forms we give conditions for the appearance of homoclinic solutions, heteroclinic solutions and periodic solutions, which correspondingly give conditions of existence for solitary waves, kink waves and periodic waves, three basic types of bounded traveling waves. Furthermore, their evolutions are discussed to investigate the existence of other types of bounded traveling waves, such as the oscillatory traveling waves corresponding to connections between an equilibrium and a periodic orbit and the oscillatory kink waves corresponding to connections of saddle–focus. © 2010 Elsevier Ltd. All rights reserved.

1. Introduction Being a reduced object of the Navier–Stokes equation to the case of space dimension 1, the Burgers equation [1] ut + auux − uxx = 0, where constant a controls the nonlinearity, is regarded as an elementary partial differential equation modeling the far field wave propagation in a nonlinear dissipative system. It can be used to describe many physical phenomena, such as sound waves in viscous media, magnetohydrodynamic waves in media with finite electrical conductivity, etc. In 1986, Satsuma [2] modified it further, by adding a nonlinear reaction term, to the form ut + auux − uxx + bu(u − 1)(u − r ) = 0,

(1.1)

called the Burgers–Huxley equation (see [3]), where the constant b describes a nonlinear source. Eq. (1.1) includes several known evolution equations, e.g., the FitzHugh–Nagumo equation [4] when a = 0 and the Newell–Whitehead equation [5] when a = 0 and r = −1, and models more extensive classes of wave propagation in biological and chemical systems. To understand complicated nonlinear wave phenomena, traveling wave solutions of a PDE play important roles. In general, three basic types of bounded traveling waves can occur for a PDE; these are periodic waves, kink waves and solitary waves. They are also known, respectively, as periodic wave trains, fronts, and pulses. Among them, the soliton, as a special type of solitary wave, is of great importance because it possesses the properties of both a wave and a particle. It is used widely

✩ Supported by the National Natural Science Foundation, PR China (Grant No. 11061039) and Scientific Research Foundation of CUIT (KYTZ200910).

∗

Corresponding author at: School of Mathematics and Statistics, Yunnan University, Kunming, Yunnan 650091, PR China. E-mail address: [email protected] (Y. Zhou).

0362-546X/$ – see front matter © 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.na.2010.09.012

1048

Y. Zhou et al. / Nonlinear Analysis 74 (2011) 1047–1060

in many fields and has been studied deeply by many people (see [6–10]). The Burgers–Huxley equation (1.1) exhibits many nonlinear behaviors, that are also reflected partially in its traveling waves. Therefore, its traveling waves are frequently studied. In 1986, assuming 0 < r < 1, Satsuma [2] reduced the equation to a bilinear form and obtained kink wave solutions. Later, Wang et al. [11] applied nonlinear transformations to obtain generalized kink wave solutions of Eq. (1.1) and its generalized form. More exact solutions, including various kink wave solutions and periodic solutions with respect to spatial variable x, were obtained in [12] by a reduction with the Cole–Hopf transformation. In addition, symmetries were investigated by Estévez and Gordoa [13,14] in the 1990’s. Since Eq. (1.1) is not an exactly solvable equation, because it does not pass the Painlevé test [15], in recent years more attention has been paid to numerical methods, such as the homotopy analysis method (HAM) [16] and the spectral collocation method [17], to solve it and its generalized form. Although there have been some profound results concerning the traveling waves of Eq. (1.1) which have contributed to our understanding of nonlinear physical phenomena and wave propagation, there are still some unresolved problems because Eq. (1.1) is not an exactly solvable equation. For instance, can any other type of bounded traveling wave solutions occur for Eq. (1.1)? If they exist, what dynamical behavior do they have? If these problems can be solved, the study of both exact traveling wave solutions and numerical solutions of Eq. (1.1) will be promoted further. In fact, these problems have involved bifurcation problems. Recall that heteroclinic orbits are trajectories which have two distinct equilibria as their α and ω-limit sets, and homoclinic orbits are trajectories whose α and ω-limit sets consist of the same equilibrium. So, the three basic types of bounded traveling waves mentioned above correspond to periodic, heteroclinic, and homoclinic orbits of the traveling wave system of a PDE respectively (see [18,19]). It is just this relationship that makes the bifurcation theory of dynamical system become an effective method to investigate bifurcations of traveling waves of PDEs. In recent decades, much effort has been devoted to bifurcations of traveling waves of PDEs. In 1997 Peterhof et al. [20] investigated the persistence and continuation of exponential dichotomies for solitary wave solutions of semilinear elliptic equations on infinite cylinders so that a Lyapunov–Schmidt reduction can be applied near solitary waves. Sánchez-Garduño and Maini [21] considered the existence of one-dimensional traveling wave solutions in nonlinear diffusion degenerate Nagumo equations and employed a dynamical systems approach to prove the bifurcation of a heteroclinic cycle. Later Katzengruber et al. [18] analyzed bifurcations of traveling waves, such as the Hopf bifurcation, multiple periodic orbit bifurcation, homoclinic bifurcation and heteroclinic bifurcation, in a standard model of electrical conduction in extrinsic semiconductors, which in scaled variables is actually a singular perturbation problem of a 3-dimensional ODE system. In 2002 Constantin and Strauss [22] constructed periodic traveling waves with vorticity for the classical inviscid water wave problem under the influence of gravity, described by the Euler equation with a free surface over a flat bottom, and used global bifurcation theory to construct a connected set of such solutions, containing flat waves as well as waves that approach flows with stagnation points. In 2003 Huang et al. [23] employed the Hopf bifurcation theorem to establish the existence of traveling front solutions and small amplitude traveling wave train solutions for a reaction–diffusion system based on a predator–prey model, which are equivalent to heteroclinic orbits and small amplitude periodic orbits in R4 respectively. Besides these, many results on the bifurcations of traveling waves for the Camassa–Holm equation, KdV equation and Zhiber–Shabat equation can be found in [24–27]. Motivated by the above results and enlightened by [18], for the more general case ab ̸= 0, r ∈ R, we try to give existence conditions and investigate the evolution behaviour of bounded traveling waves of Eq. (1.1) by studying bifurcations of its traveling wave solutions. Overcoming the difficulties of nonintegrability, we exhibit bifurcations of codimension 1 and 2 of the traveling wave system of Eq. (1.1), including the Hopf bifurcation, transcritical bifurcation, Bogdanov–Takens bifurcation, homoclinic orbits bifurcation, heteroclinic orbits bifurcation and Poincaré bifurcation, by computing center manifolds, normal forms and Abelian integrals. Using these bifurcation results and the relationships [18,19] between bounded traveling wave solutions and bounded orbits of the corresponding traveling wave system of Eq. (1.1), we obtain the existence conditions of three basic types of traveling wave solutions, including kink wave, solitary wave and periodic wave solutions. Furthermore, by discussing their evolution, some other bounded traveling wave solutions are found. In summary, for the Burgers–Huxley equation, there exist many types of bounded traveling wave solutions, including a type of monotone kink wave solution corresponding to connections of saddle–saddle, a type of oscillatory kink wave solution corresponding to connections of saddle–focus, solitary wave solutions corresponding homoclinic orbits, periodic wave solutions corresponding to periodic orbits and a type of oscillatory traveling wave solution corresponding to connections between an equilibrium and a periodic orbit. 2. Hopf bifurcation It is well known that a traveling wave solution is a special solution of the form u(x, t ) = u(x − ct ), where c ̸= 0 is the wave velocity. So, we can make the transformation ξ = x − ct to change Eq. (1.1) into the traveling wave system

−u′′ − cu′ + auu′ + bu(u − 1)(u − r ) = 0, which has the equivalent form



u′ = v, v ′ = −c v + auv + bu(u − 1)(u − r ),

where ′ denotes d/dξ . All equilibria of system (2.2) are of the form (u0 , 0), where u0 satisfies u0 (u0 − 1)(u0 − r ) = 0.

(2.2)

Y. Zhou et al. / Nonlinear Analysis 74 (2011) 1047–1060

1049

In the case that r ̸= 0 and r ̸= 1, system (2.2) has exactly three equilibria A(0, 0), B(1, 0) and C (r , 0). A is hyperbolic since bc ̸= 0. B is of center type when c = a and either b < 0, r < 1 or b > 0, r > 1; otherwise, B is hyperbolic. C is of center type when c = ar and either b < 0, r ∈ (−∞, 0) ∪ (1, +∞) or b > 0, r ∈ (0, 1); otherwise, C is hyperbolic. In the case that r = 0, system (2.2) has two equilibria A(0, 0) and B(1, 0). B is of center type when c = a, b < 0, which can be regarded as a special case in the last case; otherwise, B is hyperbolic. A is a degenerate equilibrium with the singular Jacobian matrix

[ J =

0 0

]

1 . −c

In the case that r = 1, system (2.2) has two equilibria A(0, 0) and B(1, 0). Equilibrium A is hyperbolic but B is degenerate with the singular Jacobian matrix

[ J =

0 0

1

]

a−c

.

In particular, when a = c, at the equilibrium B the Jacobian matrix becomes nilpotent and has double zero eigenvalues. In the non-degenerate cases we need to identify the focus from the center for equilibria of center type. Theorem 1. When a = c the equilibrium B of system (2.2) is of center type. In addition, it is either a center if r = 2 or a weak focus of multiplicity 1 if r ̸= 2. The weak focus is stable (or unstable) if bc (r − 2) < 0 (or >0) and a stable (or unstable) limit cycle arises from the Hopf bifurcation near B when parameters a, b, c , r vary from {(a, b, c , r )|a = c , bc (2 − r ) < 0} (or {(a, b, c , r )|a = c , bc (2 − r ) > 0}) to {(a, b, c , r )|a > c , bc (2 − r ) < 0, a − c small} (or {(a, b, c , r )|a < c , bc (2 − r ) > 0, c − a small}). In contrast, when c = ar equilibrium C of system (2.2) is of center type. In addition, it is either a center if r = 1/2 or a weak focus of multiplicity 1 if r ̸= 1/2. The weak focus is stable (or unstable) if ab(2r − 1) < 0 (or >0) and a stable (or unstable) limit cycle arises from the Hopf bifurcation near C when parameters a, b, c , r vary from {(a, b, c , r )|c = ar , ab(2r − 1) < 0} (or {(a, b, c , r )|c = ar , ab(2r − 1) > 0}) to {(a, b, c , r )|c < ar , ab(2r − 1) < 0, ar − c small} (or {(a, b, c , r )|c > ar , ab(2r − 1) > 0, c − ar small}). Proof. When a = c, expanding system (2.2) at B : (1, 0), we get



u′ = v, v ′ = (b − br )u + (2b − br )u2 + cuv + bu3 .

(2.3)

If r = 2 additionally, (2.3) takes the following form



u′ = v := U (u, v), v ′ = −bu + cuv + bu3 := V (u, v),

(2.4)

implying that B is a center because U (−u, v) = U (u, v), V (−u, v) = −V (u, v). √ When a = c , r ̸= 2, it is obvious that br − b > 0 when B is of center type. Let η1 = br − b. With the change of variables u = X , v = −η1 Y system (2.3) can be normalized as

 ′ X = −η1 Y , Y ′ = η 1 X −

1

η1

((2b − br )X 2 − η1 cXY + bX 3 ).

One can compute the first Lyapunov number LB (1) = bc (2 − r )/8η13 ̸= 0, implying that B is a weak focus of multiplicity 1 and at most one limit cycle arises from a Hopf bifurcation. Similarly, when c = ar, expanding system (2.2) at C : (r , 0), we get



u′ = v, v ′ = br (r − 1)u + auv + (2br − b)u2 + bu3 .

(2.5)

If r = 1/2 additionally, (2.5) is of the form

 ′ u = v,

b

v ′ = − u + auv + bu3 , 4

which is time-reversible as we saw in (2.4). This proves that C is a center. √ When c = ar , r ̸= 1/2, it is obvious that br − br 2 > 0 when C is of center type. Let η2 = br − br 2 . With the change of variables u = X , v = −η2 Y , system (2.5) can be normalized as

 ′ X = −η2 Y , Y ′ = η 2 X −

1

η2

(−aη2 XY + (2br − b)X 2 + bX 3 ).

1050

Y. Zhou et al. / Nonlinear Analysis 74 (2011) 1047–1060

Similarly, we can compute the first Lyapunov number LC (1) = ab(2r − 1)/8η23 ̸= 0, implying that C is also a weak focus of multiplicity 1. Other results are obtained directly from the Hopf bifurcation theory [28,29].  In fact, the existence conditions of limit cycles in Theorem 1 correspond to the local existence conditions of periodic wave solutions of system (1.1). 3. Degeneracy of codimension 1 In this section we identify the types of degenerate equilibria listed in Section 2 at each of which the Jacobian matrix has a single zero eigenvalue. Then we discuss their local bifurcations. Theorem 2. When r = 0 the degenerate equilibrium A : (0, 0) of system (2.2) is a saddle–node, at which a transcritical bifurcation occurs as r varies near 0. When r = 1 and a ̸= c, the degenerate equilibrium B : (1, 0) of system (2.2) is a saddle–node, at which a transcritical bifurcation occurs as r varies near 1. Proof. When r = 0, the expanded system of system (2.2) at the degenerate equilibrium A : (0, 0)



u′ = v, v ′ = −c v + auv − bu2 + bu3 ,

(3.6)

can be transformed by the change of variables X =u+

v c

v

,

Y =− , c

τ = −c ξ ,

as a polynomial system with the 2-jet

 b (ac + 2b) (ac + b) 2  X ′ = 2 X 2 + XY + Y := Φ1 (X , Y ), 2 2 c

c

c

(3.7)

 Y ′ = Y − b X 2 − (ac + 2b) XY − (ac + b) Y 2 := Ψ1 (X , Y ), 2 2 2 c

c

c

where ′ denotes d/dτ . From Ψ1 (X , Y ) = 0 we can obtain the implicit function Y = g (X ) :=

b c2

X2 +

(ac + 2b)b c4

X 3 + O(X 4 ).

(3.8)

From (3.7) and (3.8), X ′ = Φ1 (X , g (X )) =

b c2

X2 +

(ac + 2b)b c4

X 3 + O(X 4 ),

where the nonzero term of lowest degree is of even degree, implying by [29, Theorem 7.1, Chapter 2] that A : (0, 0) is a saddle–node. In more detail the parabolic sector of the saddle–node lies in the left-half (resp. right-half) plane in the (X , Y )coordinates for b < 0 (resp. b > 0). In order to discuss local bifurcations at the saddle–node A, we equivalently consider the suspended system of (2.2), that is,

 u′ = v, v ′ = −c v + auv + bu(u − 1)(u − r ), r ′ = 0. which can be normalized by the transformation X = u + vc , Y = −v , r = r as c

 X ′ = F (X , Y , r ), Y ′ = −cY − F (X , Y , r ), r ′ = 0 ,

(3.9)

where b F (X , Y , r ) := − X 2 − c



ac + b



c

b

b

3b

c

c

c

+ X3 + Y 3 +

Y2 −



X 2Y −

ac + 2b

 XY +

c b c

X 2r +

3b c

b c

XY 2 −

Xr + b c

b c

Yr

Y 2r −

2b c

XYr .

By [30] system (3.9) has a 2-dimensional center manifold Y = h(X , r ) near the origin, which has the expansion h(X , r ) = b X 2 − cb2 Xr + O(|X |3 + |r |3 ) because the second order approximate h2 (X , r ) := cb2 X 2 − cb2 Xr satisfies that (M h2 )(X , r ) = c2

Y. Zhou et al. / Nonlinear Analysis 74 (2011) 1047–1060

1051

O(|X |3 + |r |3 ), where

(Mh2 )(X , r ) := (∂/∂ X )h2 (X , r )F (X , h2 (X , r ), r ) + ch2 (X , r ) + F (X , h2 (X , r ), r ). Substituting the expansion of h(X , r ) in the system, we get X ′ = µA1 (r )X + µA2 (r )X 2 + O(|X |3 + |r |3 ),

(r ′ = 0),

where µA1 (r ) and µA2 (r ) are given in Appendix A. Rescaling with dτ = µA2 (r )dt, we get X ′ = µA (r )X + X 2 + O(|X |3 + |r |3 ),

(r ′ = 0),

where

µA (r ) :=

µA1 (r ) r (c 2 − br )c 2 = − . c 4 + (c 4 − c 3 a − 3bc 2 )r + (abc + b2 − 2c 2 b)r 2 + b2 r 3 µA2 (r )

Obviously, near 0 we see that µA (r ) < 0 as r > 0, = 0 as r = 0 and >0 as r < 0. It follows from [31, Section 3.4] that a transcritical bifurcation occurs at the saddle–node A. In contrast, when r = 1 and a ̸= c, system (2.2) can be expanded at the degenerate equilibrium B : (1, 0) as



u′ = v, v ′ = (a − c )v + auv + bu2 + bu3 .

(3.10)

With the homeomorphism u = −U , v = V and the time reversion ξ = −τ we can change (3.10) into



U′ = V , V ′ = (c − a)V + aUV − bU 2 + bU 3 ,

the same form as (3.6), where U ′ , V ′ present dU /dτ , dV /dτ respectively. Thus we assert that B is also a saddle–node and a transcritical bifurcation occurs as r crosses 1.  4. Degeneracy of codimension 2 In this section we discuss the bifurcation of codimension 2 at the degenerate equilibrium B : (1, 0) at which the Jacobian matrix is nilpotent. Theorem 3. When r = 1 and a = c, the degenerate equilibrium B : (1, 0) is a cusp. Moreover, in a neighborhood U of the point (c , 1) in the (a, r )-space of parameters there are four curves: SN + := {(a, r ) ∈ U |r = 1, ab(a − c ) > 0}, SN − := {(a, r ) ∈ U |r = 1, ab(a − c ) < 0}, H := {(a, r ) ∈ U |r = 2 − a/c , bc (a − c ) > 0}, HL := {(a, r ) ∈ U |r = (8a − 14c )/(a − 7c ), (a − c )(7c − a)b > 0 or r = (13a − 14c )/(6a − 7c ), (a − c )(7c − 6a)b < 0}, such that system (2.2) produces a saddle–node bifurcation at B as (a, r ) crossing SN + ∪ SN − , a Hopf bifurcation at B as (a, r ) crossing H, and a homoclinic bifurcation at B as (a, r ) crossing HL. Proof. When r = 1 and a = c, system (2.2) can be expanded at the equilibrium B : (1, 0) as the form



u′ = v, v ′ = bu2 (1 + u) + auv,

implying that B is a cusp shown by Theorem 7.3 in [29, Chapter 2]. In order to discuss bifurcations at the cusp B, for (a, r ) close to (c , 1), consider small parameters

ϵ1 := r − 1,

ϵ2 := a − c

(4.11)

and the 2-jet of the expansion of (2.2) at the cusp B, i.e.,



u′ = v, v ′ = −bϵ1 u + ϵ2 v + (b − bϵ1 )u2 + (c + ϵ2 )uv.

(4.12)

In order to reduce (4.12) to the standard form of versal unfolding of a cusp [32, Section 13.2], we first rescale (4.12) to normalize the terms of u2 and uv with the transversality, i.e., with the change of variables u=

b − bϵ1

(c + ϵ2 )

2

X1 ,

v=

(b − bϵ2 )2 Y1 , (c + ϵ2 )3

ξ=

c + ϵ2 b − bϵ 1

τ,

1052

Y. Zhou et al. / Nonlinear Analysis 74 (2011) 1047–1060

2

H

HL

II

1.5 I

III

r

SN-

SN+

1 III I 0.5 H HL II 0 0

0.5

1

1.5

2

a

Fig. 1. Bifurcation curves for the cusp B when b = 1, c = 1.

we transform (4.12) into the form

 ′ X1 = Y1 , Y1′ = −

b(c + ϵ2 )2 ϵ1

(b − bϵ1 )2

X1 +

ϵ2 (c + ϵ2 ) Y1 + X12 + X1 Y1 . b − bϵ1

(4.13)

Then, we eliminate the term of X1 from (4.13) with the transformation X1 = X2 +

b(c + ϵ2 )2 ϵ1 2(b − bϵ1 )2

,

Y1 = Y2

and obtain



X2′ = Y2 , Y2′ = µ1 + µ2 Y2 + X22 + X2 Y2 ,

(4.14)

where

µ1 = −

(c + ϵ2 )4 ϵ12 , 4b2 (ϵ1 − 1)4

µ2 =

(c + ϵ2 )(2ϵ2 − ϵ1 ϵ2 + ϵ1 c ) . 2b(ϵ1 − 1)2

(4.15)

As shown in [32, Section 13.2] and in [31, Section 7.3], the standard form (4.14) of versal unfolding of a cusp displays the Bogdanov–Takens bifurcation, including a saddle–node bifurcation as (µ1 , µ2 ) crossing SN + ∪ SN − , where SN + := {(µ1 , µ2 )|µ1 = 0, µ2 > 0} and SN − := {(µ1 , µ2 )|µ1 = 0, µ2 < 0}, a Hopf bifurcation as (µ1 , µ2 ) crossing H := 5

{(µ1 , µ2 )|µ1 = −µ22 + O(|µ2 | 2 ), µ2 > 0}, and a homoclinic bifurcation as (µ1 , µ2 ) crossing H L := {(µ1 , µ2 )|µ1 = − 49 µ2 + O(|µ2 |3 ), µ2 > 0}. 25 2 In what follows we give conditions of bifurcations for the original parameters a, r. The relations between those parameters a, b, c , r and µ1 and µ2 are given by (4.11) and (4.15). Such relations enable us to rewrite the bifurcation curves SN + , SN − , H and H L equivalently as those corresponding ones SN + , SN − , H and HL stated in the theorem.  Choosing b = 1, c = 1, we can plot those bifurcation curves SN + , SN − , H and HL as in Fig. 1. From this figure we clearly see that a homoclinic orbit occurs on homoclinic bifurcation curve HL. When parameter (a, r ) crosses the curve HL and enters the region II, the homoclinic orbit is broken, with a periodic orbit arising. The periodic orbit exists until (a, r ) reaches the Hopf bifurcation curve H. But, if (a, r ) crosses the curve HL and enters the region I, the homoclinic orbit is broken with no periodic orbit arising. In fact, our Theorem 3 and its proof clearly show that in a small enough neighborhood of (c , 1) a solitary wave arises in system (1.1) when parameter (a, r ) ∈ is broken with a periodic arising when (a, r ) crosses the curve  HL, which   14cwave   HL −13a −8a 2c −a and enters the one of eight regions: r < 14c , a > c > 0 , b < 0 , < r < , a > c > 0 , b > 0 , r > 7c −6a 7c −a c

      −8a −8a −13a , 0 < a < c , b > 0 , 2cc−a < r < 14c , 0 < a < c , b < 0 , r > 14c , c < a < 0, b > 0 , 14c 0 , r < 7c −6a , a < c < 0, b < 0 and c < r < 7c −a , c < a < 0, b < 0 . The periodic wave persists until (a, r ) reaches the Hopf bifurcation curve H. 14c −13a 7c −6a 2c −a a c

Y. Zhou et al. / Nonlinear Analysis 74 (2011) 1047–1060

1053

5. Non-local bifurcations Since there possibly exist homoclinic orbits and periodic orbits in a large range for system (2.2), we need to consider non-local bifurcations of system (2.2) with r ̸= 0 and r ̸= 1. Although it is difficult to study the non-lacal bifurcations, using the Melnikov method and some skill in dealing with the Abelian integrals, we discuss the homoclinic bifurcation, heteroclinic bifurcation and Poincaré bifurcation of system (2.2) so as to give conditions for the existence of homoclinic orbits, heteroclinic orbits and periodic orbits which correspond to the conditions for existence of solitary waves, kink waves and periodic waves respectively. Rescaling parameters a and c by a small ϵ > 0, i.e., a = ϵ a¯ , c = ϵ c¯ , we rewrite system (2.2) as



u′ = v, v ′ = bu(u − 1)(u − r ) + ϵ(auv − c v),

(5.16)

where we drop the bar from a and c for convenience. Consider the auxiliary system



u′ = v, v ′ = bu(u − 1)(u − r ),

(5.17)

which has the Hamiltonian energy function H (u, v) =

1 2

1

1

1

2

3

4

v 2 − bru2 + (b + br )u3 − bu4 .

(5.18)

Obviously, system (5.16) and system (5.17) have exactly the same three equilibria A : (0, 0), B : (1, 0) and C : (r , 0). We can calculate the energies hA , hB , hC at equilibria A, B and C : hA := H (0, 0) = 0,

hB := H (1, 0) =

1 12

1

b−

6

br ,

hC := H (r , 0) =

1 12

br 4 −

1 6

br 3 .

The dynamical behaviors of these equilibria of system (5.17) can be described in the three cases: Case (I): r > 1. In this case A, B and C are center, saddle and center respectively when b < 0 but are saddle, center and saddle respectively when b > 0. Case (II): 0 < r < 1. In this case A, C and B are center, saddle and center respectively when b < 0 but are saddle, center and saddle respectively when b > 0. Case (III): r < 0. In this case C , A and B are center, saddle and center respectively when b < 0 but are saddle, center and saddle respectively when b > 0. In case (I), system (5.17) has the following properties: (I-1): b < 0. In this subcase there are two homoclinic orbits Υ A and Υ C , which both connect the saddle B, and centers A and C are surrounded respectively by the families of periodic orbits

Γ1A

(h) : H (u, v) = h,

Γ1C (h) : H (u, v) = h,



1

0,

h∈

 h∈

12

1 12

b−

br 4 −

1 6 1 6

 br

br 3 ,

, 1 12

b−

1 6

 br

.

1 Moreover, Γ1A (h) (resp. Γ1C (h)) tends to A (resp. C ) as h → 0 resp. 12 br 4 − 16 br 3 and tends to Υ A (resp. Υ C ) as





− 16 br, as shown in Fig. 2(a). (I-2): b > 0 and 1 < r < 2. In this subcase there is a homoclinic orbit Υ1B connecting the saddle C and center B is surrounded h→

1 b 12

by the family of periodic orbits

Γ2B (h) : H (u, v) = h,

 h∈

1 12

b−

Moreover, Γ2B (h) tends to center B as h →

1 6

br ,

1 12

br 4 −

1 6

br 3



.

1 − 16 br and tends to Υ1B as h → 12 br 4 − 16 br 3 , as shown in Fig. 2(b). B (I-3): b > 0 and r = 2. In this subcase there are two heteroclinic orbits ΥU (the upper one) and ΥLB (the lower one), which 1 b 12

both connect the saddles A and C , and center B is surrounded by the family of periodic orbits

Γ3B (h) : H (u, v) = h, Moreover,

Γ3B

 h∈

 b − ,0 . 4

(h) tends to center B as h → − 4b and tends to ΥUB and ΥLB as h → 0, as shown in Fig. 2(c).

1054

Y. Zhou et al. / Nonlinear Analysis 74 (2011) 1047–1060

(a) b = −1, r = 2.

(b) b = 1, r = 1.5.

(c) b = 1, r = 2.

(d) b = 1, r = 2.5.

Fig. 2. Phase portraits of system (5.17) for r > 1.

(I-4): b > 0 and r > 2. In this subcase there is a homoclinic orbit Υ2B connecting the saddle A and center B is surrounded by the family of periodic orbits

Γ4B (h) : H (u, v) = h, Moreover,

Γ4B

 h∈

1 12

b−

(h) tends to center B as h →

1 6



br , 0 .

1 b 12

− 16 br and tends to Υ2B as h → 0, as shown in Fig. 2(d).

Let P1∗ (r ) :=

(r + 1)(f2 (r ) − f4 (r )) , √ 3 2(f1 (r ) − f3 (r )) √

P1∗∗ (r ) :=

(r + 1)(f2 (r ) + f4 (r )) , √ 3 2(f1 (r ) + f3 (r ))

√ (7r − 5)(r − 2)f5 (r ) + r (r − 1)(r + 1)(10r 2 − 19r + 10) P2 (r ) := , √ √ 2 2(2r − 1)f5 (r ) + 12 r (r − 1)(r 2 − r + 1) √ √ 2 (10r 2 + 11r + 10)f6 (r ) + r (r + 1)(10r 2 − 19r + 10) 6 ∗ P3 (r ) := , √ √ 2 2(r + 1)f6 (r ) + 12 r (r 2 − r + 1) ∗

−

2 6

where the fj s are functions of r defined in Appendix A. Theorem 4. In subcase (I-1), for sufficiently small ϵ , system (5.16) has a homoclinic orbit if either c /a = P1∗ (r ) or c /a = P1∗∗ (r ), and at least a limit cycle if either 0 < c /a < P1∗ (r ) or P1∗∗ (r ) < c /a < r. In subcase (I-2), for sufficiently small ϵ , system (5.16) has a homoclinic orbit if c /a = P2∗ (r ) and a unique limit cycle if 1 < c /a < P2∗ (r ). In subcase (I-3), for sufficiently small ϵ , system (5.16) has not only two heteroclinic orbits but also uncountably many periodic orbits inside the heteroclinic loops if a = c. In subcase (I-4), for sufficiently small ϵ , system (5.16) has a homoclinic orbit if c = aP3∗ (r ) and a unique limit cycle if P3∗ (r ) < c /a < 1. Proof. First, we prove those results in subcase (I-1). Obviously, homoclinic√orbits Υ A and Υ C intersect the u-axis at ∗ ∗ points (u∗1 , 0) and (u∗∗ −2 + 2r + 4r 2 )/3 and u∗∗ = u∗∗ 1 , 0) respectively, where u1 = u1 (r ) := (2r − 1 − 1 1 (r ) :=

Y. Zhou et al. / Nonlinear Analysis 74 (2011) 1047–1060

√ (2r − 1 + −2 + 2r + 4r 2 )/3. Consider the Melnikov function  c  (auv − c v)du = aI1 − cI0 = −aI0 M (h, a, c ) = − P ( h) , a

Γ

1055

(5.19)

where Ik := Γ uk v du, k = 0, 1, P (h) := I1 /I0 and Γ is either homoclinic loops Υ A ∪ {B}, Υ C ∪ {B} or periodic orbits Γ1A (h), Γ1C (h). First, we claim that



lim P (h) = 0,

lim P (h) = r ,

h→hA

lim P (h) =



h→hB

(5.20)

h→hC

P1∗ (r ), P1∗∗ (r ),

as u ∈ (u∗1 , 1), as u ∈ (1, u∗∗ 1 ).

(5.21)

In fact, let the periodic orbit Γ1A (h) (resp. Γ1C (h)) intersect the u-axis at two points (a1 (h), 0) and (b1 (h), 0) (resp. (a2 (h), 0) and (b2 (h), 0)), where a1 (h) < b1 (h), a2 (h) < b2 (h). Then I1 (h)

lim P (h) = lim

h→hA I0

h→hA

I1 (h)

lim P (h) = lim

I0 (h)

h→hC

h→hC

(h)

 b1 (h) a (h)

uv du

= lim  1b (h) 1 h→hA

v du

a1 (h)

 b2 (h) a (h)

uv du

= lim  2b (h) 2 h→hC

a2 (h)

v du

= lim u(h) = 0, h→hA

= lim u(h) = r , h→hC

implying (5.20). In order to prove (5.21), we note that the periodic orbit Γ1A (h) (resp. Γ1C (h)) approaches the homoclinic orbit Υ A (resp. Υ C ) as h → hB . Thus limh→hB P (h) = limh→hB I1 (h)/I0 (h) = I1∗ /I0∗ , where I0∗ := I0

 −

b(2r − 1)



 =

12

ΓhA ∪{B}

v du = 2

B

I1∗ := I1

 −

b(2r − 1)



12

 = ΓhA ∪{B}

1

∫

uv du = 2

u∗ 1

1

∫

B

V1 (u)du,

u∗ 1

uV1 (u)du,

and the function

 v = V1 (u) :=

−b 6

 (1 − u) −3u2 − 2u + 4ru − 1 + 2r ,

u ∈ (u∗1 , 1),

is the explicit expression of the curve Υ A . Direct computation implies that I1∗ /I0∗ = P1∗ (r ), a function defined just before Theorem 4. Similarly, considering the homoclinic orbit Υ C , we can calculate lim P (h) = lim

I1 ( h )

h→hB I0

h→hB

( h)

=

I1 (−b(2r − 1)/12) I0 (−b(2r − 1)/12)

= P1∗∗ (r ).

Thus (5.21) is proved. By (5.21) we see that M (hB , a, c ) has a unique zero c = aP1∗ (r ) (resp. c = aP1∗∗ (r )) as u ∈ (u∗1 , 1) (resp. u ∈ (1, u∗∗ 1 )). Therefore, system (5.16) has a homoclinic orbit, for sufficiently small ϵ , if either c /a = P1∗ (r ) or c /a = P1∗∗ (r ). In order to discuss limit cycles, we use P1 (h) and P2 (h) to distinguish the function P (h) on curves Γ1A (h) and Γ1C (h) respectively. We claim that P1 (h) is strictly increasing in h ∈ (0, hB ) and P2 (h) is strictly decreasing in h ∈ (hC , hB ). In fact, the energy function (5.18) is in the form of variable separation, the same as (17) in [33]. In order to prove the monotonicity of P1 (h), by Theorem 2 in [33] (see Appendix B) it suffices to verify the following conditions: (LZ1): G′ (u)(u − 0) > 0 (or < 0) for u ∈ (u∗1 , 1) \ {0}, (LZ2): f1 (u)f1 (˜u) > 0 for u ∈ (u∗1 , 0), (LZ3): ζ ′ (u) < 0 for u ∈ (u∗1 , 0), where G(u) := − 21 bru2 + 13 (b + br )u3 − 41 bu4 , u˜ = u˜ (u) is defined by G(u) = G(˜u), fi (u) = ui−1 , i = 1, 2, and

ζ (u) :=

f2 (u)G′ (˜u) − f2 (˜u)G′ (u) f1 (u)G′ (˜u) − f1 (˜u)G′ (u)

.

(5.22)

Condition (LZ1) can be verified, since G′ (u)u = −bu2 (u − 1)(u − r ) > 0 for u ∈ (u∗1 , 1) \ {0}. Condition (LZ2) holds naturally. Furthermore, from (5.22) one can calculate

ζ ′ (u) =

M (˜u, u) + M (u, u˜ )du˜ /du

[

u2

+ u˜ 2 + uu˜ − (1 + r )(u + u˜ ) + r ]2

,

(5.23)

1056

Y. Zhou et al. / Nonlinear Analysis 74 (2011) 1047–1060

where M (u, u˜ ) := u(u − 1)(r − u)(−u + 1 + r − 2u˜ ). Note that M (˜u, u) < 0, M (u, u˜ ) > 0 and du˜ du

G′ (u)

=

=−

G′ (˜u)

bu(u − 1)(u − r ) bu˜ (˜u − 1)(r − u˜ )

< 0,

since u∗1 < u < 0 < u˜ < 1 < r. It follows from (5.23) that ζ ′ (u) < 0 for u ∈ (u∗1 , 0). This concludes by Theorem 2 in [33] that P1′ (h) > 0 for 0 < h < hB . The claimed monotonicity of P2 (h) can be proved similarly. The above claimed monotonicity implies that if 0 < c /a < P1∗ (r ) (resp. P1∗∗ (r ) < c /a < r) then there exists a ∗ ∗∗ ′ ∗ unique h∗1 ∈ (0, hB ) (resp. h∗∗ 1 ∈ (hC , hB )) such that M (h1 , a, c ) = 0 (resp. M (h1 , a, c ) = 0). Moreover, M (h1 , a, c ) = aI0 (h∗1 )P1′ (h∗1 ) − aI0′ (h∗1 )(c /a − P1 (h∗1 )) = aI0 (h∗1 )P1′ (h∗1 ) ̸= 0 because I0 (h∗1 ), being the area of the region surrounded by Γ1A (h∗1 ), is positive. For the same reason, M ′ (h∗∗ 1 , a, c ) ̸= 0. By the Poincaré bifurcation theory [34] we obtain the results of Theorem 4 in case (I-1). Next, we prove those results in subcase (I-2). In this case the auxiliary system (5.17) has a homoclinic orbit Υ1B , which corresponds to the level curve

 v=±

b 6

 (r − u) 3u2 − 4u + 2ru − 2r + r 2 ,

u ∈ (u∗2 , r ).

√

One can easily compute the intersection point (u∗2 , 0) of Υ1B and the u-axis, where u∗2 = u∗2 (r ) := (2 − r + 4 + 2r − 2r 2 )/3. Similarly to subcase (I-1), we need to discuss zeros of the Melnikov function M (h, a, c ) defined in (5.19), where we choose Γ to be the homoclinic loop Υ1B ∪ {C } or the periodic orbits Γ2B (h). Using the same arguments as in subcase (I-1), we obtain limh→hB P (h) = 1 and limh→hC P (h) = P2∗ (r ), implying that M (hC , a, c ) has a unique zero c = aP2∗ (r ). Thus the result of homoclinic orbit in subcase (I-2) is proved. 1  1 In order to find limit cycles of system (5.16), we claim that P (h) is strictly increasing for h ∈ 12 b − 16 br , 12 br 4 − 16 br 3 . Our strategy remains to use Theorem 2 in [33]. In this subcase conditions (LZ1) and (LZ2) can be checked easily but it is too complicated to verify condition (LZ3) as done in subcase (I-1). In order to overcome the difficulty, we note that du˜ /du = G′ (u)/G′ (˜u) and G′ (˜u) = −bu˜ (˜u − 1)(˜u − r ) > 0 for u˜ ∈ (1, r ). It implies that the sign of ζ ′ (u) is the same as the function

Φ (u, u˜ ) := M (˜u, u)G′ (˜u) + M (u, u˜ )G′ (u). Applying the change of variables x = u + u˜ , y = uu˜ , we obtain

 (x, y) := Φ (u(x, y), u˜ (x, y)) Φ = −bx7 + (3rb + 3b)x6 + (−3b + 5by − 3r 2 b − 8rb)x5 + ((−14b − 14rb)y + 7r 2 b + r 3 b + b + 7rb)x4 + (−4by2 + (32rb + 13b + 13r 2 b)y − 5r 2 b − 2r 3 b − 2rb)x3 + ((11b + 11rb)y2 + (−24rb − 24r 2 b − 4r 3 b − 4b)y + r 2 b + r 3 b)x2 + (−3by3 + (−9b − 16rb − 9r 2 b)y2 + (6r 3 b + 6rb + 13r 2 b)y)x + (2b + 2rb)y3 + (2r 3 b + 6r 2 b + 6rb + 2b)y2 + (−2r 3 b − 2r 2 b)y.

(5.24)

Note that 1 + u2 < x < 1 + r since 0 < u2 < u < 1 < u˜ < r < 2. On the other hand, from the equality G(u) = G(˜u) we get ∗

∗

y = y(x) :=

x(3x2 − 4(r + 1)x + 6r ) 2(3x − 2r − 2)

.

It follows from (5.24) that

 (x, y(x)) = Φ

bx(x − 2)(x − 2r ) 8(3x − 2r − 2)2

Ψ (x, r ),

(5.25)

where

Ψ (x, r ) := 9x6 − 36x5 r + 66x4 r 2 − 56x3 r 3 + 16x2 r 4 − 36x5 + 114x4 r − 180x3 r 2 + 136x2 r 3 − 32xr 4 + 66x4 − 180x3 r + 240x2 r 2 − 144xr 3 + 24r 4 − 56x3 + 136x2 r − 144xr 2 + 48r 3 + 16x2 − 32xr + 24r 2 . Noting that G(u) − G(2 − u) = 23 b(u − 1)3 (−2 + r ) > 0 for u ∈ (u∗2 , 1), we have G(˜u) = G(u) > G(2 − u) which implies that u˜ > 2 − u, i.e., 2 < x < 1 + r since G′ (˜u) > 0 for u˜ ∈ (1, r ). Therefore, the sign of ζ ′ (u) is opposite to the sign of Ψ (x, r ) because the numerator bx(x − 2)(x − 2r ) of (5.25) is negative for 2 < x < 1 + r. In order to determine the sign of Ψ (x, r ), we need to estimate the range of x more precisely. From the fact G′ (u)+ G′ (˜u) = 1 b(x − 2)(x − 2r )x < 0 for all x ∈ (2, 1 + r ), 2 x′u = 1 +

du˜ du

=

G′ (˜u) + G′ (u) G′ (˜u)

< 0,

implying that x decreases as u increases from u∗2 to 1. Hence 2 < x < r + u∗2 .

Y. Zhou et al. / Nonlinear Analysis 74 (2011) 1047–1060

1057

Now we claim

Ψ (x, r ) > 0

(5.26)

in D := {(x, r ) ∈ R : 2 < x < r + u2 (r ), 1 < r < 2}. In fact, from the system of polynomial equations (∂/∂ x)Ψ (x, r ) = 0 and (∂/∂ r )Ψ (x, r ) = 0 we can solve with the Maple software five critical points of Φ (x, r ): ∗

2

S1 : (2, 1),

S2 : (2, 2),

S3 : (1, 1/2),

S4 : (0, 0),

S5 : (0, −1),

none of which is however included in the region D. We further investigate Ψ (x, r ) on the boundary ∂ D = C1 ∪ C2 , where C1 := {(x, r ) ∈ R2 : x = 2, 1 ≤ r ≤ 2} and C2 := {(x, r ) ∈ R2 : x = r + u∗2 (r ), 1 ≤ r ≤ 2}. On C1 we have Ψ (x, r ) = 24(r − 1)2 (r − 2)2 ≥ 0 but on C2 we have Ψ (x, r ) = (4/81)(1 + r )(2 − r )h(r ), where h(r ) = 4r 4 − 8r 3 − 12r 2 + 16r + (5 − 10r )(2(r + 1)(2 − r ))3/2

 + (24r 3 − 36r 2 + 36r − 12) 2(r + 1)(2 − r ) + 16. With the Maple software we find that h(r ) has a unique critical point r = 1.8363465 on the interval [1, 2], which is a maximal point. On the other hand, h(1) = h(2) = 0. It concludes that h(r ) ≥ 0 on [1, 2] and therefore Ψ (x, r ) ≥ 0 on C2 , completing the proof of the claimed (5.26). Thus we see that ζ ′ (u) < 0 for u∗2 < u < 1, i.e., condition (LZ3) is verified. By 1  1 Theorem 2 in [33] we have P ′ (h) > 0 for h ∈ 12 b − 16 br , 12 br 4 − 16 br 3 . Similarly to the discussion in subcase (I-1) we can 1 1 prove that M (h, a, c ) has a unique zero for h ∈ 12 b − 16 br , 12 br 4 − 16 br 3 when 1 < c /a < P2∗ (r ), showing the result of the limit cycle in subcase (I-2). In order to prove the results in subcase (I-3), we note that the auxiliary √ system (5.17) has two heteroclinic orbits connecting saddles A and C , which correspond to the level curves v = ± b/2u(2 − u) for 0 < u < 2. We need to discuss zeros of the Melnikov function M (h, a, c ) in (5.19), where we choose Γ to be the herteroclinic orbits ΥUB , ΥLB or a periodic orbit Γ3B . As in subcase (I-1), we first compute that limh→hB P (h) = 1 and limh→0 P (h) = 1. It means that either M (hA , a, c ) or M (hC , a, c ) has a unique zero c = a and the result of heteroclinic orbits are obtained. Furthermore, by our Theorem 1, equilibrium B is just a center in this subcase. Thus, the result of periodic orbits in subcase (I-3) is obtained. In subcase (I-4) the auxiliary system (5.17) has a homoclinic orbit Υ2B , which corresponds to the level curve



 v=±

b  2 u 3u − 4u − 4ru − 2r + 6r , 6



u ∈ (0, u∗4 ).

√

We can calculate the intersection point (u∗4 , 0) of Υ2B and the u-axis, where u∗4 = u∗4 (r ) := (2 + r − 4 − 10r + 4r 2 )/3. Then we discuss zeros of the Melnikov function M (h, a, c ) in (5.19), where we choose Γ to be either the homoclinic loop Υ2B ∪ {A} or a periodic orbit Γ4B . Using the same arguments as in subcase (I-2), we obtain that limh→hB P (h) = 1, limh→hA P (h) = P3∗ (r ) 1  and P ′ (h) < 0 for h ∈ 12 b − 16 br , 0 . It implies that M (hA , a, c ) has a unique zero c = aP3∗ (r ) and M (h, a, c ) has a unique zero for h ∈



1 b 12

 − 16 br , 0 if P3∗ (r ) < c /a < 1. Thus results in subcase (I-4) are proved. This completes the proof.



We similarly discuss the cases (II) and (III). First, the dynamical behaviors of the auxiliary system (5.17) are described as follows: (II-1) (or (III-1)): b < 0. In this subcase there are two homoclinic orbits connecting the same saddle C (or A) and centers A (or C ) and B are surrounded by a family of periodic orbits respectively. (II-2) (or (III-2)): b > 0 and 12 < r < 1 (or r < −1). In this subcase there is a homoclinic orbit connecting the saddle B and the center C (or A) is surrounded by a family of periodic orbits. (II-3) (or (III-3)): b > 0 and r = 12 (or r = −1). In this subcase there are two heteroclinic orbits connecting the saddles A (or C ) and B and the center C (or A) is surrounded by a family of periodic orbits. (II-4) (or (III-4)): b > 0 and 0 < r < 21 (or −1 < r < 0). In this subcase there is a homoclinic orbit connecting the saddle A (or C ) and the center C (or A) is surrounded by a family of periodic orbits. Let P4∗ (r ) :=

√

√ P5 (r ) = ∗

f9 (r ) + f10 (r )

36 2 − r (f7 (r ) + f8 (r ))

,

P4∗∗ (r ) :=

√

f9 (r ) − f10 (r )

36 2 − r (f7 (r ) − f8 (r ))

,

√

2(r + 1)[(−20r 3 + 48r 2 − 33r + 7)f11 (r ) + 6 2 − 2r (10r 2 − 19r + 10)]

√

12[(−2r 3 + 3r 2 + 3r − 2)f11 (r ) + 12 1 − r (−r 2 − r + 1)] f13 (r ) − f14 (r ) f13 (r ) + f14 (r ) , P6∗∗ (r ) := , P6∗ (r ) := 12(f12 (r ) − f3 (r )) 12(f12 (r ) + f3 (r )) √

P7 (r ) := ∗

2 12

√ √ (r + 1)[(5r − 7)(2r − 1)2 f11 (r ) − 6 2 1 − r (10r 2 − 19r + 10)] , √ √ 2(r − 2)(2r − 1)(r + 1)f11 (r ) − 12 1 − r (r 2 − r + 1)

,

1058

Y. Zhou et al. / Nonlinear Analysis 74 (2011) 1047–1060

√

√ − 122 (7r − 5)(r − 2)f15 (r ) + r (r − 1)(r + 1)(10r 2 − 19r + 10) P8 (r ) := , √ √ 2(2r − 1)f15 (r ) + 12 r (r − 1)(r 2 − r + 1) where the fj (r )s are functions of r defined in Appendix A. Then, using the same arguments as in case (I), we obtain the ∗

following results for system (5.16). Theorem 5. In subcase (II-1), for a sufficiently small ϵ , system (5.16) has a homoclinic orbit if either c = aP4∗ (r ) or c = aP4∗∗ (r ) and at least a limit cycle if either 0 < c /a < P4∗ (r ) or P4∗∗ (r ) < c /a < 1. In subcase (II-2), for sufficiently small ϵ , system (5.16) has a homoclinic orbit if c = aP5∗ (r ) and a unique limit cycle if r < c /a < P5∗ (r ). In subcase (II-3), for sufficiently small ϵ , system (5.16) has not only two heteroclinic orbits but also uncountably many periodic orbits inside the heteroclinic loops if c = a/2. In subcase (II-4), for sufficiently small ϵ , system (5.16) has a homoclinic orbit if c = aP3∗ (r ) and a unique limit cycle if P3∗ (r ) < c /a < r. Theorem 6. In subcase (III-1), for sufficiently small ϵ , system (5.16) has a homoclinic orbit if either c = aP6∗ (r ) or c = aP6∗∗ (r ) and at least a limit cycle if either r < c /a < P6∗ (r ) or P6∗∗ (r ) < c /a < 1. In subcase (III-2), for sufficiently small ϵ , system (5.16) has a homoclinic orbit if c = aP7∗ (r ) and a unique limit cycle if 0 < c /a < P7∗ (r ). In subcase (III-4), for sufficiently small ϵ , system (5.16) has a homoclinic orbit if c = aP8∗ (r ) and a unique limit cycle if P8∗ (r ) < c /a < 0. 6. Discussion By studying bifurcations of traveling wave system of Burgers–Huxley equation, we obtain the existence conditions for periodic orbits, homoclinic orbits and heteroclinic orbits which correspond to the existence conditions of periodic waves, solitary waves and kink waves of the Burgers–Huxley equation. From these results, the evolutions of these bounded traveling waves can be discussed further and some other types of bounded traveling wave can be found. Take the subcase (I-4) in Theorem 4 for example. When c /a = P3∗ (r ) and ϵ is small, Eq. (5.16) has a homoclinic orbit. Correspondingly, Eq. (1.1) has a solitary wave. When bifurcation parameter c /a varies from c /a = P3∗ (r ) to c /a < P3∗ (r ), the homoclinic orbit breaks, with the stable manifold of saddle A(0, 0) tending to the focus B(1, 0) when ξ → −∞, as shown in Fig. 3(a) and (b). It means, in this case, the solitary wave of Eq. (1.1) evolves into an oscillatory kink wave with amplitude increasing when ξ → +∞. When bifurcation parameter c /a varies from c /a = P3∗ (r ) to P3∗ (r ) < c /a < 1, the homoclinic orbit breaks, with a periodic orbit arising. Correspondingly, it means the solitary wave of Eq. (1.1) breaks, with a periodic wave arising. In addition, the unstable manifold of saddle A(0, 0) tends to the periodic orbit when ξ → +∞. Simultaneously, inside the periodic orbit, there exist a family of bounded orbits which tend to the periodic orbit (the focus B(1, 0)) when ξ → +∞ (ξ → −∞), as shown in Fig. 3(c) and (d). The two types of bounded orbits correspond to bounded oscillatory traveling waves of Eq. (1.1), which indicate periodicity when ξ → +∞. When the bifurcation parameter c /a varies from P3∗ (r ) < c /a < 1 to c /a ≥ 1, the periodic orbit disappears. The unstable manifold of saddle A(0, 0) tends to the equilibrium B(1, 0) when ξ → +∞, as shown in Fig. 3(e) and (f). Correspondingly, it means the periodic wave of Eq. (1.1) disappears and the bounded oscillatory traveling wave corresponding to connection between the saddle A(0, 0) and the periodic orbit evolves to an oscillatory kink wave. Other cases can be discussed similarly. Appendix A

(−ac − b) b2

   b2 b b (−ac − 2b) b b2 2 + 2 r + − − + r− c3 c5 c5 c3 c c3 c3 c  √  √ √ 2(r − 2) f1 (r ) = 2 2(r − 2)(2r − 1)(r + 1) arcsin √ + 12 r − 1(r 2 − r + 1), 2 2r + r − 1  √  √ 2(r − 2) 2 f2 (r ) = (5r − 7)(2r − 1) arcsin √ + 2r − 2(30r 2 − 57r + 30), 2 2r + r − 1 √ π f3 (r ) = 2π (r − 2)(2r − 1)(r + 1), f4 (r ) = (5r − 7)(2r − 1)2 , 2 √ √ √    4r − 2 + 3 2r (r − 1) 2(1 + r ) − 3 r f5 (r ) = (r + 1)(r − 2) ln , f6 (r ) = (2r − 1)(r − 2) ln √ , √ 2(r + 1)(2 − r ) (2r − 1)(r − 2)  √   √ 2(2r − 1) f7 (r ) = −2 2(r − 2)(2r − 1)(r + 1) arcsin √ + 12 r (1 − r )(r 2 − r + 1), 2 −r + r + 2 √ 3 2 f8 (r ) = π 2(−r + 3r + 3r − 2), 

µA1 =

br c 2 − br



,

µA2 = −

b3 r 3



+

Y. Zhou et al. / Nonlinear Analysis 74 (2011) 1047–1060

1059

1

0.5

0

-0.5

-1

-1.5

(a) Simulation by Maple for c /a < P3 (r ).

0

0.5

1

1.5

2

2.5

(b) Simulation by Matlab for c /a < P3∗ (r ).

∗

(c) Simulation by Maple for c /a ∈ (P3∗ (r ), 1).

0.8

0.4

0.6

0.3

0.4

0.2

0.2

0.1 0

0

-0.1

-0.2

-0.2 -0.4

-0.3

-0.6

-0.4

-0.8

-0.5

-1 0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

-0.6

1.8

(d) Simulation by Matlab for c /a ∈ (P3∗ (r ), 1).

(e) Simulation by Maple for c /a > 1.

0.7

0.8

0.9

1

1.1

1.2

1.3

1.4

1.5

(f) Simulation by Matlab for c /a > 1.

Fig. 3. Phase portraits of system (5.16) for ϵ = 0.5, b = 1, r = 2.5.

√

f9 (r ) = 6 2 − r |r + 1|

 √

 √

2(7r − 5)(r − 2) arcsin

2(2r − 1)



 

+ 6 r (1 − r )(10r − 19r + 10) , √ −r 2 + r + 2   √ √ 2(2r − 1)(r + 1) 2 f10 (r ) = 3 2π |r + 1| 2 − r (7r − 5)(r − 2) , f11 (r ) = ln , √ (4 − 2r + 3 2 − 2r )2  √  √ √ 2(r + 1) f12 (r ) = 2 2(r − 2)(2r − 1)(r + 1) arcsin √ + 12 −r (r 2 − r + 1), 2 2r − 5r + 2  √  √ √ 2(r + 1) 2 f13 (r ) = 2 2(r − 2)(2r − 1)(10r + 11r + 10) arcsin √ + −r (120r 3 − 108r 2 − 108r + 120), 2r 2 − 5r + 2   √ 2(r + 1)(2 − r ) 4 3 2 f14 (r ) = π 2(20r − 28r − 15r − 28r + 20), f15 (r ) = (r + 1)(r − 2) ln . √ (2 − 4r + 3 2r 2 − 2r )2 2

2

Appendix B Theorem 2 in [33]. Let the Hamiltonian function H (x, y) := φ(x)y2 + Φ (x) be defined on [α, A] × [β, B ], where φ, Φ ∈ C 1 and φ(x) has fixed sign, and Γh be the compact component of the curve H (x, y) = h, h1 < h < h2 . Suppose that (i) there is a ∈ (α, A) such that Φ ′ (x)(x − a) > 0 (or <0) for x ∈ (α, A) \ {a}, and (ii) f1 (x), f2 (x) are continuous and f1 (x)f1 (˜x) > 0 for x ∈ (α, a), where x˜ := x˜ (x) is defined by Φ (x) = Φ (˜x) for x ≤ a ≤ x˜ . Then the ratio of two Abelian integrals

w(h) :=

∫ Γh

f2 (x)ydx

∫ Γh

f1 (x)ydx

is decreasing (resp. increasing) for h1 < h < h2 if the function

ζ (x) :=

f2 (x)Φ ′ (˜x) − f2 (˜x)Φ ′ (x) f1 (x)Φ ′ (˜x) − f1 (˜x)Φ ′ (x)

is increasing (resp. decreasing) for α < x < a.

1060

Y. Zhou et al. / Nonlinear Analysis 74 (2011) 1047–1060

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34]

J.M. Burgers, A mathematical model illustrating the theory of turbulence, Adv. Appl. Mech. 1 (1948) 171–199. J. Satsuma, Topics in Soliton Theory and Exactly Solvable Nonlinear Equations, World Scientific, Singapore, 1987. L. Debnath, Nonlinear Partial Differential Equations, Birkhauser, Boston, 1997. R. Fitzhugh, Mathematical Models of Excitation and Propagation in Nerve, Biological Engineering, McGraw-Hill, New York, 1969. B.-Q. Lu, et al., Exact traveling wave solution of one class of nonlinear diffusion equation, Phys. Lett. A 175 (1993) 113–115. M.J. Ablowitz, P.A. Clarkson, Solitons, Nonlinear Evolution Equations and Inverse Scattering, Cambridge University Press, 1991. G. Eilenberger, Solitons, Mathematical Method for Physicists, Springer, Berlin, 1981. C.-H. Gu, et al., Soliton Theory and its Application, Springer, Berlin, 1995. B.-L. Guo, X.-F. Pang, Soliton, Science Press, Beijing, 1987. Y.-S. Li, Soliton and Integrable System, Shanghai Scientific and Technological Education Publishing House, Shanghai, 1999. X.-Y. Wang, Z.-S. Zhu, Y.-K. Lu, Solitary wave solutions of the generalised Burgers–Huxley equation, J. Phys. A 23 (1990) 271–274. O.Yu. Yefimova, N.A. Kudryashov, Exact solutions of the Burgers–Huxley equation, J. Appl. Math. Mech. 68 (2004) 413–420. P.G. Estévez, Non-classical symmetries and the singular manifold method: the Burgers and the Burgers–Huxley equations, J. Phys. A 27 (1994) 2113–2127. P.G. Estévez, P.R. Gordoa, Nonclassical symmetries and the singular manifold method: theory and six examples, Stud. Appl. Math. 95 (1995) 73–113. N.A. Kudryashov, Analytical Theory of Non-linear Differential Equations, Mosk Inzk-Firz Inst., Moscow, 2002. A. Molabahrami, F. Khani, The homotopy analysis method to solve the Burgers–Huxley equation, Nonlinear Anal. RWA 10 (2009) 589–600. M.T. Darvishi, et al., Spectral collocation method and Darvishis preconditionings to solve the generalized Burgers–Huxley equation, Commun. Nonlinear Sci. Numer. Simul. 13 (2008) 2091–2103. B. Katzengruber, M. Krupa, P. Szmolyan, Bifurcation of traveling waves in extrinsic semiconductors, Physica D 144 (2000) 1–19. J.-B. Li, H.-H. Dai, On the Study of Singular Nonlinear Travelling Wave Equation: Dynamical System Approach, Science Press, Beijing, 2007. D. Peterhof, B. Sandstede, A. Scheel, Exponential dichotomies for solitary-wave solutions of semilinear elliptic equations on infinite cylinders, J. Differential Equations 140 (1997) 266–308. F. Sánchez-Garduño, P.K. Maini, Traveling wave phenomena in non-linear diffusion degenerate Nagumo equations, J. Math. Biol. 35 (1997) 713–728. A. Constantin, W. Strauss, Exact periodic traveling water waves with vorticity, C. R. Acad. Sci. Paris, Ser. I 335 (2002) 797–800. J. Huang, et al., Existence of traveling wave solutions in a diffusive predator–prey model, J. Math. Biol. 46 (2003) 132–152. B. He, et al., Bifurcations of travelling wave solutions for a variant of Camassa–Holm equation, Nonlinear Anal. RWA 9 (2008) 222–232. J.-B. Li, Z.-R. Liu, Smooth and non-smooth traveling waves in a nonlinearly dispersive equation, Appl. Math. Model. 25 (2000) 41–56. Z.-R. Liu, C.-X. Yang, The application of bifurcation method to a higher-order KdV equation, J. Math. Anal. Appl. 275 (2002) 1–12. Y.-N. Tang, et al., Bifurcations of traveling wave solutions for Zhiber–Shaba equation, Nonlinear Anal. TMA 67 (2007) 648–656. N.G. Lloyd, Limit cycles of polynomial systems—some recent developments, in: T. Bedford, J. Swift (Eds.), New Directions in Dynamical Systems, Cambridge Univeristy Press, Cambridge, 1988. Z.-F. Zhang, et al., Qualitative Theory of Differential Equations, Amer. Math. Soc., Providence, 1992. J. Carr, Applications of Centre Manifold Theory, Springer, New York, 1981. J. Guckenheimer, P. Holmes, Nonlinear Oscillations, Dynamical Systems and Bifurcation of Vector Fields, Springer, New York, 1983. S.-N. Chow, J.K. Hale, Method of Bifurcation Theory, Springer, New York, 1982. C.-Z. Li, Z.-F. Zhang, A criterion for determining the monotonicity of the ratio of two Abelian integrals, J. Differential Equations 124 (1990) 407–424. S.-N. Chow, C.-Z. Li, D. Wang, Normal Forms and Bifurcation of Planar Vector Fields, Cambridge University Press, New York, 1994.