Periodic orbits in 3-dimensional systems and application to a perturbed Volterra system

Available online at www.sciencedirect.com

ScienceDirect J. Differential Equations 260 (2016) 2750–2762 www.elsevier.com/locate/jde

Periodic orbits in 3-dimensional systems and application to a perturbed Volterra system ✩ Chengzhi Li a,b,∗ , Zhien Ma b , Yicang Zhou b a LMAM and School of Mathematical Sciences, Peking University, Beijing 100871, China b School of Mathematics and Statistics, Xi’an Jiaotong University, Xi’an 710049, China

Received 15 July 2015; revised 5 October 2015 Available online 23 October 2015

Abstract We give a theorem about the existence of periodic orbits for a 3-dimensional system, then use it to study a perturbation of a 3-dimensional Volterra system, and prove the existence and uniqueness of its periodic orbits on a 2-dimensional invariant manifold. © 2015 Elsevier Inc. All rights reserved. MSC: 34C07; 34C08; 37G15 Keywords: Invariant manifold; Periodic orbit; Volterra system; Bifurcation

1. A general theorem It is well-known that to study the number of periodic orbits for a perturbation of a 2-dimensional Hamiltonian system, the famous Poincaré–Pontryagin theorem is a crucial tool, see for example Section 2.1 in Part II of [3]. In this article we want to generalize this theory to 3-dimensional systems, then use it to study a perturbation of a 3-dimensional Volterra system. ✩

This work was partially supported by grants NSFC-11171267 (C. Li and Z. Ma) and NSFC-11271027 (C. Li).

* Corresponding author.

E-mail addresses: [email protected] (C. Li), [email protected] (Z. Ma), [email protected] (Y. Zhou). http://dx.doi.org/10.1016/j.jde.2015.10.018 0022-0396/© 2015 Elsevier Inc. All rights reserved.

C. Li et al. / J. Differential Equations 260 (2016) 2750–2762

2751

Consider a 3-dimensional system ⎧ dx ⎪ = p(x, y, z) + μf (x, y, z), ⎪ ⎨ dt dy Xμ : dt = q(x, y, z) + μ g(x, y, z), ⎪ ⎪ ⎩ dz = −z + z2 r(x, y, z) + μ e(x, y, z),

(1.1)

dt

where p, q, r, f, g, e ∈ C l (D), l ≥ 1, D = {(x, y, z) | (x, y) ∈ E, |z| ≤ β}, E is a compact region, containing (x, y) = (0, 0), β is a positive constant, and μ is a small parameter. We make the following basic assumptions for X0 , i.e. the system (1.1) with μ = 0: (A) X0 has a singularity at (x, y, z) = (0, 0, 0) and has a first integral (x, y, z) = h near z = 0, where  ∈ C m (D  ), m ≥ 2, D  = {(x, y, z) | (x, y) ∈ E, |z| ≤ δ < β} for some δ > 0. Let h = {(x, y, z) ∈ D  | (x, y, z) = h}, then {h | h ∈ (a, b)} is the family of integration surfaces of X0 . (B) Let γh = {(x, y) | (x, y, 0) = h}. Assume {γh | h ∈ (a, b)} is a family of closed curves, surrounding (x, y) = (0, 0), continuously and monotonically expanding from (x, y) = (0, 0) as h increases (or decreases). Let ∪h∈(a,b) γh = E  ⊂ E. (C) Along each curve γh the surface h transversally cuts the (x, y)-plane, i.e. 

∂ ∂x

2



∂ + ∂y

2 = 0. {(x,y)∈E  ,z=0}

We expand p, q, x , y in z at z = 0 as follows p(x, y, z) = p0 (x, y) + p1 (x, y)z + o(|z|), q(x, y, z) = q0 (x, y) + q1 (x, y)z + o(|z|), ∂ = ϕ0 (x, y) + ϕ1 (x, y)z + o(|z|), ∂x ∂ = ψ0 (x, y) + ψ1 (x, y)z + o(|z|). ∂y

(1.2)

