Limit cycles for quadratic and cubic planar differential equations under polynomial perturbations of small degree

In this paper we consider planar systems of differential equations of the form \begin{document}$\left\{ \begin{array}{lcl} \dot x&=&-y+\delta p(x,y)+\varepsilon P_n(x,y),\\ \dot y&=x'>where \begin{document}$δ, \varepsilon$\end{document} are small parameters, $(p, q)$ are quadratic or cubic homogeneous polynomials such that the unperturbed system ($\varepsilon=0$) has an isochronous center at the origin and $P_n, Q_n$ are arbitrary perturbations. Estimates for the maximum number of limit cycles are provided and these estimatives are sharp for $n≤q 6$ (when $p, q$ are quadratic). When $p, q$ are cubic polynomials and $P_n, Q_n$ are linear, the problem is addressed from a numerical viewpoint and we also study the existence of limit cycles.

1. Introduction. In this paper we consider the following 2-parameters planar polynomial differential system ẋ = −y + δp(x, y) + εP n (x, y), y = x + δq(x, y) + εQ n (x, y), where (p, q) are quadratic or cubic homogeneous polynomials and P n , Q n are ndegree polynomials. The parameter ε > 0 is small, and δ ∈ R. We suppose also that (1) for ε = 0 (and any value of δ ∈ R) has an isochronous center at the origin. Our aim is to study the number of limit cycles that bifurcates from the linear and non-linear centers of (1), considering polynomials perturbations. We also compare the maximum number of limit cycles that can bifurcate from linear and non-linear centers, when perturbed with a fixed perturbation.
The number of limit cycles in systems that are perturbations of a linear nondenegerate center has been extensively studied in [3,18,17]. This problem is considered in a large number of papers (see [3] and references), including higher dimensions analogues. The averaging theory has been the main technique to attack this problem, especially in the last years.
The existence of limit cycles in perturbations of non-linear centers is also a very active research topic, and the main tools employed are the Abelian integrals [8], Melnikov functions [11] and averaging theory [21]. In [3] is proved that if we consider the special case δ = 1 in (1) and a quadratic perturbation then at most 2 limit cycles bifurcate from this center. The same result was proved in [6]. A computationalnumeric approach to this problem (using general perturbations up to some degree) can be found in [26,13]. For cubic centers in particular form, Llibre and Giné studied the number of limit cycles in [9].
Worth to mention that the interest on this kind of problem (counting the number of limit cycles) is mainly motivated by the Hilbert 16th problem [12].
To work with arbitrary perturbations of non-linear centers is a very hard task. In [13] authors consider particular perturbations up to some degree, and in [20] authors consider a very specific family of polynomial perturbations without restrictions on the degree.
The choice of the polynomials p(x, y) and q(x, y) is justified by the classification of isochronous planar centers, due to Loud in the quadratic case and to Pleshkan in the cubic case.
In [19] is proved that the origin is a isochronous center of a quadratic system if and only if the system can be brought to one of the systems (2), (3), (4) or (5) through a linear change of coordinates and a rescaling of time: ẋ = −y + 16 3 x 2 − 4 3 y 2 y = x(1 + 8 3 y) All these systems are integrable. In this paper we consider variations of systems (2) and (3), that we highlight now: ẋ = −y + δ(x 2 − y 2 ) y = x(1 + 2δy) (S1) ẋ = −y + δx 2 y = x(1 + δy) (S2) The main results we prove in this paper are the following.
Theorem 1.1. Regarding systems (S1) and (S2). a) If we consider an n-degree polynomial perturbation of system (S1), ẋ = −y + δ(x 2 − y 2 ) + εP n (x, y), y = x(1 + 2δy) + εQ n (x, y), where δ = 0, ε > 0 is small, P n , Q n are n-degree polynomials, then the maximum number of limit cycles bifurcating from the nonlinear system (S1p) which are detected by first order averaging is given by Table 1 and it is achieved. b) If we consider an n-degree polynomial perturbation of system (S2), ẋ = −y + δx 2 + εP n (x, y), y = x(1 + δy) + εQ n (x, y), where δ = 0, ε > 0 is small, P n , Q n are n-degree polynomials, then the maximum number of limit cycles bifurcating from the nonlinear system (S2p) which are detected by first order averaging is given by Table 2 and it is achieved.
Based on Tables 1 and 2, we present the following conjectures.
Conjecture 1: Consider system (S1p) for δ = 0, ε > 0 small and P n , Q n ndegree polynomials. Then the maximum number of limit cycles bifurcating from the nonlinear center (S1p) (via first order averaging) is n − 2 for n ≥ 3.
Conjecture 2: Consider system (S2p) for δ = 0, ε > 0 small and P n , Q n ndegree polynomials. Then the maximum number of limit cycles bifurcating from the nonlinear center (S2p) (via first order averaging) is given by n for n ≥ 2.
We remark that the case δ = 0 of Theorem 1.1 is a classical result discussed in [7], [15] and references therein. Theorem 1.1 can also be proved using the techniques presented in [10].
We also notice that the case δ = 1 with P n , Q n quadratic is studied in [16]. In Theorem 4 of [14], the authors obtained an upper bound of limit cycles when P n and Q n are cubic, but they didn't provide the upper bound for a perturbation of any degree.
The main objective of Theorem 1.1 is to provide the maximum number of limit cycles (by first order averaging) for a perturbation of arbitrary degree, and show that it seems to follow a certain rule with respect to the perturbation degree, and it originates (and justify) Conjectures 1 and 2.
The cubic case is studied in [22] and a classification is obtained. In particular, is proved that the origin is a isochronous center of a cubic system if and only if the system can be brought to one of the systems (6), (7), (8), (9) through a linear change of coordinates and a rescaling of time: Again, all of these systems are integrable.
Using the averaging theory, we relate the existence of limit cycles to the existence of simple zeroes of a function that involves elliptic integrals, so it is not possible to solve analytically, and we employ some numerical methods. Theorems 1.2 and 1.3 can also be proved using the techniques presented in [10]. In the cubic case, we just consider linear perturbations. Theorem 1.2. We consider the 2-parameter planar polynomial system that is a linear perturbation of (6). If δ = 0 then there is no limit cycle. If δ = 0 we present numerical evidences that there is a set Ω ⊂ [0, 1] 4 in the 4-dimensional space of the parameters a 1,0 , a 0,1 , b 1,0 , b 0,1 such that system (T1p) admits a limit cycle if the parameters are in Ω.
Remark 1. As the result above is numerical, the set Ω is discrete. We conjecture that: i) it is possible to take Ω as an convex set and ii) the reciprocal is also true, that is, system (T1p) has a limit cycle only if the parameters are in Ω, but our method does not allow to conclude this. Theorem 1.3. We consider the 2-parameter planar polynomial system that is a linear perturbation of (7). If δ = 0 then there is no limit cycle. If δ = 0 there is at most one limit cycle via first order averaging, that shrinks to the origin when δ → 0 + . This paper is organized as follows. In Section 2 we provide classifications for quadratic and cubic isochronous integrable centers. In Section 3 we describe the results of the averaging theory that we shall need to our purposes. In Section 4 we consider the case of quadratic centers and prove Theorem 1.1. In Section 5 we consider the case of cubic centers and prove Theorem 1.2. In Section 6 we provide some examples, including numerics.
2. Quadratic and cubic centers. In this section we present some classical classification results on planar isochronous centers. More details can be found on [5] and references therein.
Now we list some useful informations about the quadratic systems (2)- (5) and also about the cubic systems (6)-(9), including their expressions in (generalized) polar coordinates. For details, see [2].

