The number of limit cycles from a cubic center by the Melnikov function of any order

JID:YJDEQ AID:9972 /FLA

[m1+; v1.304; Prn:10/09/2019; 10:19] P.1 (1-32)

Available online at www.sciencedirect.com

ScienceDirect J. Differential Equations ••• (••••) •••–••• www.elsevier.com/locate/jde

The number of limit cycles from a cubic center by the Melnikov function of any order Peixing Yang, Jiang Yu ∗ School of Mathematical Sciences, Shanghai Jiaotong University, Shanghai 200240, P. R. China Received 27 January 2019; revised 19 June 2019; accepted 31 August 2019

Abstract In this paper, we consider the system x˙ = y(1 + x)2 − P (x, y), y˙ = −x(1 + x)2 + Q(x, y) where P (x, y) and Q(x, y) are arbitrary quadratic polynomials. We study the maximum number of limit cycles bifurcating from the periodic orbits by using the Melnikov function of any order. We prove that the upper bound for the number of limit cycles is 3 and reached. © 2019 Elsevier Inc. All rights reserved. MSC: 34C05; 34C07; 37G15 Keywords: Melnikov functions; Bifurcation; Limit cycles; Generators

1. Introduction and main results The development of the theory for planar ordinary differential equations has been stimulated by many applications. One of the most challenging problems is the second part of the classical Hilbert’s 16th problem which is to determine the number of limit cycles for the planar polynomial ordinary differential systems with degree n. Many methodologies have been developed, such as the averaging theory or the Melnikov function theory, which are used to study the number of limit cycles. * Corresponding author.

E-mail address: [email protected] (J. Yu). https://doi.org/10.1016/j.jde.2019.08.053 0022-0396/© 2019 Elsevier Inc. All rights reserved.

JID:YJDEQ AID:9972 /FLA

[m1+; v1.304; Prn:10/09/2019; 10:19] P.2 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

2

For the system 

x˙ = y + p(x, y), y˙ = −x + q(x, y),

where p(x, y) and q(x, y) are polynomials with degree n. The unperturbed system has a family of periodic orbits when  = 0. Giacomini et al. [11] obtained up to first order in  at most [ n−1 2 ] bifurcated limit cycles, where [·] denotes the integer part function. On the other hand, many researchers studied the perturbations of the integrable non-Hamiltonian system as follows: 

x˙ = yh(x, y) + p(x, y), y˙ = −xh(x, y) + q(x, y),

(1.1)

where h(x, y) is a polynomial in x and y with h(0, 0) = 0, while p(x, y) and q(x, y) are polynomials with degree n. We remark σ the maximum number of limit cycles which bifurcate from the periodic orbits of the unperturbed system (1.1) with  = 0. Llibre et al. [17,19] studied system (1.1) with h(x) = 1 + x and obtained at least n limit cycles bifurcating from the periodic orbits of the unperturbed system by using the averaging theory of first order, and obtained at most 2n − 1 limit cycles by using the averaging theory of second order. Authors in [12] considered the case that h(x, y) is a conic and h(0, 0) = 0, while p(x, y) and q(x, y) are arbitrary cubic polynomials. They found at least six limit cycles bifurcating from the periodic orbits of the unperturbed system by using the averaging theory of first order. Buic˘a and Llibre [4] researched system (1.1) with h(x, y) = (x + a)(y + b) where a, b ∈ R n−1 n−1 and ab = 0. They obtained that 3[ n−1 2 ] + 2 ≤ σ ≤ 3[ 2 ] + 4 when a = b and 2[ 2 ] + 1 ≤ n−1 σ ≤ 2[ 2 ] + 2 when a = b by using the Melnikov function of first order, which is also called Abelian integral. Authors in [18] used the averaging theory of second order to study system (1.1) with h(x, y) = (a1 x + a0 )(b1 y + b0 ) where a0 b0 = 0 and they obtained at most 17n + 15 limit cycles which bifurcate from the periodic orbits of the unperturbed system. System (1.1) with h(x) = 1 + Bx + Ax 2 where A = 0 was studied by Xiang and Han [22], and they obtained the maximum number of limit cycles which is at most 2n + 2 − (−1)n when B 2 − 4A = 0, and n with B 2 − 4A = 0 by using the Melnikov function of first order. Li et al. [16] gave an upper bound on the number of limit cycles for system (1.1) with h(x, y) = y 2 + ax 2 + bx + c and c = 0 by using the averaging theory of first order, which is 2n + 3. Coll et al. [5] used the averaging theory of first order to provide a boundary in the number n−1 of limit cycles, which is 4[ n−1 2 ] + 4 ≤ σ ≤ 5[ 2 ] + 14 for system (1.1) with h(x, y) = (x + a)(y + b)(x + c) with abc = 0. System (1.1) with h(x) = (1 − x)m has been studied by Xiang and Han [21]. They proposed that the maximum number of limit cycles is n − m + 1 by using the Melnikov function of first order, and when m = 1, the upper bound can be reached, namely, the system has precisely n limit cycles. Then Gasull et al. [8] also used the Melnikov function of first order to study the similar system (1.1) with h(y) = (1 − y)m and obtained that σ ≤ [ m+n 2 ] − 1 when n < m − 1 and σ ≤ n when n ≥ m − 1.

JID:YJDEQ AID:9972 /FLA

[m1+; v1.304; Prn:10/09/2019; 10:19] P.3 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

3

The works [1,2] studied system (1.1) with h(x, y) = (x 2 − a 2 )(y 2 − b2 ) where a, b ∈ R − {0}, l  aj x j where l = 2n + 2 or 2n + 3. and h(x) = 1 + x 4 with p(x, y) = 0 and q(x, y) = y 2m−1 Gasull et al. [7,9] considered system (1.1) with h(x, y) = K 

j =0 K 1

K 2

j =1

l=1

(x − aj )

(y − bl ) and h(x, y) =

((x − aj )2 + (y − bj )2 ). Xiong [23] researched system (1.1) with h(x, y) =

j =1 y 2 − bi2 ]Ki

m 

[(x − ai )2 +

i=1

where ai = 0, bi > 0 and ai2 = bi2 by using the Melnikov function of first order. As we can see, most of the papers are studied through the first order Melnikov function because of the complexity of calculating and studying the higher order Melnikov function. To the best of our knowledge, there are few papers studying the perturbations of integrable systems by the Melnikov function of higher order, referenced in [3,6,10,13–15,20]. As we all know, the number of limit cycles of system can be studied through the first nonidentically-zero Melnikov function Mk (h), namely, the k−th order Melnikov function. Li and Zhang [15] derived an explicit formula of the Melnikov function of second order with some restrictions which here we denote by M˜ 2 (h). Jiang and Han [14] gave a explicit form of M2 (h) without any restriction and proved M2 (h) = M˜ 2 (h). Françoise [6] provided a algorithm to calculate the Melnikov function of high order, and Iliev [13] gave a method for calculating the Melnikov function of second order with elliptic or hyperellipitic Hamiltonians under polynomial perturbations, and it also gave an explicit expression of M2 (h). Recently Gavrilov and Iliev [10] studied a cubic system in the plane which is a polynomial Hamiltonian system with degree 4 under cubic polynomial perturbations by using the Melnikov function of high order. Among the above papers, there are few papers that studied perturbations of polynomial systems with a straight line of critical points by using the Melnikov function of high order. Buic˘a et al. [3] considered system (1.1) with h(x, y) = 1 + x under quadratic perturbations and they proved the maximum number of limit cycles is 3 by the Melnikov function of the first three order and it has precisely 3 limit cycles. In this paper, we consider system (1.1) with h(x, y) = (1 + ax + by)2 and b = 0, making the linear change of variable as    x cosθ = y sinθ

−sinθ cosθ

  u , v

where tanθ = ab , and a, b ∈ R. Then we make a scaling of cu → U, cv → V and denote U and V by x and y. After the transformation, system (1.1) with h(x, y) = (1 + ax + by)2 passes to another system (1.1) with h(x, y) = (1 + x)2 . So in the following, we only consider the case that h(x) = (1 + x)2 by using the Melnikov function of any order based on the algorithm of Françoise [6] and Iliev [13]. More precisely, we consider a perturbed system as 

x˙ = y(1 + x)2 − P (x, y), y˙ = −x(1 + x)2 + Q(x, y),

(1.2)

where  > 0 is a small parameter and P (x, y) and Q(x, y) are arbitrary quadratic polynomials in x, y given by

JID:YJDEQ AID:9972 /FLA

4

[m1+; v1.304; Prn:10/09/2019; 10:19] P.4 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

P (x, y) = a00 + a10 x + a01 y + a20 x 2 + a11 xy + a02 y 2 , Q(x, y) = b00 + b10 x + b01 y + b20 x 2 + b11 xy + b02 y 2 . For the case h(x) = (1 + x)m , Buic˘a et al. [3] studied the case m = 1 under the quadratic polynomials perturbation through the Melnikov function of the first three order. Here, we consider the case h(x) = (1 + x)2 under the same quadratic polynomial perturbation in order to study the change of the Melnikov function of any order and the number of limit cycles from system (1.2) as m grows. Further more, Gasull et al. [8] considered system (1.1) with h(y) = (1 − y)m under the perturbation of polynomials in degree n. They obtained that when n ≥ m − 1, the maximum number of limit cycles is n, which is irrelevant to m, by using the Melnikov function of first order. Here we guess when n ≥ m − 1, the maximum number of limit cycles is related to the degree m by using the Melnikov function of higher order. Suppose  = 0, there is a family of periodic orbits h surrounding the orgin, where h is defined by H (x, y) = 12 (x 2 + y 2 ) = h in the region 0 < x 2 + y 2 < 1. As we all know, if we fix a transversal segment to the flow in (1.2) and use the energy level h to parameterize it, we can get the corresponding displacement function as follows, 1 d(h, ) = M1 (h) +  2 M2 (h) +  3 M3 (h) + · · · , h ∈ (0, ), 2

(1.3)

where Mk (h) is called the k−th order Melnikov function of system (1.2). It’s well known that the Melnikov function of first order of system (1.2) has the form  M1 (h) = (h)

Q(x, y) P (x, y) dx + dy, (1 + x)2 (1 + x)2

if M1 (h) ≡ 0, each simple zero h0 ∈ (0, 1/2) of M1 (h) corresponds to a limit cycle of the perturbed system (1.2) near h0 . If M1 (h) ≡ 0, the next term M2 (h) is important to study the limit cycles of the perturbed system. Similarly, if M2 (h) ≡ 0, we should use M3 (h) for studying the limit cycles and the procedure goes on, for more results about the Melnikov functions, refer to [3,6,10,13–15,20]. Our main results are: Theorem 1.1. For k ≥ 1, let Mk (h) be the Melnikov functions associated to system (1.2). Then M1 (h) has at most 2 zeros, taking into account their multiplicities. If M1 (h) ≡ 0, then M2 (h) has also at most 2 zeros, taking into account their multiplicities. If M1 (h) ≡ M2 (h) ≡ 0, then M3 (h) has also at most 2 zeros, taking into account their multiplicities. If M1 (h) ≡ M2 (h) ≡ M3 (h) ≡ 0, then M4 (h) has at most 3 zeros, taking into account their multiplicities. If M1 (h) ≡ M2 (h) ≡ M3 (h) ≡ M4 (h) ≡ 0, then Mk (h) ≡ 0 for all k ≥ 5, namely, the perturbation is integrable. Therefore, taking into account the expansion of the displacement map (1.3) of any order in , the system (1.2) has at most three limit cycles which bifurcate from the period annulus of the unperturbed system and this upper bound for the number of limit cycles is reached. Theorem 1.2. The perturbation (1.2) is integrable if and only if at least one set of the following conditions is satisfied:

JID:YJDEQ AID:9972 /FLA

[m1+; v1.304; Prn:10/09/2019; 10:19] P.5 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

5

1) a00 = a02 = a10 = a20 = b01 = b11 = 0; 2) a00 = a10 = −b01 , a01 = a11 = −2b02 , a02 = a20 = b11 = 0; 3) a00 = a20 − b01 , a01 = a11 , a10 = 2a20 − b01 , b02 = − 12 a11 , a02 = b11 = 0; 4) a00 = a10 = a20 = b00 = b01 = b02 = 0, a02 = b11 , b10 = a01 , b20 = a11 . Remark 1.1. Consider the system (1.1) with h(x) = (1 + x)m , Buic˘a et al. [3] studied system (1.1) with h(x) = 1 + x under the perturbations of arbitrary polynomials of degree n = 2 by using the Melnikov functions up to the third order. They didn’t study the Melnikov function of fourth order further more because the expression was too complex. And they obtained three limit cycles by the Melnikov function of third order. But with the increase of m, in the present paper, we consider system (1.1) with h(x) = (1 + x)2 with arbitrary polynomials of degree n = 2 by the Melnikov function of any high order. We obtain three limit cycles by the Melnikov function of fourth order. Moreover if M1 (h) = M2 (h) = M3 (h) = M4 (h) ≡ 0, Mk (h) ≡ 0 for k ≥ 5. With regard to the results of the two articles, we can clearly observe the differences between them by the following table: The upper bounds h(x, y) (1 + x) (1 + x)2