From hypotheses (A) and (B) we know that system X0 |z=0 is integrable on its invariant (x, y)-plane. Let λ = λ(x, y) = 0 be the integrating factor for (x, y) ∈ E  , such that ϕ0 (x, y) = λ(x, y)q0 (x, y),

ψ0 (x, y) = −λ(x, y)p0 (x, y).

(1.3)

Finally, we let U (x, y) = λ(x, y)f (x, y, 0) + [λ(x, y)p1 (x, y) + ψ1 (x, y)]e(x, y, 0), V (x, y) = λ(x, y)g(x, y, 0) + [λ(x, y)q1 (x, y) − ϕ1 (x, y)]e(x, y, 0),

(1.4)

2752

C. Li et al. / J. Differential Equations 260 (2016) 2750–2762

where f , g, e are given in (1.1), and define the integral

I (h) =

U (x, y)dy − V (x, y)dx.

(1.5)

γh

Definition 1.1. If there exist an h∗ ∈ (a, b) and a σ > 0 such that system Xμ has a periodic orbit μ in R3 for 0 < |μ| < σ , and μ tends to γh∗ (in Hausdorff distance) as μ → 0, then we say that μ bifurcates from γh∗ . We say that a periodic orbit of Xμ bifurcates from the annulus E  = ∪h∈(a,b) γh , if there is an h ∈ (a, b) such that bifurcates from γh . We will show that such periodic orbit μ , bifurcating from a γh , stays on a 2-dimensional invariant manifold Mμ of system Xμ . Theorem 1.1. Suppose that the hypotheses (A), (B) and (C) are satisfied, and the integral I (h), defined in (1.5), is not identically zero for h ∈ (a, b), then the following statements hold: (1) If Xμ has a periodic orbit bifurcating from γh∗ , then I (h∗ ) = 0. (2) If there exists an h∗ ∈ (a, b) such that I (h∗ ) = 0 and I  (h∗ ) = 0, then Xμ has a unique periodic orbit μ bifurcating from γh∗ . Moreover, if I  (h∗ ) < 0 (resp. I  (h∗ ) > 0) then μ is attracting (resp. repelling) on Mμ if μ > 0, and μ is repelling (resp. attracting) on Mμ if μ < 0. (3) If there exists an h∗ ∈ (a, b) such that I (h∗ ) = I  (h∗ ) = · · · = I (k−1) (h∗ ) = 0, and I (k) (h∗ ) = 0, then Xμ has at most k periodic orbits bifurcating from the same γh∗ . (4) The total number of the periodic orbits of Xμ , bifurcating from the annulus ∪h∈(a,b) γh , is bounded by the maximum number of isolated zeros (taking into account their multiplicities) of the integral I (h) for h ∈ (a, b), if the later number is finite. Proof. By the assumptions (A) and (B) and the third equation of system Xμ , given in (1.1), we know that the linearized system of X0 at the singularity (0, 0, 0) has the center subspace {(x, y, z) | z = 0} and the stable subspace {(x, y, z) | x = 0, y = 0}. Besides, the (x, y)-plane is invariant for non-linear system X0 , hence it is the center manifold of X0 . We add dμ dt = 0 to system (1.1), then by using the center manifold theorem to this expanded 4-dimensional system and by the compactness of E, there exists a μ ¯ > 0 such that in region G = {(x, y, z, μ) | (x, y, z) ∈ D, |μ| ≤ μ} ¯ the expanded system has a center manifold W c , given by a function z = η(x, y, μ) for (x, y) ∈ E and |μ| ≤ μ, ¯ satisfying η(0, 0, 0) = 0, Dη(0, 0, 0) = (0, 0, 0); moreover, η → 0 as μ → 0, uniformly with respect to (x, y) ∈ E (see, for example, Chapter 1 of [2]). Note that for each μ0 , 0 < |μ0 | < μ, ¯ the intersection W c ∩ {μ = μ0 }, denoted by Mμ0 , is invariant with respect to the expanded system, and hence is also invariant with respect to system Xμ . Hence, for such a μ0 , Mμ0 is a 2-dimensional surface, given by z = η(x, y, μ0 ) for (x, y) ∈ E. By assumptions (A)–(C), if we let μ¯ small enough then for each μ, 0 < |μ| < μ, ¯ the surface h is also transversal to Mμ . Let γ˜h = h ∩ Mμ , then the set {γ˜h | h ∈ (a, b)} on Mμ is a family of ovals, surrounding the origin and continuously and monotonically depending on h. For any fixed h ∈ (a, b) and fixed μ, 0 < |μ| < μ, ¯ we take a point P0 on γ˜h , and take a curved segment  on Mμ , transversal to γ˜h at P0 . If we choose the function  itself to parameterize ,