System
A strong first integral for this system is H(x, y) = x 2 + y 2 (1 + y) 2 and a reciprocal integrating factor is V = (1 + y) 3 .
A strong first integral for this system is H(x, y) = x 2 + y 2 1 + 2xy and a reciprocal integrating factor is V = (1 + 2xy) 2 .
A strong first integral for this system is H(x, y) = 3. Averaging theory. Here we briefly describe some results on periodic averaging of first order. This is the simplest form of averaging, and is concerned with approximating solutions of a non-autonomous differential equation by solutions of a autonomous one. In particular, the first order averaging method is equivalent to the study of the first order Melnikov function (both are equivalent to the study of the displacement function). For more references on this, see [11]. The next theorem is the classical averaging theorem for periodic differential system.
Systems in the form of (10) are said to be in the standard form. We remark that the regularity condition that f, g are of class C 2 is not really necessary, but simplify the statement. For a proof and more comments, see chapter 6 of [21]. Theorem 3.1 is about the existence of periodic solutions for non-autonomous systems. The following construction allows us to use this theorem for proving the existence of limit cycles. Consider the planar system where p, q, P, Q : R 2 → R are continuous functions. Let us first consider the case p(x, y) = −y and q(x, y) = x. Then, for ε = 0, system (12) is the linear center.
By Theorem 3.1, each simple zero of G is associated with a periodic solution of (15), and then, with a limit cycle of (12). Now we turn to the general case, where p(x, y) and q(x, y) are general polynomials. Our aim is to put (12) in the standard form to apply the Averaging Theorem 3.1.
From now, we assume that (12), for ε = 0, has a first integral H, and a continuous family of closed orbits This allow us to find a change of coordinates that put (12) in the standard form. Assume that xq(x, y) − yp(x, y) = 0 for all (x, y) in Γ h . Let be a continuous function such that and all φ ∈ [0, 2π). The change of coordinates x = ρ(R, φ) cos(φ), y = ρ(R, φ) sin(φ) applied to system (12) give us Ṙ = εL(R, φ), for some smooth functions L, S. Now, using the same argument as in the linear center, we can obtain a non-autonomous system in the form where the prime is the derivative with respect to φ. Note that (17) is in the standard form.
Applying the Averaging Theorem 3.1 to (17) and writing the expression of L, we obtain the following result: , Theorem 5.2). Consider that system (12) has a first integral H for ε = 0, and a continuous family of closed orbits Let µ(x, y) be an integrating factor for system (12). If where µ, P, p, Q, q depends on x = ρ(R, φ) cos(φ) and y = ρ(R, φ) sin(φ), then each simple zero of F (R) give us a limit cycle of (12).
Remark 2. The integrand of F (R) is exactly the function L(R, φ) in the first line of system (16).
4. Proof of Theorem 1.1. Here we consider the following 2-parameter systems where P n , Q n are n-degree polynomials. The parameter ε > 0 is small, and δ ∈ R.
To study the number of limit cycles in these systems for ε > 0 small, by Theorem 3.2, we have to determine the isolated zeroes of the function 18. We divide the proof in the following cases: 4.1. Non-linear quadratic center (S1) (δ = 0). Consider system (S1p) with is a first integral for this system.
In order to apply Theorem 3.2, we must find ρ such that: Solving this equation for this system, we obtain: Notice that we must assume (δR) 2 ≤ 1 to get ρ well-defined. Therefore, the domain of R is (−1/|δ|, 1/|δ|) and we have that h 2 = 1/δ 2 .
Since h 1 = 0, it follows that the period annulus occurs when R ∈ 0, 1 |δ| . The integrating factor associated to the first integral H 1 δ is given by: Now, we use these informations to compute the Melnikov function F without particularize the degree n of the polynomials P n , Q n , therefore we must fix n between 1 and 7.
As we mention, the period annulus occurs when R ∈ 0, 1 |δ| , so we seek zeroes of the Melnikov functions in this domain.
For simplicity, we introduce a new variable Z given by R = After solving the corresponding integrals, we obtain the expressions for the functions F 1 (Z), where ρ n 1 is a C ∞ function without real zeroes and f  Table 3.  Table 3. Degree of the polynomial part of F  To conclude the proof, now we present an example for every degree, proving that the maximum number of limit cycles can be achieved.
Consider the family of differential equations where P n , Q n are n-degree polynomials, for 1 ≤ n ≤ 7.
has only one zero in (0,1), so the system has one limit cycle.
and Q 4 (x, y) = x 4 + y 4 + y, then f has two zeroes in (0,1), so the system has two limit cycles.
In Appendix C we provide general expressions on the coefficients of f [n] 1 for 1 ≤ n ≤ 7, to obtain the maximum number of limit cycles.
In Table 4 we summarize the relations between the degree of the perturbation, the degree of the polynomial part of F  Table 4. Degree of the perturbation in (19), degree of the polynomial part of F 4.2. Non-linear quadratic center (S2) (δ = 0). Now, consider system (S2p) with δ = 0, P n (x, y) = n i+j=1 a i,j x i y j and Q n (x, y) = n i+j=1 b i,j x i y j . Notice that, for ε = 0, the system is integrable for each δ ∈ R.
In fact H 2 δ (x, y) = x 2 + y 2 (1 + δy) 2 is a first integral for this system. In order to apply Theorem 3.2, we must find ρ such that: Solving this equation for this system, we obtain: .
Since P n (x, y) = n i+j=1 a i,j x i y j and Q n (x, y) = n i+j=1 b i,j x i y j , (20) writes as: (21) Also, it is not possible to obtain the expression of F without particularize the degree n of the polynomials P n , Q n , therefore we must fix n between 1 and 6.
As we mention, the period annulus occurs when R ∈ 0, 1 |δ| , so we seek zeroes of the Melnikov functions in this domain.
For simplicity, we introduce a new variable Z given by R = √ 1−Z 2 |δ| . We observe that R ∈ 0 , 1 |δ| is a zero of the Melnikov function F 1 , for n = 1, . . . , 6 (see Appendix B). As in the last, we can write F 2 (Z), where ρ n 2 is a C ∞ function without real zeroes in (0, 1) and f   Table 5.
Again, to conclude the proof, we present an example of every degree, proving that the maximum number of limit cycles can be achieved also in this case.
To conclude the proof, now we present an example for every degree, proving that the maximum number of limit cycles can be achieved.
Consider the family of differential equations ẋ = −y + δx 2 + εP n (x, y), where P n , Q n are n-degree polynomials, for 1 ≤ n ≤ 6. For n = 1, the only f [1] 1 (Z) has no non-zero roots, so the system has no limit cycles.
has three zeroes in (0,1), so the system has three limit cycles. has four zeroes in (0,1), so the system has four limit cycles. has five zeroes in (0,1), so the system has five limit cycles.
has six zeroes in (0,1), so the system has six limit cycles. In Appendix D we provide general expressions on the coefficients of f 2 , for 1 ≤ n ≤ 6, to obtain the maximum number of limit cycles.
In Table 6 we summarize the relations between the degree of the perturbation, the degree of the polynomial part of F  Table 6. Degree of the perturbation in (25), degree of the polynomial part of F and note that H 3 δ (x, y) = x 2 + y 2 1 + 2δxy is a first integral for this system, with integrating factor µ = 2/(1 + 2δxy) 2 . To study the number of limit cycles in this system for ε > 0 small and δ = 0, by Theorem 3.2, we have to construct function (18) and determine its simple zeroes.
In this case, (18) is given by Considering δ > 0, note that F δ (R) = 0 has a positive and simple solution if and only if a 0,1 + b 1,0 = 0, and the solution is given by 0 + a 0,1 ) .