Mk (h)

M1 (h)

M2 (h)

M3 (h)

M4 (h)

M5 (h)

2 2

2 2

3 2

– 3

– 0

It is worth mentioning that the maximum number of limit cycles is obtained by using the k−th order Melnikov function with the conditions that Mj ≡ 0, j = 1, · · · , k − 1. Tigan[20] studied system (1.1) with h(y) = 1 + y under quadratic perturbations by using the Melnikov function of any order, and obtained the maximum number of limit cycles by the Melnikov function of any order. Then it has shown that the upper bounds for the number of limit cycles are 2, 3, 2 and 0 by M1 (h), M2 (h), M3 (h) and M4 (h), respectively. It’s contradictory to the results obtained by Buic˘a et al. [3]. Hence there are some mistakes in the proving process. The author in [20] ignored the independence of the generators. In fact, in [20], the generators of its Melnikov function of first order are not independent. Remark 1.2. Gasull et al. [8] considered system (1.1) with h(x, y) = (1 + x)m under perturbations of polynomials in degree n, and proposed that the maximum number of limit cycles is n when n ≥ m − 1 by the Melnikov function of first order. It means that the value m in (1 + x)m has no influence on the number of limit cycles by the Melnikov function of first order when n ≥ m − 1. While by using the Melnikov function of higher order, the upper bound for the number of limit cycles is related to the degree m. For m = 2 and n = 2, in the present paper we improve the maximum number of limit cycles to 3 by the Melnikov function of any order, that is, the growing of m has an impact on the maximum number of limit cycles when n ≥ m − 1 by the Melnikov function of higher order. 2. Preliminaries We list the following Lemmas and Remarks that shall be used to prove the main theorems. The algorithm of calculating Mk (h) is given by Françoise [6] and Iliev [13]. We rewrite the Assertion 2.1 of [13] in the following Lemma.

JID:YJDEQ AID:9972 /FLA

6

[m1+; v1.304; Prn:10/09/2019; 10:19] P.6 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

Lemma 2.1. Assume h are the period annulus defined by H (x, y) = h, the polynomial function  H (x, y) and the 1-form ω satisfy h ω ≡ 0, if and only if there are two analytic functions q(x, y) and Q(x, y) in a neighborhood of h such that ω = qdH + dQ. In fact, the proof of Lemma 2.1 uses only the fact that they are uniformly Lipschitz continuous there. In the present paper, we consider system (1.2) up to first order in , and we rewrite ω as ω=

Q(x, y) P (x, y) dx + dy, (1 + x)2 (1 + x)2

then the system (1.2) can be written in a Pfaffian form dH = ω. Following Françoise’s algorithm [6] and Lemma 2.1, we propose the following lemma:  Lemma 2.2. The Melnikov function of first order is given by M1 (h) = h ω, rewrite ω in the   ¯ 0 + N0 , then M1 (h) = form ω = q¯0 dH + d Q H =h ω = H =h N0 ;  ˜ 0. If M1 (h) ≡ 0, that is h N0 ≡ 0, according to Lemma 2.1, we have N0 = q˜0 dH + d Q ¯ 0 + Q˜ 0 . In this way, ω = q0 dH + dQ0 , q0 ω is written in the Denote by q0 = q¯0 + q˜0 , Q0 = Q   ¯ 1 + N1 , then M2 (h) = decomposed form q0 ω = q¯1 dH + d Q q ω = 0  h N1 ; h  ˜ 1 . In this If M2 (h) ≡ 0, that is h N1 ≡ 0, according to Lemma 2.1, we have N1 = q˜1 dH + d Q ¯ ˜ way, q0 ω = q1 dH + dQ1 , where q1 = q¯1 + q˜1 , Q  1 = Q1 +Q1 , q1 ω is written in the decomposed ¯ 2 + N2 , then M3 (h) = form q1 ω = q¯2 dH + d Q h q1 ω = h N2 ;  ˜ 2 . In this N2 ≡ 0, according to Lemma 2.1, we have N2 = q˜2 dH + d Q If M3 (h) ≡ 0, that is h

¯ ˜ way, q1 ω = q2 dH + dQ2 , where q2 = q¯2 + q˜2 , Q  2 = Q2 +Q2 , q2 ω is written in the decomposed ¯ 3 + N3 , then M4 (h) = form q2 ω = q¯3 dH + d Q q ω = h 2 h N3 ; Follow this way, if M1 (h) = M2 (h) = · · · = Mi−1 (h) ≡ 0, according to Lemma 2.1, qj ω, j ≤ i − 2 can be written  as qj ω = qj +1 dH + dQj +1 , and when j = −1, we assume q−1 = 1. Then ¯ i−1 + Ni−1 . we have Mi (h) = h qi−2 ω = h Ni−1 , where qi−2 ω = q¯i−1 dH + d Q   Remark 2.1. In Lemma 2.2, we know Mi (h) = h qi−2 ω = h qi−2 q0 dH + qi−2 dQ0 . It is apparent that qi−2 is important to calculate Mi (h), but in the decomposition qi ω = qi+1 dH + dQi+1 , dQi is not used in the subsequent calculative process when i > 0. Therefore we don’t show the specific form of dQi for i ≥ 1 when we give the decomposition of qi ω in the following sections. In fact, we only need consider Q0 in the decomposition of 1-form ω among all the functions Qi , i ∈ Z by qi ω = qi (q0 dH + dQ0 ). In Lemma 2.2 and Remark 2.1, the decomposition of qi ω is not unique. The complex of expressions for Mk (h) is related to the decomposition of qi ω. From Remark 2.1, Qi for i ≥ 1 are not used in the subsequent calculations. So the principle of our decomposition method is to simplify the expressions of qi . In other words, the complexity is reflected in the expressions of Qi to some extent. For example, in [3], Mk (h) for k ≥ 4 have not been further studied because the expression of M4 (h) is too complex. In the next lemma, we give an algorithm to decompose qi ω

JID:YJDEQ AID:9972 /FLA

[m1+; v1.304; Prn:10/09/2019; 10:19] P.7 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

7

which can simplify the expression of Mk (h). Under the decomposition, the calculating process for Mi (h) can be relatively simple. Denote by k ωij =

xi yj dx, (1 + x)k

δijk =

xi yj dy, (1 + x)k

0 ≤ i + j ≤ k,

k≥1

and  Jk (h) =

k δ00 ,

k ≥ 1.

(2.1)

H (x, y)=h k , δk Then we claim that Jk (h) are the generators of Mi (h). In the following we decompose ωij ij k k k k into the combination of δ00 , hij (x, H )dH and dQij (x, H ), where dQij are the perfect differential functions in x and H .

Lemma 2.3. All the 1-forms ωij and δij can be expressed as follows: (i) For 0 ≤ i + j ≤ 2 with k = 2, we have 1 ), 1+x

y 1 , ) + δ00 1+x 2 2H 1 2 ω02 = dH + d(− −x + + 2ln(1 + x)), 1+x 1+x 1+x y 2 2 2 1 = dH − 2H δ00 + δ00 − 2δ00 + dy, ω11 (1 + x)2