C. Li et al. / J. Differential Equations 260 (2016) 2750–2762

2753

and denote by γ (h, μ) a piece of the orbit of Xμ on its invariant surface Mμ between the starting point P0 (measured by  = h) and the next intersection point of the orbit with , measured by  = P (h, μ). In this way we construct the Poincaré map for system Xμ along  on Mμ . Now let the difference d(h, μ) = P (h, μ) − h, then the orbit of Xμ , passing through the point P0 , is closed if and only if d(h, μ) = 0 for some (h, μ). It is clear that tP(h,μ) 

P (h,μ)

d(h, μ) =

d = th

h

 ∂ dx ∂ dy ∂ dz dt, + + ∂x dt ∂y dt ∂z dt γ (h,μ)

(1.6)

where th and tP (h,μ) are the moments of time t on orbit γ (h, μ), corresponding to the starting point P0 with  = h and the end point with  = P (h, μ) respectively. Note that h is the integration surface of X0 , hence for (x, y) ∈ E and |z| < δ we have ∂ ∂ ∂ p(x, y, z) + q(x, y, z) + (−z + z2 r(x, y, z)) ≡ 0. ∂x ∂y ∂z

(1.7)

Substituting (1.2) into above equality and using (1.3), we obtain ∂ = ϕ0 p1 + ϕ1 p0 + ψ0 q1 + ψ1 q0 . ∂z z=0

(1.8)

Since γ (h, μ) is a piece of orbit of Xμ , we can use equation (1.1) to compute (1.6). By using (1.7), we find tP(h,μ) 

d(h, μ) = μ th

 ∂ ∂ ∂ dt. f+ g+ e ∂x ∂y ∂z γ (h,μ)

(1.9)

We expand the integral in (1.9) with respect to μ at μ = 0, and use the fact that γ (h, 0) = γh stays on {z = 0}, because Mμ is given by z = η(x, y, μ), and η(x, y, μ) → 0 as μ → 0. By using (1.2), (1.3), (1.8), (1.4), (1.5) and the fact that along γh we have p0 dt = dx, q0 dt = dy, we finally obtain d(h, μ) = μ[I (h) + μ ξ(h, μ)],

(1.10)

where the function ξ(h, μ) is smooth (depending on the smoothness of f , g, e and ), and is uniformly bounded for (h, μ) in a compact region near (h, 0). Now we prove the statements (1)–(4) by using (1.10). (1) Suppose that a periodic orbit μ of Xμ bifurcates from γh∗ . That is, there exist a σ > 0 and hμ , such that hμ → h∗ as μ → 0, and d(hμ , μ) = μ[I (hμ ) + μ ξ(hμ , μ)] ≡ 0,

0 < |μ| < σ.

Dividing by μ on both sides, and taking limit as μ → 0, we obtain I (h∗ ) = 0.

2754

C. Li et al. / J. Differential Equations 260 (2016) 2750–2762