RICARDO M. MARTINS AND OTÁVIO M. L. GOMIDE
Then we have at most one limit cycle, thus proving Theorem 1.3.
Corollary 1. A differential system with the form with ε > 0 small, δ > 0 and a i,j , b i,j ∈ R has a limit cycle if and only if a 0,1 + b 1,0 = 0.
Now we proceed to the proof of Theorem 1.2. Recall system (T1p) and note that H 4 δ (x, y) = (x 2 + y 2 ) 2 1 + 4δxy is a first integral for this system, with integrating factor µ = 4(x 2 + y 2 )/(1 + 4δxy) 2 . Again, we construct function (18) to this case and determine its simple zeroes. In the case of system (T1p), the function (18) is given by where D(R, ψ), A i,j (R, ψ), B i,j (R, ψ) are given in Appendix E. This integral can be rewritten with respect to elliptic integrals of the first and second kinds: where where EllipticF(z, k),EllipticK(k) are the incomplete and complete integrals of first kind and EllipticE(z, k) is the (complete and incomplete) elliptic integral of second kind, given by More details on these functions can be found on [25,23]. It is not possible to obtain the real zeros of F δ (R) by using algebraic or analytic methods, due to the appearance of the elliptic integrals. Our strategy is to employ a numerical procedure to obtain the conclusions of Theorem 1.2. We denote by Ω ⊂ M the subset of the admissible vertices. If m ∈ Ω then the equation F m δ (R) = 0 has a positive solution, so there is a limit cycle for system (T1p), with the parameters (a 1,0 , a 0,1 , b 1,0 , b 0,1 ) = m.