2 ω00 = d(−

2 = d(− ω01

2 = d(ln(1 + x) − ω10

x ), 1+x

2 = d(x − ω20

1 − 2ln(1 + x)), 1+x

2 = δ01

1 x dH − d(ln(1 + x) − ), 1+x (1 + x)2

2 = δ11

x 1 dH − d(x − − 2ln(1 + x)), 1+x (1 + x)2

2 2 1 = δ00 − 2δ00 + dy, δ20

2 1 2 = δ00 − δ00 , δ10

2 2 2 1 δ02 = 2H δ00 − δ00 + 2δ00 − dy.

(ii) For 0 ≤ i + j < k with k > 2, we have if j is even, k ωij = hkij (x, H )dH + d(Qkij (x, H )),

δijk =

j/2 k−i−j +2r k (2H )r L(δ00 , · · · , δ00 ), r=0

if j is odd,

(2.2) (2.3)

JID:YJDEQ AID:9972 /FLA

[m1+; v1.304; Prn:10/09/2019; 10:19] P.8 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

8

k ωij = d(Sijk (x, y)) +

(j −1)/2

k−i−j +2r

(2H )r L(δ00

k−1 , · · · , δ00 ),

(2.4)

r=0

δijk = hkij (x, H )dH + d(Qkij (x, H )),

(2.5)

where hkij (x, H ) and Qkij (x, H ) are functions of x and H , while Sijk (x, y) are functions of x and 1 , · · · , δ k ) is the linear combination of δ l (l = 1, 2, · · · , k) y. It’s worth mentioning that L(δ00 00 00 k  l with a ∈ Z. 1 , · · · , δk ) = in Z, namely, L(δ00 al δ00 l 00 l=1

2 and δ 2 by partial integration and the laurent expansion method. Proof. (i) We decompose ωij ij Here we omit some formulas whose decompositions are obvious. First of all, we have by definition 2 ω01 =

y 1 y 1 dx = yd(− , ) = d(− ) + δ00 2 (1 + x) 1+x 1+x

2H − x 2 dx (1 + x)2 1 2H dx − d(x − − 2ln(1 + x)) = 2 1+x (1 + x) 1 1 = 2H d(− ) − d(x − − 2ln(1 + x)) 1+x 1+x 2H 2H 1 = dx + d(− −x + + 2ln(1 + x)) 2 1+x 1+x (1 + x)

2 = ω02

and 2 ω11 =

xy y dx = dx 2 2 (1 + x) 2(1 + x)2 y = d(2H − y 2 ) 2(1 + x)2 =

y2 y dH − dy 2 (1 + x) (1 + x)2

2H − x 2 y dH − dy (1 + x)2 (1 + x)2 y 1 1 1 = dH − 2H dy + dy − 2 dy + dy 2 2 2 (1 + x) (1 + x) (1 + x) (1 + x) y 2 2 1 dH − 2H δ00 + δ00 − 2δ00 + dy. = (1 + x)2 =

For δij2 (0 ≤ i + j ≤ 2), with laurent expansion at x = −1, we have 2 1 2 δ10 = δ00 − δ00 ,

2 2 1 δ20 = δ00 − 2δ00 + dy,

JID:YJDEQ AID:9972 /FLA

[m1+; v1.304; Prn:10/09/2019; 10:19] P.9 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

9

and 2 δ02 =

y2 2H − x 2 2 2 1 dy = dy = 2H δ00 − δ00 + 2δ00 − dy. (1 + x)2 (1 + x)2

2 and δ 2 , we obtain For δ01 11 2 δ01 =

y 1 dy = dy 2 2 (1 + x) 2(1 + x)2 1 = d(2H − x 2 ) 2(1 + x)2 1 x = dH − dx (1 + x)2 (1 + x)2 1 x = dH − d(ln(1 + x) − ) (1 + x)2 1+x

and 2 δ11 =

xy x dy = dy 2 (1 + x)2 2(1 + x)2 x = d(2H − x 2 ) 2(1 + x)2 x2 x dH − dx (1 + x)2 (1 + x)2 x 1 = dH − d(x − − 2ln(1 + x)). (1 + x)2 1+x =

(ii) For j is even, taking j = 2m, m ∈ Z, we have k ωij =

xi yj x i y 2m dx = dx (1 + x)k (1 + x)k =

=

x i (2H − x 2 )m dx (1 + x)k m  xi (2H )r (−1)m−r x 2(m−r) r=0

(1 + x)k

dx

m (−1)m−r x i+2(m−r) (2H )r dx. = (1 + x)k r=0

Denote by

x Qkijr

= 0

(−1)m−r s i+2(m−r) ds, (1 + s)k

(2.6)

JID:YJDEQ AID:9972 /FLA

[m1+; v1.304; Prn:10/09/2019; 10:19] P.10 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

10

then m

(2H )r

r=0

=

(−1)m−r x i+2(m−r) dx (1 + x)k

m (2H )r d(Qkijr ) r=0

= d(Qkij0 ) +

m [d((2H )r Qkijr ) − 2r Qkijr H r−1 dH ], r=1

that is, k ωij = hkij (x, H )dH + d(Qkij (x, H )),

(2.7)

where hkij (x,

H) =

m

−2

r

Qkijr H r−1 ,

Qkij (x,

H) =

r=1

m

(2H )r (Qkijr ).

r=0

Next we consider δijk , for j even. We get δijk =

xi yj x i (y 2 )m dy = dy (1 + x)k (1 + x)k = =

x i (2H − x 2 )m dy (1 + x)k m x i (2H )r (−1)m−r (x 2 )m−r

(1 + x)k

r=0

=

m

(2H )r (−1)m−r

x i+2m−2r dy (1 + x)k

(2H )r (−1)m−r

(x + 1 − 1)i+2m−2r dy (1 + x)k

r=0

=

m r=0

=

m i+2m−2r r=0

=

j/2

dy

k−n (2H )r ar,k−n δ00

n=0 k−i−j +2r

(2H )r L(δ00

k , · · · , δ00 ),

r=0

where n ar, k−n = (−1)m−r Ci+2m−2r (−1)i+2m−2r−n ,

(2.8)

JID:YJDEQ AID:9972 /FLA

[m1+; v1.304; Prn:10/09/2019; 10:19] P.11 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

11

then we get formula (2.3). For j which is odd, taking j = 2m + 1, it follows k ωij =

=

x i y 2m+1 dx (1 + x)k (x + 1 − 1)i y 2m+1 dx (1 + x)k i 

= =

a=0

Cia (1 + x)a (−1)i−a y 2m+1 (1 + x)k

i

dx

Cia (1 + x)a−k (−1)i−a y 2m+1 dx

a=0

=

i

k−a Cia (−1)i−a ω0j .

a=0 k−a For ω0j (a = 0, · · · , i), by definition, we obtain

k−a ω0j =

yj dx (1 + x)k−a

= y j d(−

1 ) (k − a − 1)(1 + x)k−a−1

jy j −1 yj ) + dy = d(− (k − a − 1)(1 + x)k−a−1 (k − a − 1)(1 + x)k−a−1 = d(−

(2.9)

j yj )+ δ k−a−1 . (k − a − 1)(1 + x)k−a−1 k − a − 1 0, j −1

Combining the formula (2.8)–(2.9), we get (2.4). Finally, we consider δijk when j is odd, a similar calculation with 2.6 gives to δijk =

x i y 2m+1 1 x i y 2m dy = dy 2 k (1 + x) 2 (1 + x)k =

1 x i y 2m d(2H − x 2 ) 2 (1 + x)k

x i+1 y j −1 x i (2H − x 2 )m dH − dx = (1 + x)k (1 + x)k =

x i (2H − x 2 )m k dH − ωi+1, j −1 . (1 + x)k

Noticing that (2.7) and (2.10), we get the formula (2.5). This ends the proof. 2

(2.10)

JID:YJDEQ AID:9972 /FLA

[m1+; v1.304; Prn:10/09/2019; 10:19] P.12 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

12

k and δ k with 0 ≤ i + j ≤ 2 are used for the Melnikov Remark 2.2. The decomposition of ωij ij function of first order. When we consider the Melnikov function of any high order of system (1.2) k and δ k with i + j = k for k > 2 don’t appear. under the quadratic perturbation, the terms ωij ij

Lemma 2.4. The generators Jk (h) given by (2.1) can be easily obtained via Maple software as  J1 (h) = h



J2 (h) = h



J3 (h) = h



J4 (h) = h



J5 (h) = h



J6 (h) = h



J7 (h) = h

1 1 dy = 2π(1 − √ ), 1+x 1 − 2h 1 4πh dy = − √ , 2 (1 + x) ( 1 − 2h)3 1 6πh dy = − √ , (1 + x)3 ( 1 − 2h)5 1 4πh(h + 2) dy = − √ , 4 (1 + x) ( 1 − 2h)7

(2.11)

1 5πh(3h + 2) dy = − √ , (1 + x)5 ( 1 − 2h)9 1 6π(h2 + 6h + 2)h dy = − , √ 6 (1 + x) ( 1 − 2h)11 1 7π(5h2 + 10h + 2)h dy = − . √ (1 + x)7 ( 1 − 2h)13

3. The calculation of Mk (h) With the above notations, ω can be written as the following form, 2 2 2 2 2 2 ω = a00 δ00 + a10 δ10 + a01 δ01 + a20 δ20 + a11 δ11 + a02 δ02 2 2 2 2 2 2 + b00 ω00 + b10 ω10 + b01 ω01 + b20 ω20 + b11 ω11 + b02 ω02 .

(3.1)

3.1. The Melnikov function of first order M1 (h) We begin with the decomposition of 1-form ω. Lemma 3.1. The 1-form ω given by (3.1) can be expressed in the following way ¯ 0 + N0 , ω = q¯0 dH + d Q where

(3.2)

JID:YJDEQ AID:9972 /FLA

[m1+; v1.304; Prn:10/09/2019; 10:19] P.13 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

q¯0 = q¯0 (x, y) =

b11 y 2b02 + a11 a01 − a11 + , + 2 (1 + x) 1+x (1 + x)2

¯ 0 (x, y, H ) = ¯ 0 = dQ dQ +

13

b00 (b10 − a01 )x (b20 − a11 )x 2 dx + dx + dx (1 + x)2 (1 + x)2 (1 + x)2

b01 2b02 b02 y 2 b01 y dx − dx + (b11 + a20 − a02 )dy − dy + dH, 2 (1 + x) 1+x (1 + x)2 1+x

1 2 + (b11 + a00 − a10 + a20 − a02 )δ00 N0 = (b01 − 2b11 + a10 − 2a20 + 2a02 )δ00 2 + (2a02 − 2b11 )H δ00 . 2 and δ 2 with 0 ≤ i + j ≤ 2. Proof. Using Lemma 2.3 (i), we can obtain the decompositions of ωij ij A routine computation gives rise to the formula (3.2). 2

According to Lemma 2.2, it follows that  M1 (h) =

 ω=

h

¯ 0 + N0 . q¯0 dH + d Q

H (x, y)=h

As we have defined by (2.1), let Jk (h) = easy to see that



k H (x, y)=h δ00 ,



k ≥ 1 as the generators of Mi (h). It is

 q¯0 dH = 0,

H (x, y)=h

¯ 0 = 0, dQ

H (x, y)=h

so  M1 (h) =

N0 = A1 J1 (h) + A2 J2 (h) + A3 hJ2 (h),

(3.3)

H (x, y)

where A1 = b01 − 2b11 + a10 − 2a20 + 2a02 , A2 = b11 + a00 − a10 + a20 − a02 , A3 = 2a02 − 2b11 . Theorem 3.1. M1 (h) has at most two zeros, namely, system (1.2) has at most two limit cycles by the Melnikov function of first order and the upper bound for the number of limit cycles is reached. Proof. Taking z = into (3.3) gives to M1 (

√

1 − 2h, z ∈ (0, 1). Substituting (2.11) from Lemma 2.4 and z =

1 − z2 A3 2π 1 ) = − 3 ( A3 z4 − A1 z3 + (A1 − A2 − A3 )z2 + A2 + ). 2 z 2 2

√

1 − 2h

JID:YJDEQ AID:9972 /FLA

[m1+; v1.304; Prn:10/09/2019; 10:19] P.14 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

14

Then, the zeros of M1 (h), h ∈ (0, 12 ) is equivalent to the zeros of 1 A3 M1 (z) = A3 z4 − A1 z3 + (A1 − A2 − A3 )z2 + A2 + , z ∈ (0, 1). 2 2 It is easy to show that {1, z2 , z3 , z4 } is a Chebyshev system for z > 0. So the equation above has at most three zeros for z > 0. Clearly, z = 1 is one zero of M1 (z) (that is, h = 0), which corresponds to the center of system (1.2). Hence M1 (h) has at most two zeros in h ∈ (0, 12 ), taking into account their multiplicities. Moreover, it is easy to check that there are enough coefficients such that M1 (h) has exactly two zeros. This completes the proof of Theorem 3.1. 2 3.2. The Melnikov function of second order M2 (h) If M1 (h) ≡ 0, then A1 = A2 = A3 = 0, that is, a00 = a20 − b01 ,

a10 = 2a20 − b01 ,

a02 = b11 .

(3.4)

From now on we assume that M1 (h) ≡ 0. Then the calculation of M2 (h) makes sense. By sub˜ 0 = 0 by Lemma 2.2. It is stituting (3.4) into N0 , we obtain N0 = 0, which follows q˜0 = 0, d Q obvious that q0 = q¯0 =

b11 y 2b02 + a11 a01 − a11 + , + 2 1+x (1 + x) (1 + x)2

¯0 = dQ0 = d Q +

b00 (b10 − a01 )x (b20 − a11 )x 2 dx + dx + dx 2 2 (1 + x) (1 + x) (1 + x)2

(3.5)

b01 2b02 b02 y 2 b01 y dx − dx + a20 dy − dy + dH 2 2 (1 + x) 1+x (1 + x) 1+x

and ω = q0 dH + dQ0 . Lemma 3.2. The 1-form q0 ω can be decomposed as follows: ¯ 1 + N1 , q0 ω = q¯1 dH + d Q

(3.6)

where 2b02 1 2 + [ (a01 b02 − a11 b02 + b01 b11 ) − b01 b11 ] 1+x 3 (1 + x)3 1 + [b02 (2b02 + a11 ) + a20 b11 ] , (1 + x)2

q¯1 = q02 − q0

1 2 3 3 N1 = B1 δ00 + B2 δ00 + B3 δ00 + B4 H δ00 ,

with

(3.7)

JID:YJDEQ AID:9972 /FLA

[m1+; v1.304; Prn:10/09/2019; 10:19] P.15 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

15

B1 = −b02 b11 − b11 (a11 − b20 ) + a20 (2b02 + a11 ), 1 1 1 B2 = 2b02 b11 − a01 b11 − a11 b01 + a11 b11 − b01 b02 + b10 b11 2 2 2 − b11 b20 + a01 a20 − a11 a20 , 1 B3 = (a01 b01 + a01 b11 − a11 b01 − a11 b11 + b00 b11 − b10 b11 + b11 b20 ) 3 − b02 b11 + b01 a11 − b01 a01 , B4 = 2b02 b11 . ¯ 1 is a function in (x, y, H ). In the decomposition of the formula (3.6), Q Proof. We rewrite q0 (x, y) and dQ0 (x, y, H ) as the following form q0 (x, y) = f0 (x) + g0 (x, y), dQ0 (x, y, H ) = 1 (x) + 2 (x, y) + h(x)dH, where f0 (x) =

2b02 + a11 a01 − a11 , + 1+x (1 + x)2

1 (x) =

b00 (b10 − a01 )x (b20 − a11 )x 2 dx + dx + dx, 2 2 (1 + x) (1 + x) (1 + x)2

g0 (x, y) =

b11 y , (1 + x)2

b01 y b01 b02 y 2 dx − dx + a20 dy, dy + 2 (x, y) = (1 + x)2 1+x (1 + x)2 2b02 . h(x) = − 1+x

(3.8)

We write q0 ω as q0 ω = q0 (q0 dH + dQ0 ) = q02 dH + q0 dQ0 .

(3.9)

In the above formula, we first consider how to decompose q0 dQ0 . A routine computation gives rise to q0 dQ0 = (f0 (x) + g0 (x, y))( 1 (x) + 2 (x, y) + h(x)dH ) = f0 (x) 1 (x) + f0 (x) 2 (x, y) + g0 (x, y) 1 (x) + g0 (x, y) 2 (x, y) + q0 h(x)dH 4 = f0 (x) 1 (x) + q0 h(x)dH + α03

y3 y2 y 4 4 dx + α dx + α01 dx 02 (1 + x)4 (1 + x)4 (1 + x)4

y2 y y y 3 2 3 dx + α01 dx + α01 dx + β01 dy 3 3 2 (1 + x) (1 + x) (1 + x) (1 + x)3 1 y 1 y 3 2 2 1 dy + β01 dy + β00 dy + β00 dy + β00 3 2 2 1+x (1 + x) (1 + x) (1 + x)

3 + α02

(3.10)

JID:YJDEQ AID:9972 /FLA

[m1+; v1.304; Prn:10/09/2019; 10:19] P.16 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

16

4 4 4 4 4 4 3 3 = f0 (x) 1 (x) + q0 h(x)dH + α03 ω03 + α02 ω02 + α01 ω01 + α02 ω02 3 3 2 2 3 3 3 3 2 2 2 2 1 1 + α01 ω01 + α01 ω01 + β01 δ01 + β00 δ00 + β01 δ01 + β00 δ00 + β00 δ00 , k and δ k respectively. In (3.10), α k and β k are where αijk and βijk are the coefficients of ωij ij ij ij expressed as follows, 4 α03 = b02 b11 ,

4 α02 = a01 b02 − a11 b02 + b01 b11 ,

4 α01 = a01 b01 + a01 b11 − a11 b01 − a11 b11 + b00 b11 − b10 b11 + b11 b20 , 3 = −a01 b11 + a11 b01 + 2a11 b11 + 2b01 b02 + b10 b11 − 2b11 b20 , α01 3 α02 = b02 (2b02 + a11 ),

2 α01 = −b11 (a11 − b20 ),

3 β00 = b01 (a11 − a01 ),

2 β01 = a20 b11 ,

3 β01 = −b01 b11 ,

1 β00 = a20 (2b02 + a11 ),

2 = a01 a20 − a11 a20 − a11 b01 − 2b01 b02 . β00 4 , ω3 , δ 3 and δ 2 can be decomposed into hk dH + From Lemma 2.3, we know that ω02 ij 02 01 01 4 , ω4 , ω3 , ω2 , δ 3 , δ 2 and δ 1 can be decomposed into generators, the dQkij , similarly, ω03 01 01 01 00 00 00 decompositions have shown as follows,

2 1 1 2H 1 1 dH + d( + − ), − 3 2 3 3 (1 + x) 1 + x (1 + x) 3(1 + x) 3(1 + x)3 1 2 H 1 3 ω02 = dH − d(ln(1 + x) + + ), − 2 2 1 + x 2(1 + x) (1 + x) (1 + x)2 1 1 1 3 δ01 = dH + d( ), − 1 + x 2(1 + x)2 (1 + x)3 1 x 2 δ01 = dH − d(ln(1 + x) − ), (1 + x)2 1+x 4 ω02 =

4 3 1 2 3 = 2H δ00 − δ00 + 2δ00 − δ00 + d(− ω03

1 3 y 4 ω01 = δ00 − d( ), 3 3(1 + x)3 y 2 1 = δ00 − d( ). ω01 1+x

y3 ), 3(1 + x)3

1 2 y 3 ω01 = δ00 − d( ), 2 2(1 + x)2 (3.11)

Combining (3.9)–(3.11), we obtain (3.6). According to Remark 2.1, we didn’t give the ex¯ 1 . Hence the proof has been completed. 2 pression of Q From Lemma 2.2, we have  M2 (h) = H (x, y)=h

which follows,

 q0 ω = H (x, y)=h

N1 ,

JID:YJDEQ AID:9972 /FLA

[m1+; v1.304; Prn:10/09/2019; 10:19] P.17 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

17

 M2 (h) =

N1 = B1 J1 (h) + B2 J2 (h) + B3 J3 (h) + B4 hJ3 (h),

(3.12)

H (x, y)=h

where Jk (h)(k = 1, 2, 3) and Bi (i = 1, 2, 3, 4) are defined by (2.1) and (3.7). Theorem 3.2. If M1 (h) ≡ 0, M2 (h) has at most two zeros, namely, system (1.2) has at most two limit cycles by the Melnikov function of second order and the maximum number of limit cycles is reached. Proof. Taking z = into (3.12) gives to

√

1 − 2h, z ∈ (0, 1). Substituting (2.11) from Lemma 2.4 and z =

M2 (

√

1 − 2h

1 − z2 2π ) = − 5 M2 (z), 2 z

where 3 3 3 3 3 M2 (z) = −B1 z5 + (B1 − B2 + B4 )z4 + (B2 − B3 − B4 )z2 + B3 + B4 . 4 2 2 2 4 By the definition of Chebyshev system we get that {1, z2 , z4 , z5 } is a Chebyshev system for z > 0. So M2 (z) has at most three zeros for z > 0. Furthermore, through calculating, we find z = 1 is one zero of M2 (z) = 0 (that is, h = 0), which corresponds to the center of system (1.2). Therefore M2 (z) has at most two zeros for z ∈ (0, 1) and so does M2 (h) for h ∈ (0, 12 ). With some proper coefficients, M2 (h) has exactly two zeros. Hence we complete the proof. 2 3.3. The Melnikov function of third order M3 (h) If M2 (h) ≡ 0, we have 3 3 3 3 3 −B1 = B1 − B2 + B4 = B2 − B3 − B4 = B3 + B4 = 0. 4 2 2 2 4

(3.13)

With the help of Maple software, the formula (3.13) holds if and only if one of the following four cases holds. Case 1: a00 = a02 = a10 = a20 = b01 = b11 = 0.

(3.14)

Case 2: a02 = a20 = b11 = 0,

a00 = a10 = −b01 ,

a01 = a11 = −2b02 .

(3.15)

Case 3: a02 = b11 = 0, Case 4:

a00 = a20 − b01 ,

a01 = a11 ,

a10 = 2a20 − b01 ,

1 b02 = − a11 . 2

(3.16)

JID:YJDEQ AID:9972 /FLA

[m1+; v1.304; Prn:10/09/2019; 10:19] P.18 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

18

b11 = 0,

a00 = a20 − b01 ,

a02 = b11 ,

a10 = 2a20 − b01 ,

b00 b11 = −2a01 a20 + 2a01 b01 + a11 a20 − a11 b01 − 2a20 b02 + 2b01 b02 , b10 b11 = −2a01 a20 + a01 b11 + a11 b01 − 4a20 b02 + 2b01 b02 + b02 b11 ,

(3.17)

b20 b11 = −a11 a20 + a11 b11 − 2a20 b02 + b02 b11 . We assume that M1 (h) = M2 (h) ≡ 0, in following contents we will talk about the four cases. With this assumption, M3 (h) should be calculated. First we consider case 1 under the condition that M1 (h) = M2 (h) ≡ 0. Theorem 3.3. Assume that case 1 holds, we have Mk (h) ≡ 0, (1.2) is integrable.

k ≥ 3, that is the perturbation

Proof. For case 1, M1 (h) = M2 (h) ≡ 0, by substituting (3.14) into (3.5), we have N0 = 0 and ω = q0 dH + dQ0 , where q0 = dQ0 =

2b02 + a11 a01 − a11 , + 1+x (1 + x)2 y2 b00 x x2 dx + (b10 − a01 ) dx + (b20 − a11 ) dx + b02 dx 2 2 2 (1 + x) (1 + x) (1 + x) (1 + x)2 1 dH. + 2b02 1+x

By substituting (3.14) into (3.6), we obtain ¯ 1 + N1 , q0 ω = q¯1 dH + d Q where N1 = 0 and q¯1 =

2 − 2a a + a 2 2 + 3a b + 2b2 a01 a11 01 11 11 02 11 02 + (1 + x)4 (1 + x)2

+

2 − 8a b 1 6a01 a11 + 8a01 b02 − 6a11 11 02 . 3 3 (1 + x)

˜ 1 = 0 by N1 = 0. Hence, It follows q˜1 = 0, d Q q1 = q¯1 =

2 − 2a a + a 2 2 + 3a b + 2b2 a01 a11 01 11 11 02 11 02 + (1 + x)4 (1 + x)2

+ which implies,

2 − 8a b 1 6a01 a11 + 8a01 b02 − 6a11 11 02 , 3 3 (1 + x)

JID:YJDEQ AID:9972 /FLA

[m1+; v1.304; Prn:10/09/2019; 10:19] P.19 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

q1 ω = q1 q0 dH + q1 dQ0 . We rewrite dQ0 as dQ0 = 1 (x) + 2 (x, y) + h(x)dH, where 1 (x) =

b00 x x2 dx + (b − a ) dx + (b − a ) dx, 10 01 20 11 (1 + x)2 (1 + x)2 (1 + x)2

2 (x, y) = b02

y2 dx, (1 + x)2

h(x) = 2b02

1 . 1+x

Then q1 ω can be written as q1 ω = q1 q0 dH + q1 dQ0 = (q1 q0 + q1 h(x))dH + q1 1 (x) + q1 2 (x, y). By directing calculation, it follows, 4 5 6 q1 2 (x, y) = m402 ω02 + m502 ω02 + m602 ω02 ,

where 2 2 m402 = (a11 + 3a11 b02 + 2b02 )b02 ,

1 2 − 8a11 b02 )b02 , m502 = (6a01 a11 + 8a01 b02 − 6a11 3 2 2 − 2a01 a11 + a11 )b02 . m602 = (a01 4 , ω5 and ω6 as From Lemma 2.3, we obtain the decomposition of ω02 02 02 4 ω02 =

2 2 H 1 1 1 1 dH + d(− + − + ), 3 (1 + x)3 3 (1 + x)3 3(1 + x)3 (1 + x)2 1 + x

5 = ω02

1 1 H 1 2 1 1 dH + d(− + − + ), 4 4 4 3 2 (1 + x) 2 (1 + x) 4(1 + x) 3(1 + x) 2(1 + x)2

6 ω02 =

2 2 H 1 1 1 1 dH + d(− + − + ). 4 5 5 5 5 (1 + x) 5 (1 + x) 2(1 + x) 3(1 + x)3 5(1 + x)

¯ 2 + N2 where N2 = 0 and Thus, we have q1 ω = q¯2 dH + d Q q¯2 = q1 q0 + q1 h(x) + From N2 = 0 we have

2 1 2 1 1 1 + + . 3 (1 + x)3 2 (1 + x)4 5 (1 + x)5

19

JID:YJDEQ AID:9972 /FLA

[m1+; v1.304; Prn:10/09/2019; 10:19] P.20 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

20

 M3 (h) =

N2 = 0.

H (x, y)=h

˜ 2 = 0 gives to q˜2 = 0. Hence we have Furthermore, N2 = q˜2 dH + d Q ¯ 2  q2 dH + dQ2 . q1 ω = q¯2 dH + d Q k . With the same calculation, we obtain the expression of q2 ω is a linear combination of ω02 k k k And from Lemma 2.3 we know ω02 can be decomposed into hij dH + dQij . Hence we can get  the conclusion that Ni = 0, i ∈ Z, and Mk (h) = H (x, y)=h Nk−1 = 0, k ≥ 3. This completes the proof. 2

We talk about case 2 and case 3 in the next theorem. Theorem 3.4. Assume that case 2 or case 3 holds, we have Mk (h) ≡ 0, k ≥ 3, that is the perturbation (1.2) is integrable. Proof. Consider the case 2, by substituting (3.15) into (3.5), we have q0 = 0 and N0 = 0. Similarly, by substituting (3.15) into (3.6), we obtain N1 = 0. According to Lemma 2.2, we have q0 ω = q1 dH + dQ1 = 0, that is q1 = 0. Hence,  M3 (h) =

 q1 ω =

H (x, y)=h

N2 = 0.

H (x, y)=h

Performing the same procedure for i from 1 to k, we have qi = 0 and Ni = 0, then 



Mk (h) =

qk−2 ω =

H (x, y)=h

Nk−1 = 0,

k ≥ 3,

H (x, y)=h

which means the perturbation (1.2) is integrable. For case 3, we substitute (3.16) into (3.5) and (3.6), we have q0 = 0, N0 = 0 and N1 = 0. Similarly, we can prove that with the case 3, the perturbation (1.2) is also integrable. Then we complete the proof. 2 Finally we talk about case 4. Lemma 3.3. Assume that the case 4 holds, we get ¯ 2 + N2 , q1 ω = q¯2 dH + d Q where

(3.18)

JID:YJDEQ AID:9972 /FLA

[m1+; v1.304; Prn:10/09/2019; 10:19] P.21 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

21

2b02 1 2 + [ (a01 b02 − a11 b02 + b01 b11 ) − b01 b11 ] 1+x 3 (1 + x)3 y 1 + [b02 (2b02 + a11 ) + a20 b11 ] + b02 b11 . 2 (1 + x) (1 + x)3

q1 = q02 − q0

q¯2 = q1 q0 + q1 h(x) + τ1

1 1 1 y2 + τ + τ + τ , 2 3 4 (1 + x)3 (1 + x)4 (1 + x)5 (1 + x)5

2 3 4 4 5 5 + C2 δ00 + C3 δ00 + C4 H δ00 + C5 δ00 + C6 H δ00 N2 = C1 δ00

with 1 2 2 3 2 b02 + 20a11 a20 b11 + 30a11 b02 + 35a20 b02 b11 + 20b02 + 2b02 b11 ), (10a11 15 1 2 2 τ2 = (30a01 a11 b02 + 30a01 a20 b11 + 40a01 b02 − 30a11 b02 − 30a11 a20 b11 , 30

τ1 =

2 2 − 20b01 b02 b11 − 3b02 b11 ), − 15a11 b01 b11 − 40a11 b02

1 2 2 b02 − 4a01 a11 b02 − 2a01 b01 b11 + 2a11 b02 + 2a11 b01 b11 ), τ3 = (2a01 5 4 2 τ4 = b02 b11 , 5 7 3 2 2 2 a20 − a11 a20 b02 + a11 b02 b11 − 3a20 b02 + b02 b11 , C1 = −a11 2 2 3 2 2 2 − b02 b11 (2a11 + 3b02 ) + a20 (a11 + 3a11 b02 + a20 b11 + 2b02 − b11 ), 4 2 2 2 1 C2 = − a01 a20 b02 + a01 b02 b11 + a11 a20 b02 + a11 b01 b02 , 3 3 3 3 2 3 4 2 2 − b02 b11 + b02 b11 (2a11 + 3b02 ), − a11 b02 b11 − a20 b01 b11 + b01 b02 3 3 2 3 2 2 + 2a20 b11 − b11 (2a01 b02 − 2a11 b02 + b01 b11 ) + b01 b11 , 5 1 1 1 2 2 2 b11 − 2b01 b11 − a20 b11 C3 = (a01 − a11 )b01 b02 + a11 b02 b11 + b01 2 2 4 3 2 6 − b02 b11 ( a01 + 2a11 + 3b02 ) + b11 (2a01 b02 − 2a11 b02 + b01 b11 ), 4 3 5 3 2 , (3.19) C4 = b02 b11 (2a11 + 3b02 ) + 2a20 b11 2 3 2 C5 = − b11 (2a01 b02 − 2a11 b02 + b01 b11 ) + b01 b11 , 5 6 2 . C6 = b11 (2a01 b02 − 2a11 b02 + b01 b11 ) − 2b01 b11 5   Proof. With the case 4, we have M1 (h) = M2 (h) ≡ 0, which means H (x, y)=0 N0 = H (x, y)=0 N1 = 0.

JID:YJDEQ AID:9972 /FLA

[m1+; v1.304; Prn:10/09/2019; 10:19] P.22 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

22

By substituting (3.17) into N1 , we have N1 = b02 b11 (2H

1 3 1 1 dy + dy − dy). 2 (1 + x)2 (1 + x)3 (1 + x)3

˜ 1 by From Lemma 2.1, N1 can be rewritten as N1 = q˜1 dH + d Q show how to rewrite N1 . From the formula (3.20) we have N1 = b02 b11 (2H



H (x, y)=0 N1

(3.20) ≡ 0. Next we

1 3 1 1 dy + dy − dy) 3 2 2 (1 + x) (1 + x) (1 + x)3

x2 + y2 − 1 3 1 dy + dy) (1 + x)3 2 (1 + x)2 x −1 y 3 1 dy + dy 2 + dy) = b02 b11 ( 2 (1 + x)2 (1 + x)2 2(1 + x)3

= b02 b11 (

x + 12 y xy = b02 b11 ( dy + dH − dx) (1 + x)2 (1 + x)3 (1 + x)3 1 y xy 1 = b02 b11 ( dy + dH − dx) dy − 1+x 2(1 + x)2 (1 + x)3 (1 + x)3 y y y ) + b02 b11 dH. − = b02 b11 d( 1 + x (1 + x)2 (1 + x)3

(3.21)

y So we obtain q˜1 = b02 b11 (1+x) 3 . According to Lemma 2.2, combining (3.7) and (3.17) gives to

q1 = q¯1 + q˜1 = q02 − q0

2b02 1 2 + [ (a01 b02 − a11 b02 + b01 b11 ) − b01 b11 ] 1+x 3 (1 + x)3

+ [b02 (2b02 + a11 ) + a20 b11 ]

y 1 + b02 b11 . (1 + x)2 (1 + x)3

Then we have q1 ω = q1 (q0 dH + dQ0 ) = q1 q0 dH + q1 dQ0 . From the formula (3.8), a direct computation gives rise to q1 dQ0 = q1 ( 1 (x) + 2 (x, y)) + q1 h(x)dH = q1 1 (x) + q1 h(x)dH +

k−2 6 k=3 j =1

where

k k ξ0j ω0j +

k−2 5 k=2 j =0

k k η0j δ0j ,

(3.22)

JID:YJDEQ AID:9972 /FLA

[m1+; v1.304; Prn:10/09/2019; 10:19] P.23 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

23

3 2 2 2 ξ01 = −2a11 a20 − 7a11 a20 b02 + 2a11 b02 b11 − 6a20 b02 + 3b02 b11 , 4 2 2 = −6a01 a11 a20 − 10a01 a20 b02 + 2a01 b02 b11 + 6a11 a20 + 3a11 b01 ξ01 2 2 + 10a11 a20 b02 + 10a11 b01b02 − 4a11 b02 b11 + a20 b01b11 + 8b01b02 − 3b02 b11 , 5 2 = −4a01 a20 + 8a01 a11 a20 + 8a01 a11 b01 + 38a01 b01 b02 /3 − 2a01 b02 b11 ξ01 2 2 2 − 4a11 a20 − 8a11 b01 − 38a11 b01 b02 /3 + 2a11 b02 b11 − b01 b11 /3, 6 2 2 = 5b01 (a01 − 2a01 a11 + a11 ), ξ01 4 2 2 3 2 = a11 b02 − a11 a20 b11 + 3a11 b02 − a20 b02 b11 + 2b02 + b02 b11 , ξ02 5 2 2 = 2a01 a11 b02 − 2a01 a20 b11 + 8a01 b02 /3 − 2a11 b02 + 2a11 a20 b11 ξ02 2 2 + 3a11 b01 b11 − 8a11 b02 /3 + 14b01 b02 b11 /3 − b02 b11 , 6 2 2 = a01 b02 − 2a01 a11 b02 + 4a01 b01 b11 + a11 b02 − 4a11 b01 b11 , ξ02 5 = b02 b11 (2a11 + 3b02 ), ξ03 6 = b11 (2a01 b02 − 2a11 b02 + b01 b11 ), ξ03

6 2 ξ04 = b02 b11 ,

2 2 2 = a20 (a11 + 3a11 b02 + a20 b11 + 2b02 ), η00 3 2 2 = 2a01 a11 a20 + 8a01 a20 b02 /3 − 2a11 a20 − a11 b01 − 8a11 a20 b02 /3 η00 2 , − 3a11 b01 b02 − 4a20 b01 b11 /3 − 2b01 b02 4 2 2 2 = a01 a20 − 2a01 a11 a20 − 2a01 a11 b01 − 8a01 b01 b02 /3 + a11 a20 + 2a11 b01 η00 2 + 8a11 b01 b02 /3 + b01 b11 /3, 3 = a20 b11 (2a11 + 3b02 ), η01

5 2 2 η00 = −b01 (a01 − 2a01 a11 + a11 ),

5 η01 = −2b11 b01 (a01 − a11 ),

4 = b11 (2a01 a20 − 2a11 a20 − 2a11 b01 − 3b01 b02 ), η01 4 2 = a20 b11 , η02

5 2 η02 = −b01 b11 ,

5 η03 = 0.

k , ωk and δ k can be decomposed into hk dH + dQk , we give From Lemma 2.3, we know ω02 ij ij 04 01 the decompositions as

4 ω02 =

2 2H x2 dH + d(− ) − dx, 3(1 + x)3 3(1 + x)3 (1 + x)4

5 ω02 =

1 H x2 dH + d(− ) − dx, 2(1 + x)4 2(1 + x)4 (1 + x)5

6 ω02 =

2 2H x2 dH + d(− )− dx, 5 5 5(1 + x) 5(1 + x) (1 + x)6

6 =( ω04

8H 4 2 4 − + − )dH 5(1 + x)5 3(1 + x)3 (1 + x)4 5(1 + x)5

+ d(−

4H 2H 4H x4 4H 2 + − + )+ dx, 3 4 5 5 3(1 + x) (1 + x) 5(1 + x) 5(1 + x) (1 + x)6

JID:YJDEQ AID:9972 /FLA

[m1+; v1.304; Prn:10/09/2019; 10:19] P.24 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

24 3 δ01 =

1 x dH − dx, 3 (1 + x) (1 + x)3

4 = δ01

1 x dH − dx, (1 + x)4 (1 + x)4

5 = δ01

1 x dH − dx. 5 (1 + x) (1 + x)5

k , ωk and δ k can be decomposed as follows, As the same, ω01 03 02

1 2 y 3 ω01 = δ00 + d(− ), 2 2(1 + x)2 1 3 y 4 = δ00 + d(− ), ω01 3 3(1 + x)3 1 4 y 5 = δ00 + d(− ), ω01 4 4(1 + x)4 1 5 y 6 = δ00 + d(− ), ω01 5 5(1 + x)5 3 y3 5 4 4 3 2 = (2H δ00 − δ00 + 2δ00 − δ00 ) + d(− ), ω03 4 4(1 + x)4 3 y3 6 5 5 4 3 = (2H δ00 − δ00 + 2δ00 − δ00 ) + d(− ), ω03 5 5(1 + x)5 4 4 4 3 2 δ02 = 2H δ00 − δ00 + 2δ00 − δ00 , 5 5 5 4 3 = 2H δ00 − δ00 + 2δ00 − δ00 . δ02

(3.23)

¯ 2 is a function in x, y and H . Finally, combining (3.22)–(3.23), we obtain (3.18). Here d Q ¯ Through Remark 2.1, we didn’t give the specific form of d Q2 . Hence we get the proof completed. 2 It is easy to get  M3 (h) =

N2

H (x, y)=h

(3.24)

= C1 J2 (h) + C2 J3 (h) + C3 J4 (h) + C4 hJ4 (h) + C5 J5 (h) + C6 hJ5 (h), where Jk (h), k = 2, 3, 4, 5 and Ci , i = 1, · · · , 6 are defined by (2.1) and (3.19) respectively. Theorem 3.5. If M1 (h) = M2 (h) ≡ 0, and the case 4 holds, M3 (h) has at most two zeros, namely, system (1.2) has at most two limit cycles by the Melnikov function of third order and the upper bound of the number of limit cycles is reached.

JID:YJDEQ AID:9972 /FLA

[m1+; v1.304; Prn:10/09/2019; 10:19] P.25 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

25

Proof. Substituting (2.11) from Lemma 2.4 into (3.24) gives to π M3 (h) = − √ h(2h − 1)(c2 h2 + c1 h + c0 ), ( 1 − 2h)9

(3.25)

where 2 2 2 2 c0 = −8a11 a20 b02 + 4a11 b02 b11 + 16a20 b11 − 16a20 b02 − 8a20 b11 + 6b02 b11 , 2 c1 = 8a01 a20 b02 + 2a01 b01 b02 − 4a01 b02 b11 − 6a11 b01 b02 + 2a11 b02 b11 − 16a20 b11 2 2 2 2 2 2 + 12a20 b01 b11 + 16a20 b02 + 4a20 b11 + b01 b11 − 8b01 b02 − 2b01 b11 − 3b02 b11 , 2 c2 = −4a01 a20 b02 + 4a01 b01 b02 + 2a11 a20 b02 − 2a11 b01 b02 + 4a20 b11 − 6a20 b01 b11 2 2 2 − 4a20 b02 + 2b01 b11 + 4b01 b02 .

From the above formula (3.25) we obtain that M3 (h) has at most four zeros, but h = 0 and h = 12 are corresponding to the center and the homoclinic orbit of system (1.2) respectively. So system (1.2) has at most two zeros in h ∈ (0, 12 ). It’s easy to check that there are enough coefficients such that M3 (h) has exactly 2 simple zeros. The proof has been completed. 2 3.4. The Melnikov function of fourth order M4 (h) Now if M1 (h) = M2 (h) = M3 (h) ≡ 0, M4 (h) is necessary to be considered. With the condition of M1 (h) = M2 (h) = M3 (h) ≡ 0, we have c0 = c1 = c2 = 0. With the help of Maple software, we get four cases as Case 4a: b11 = 0, a00 = a10 = a20 = b00 = b01 = b02 = 0, a02 = b11 , b10 = a01 , b20 = a11 .

(3.26)

Case 4b: b11 = 0, a00 = b00 = 0, a02 = b11 , a10 = b01 = a20 ,

Case 4c:

a01 =

2 b − 8a b2 − 3a b2 + 3b2 b 1 6a20 11 20 02 20 11 02 11 , 2 b02 (2a20 − b11 )

a11 =

2 b − 8a b2 − 4a b2 + 3b2 b 1 8a20 11 20 02 20 11 02 11 , 2 b02 (2a20 − b11 )

b10 = −

3 − 8a 2 b + 3a b2 + 3a b2 − b2 b 1 4a20 20 02 20 11 20 11 02 11 , 2 b02 (2a20 − b11 )

b20 = −

3 − 12a 2 b + 3a b2 + 4a b2 − b2 b 1 8a20 20 02 20 11 20 11 02 11 . 2 b02 (2a20 − b11 )

(3.27)

JID:YJDEQ AID:9972 /FLA

[m1+; v1.304; Prn:10/09/2019; 10:19] P.26 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

26

2 b11 = 0, a00 = b00 = 0, a02 = b11 , a10 = b01 = a20 = b11 , 5 2 + 8b2 2 + 24b2 25b 1 5b02 1 11 02 11 a11 = , b20 = , 10 b02 50 b02 b10 =

(3.28)

2 + 8b2 1 5a01 b02 + 10b02 11 . 25 b02

Case 4d: 2 2 b11 = 0, a00 = a20 − b11 , a02 = b11 , a10 = 2a20 − b11 , 5 5 a01 =

2 b − 80a b2 − 48a b2 + 35b2 b + 4b3 1 80a20 11 20 02 20 11 02 11 11 , 20 b02 (2a20 − b11 )

a11 =

2 b − 8a b2 − 4a b2 + 3b2 b 1 8a20 11 20 02 20 11 02 11 , 2 b02 (2a20 − b11 )

(3.29)

b00 = −

2 − 15a b + 2b2 2 25a20 20 11 11 , 25 b02

b10 = −

3 − 208a 2 b + 30a b2 + 72a b2 − 11b2 b − 4b3 1 160a20 20 02 20 11 11 20 11 02 11 , 20 b02 (2a20 − b11 )

b20 = −

3 − 12a 2 b + 3a b2 + 4a b2 − b2 b 1 8a20 20 02 20 11 20 11 02 11 . 2 b02 (2a20 − b11 )

2 b01 = b11 , 5

First we consider the case 4a. Theorem 3.6. Assume that the case 4a holds, which implies M1 (h) = M2 (h) = M3 (h) ≡ 0, we have Mk (h) ≡ 0, k ≥ 4, that is the perturbation (1.2) is integrable. Proof. For the case 4a, substituting (3.26) into (3.5) and (3.2), we get dQ0 = 0, N0 = 0. It follows ω = q0 dH and qi = q0i+1 by Lemma 2.2. Thus  Mk (h) =

 qk−2 ω =

H (x, y)=h

q0k dH = 0

k ≥ 2.

H (x, y)=h

Hence the theorem has been proved. 2 Before talking about the next three cases, we first give the general form of M4(h). If M1 (h) = M2 (h) = M3 (h) ≡ 0, we have 

 N0 =

H (x, y)=h

Solving for



H (x, y)=h N2

H (x, y)=h

 N1 =

N2 ≡ 0.

H (x, y)=h

≡ 0, in terms of Ci , (i = 1, · · · , 6), which defined in (3.19),

JID:YJDEQ AID:9972 /FLA

[m1+; v1.304; Prn:10/09/2019; 10:19] P.27 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

27

1 7 5 9 7 C1 = C3 + C5 , C2 = − C3 − C5 , C4 = −2C3 − C5 , C6 = −2C5 . 2 8 3 4 2 Hence using the same method for obtaining formula (3.21), it yields 2 3 4 4 5 5 N2 = C1 δ00 + C2 δ00 + C3 δ00 + C4 H δ00 + C5 δ00 + C6 H δ00

1 2 5 3 7 2 9 3 7 4 4 4 5 5 = C3 ( δ00 + δ00 − δ00 − 2hδ00 ) + C5 ( δ00 − δ00 − hδ00 + δ00 − 2hδ00 ) 2 3 8 4 2 = q˜2 dH + d Q˜ 2 , where y 7 y q˜2 = (−C3 − C5 ) − C5 4 4 (1 + x) (1 + x)5

(3.30)

˜ 2 is a continuous function in x, y. Then we have q2 = q¯2 + q˜2 , where q¯2 and q˜2 are given and Q by (3.18) and (3.30). An application to Lemma 2.2 gives to    M4 = q2 ω = q2 (q0 dH + dQ0 ) = q2 dQ0 . H (x, y)=h

H (x, y)=h

H (x, y)=h

After a tedious computation, we obtain  M4 (h) =

k−3 8

H (x, y)=h k=4 j =1

k λk0j ω0j +

k−3 7

k μk0j δ0j ,

(3.31)

k=3 j =0

where λk0j and μk0j are given in the supplementary material.1 k , ωk , δ k and δ k can be decomposed into hk dH + dQk . According to Lemma 2.3, ω02 ij ij 04 01 03 k , ωk , ωk , δ k and δ k can be decomposed into generators. For studying the zeros Similarly ω01 03 02 05 00 of M4 (h), we only need to give the decompositions which can be decomposed into generators. k , ωk , ωk , δ k and δ k , hence here The formula (3.23) has given some decompositions about ω01 03 02 05 00 we only give the rest decompositions as follows, 1 6 y 7 ω01 = δ00 + d(− ), 6 6(1 + x)6 1 7 y 8 ω01 = δ00 + d(− ), 7 7(1 + x)7 1 y3 7 6 6 5 4 ω03 = (2H δ00 − δ00 + 2δ00 − δ00 ) + d(− ), 2 6(1 + x)6 3 y3 8 7 7 6 5 = (2H δ00 − δ00 + 2δ00 − δ00 ) + d(− ), ω03 7 7(1 + x)7

(3.32)

1 The supplementary material can be found at the URL: http://www.math.sjtu.edu.cn/upload/teachers/8137/Melnikov. pdf.

JID:YJDEQ AID:9972 /FLA

[m1+; v1.304; Prn:10/09/2019; 10:19] P.28 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

28

5 y5 8 7 6 5 4 3 ω05 = ((2H − 1)2 δ00 + (8H − 4)δ00 + (6 − 4H )δ00 − 4δ00 + δ00 ) + d(− ), 7 7(1 + x)7 6 6 6 5 4 = 2H δ00 − δ00 + 2δ00 − δ00 , δ02 7 7 7 6 5 = 2H δ00 − δ00 + 2δ00 − δ00 . δ02

Combining (3.23), (3.31) and (3.32), we deduce the fourth Melnikov function through Lemma 2.3 M4 (h) = D1 J3 + D2 J4 + D3 J5 + D4 hJ5 + D5 J6 + D6 hJ6 + D7 J7 + D8 hJ7 + D9 h2 J7 , (3.33) where Di , i = 1, · · · , 9 are given in the supplementary material2 and Jk , k = 3, · · · , 7 are defined by (2.1). Next we consider case 4b ∼ 4d together and give the following theorem. Theorem 3.7. Assume M1 (h) = M2 (h) = M3 (h) ≡ 0, then (i) in the case 4b, M4 (h) has at most one zero, (ii) in the case 4c, M4 (h) has at most two zeros, (iii) in the case 4d, M4 (h) has at most three zeros. In short, M4 (h) has at most three zeros, namely, system (1.2) has at most three limit cycles by the Melnikov function of fourth order and the upper bound for the number of limit cycles is reached. Moreover, in the three cases, M4 (h) ≡ 0. Proof. Substituting (3.27) and (2.11) from Lemma 2.4 into (3.33) gives to π 1 b02 b11 M4 (h) = − √ (2h − 1)2 h2 [− (dˆ1 h + dˆ0 )], 13 12 (2a20 − b11 )2 ( 1 − 2h) where 3 2 2 2 2 2 2 2 b11 − 45a20 b02 − 6a20 b11 + 33a20 b02 b11 − 6b02 b11 , dˆ0 = 12a20 4 3 2 2 2 2 2 dˆ1 = 840a20 − 1424a20 b11 + 90a20 b02 + 878a20 b11 − 66a20 b02 b11 3 2 2 4 − 236a20 b11 + 12b02 b11 + 24b11 .

So M4 (h) has at most one zero in h ∈ (0, 12 ) under the case 4b. If M4 (h) ≡ 0 for the case 4b, we get dˆ0 = dˆ1 . By solving the equation, we obtain two solutions {a20 = 12 b11 , b02 = 0, b11 = b11 } and {a20 = 0, b02 = b02 , b11 = 0}. For case 4b, b11 = 2a20 and b11 = 0, this is a contradiction to the above solutions. Thus in case 4b, M4(h) ≡ 0. In the case 4c, substituting (3.28) and (2.11) from Lemma 2.4 into (3.33) gives to π 1 b11 ˜ 2 M4 (h) = − √ (2h − 1)h2 [− (d2 h + d˜1 h + d˜0 )], 13 1875 b02 ( 1 − 2h) 2 The supplementary material can be found at the URL: http://www.math.sjtu.edu.cn/upload/teachers/8137/Melnikov. pdf.

JID:YJDEQ AID:9972 /FLA

[m1+; v1.304; Prn:10/09/2019; 10:19] P.29 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

29

where 2 2 3 2 4 2 2 4 b02 + 150a01 b02 − 120a01 b02 b11 − 225b02 + 360b02 b11 − 144b11 , d˜0 = 600a01 2 2 3 2 4 2 2 4 b02 − 1850a01 b02 − 2040a01 b02 b11 + 775b02 − 1480b02 b11 + 1008b11 , d˜1 = 600a01 3 2 4 2 2 4 − 480a01 b02 b11 + 250b02 + 4040b02 b11 + 288b11 . d˜2 = −500a01 b02

So M4 (h) has at most two zeros in h ∈ (0, 12 ) in the case 4c, namely, system (1.2) has at most two limit cycles. In this way, if M4 (h) ≡ 0, that is, d˜2 = d˜1 = d˜0 = 0, whose solutions are {b02 = 0, b11 = 0, a01 = a01 }

(3.34)

1 {a01 = b02 , b02 = b02 , b11 = 0}. 2

(3.35)

and

Here (3.34) and (3.35) are contradictory with b11 = 0, hence M4 ≡ 0 for case 4c. Finally, we consider case 4d. In this case, substituting (3.29) and (2.11) from Lemma 2.4 into (3.33) gives to π 1 b11 b02 M4 (h) = − √ h(2h − 1)[− (d¯3 h3 + d¯2 h2 + d¯1 h + d¯0 )], (3.36) 1500 (2a20 − b11 )2 ( 1 − 2h)13 where 4 3 2 2 3 4 + 30600a20 b11 − 19440a20 b11 + 5472a20 b11 − 576b11 , d¯0 = −18000a20 4 3 2 2 2 2 − 183600a20 b11 + 5625a20 b02 + 115920a20 b11 d¯1 = 108000a20 2 3 2 2 4 − 4950a20 b02 b11 − 32184a20 b11 + 1080b02 b11 + 3312b11 , 4 3 2 2 2 2 + 367200a20 b11 − 22500a20 b02 − 228680a20 b11 d¯2 = −216000a20 2 3 2 2 4 + 17875a20 b02 b11 + 62036a20 b11 − 3470b02 b11 − 6248b11 , 4 3 2 2 2 2 − 244800a20 b11 + 22500a20 b02 + 149200a20 b11 d¯3 = 144000a20 2 3 2 2 4 − 15950a20 b02 b11 − 39400a20 b11 + 2800b02 b11 + 4000b11 .

In this case, M4 (h) has at most three zeros in h ∈ (0, 12 ), which implies that system (1.2) has at most three limit cycles. In the following, we consider how to find the three zeros of (3.36). As we know, in order to find the zeros of (3.36), we only need to find the zeros of the following equation, m4 (h) = d¯3 h3 + d¯2 h2 + d¯1 h + d¯0 .

(3.37)

The coefficients d¯0 , d¯1 , d¯2 , d¯3 depend on these variables (a20 , b02 , b11 ). It is unlucky that the ¯ ¯ ¯ Jacobian determinant of ( dd¯0 , dd¯1 , dd¯2 ) with respect to the variables (a20 , b02 , b11 ) is equal to 0. 3

3

3

JID:YJDEQ AID:9972 /FLA

30

[m1+; v1.304; Prn:10/09/2019; 10:19] P.30 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

However we claim that (3.37) has still three real zeros. Then we assume h = 0 is a double root of m4 (h) = 0, thus it should satisfy m4 (h) = m 4 (h) = 0. It is easy to know that if a20 = 25 b11 , m4 (h) = m 4 (h) = 0. In this condition, we have m

4 (0) =

32 2 b (5b02 − 4b11 )(5b02 + 4b11 ) 5 11

and 1 18 4 45 2 2 m4 ( ) = b02 b11 + b11 > 0. 2 2 5 We want to find two positive zeros near h = 0 and a transversal zero in h ∈ (0, 12 ) by perturbing a20 , b02 , b11 . Assume m

4 (0) < 0, namely, (5b02 − 4b11 )(5b02 + 4b11 ) < 0, and let b11 = 1, a20 = 25 , b02 = 35 . In this condition, h = 0 is a double root, and there is a transversal zero in h ∈ (0, 12 ). Perturbing a20 , b02 , b11 as follows, 2 3 a20 = b11 − , b02 = b11 , b11 = 1 + . 5 5 Assume h() is one zero of m4 (h, ) = 0 near h = 0, by expanding h() at h = 0, we have h() = n1  + n2  2 + n3  3 + o( 3 ).

(3.38)

Substituting (3.38) into (3.37) gives to m4 (h, ) = f1 (n1 ) 2 + f2 (n1 , n2 ) 3 + f3 (n1 , n2 , n3 ) 4 + o( 4 ), where 2 f1 (n1 ) = − n1 (−225 + 56n1 ), 5 483 2 224 692 3 f2 (n1 , n2 ) = −1800 + 1575n1 − n + 90n2 − n 1 n2 + n , 5 1 5 5 1 18277 2 966 2478 3 f3 (n1 , n2 , n3 ) = −19800 + 13680n1 − n1 + 1575n2 − n1 n2 + n 5 5 5 1 224 112 2 2076 2 + 90n3 − n1 n3 − n + n n2 . 5 5 2 5 1 By m4 (h, ) = 0, we have f1 (n1 ) = 0, f2 (n1 , n2 ) = 0, f3 (n1 , n2 , n3 ) = 0. We obtain two solutions as follows by solving the three equations

JID:YJDEQ AID:9972 /FLA

[m1+; v1.304; Prn:10/09/2019; 10:19] P.31 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

n1 = 0, n2 = 20, n3 = −

31

274 9

and n1 =

225 11654515 131084762279 , n2 = , n3 = . 56 87808 19361664

Then we get two zeros h1 = 20 2 + o( 2 ),

h2 =

225  + o(). 56

When  is small enough, the sighs of h1 and h2 are decided by the lower degree of , hence h1 1 and h2 are both positive. Thus there are three zeros of m(h) = 0 in h ∈ (0, 12 ). Let  = 1000 , that 1997 3003 1001 is a20 = 5000 , b02 = 5000 , b11 = 1000 , we have three zeros as follows: h1 = 0.00001996717376, h2 = 0.004157828578, h3 = 0.1578135609. If M4 (h) ≡ 0, namely, d¯0 = d¯1 = d¯2 = d¯3 = 0. By solving the equations, we obtain two solutions {a20 = 0, b11 = 0} and {a20 = a20 , b11 = 2a20 , b02 = 0}. In case 4d, we know b11 = 0 and b11 = 2a20 , this leads to a contradiction. Hence M4 (h) ≡ 0 in the case 4d. Thus we have completed the proof of Theorem 3.7. 2 4. The proofs of main theorems In this section, we will give the proofs of the main results: The proof of Theorem 1.1. Theorem 1.1 can be proved by combining Theorem 3.1, Theorem 3.2, Theorem 3.5 and Theorem 3.7. 2 The proof of Theorem 1.2. Theorem 1.2 can be demonstrated according to Theorem 3.3, Theorem 3.4 and Theorem 3.6. 2 Acknowledgments The corresponding author is supported by the state key program of NSF of China grant number 11431008 and 11771282, NSF of Shanghai grant number 15ZR1423700. References [1] A. Atabaigi, N. Nyamoradi, H.R.Z. Zangeneh, The number of limit cycles of a quintic polynomial system, Comput. Math. Appl. 57 (4) (2009) 677–684. [2] A. Atabaigi, N. Nyamoradi, H.R.Z. Zangeneh, The number of limit cycles of a quintic polynomial system with center, Nonlinear Anal. 71 (7–8) (2009) 3008–3017. [3] A. Buic˘a, A. Gasull, J. Yang, The third order Melnikov function of a quadratic center under quadratic perturbations, J. Math. Anal. Appl. 331 (1) (2007) 443–454. [4] A. Buic˘a, J. Llibre, Limit cycles of a perturbed cubic polynomial differential center, Chaos Solitons Fractals 32 (3) (2007) 1059–1069.

JID:YJDEQ AID:9972 /FLA

32

[m1+; v1.304; Prn:10/09/2019; 10:19] P.32 (1-32)

P. Yang, J. Yu / J. Differential Equations ••• (••••) •••–•••

[5] B. Coll, J. Llibre, R. Prohens, Limit cycles bifurcating from a perturbed quartic center, Chaos Solitons Fractals 44 (4–5) (2011) 317–334. [6] J.P. Francoise, Successive derivatives of a first return map, application to the study of quadratic vector fields, Ergod. Theory Dyn. Syst. 16 (1) (1996) 87–96. [7] A. Gasull, J.T. Lázaro, J. Torregrosa, Upper bounds for the number of zeroes for some Abelian integrals, Nonlinear Anal. 75 (13) (2012) 5169–5179. [8] A. Gasull, C. Li, J. Torregrosa, Limit cycles appearing from the perturbation of a system with a multiple line of critical points, Nonlinear Anal. 75 (1) (2012) 278–285. [9] A. Gasull, R. Prohens, J. Torregrosa, Bifurcation of limit cycles from a polynomial non-global center, J. Dyn. Differ. Equ. 20 (4) (2008) 945–960. [10] L. Gavrilov, I.D. Iliev, Cubic perturbations of elliptic Hamiltonian vector fields of degree three, J. Differ. Equ. 260 (5) (2016) 3963–3990. [11] H. Giacomini, J. Llibre, M. Viano, On the nonexistence, existence and uniqueness of limit cycles, Nonlinearity 9 (2) (1996) 501–516. [12] J. Giné, J. Llibre, Limit cycles of cubic polynomial vector fields via the averaging theory, Nonlinear Anal. 66 (8) (2007) 1707–1721. [13] I.D. Iliev, On second order bifurcations of limit cycles, J. Lond. Math. Soc. (2) 58 (2) (1998) 353–366. [14] Q. Jiang, M. Han, Melnikov functions and perturbation of a planar Hamiltonian system, Chin. Ann. Math., Ser. B 20 (2) (1999) 233–246. [15] B.Y. Li, Z.F. Zhang, A note on a result of G.S. Petrov about the weakened 16th Hilbert problem, J. Math. Anal. Appl. 190 (2) (1995) 489–516. [16] S. Li, Y. Zhao, J. Li, On the number of limit cycles of a perturbed cubic polynomial differential center, J. Math. Anal. Appl. 404 (2) (2013) 212–220. [17] J. Llibre, J.S. Pérez del Río, J.A. Rodríguez, Averaging analysis of a perturbated quadratic center, Nonlinear Anal., Theory Methods Appl. 46 (1) (2001) 45–51. [18] J. Llibre, H. Wu, J. Yu, Linear estimate for the number of limit cycles of a perturbed cubic polynomial differential system, Nonlinear Anal. 70 (1) (2009) 419–432. [19] J. Llibre, J. Yu, On the upper bound of the number of limit cycles obtained by the second order averaging method, Dyn. Contin. Discrete Impuls. Syst., Ser. B, Appl. Algorithms 14 (6) (2007) 841–873. [20] G. Tigan, Using Melnikov functions of any order for studying limit cycles, J. Math. Anal. Appl. 448 (1) (2017) 409–420. [21] G. Xiang, M. Han, Global bifurcation of limit cycles in a family of multiparameter system, Int. J. Bifurc. Chaos Appl. Sci. Eng. 14 (9) (2004) 3325–3335. [22] G. Xiang, M. Han, Global bifurcation of limit cycles in a family of polynomial systems, J. Math. Anal. Appl. 295 (2) (2004) 633–644. [23] Y. Xiong, The number of limit cycles in perturbations of polynomial systems with multiple circles of critical points, J. Math. Anal. Appl. 440 (1) (2016) 220–239.