(2) Suppose that there exists an h∗ ∈ (a, b) such that I (h∗ ) = 0 and I  (h∗ ) = 0. Since we consider periodic orbit for small μ and μ = 0, instead of the function d(h, μ) we may study the ˜ μ) = d(h, μ)/μ. By (1.10) we have zeros of d(h, ˜ μ) = I (h) + μ ξ(h, μ), d(h, where ξ is smooth and uniformly bounded in a compact region near (h∗, 0). Since ˜ ∗ , 0) = I (h∗ ) = 0, d(h

d˜h (h∗ , 0) = I  (h∗ ) = 0,

by the Implicit Function Theorem, we find a σ > 0, an  > 0 and a unique function h = h(μ) ˜ defined in U ∗ = {(h, μ) : |h − h∗ | ≤ , |μ| ≤ σ }, such that h(0) = h∗ and d(h(μ), μ) ≡ 0 for (h, μ) ∈ U ∗ . Hence, the unique h(μ) gives a unique periodic orbit μ , bifurcating from γh∗ . Now we assume that there exists an h∗ ∈ (a, b) such that I (h∗ ) = 0 and I  (h∗ ) = 0, then we expand d(h, μ) for h near h∗ as follows d(h, μ) ≡ P (h, μ) − h = μ [I  (h∗ )(h − h∗ ) + o((h − h∗ )) + o(μ)]. Thus, for 0 < |μ| 1 and 0 < |h − h∗ | 1 we have (P (h, μ) − h)(h − h∗ ) < 0 if μI  (h∗ ) < 0 and (P (h, μ) − h)(h − h∗ ) > 0 if μI  (h∗ ) > 0, this gives the desired results, no matter h is monotonically increasing or decreasing with respect to increasing of h. (3) Assume that there exists an h∗ ∈ (a, b) such that I (h∗ ) = I  (h∗ ) = · · · = I (k−1) (h∗ ) = 0, and I (k) (h∗ ) = 0. We need to show that there exist a σ > 0 and an  > 0, such that for any (h, μ) ∈ U = {|h − h∗ | < , |μ| < σ }, the function d(h, μ) has at most k zeros in h. Suppose the contrary, then for any integer j there exist μj > 0 and j > 0, μj → 0 and j → 0 as j → ∞, such that for any μj the function d(h, μj )/μj has at least k + 1 zeros for |h − h∗ | < j . By using the Rolle Theorem k times we find an hj such that |hj − h∗ | < j and I (k) (hj ) + μj

∂k ξ(hj , μj ) = 0, ∂kh

which implies I (k) (h∗ ) = 0 by taking limit as j → ∞, leading to a contradiction. Of course, here we need to assume that the functions in (1.1) have enough smoothness, such that I (k) (h) is k continuous and ∂∂k h ξ(h, μ) is uniformly bounded for (h, μ) near (h∗ , 0). (4) This statement is a consequence of the first three statements. 2 Remark 1.1. If the system depends on some parameter(s), then the zero h∗ of I (h, μ) may also depend on the parameter(s), and h∗ may tend to a or b, the end-points of the interval (a, b) as the parameter approaches to some value, but a or b may correspond to some special orbits of X0 |z=0 , for example singular point, homoclinic orbit or heteroclinic orbits, or some orbits tending to infinity. This gives a difficulty to find a bound σ for μ uniformly with respect to the parameter(s). For more detailed explanation see, for example, Remark 2.5 in Part II (page 115) of [3], and the references in this book.

C. Li et al. / J. Differential Equations 260 (2016) 2750–2762

2755

2. A perturbed 3-dimensional Volterra system The Volterra system is an important model in biomathematics, see for example [7,8]. As an application of Theorem 1.1, we will prove the existence and uniqueness of periodic orbits for the following 3-dimensional system, which is a perturbation of a Volterra system. ⎧ du = u(1 − u − 3v + w) + μ (β + u2 ), ⎪ ⎨ dt dv 2 dt = v(1 + u − v − 3w) + μ (β + v ), ⎪ ⎩ dw 2 dt = w(1 − 3u + v − w) + μ (β + w ),

(2.1)

where β and μ are parameters and |μ| 1. When μ = 0, the Volterra system has 4 singularities at (0, 0, 0), (1, 0, 0), (0, 1, 0) and ( 13 , 13 , 13 ); and this system has a family of integration surfaces (cones) (u + v + w)3 = C, uvw

C > 27,

and an invariant plane u + v + w = 1. To transform system (2.1) to the form (1.1), we make the non-degenerate linear change u=

x +y 1 + , 2 3

v=

x −y 1 + , 2 3

w=

1 − x − z, 3

(2.2)

then system (2.1) is transformed to ⎧ dx ⎪ = 2y( 13 − x) + (x − 2y + 23 )z + μ [2β + 12 (x 2 + y 2 ) + 23 x + 29 ], ⎪ ⎨ dt dy 4 2 2 2 Xμ : dt = −2x − 3x + y + (−2x + y − 3 )z + μ [y(x + 3 )], ⎪ ⎪ ⎩ dz = −z + z2 + μ [−3β − 1 (x 2 + y 2 ) − (x + z)2 + 1 (2z − 1)]. dt

2

(2.3)

3

When μ = 0 system X0 has a first integral (x, y, z) =

(x + y + 23 )(x − y + 23 )( 13 − x − z) 4 = = α, 3 C (1 − z)

Since (x, y, 0) = (x − 13 )y 2 − x 2 (x + 1) +

4 27 ,

0<α<

4 . 27

we can write (x, y, 0) = α as

  1 2 γh : H (x, y) = x − y − x 2 (x + 1) = h, 3

−

4 < h < 0. 27

(2.4)

4 The level curves of {γh | h ∈ (− 27 , 0)} are shown in Fig. 1, which can be seen as the phase portraits of X0 |z=0 . γh shrinks to (x, y) = (0, 0) as h → 0− , and γh expands to the heteroclinic 4 + . loop, related to the 3 saddles and formed by the lines {y = ±(x + 23 )} and {x = 13 }, as h → − 27

2756

C. Li et al. / J. Differential Equations 260 (2016) 2750–2762

Fig. 1. The level curves of {γh |h ∈ (−4/27, 0)}. 4 It is easy to check that for (x, y) ∈ γh , h ∈ (− 27 , 0) we have



∂ ∂x



2

∂ + ∂y

2 {z=0}

 2 1 = (y 2 − 3x 2 − 2x)2 + 2y = 0. −x 3

As such, the assumptions (A)–(C) for system (2.3) are satisfied. Moreover, X0 |z=0 is a Hamiltonian system, i.e. λ(x, y) = 1 in (1.3). For system (2.3) the functions in (1.4) can be computed as 3 1 1 1 U (x, y) = −(9y + )x 3 + (3y − )x 2 − (3y 3 + y 2 + 2y + 18βy + 3β − )x 2 2 2 3 1 + y(6y 2 + y + 36β + 4), 6 3 1 3 4 1 3 V (x, y) = y − y − (3x − 9β − 1)y 2 − ( x 2 − x + 3β − )y 2 2 2 3 1 − x(3x + 2)(9x 2 + 18β + 2). 2

(2.5)

Since γh ⊂ H −1 (h), by (2.4) it is easy to see that for any polynomials P (x), Q(y) we have



P (x)dx =

γh

P (x)y dx =

γh

P (x)y dx = 0,

2

Q(y)dy = 0.

4

γh

γh

On the other hand,

y2 P (x)y dy =

γh

y1 P (x(y, h))y dy +

k

P (x(y, ˜ h))y k dy,

k

y1

y2

C. Li et al. / J. Differential Equations 260 (2016) 2750–2762

2757

Fig. 2. The two branches of γh .

where yi = yi (h); x(y, h) ≤ x(y, ˜ h) are two branches of γh , shown in Fig. 2. By integrating by ˜ i , h) for i = 1, 2, we find for an integer k ≥ 0 parts and noting x(yi , h) = x(y

1 P (x)y dy = − k+1



k

γh

P  (x)y k+1 dx.

(2.6)

γh

Substituting (2.5) into (1.5), and taking into account the above equalities, we find 2 I (h) = 3

[y 3 + (9x 2 + 9β − 1)y] dx.

(2.7)

γh

Multiplying both sides of (2.4) by dy then integrating it over γh and using (2.6), we obtain



y 3 dx =

γh

(9x 2 + 6x)y dx. γh

Substituting above equality into (2.7), we finally obtain 2 I (h) = 3

(18x 2 + 6x + 9β − 1)y dx.

(2.8)

γh

Note that A(h) = 4 , 0), h ∈ (− 27

 γh

ydx is the area of the region surrounded by γh , which is nonzero for

hence I (h) = 0 is equivalent to 1 − 9β =

6



γh (3x

2



γh

+ x)y dx

y dx

= F (h).

4 Since the right-hand side functions in (2.3) are all polynomials, F (h) ∈ C ∞ (− 27 , 0).

(2.9)

2758

C. Li et al. / J. Differential Equations 260 (2016) 2750–2762

Lemma 2.1. We have that lim F (h) = 0,

lim

h→0−

4 h→− 27

+

F (h) = 1.

Proof. By using the symmetry of γh with respect to x-axis and the mean-value theorem for integration, we have

x2 x2 2 2 (3x + x)y dx = 2 (3x + x)y(x, h) dx = 2 (3 θ + θ ) y(x, h) dx, 2

γh

x1

x1

where y = y(h, x) ≥ 0 is given by H (x, y) = h, x1 < θ < x2 with xi = xi (h) → 0 as h → 0− , see Fig. 2. Thus we immediately get F (h) → 0 as h → 0− . 4 + On the other hand, when h → − 27 , the oval γh expands to the heteroclinic loop, formed by  the lines {y = ±(x + 23 )} and {x = 13 }, hence γh y dx tends to the area of the triangle region, which is 1, and 1



3 (3x 2 + x)y dx → 2

γh

  2 1 (3x 2 + x) x + dx = , 3 6

as h → −

4+ . 27

− 23 +

4 This implies F (h) → 1 as h → − 27 .

2

4 Next, we prove that F  (h) < 0 for h ∈ (− 27 , 0). This fact is not obvious, we need to use a general result in [5], and to do some preparation. Consider a function

S(x, y) = ψ(x)y 2 + (x),

(2.10)

where ψ,  ∈ C 2 (α − , α + ) with α − < 0 < α + , ψ(x) > 0. We suppose that the continuous ovals {γs } are given by the level curves {(x, y) | S(x, y) = s, s ∈ (a, b)}, and we want to study the monotonicity property of the function  γ G(s) =  s γs

f2 (x)y dx f1 (x)y dx

,

(2.11)

where f1 and f2 are continuous in (α − , α + ). First, we make a hypothesis as follows: (H1) x   (x) > 0 (or < 0) for x ∈ (α − , α + ) \ {0}. By this condition the graph of y = (x) is U -shaped with a unique minimum (or ∩-shaped with a unique maximum) at x = 0, hence we may define a function x˜ = x(x) ˜ by (x) = (x), ˜

α − < x < 0 < x˜ < α + .

(2.12)

C. Li et al. / J. Differential Equations 260 (2016) 2750–2762

2759

Now we make the second hypothesis as follows: ˜ > 0 for x ∈ (α − , 0). (H2) f1 (x) f1 (x(x)) Lemma 2.2. (See Theorem 2 of [5].) Assume the hypotheses (H1) and (H2), let  √ ˜ ψ(x) ˜ − f2 (x) ˜  (x) ψ(x) f2 (x)  (x)  ν(x) = , √ f1 (x)  (x) ˜ ψ(x) ˜ − f1 (x) ˜  (x) ψ(x) x= ˜ x(x) ˜

(2.13)

then G (s) > 0 (resp. < 0) if ν  (x) < 0 (resp. > 0) for x ∈ (α − , 0). 4 Lemma 2.3. We have that F  (h) < 0 for h ∈ (− 27 , 0), where F (h) is give in (2.9).

Proof. To use Lemma 2.2 we need to change the sign of (2.4), such that to satisfy the condition ψ(x) > 0, that is change (2.4) to  γs : S(x, y) =

 1 − x y 2 + x 2 (x + 1) = s, 3

0
4 . 27

(2.14)

Comparing it with (2.10) we have that ψ(x) = 13 − x > 0, (x) = x 2 (x + 1), where x ∈ (− 23 , 13 ). We also have f1 (x) = 1, f2 (x) = 3x 2 + x. Hence the hypotheses (H1) and (H2) are satisfied. 4 4 It is clear that F  (h) < 0 for h ∈ (− 27 , 0) is equivalent to G (s) > 0 for s ∈ (0, 27 ). The function ν(x) can be written as √ √ ˜ x˜ + 1)(3x + 2) 1 − 3x x x(3x ˜ + 1)(3x˜ + 2) 1 − 3x˜ − xx(3 ν(x) = . (2.15) √ √ x(3 ˜ x˜ + 2) 1 − 3x˜ − x(3x + 2) 1 − 3x x= ˜ x(x) ˜ To simplify the calculations, we let τ= then − 23 < x < 0 < x˜ <

1 3

√

1 − 3x,

τ˜ =



1 − 3x, ˜

(2.16)

implies 0 < τ˜ < 1 < τ <

√

3.

(2.17)

2

dx 2   Since x = x(τ ) = 1−τ 3 , dτ = − 3 τ < 0, to prove ν (x) < 0 it is enough to prove ν¯ (τ ) > 0, where ν¯ (τ ) = ν(x(τ )). Substituting (2.16) into (2.15) we have

(τ 2 − 1)(τ˜ 2 − 1)[(τ τ˜ )2 − 2(τ + τ˜ )2 + 5τ τ˜ + 6] ν¯ (τ ) = , 2 2 2 2 2 2 3[(τ + τ˜ ) + τ τ˜ (τ + τ˜ − τ τ˜ + 4) − 4(τ + τ˜ ) + 3] τ˜ =τ˜ (τ ) where τ˜ = τ˜ (τ ) is given by substituting (2.16) into (2.12), i.e.  

1 − τ2 3



 =

 1 − τ˜ 2 , 3

(2.18)

2760

C. Li et al. / J. Differential Equations 260 (2016) 2750–2762

with (x) = x 2 (x + 1). Calculation gives (τ − τ˜ )(τ + τ˜ )(τ 2 + τ˜ 2 − τ τ˜ − 3)(τ 2 + τ˜ 2 + τ τ˜ − 3) = 0. From condition (2.17) we know that the first 3 factors above have fixed signs, hence τ˜ = τ˜ (τ ) is determined by τ 2 + τ˜ 2 + τ τ˜ − 3 = (τ + τ˜ )2 − τ τ˜ − 3 = 0.

(2.19)

It is convenient to make one more change ω = τ + τ˜ , ς = τ τ˜ , then from (2.19) we have ς = ω2 − 3. Using these new relations to (2.18), we find ν˜ (ω) = ν(τ ¯ (ω)) as follows ν˜ (ω) = =

(ς 2 − ω2 + 2ς + 1)(ς 2 − 2ω2 + 5ς + 6) 3[(ω2 − 2ς)2 + ς(ω2 − 3ς + 4) − 4ω2 + 3] ς=ω2 −3 1 2 (ω − 1)(4 − ω2 ). 3

(2.20)

  We first show dω dτ < 0, hence to prove ν¯ (τ ) > 0 it is enough to prove ν˜ (ω) < 0. In fact, dτ˜ 2τ +τ˜ from (2.19) we find dτ = − 2τ˜ +τ , hence

dω dτ˜ τ˜ − τ =1+ = < 0, dτ dτ 2τ˜ + τ √ because τ and τ˜ satisfy (2.17). Moreover, we obtain ω = τ + τ˜ ∈ ( 3, 2), since τ˜ (1) = 1 and √ τ˜ ( 3) = 0. Finally, by (2.20) we find ν˜  (ω) =

2 ω (5 − 2 ω2 ) < 0, 3

4 because ω2 ∈ (3, 4). Thus, by Lemma 2.2 G (s) > 0 for s ∈ (0, 27 ), and this implies F  (h) < 0 4 for h ∈ (− 27 , 0). 2

Theorem 2.1. For a given β ∈ (0, 19 ) there is a σ > 0, such that for any μ, 0 < |μ| < σ , system Xμ , given by (2.3), has a unique periodic orbit μ , tending to γh (in Hausdorff distance) as 4 , 0). Moreover, μ stays on a 2-dimensional invariant μ → 0, where h = F −1 (1 − 9β) ∈ (− 27 manifold of Xμ , and μ is attracting to nearby orbits in R3 if μ > 0. Proof. From (2.8) and (2.9) we have I (h) = where A(h) = have

(2.21)

γh y dx is the area of the region surrounded by γh . From (2.4) we see along γh we 1 = , by using the first equation of X0 |z=0 , given in (2.3), we obtain A (h) = T (h), 2y(x− 13 ) 4 is the period of γh , hence T (h) > 0 for h ∈ (− 27 , 0). Thus we have

∂y ∂h

which



2 A(h)[F (h) − (1 − 9β)], 3

C. Li et al. / J. Differential Equations 260 (2016) 2750–2762

2761

Fig. 3. The periodic solution of system (2.1) for β = 0.1 and μ = 0.001.

2 I  (h) = {T (h)[F (h) − (1 − 9β)] + A(h)F  (h)}. 3

(2.22)

4 By Lemmas 2.1 and 2.3 for any given β ∈ (0, 19 ) there is a unique h∗ = F −1 (1 − 9β) ∈ (− 27 , 0), −1 ∗ where F is the inverse function of F . Hence, by (2.21) and (2.22), we obtain I (h ) = 0 and I  (h∗ ) < 0. By the statement (2) of Theorem 1.1, there is a σ > 0, such that system Xμ has a unique periodic orbit μ on the 2-dimensional invariant surface Mμ for 0 < |μ| < σ . Moreover, μ is bifurcated from γh , and it is attracting to nearby orbits on Mμ for μ > 0. Besides, the third equation of (2.3) shows that the system is also attracting in z-direction outside Mμ , hence μ is attracting to nearby orbits in R3 if μ > 0. 2

Going back to the original system (2.1), we immediately obtain the following results. Theorem 2.2. For a given β ∈ (0, 19 ), there is a σ > 0 such that for any μ, 0 < |μ| < σ , system (2.1) has a unique periodic orbit μ on a 2-dimensional invariant surface, and μ is attracting to nearby orbits in R3 if μ > 0. Remark 2.1. As we mentioned in Remark 1.1, to find a σ > 0 uniformly with respect to 4 β ∈ (0, 19 ), we need to use compactness of the interval [a, b]. Now a = − 27 and b = 0. Since +

4 and A(h) → 0 as h → 0− , both give some trouble to finish the proof, T (h) → +∞ as h → − 27 see (2.21) and (2.22). Hence, a uniform result can be stated as follows: for a given ε, 0 < ε 1, there is a σε > 0, such that system Xμ has a unique periodic orbit for all β ∈ [ε, 19 − ε], where 0 < |μ| < σε . Since the singular point at (0, 0) of X0 |z=0 is a non-degenerate center, by known result we can expand β ∈ [ε, 19 − ε] to β ∈ [ε, 19 ], see for example Remark of Theorem 2 of [4], Chapter 4 of [6], or Lemma 20 of [1].

Remark 2.2. The numerical simulations show that system (2.1) has a periodic orbit μ for small as well as large μ when β ∈ (0, 19 ) (see Fig. 3 and Fig. 4).

2762

C. Li et al. / J. Differential Equations 260 (2016) 2750–2762

Fig. 4. The periodic solution of system (2.1) for β = 0.1 and μ = 0.5.

Acknowledgments The authors would like to thank the reviewer of this paper for the comments that help them to improve the earlier version of the manuscript. References [1] M. Caubergh, F. Dumortier, Hopf–Takens bifurcations and centres, J. Differential Equations 202 (2004) 1–31. [2] S.-N. Chow, C. Li, D. Wang, Normal Forms and Bifurcations of Planar Vector Fields, Cambridge University Press, 1994. [3] C. Christopher, C. Li, Limit Cycles of Differential Equations, Birkhäuser Verlag, Basel–Boston–Berlin, 2007. [4] C. Li, Abelian integrals and limit cycles, Qual. Theory Dyn. Syst. 11 (2012) 111–128. [5] C. Li, Z. Zhang, A criterion for determining the monotonicity of ratio of two Abelian integrals, J. Differential Equations 124 (1996) 407–424. [6] W. Li, Normal Form Theory and Its Applications, Science Press, Beijing, 2000 (in Chinese). [7] Zhien Ma, Mathematical Modeling and Research of Population Ecology, Anhui Education Press, 1996 (in Chinese). [8] Yanni Xiao, Yicang Zhou, Sanyi Tang, The Principle of Biomathematics, Xian Jiaotong University Press, 2012 (in Chinese).