Numerical strategy. We consider a partition
We note that in all the other cases considered in this paper, we have open conditions defining the subsets of the parameter space where there are limit cycles. We believe it is reasonable to suppose that this also holds for the current case. Then we consider that there is a set Ω ⊇ Ω open and represented by the convex hull of Ω such that for every m ∈ Ω, system (T1p) with the parameters (a 1,0 , a 0,1 , b 1,0 , b 0,1 ) = m has a limit cycle.
We applied this strategy for n = 10. In this case the set M has 11 4 = 14.641 points. Using Maple 14 TM to numerically determine if each F m δ (R) has a zero, we found that Ω consists of 8.960 admissible vertices.
Using the software Octave we calculated the convex hull H of Ω and obtained the equations of 21 limiting hyperplanes for Ω (consider Ω ⊂ [0, 1] 4 = {(x, y, z, w); 0 ≤ x, y, z, w ≤ 1}, for δ = 1 (the result is the same for every δ ∈ (0, 1]):     This complete the proof of Theorem 1.3. We summarize the results in the following corollary.
Despite our numerical method, and based on the numerical results we found, we present the follow conjecture. Conjecture 1. A differential system of the form with ε > 0 small, δ > 0 and a i,j , b i,j ∈ [0, 1] 4 has a limit cycle if and only if (a 1,0 , a 0,1 , b 1,0 , b 0,1 ) ∈ Ω, where Ω ⊂ [0, 1] 4 is a convex set that can be approximated by taking the convex hull of Ω for a large grid.
Appendix A. Full expressions of functions F are given by: