Full Terms & Conditions of access and use can be found at https://www.tandfonline.com/action/journalInformation?journalCode=gcom20 International Journal of Computer Mathematics ISSN: 0020-7160 (Print) 1029-0265 (Online) Journal homepage: https://www.tandfonline.com/loi/gcom20 Sufficient conditions for the existence of periodic solutions of the extended Duffing–Van der Pol oscillator Rodrigo D. Euzébio & Jaume Llibre To cite this article: Rodrigo D. Euzébio & Jaume Llibre (2016) Sufficient conditions for the existence of periodic solutions of the extended Duffing–Van der Pol oscillator, International Journal of Computer Mathematics, 93:8, 1358-1382, DOI: 10.1080/00207160.2015.1046847 To link to this article: https://doi.org/10.1080/00207160.2015.1046847 Accepted author version posted online: 30 Apr 2015. Published online: 09 Jun 2015. Submit your article to this journal Article views: 90 View Crossmark data https://www.tandfonline.com/action/journalInformation?journalCode=gcom20 https://www.tandfonline.com/loi/gcom20 https://www.tandfonline.com/action/showCitFormats?doi=10.1080/00207160.2015.1046847 https://doi.org/10.1080/00207160.2015.1046847 https://www.tandfonline.com/action/authorSubmission?journalCode=gcom20&show=instructions https://www.tandfonline.com/action/authorSubmission?journalCode=gcom20&show=instructions http://crossmark.crossref.org/dialog/?doi=10.1080/00207160.2015.1046847&domain=pdf&date_stamp=2015-04-30 http://crossmark.crossref.org/dialog/?doi=10.1080/00207160.2015.1046847&domain=pdf&date_stamp=2015-04-30 International Journal of Computer Mathematics, 2016 Vol. 93, No. 8, 1358–1382, http://dx.doi.org/10.1080/00207160.2015.1046847 Sufficient conditions for the existence of periodic solutions of the extended Duffing–Van der Pol oscillator Rodrigo D. Euzébioa,b and Jaume Llibrea∗ aDepartament de Matemàtiques, Universitat Autònoma de Barcelona, 08193 Bellaterra, Barcelona, Catalonia, Spain; bDepartament de Matemática, IBILCE, UNESP, Rua Cristovao Colombo 2265, Jardim Nazareth, CEP 15.054–00, Sao José de Rio Preto, SP, Brazil (Received 20 October 2014; revised version received 17 March 2015; accepted 28 March 2015) In this paper, some aspects on the periodic solutions of the extended Duffing–Van der Pol oscillator are discussed. Doing different rescaling of the variables and parameters of the system associated with the extended Duffing–Van der Pol oscillator, we show that it can bifurcate one or three periodic solutions from a two-dimensional manifold filled by periodic solutions of the referred system. For each rescaling we exhibit concrete values for which these bounds are reached. Beyond that we characterize the stability of some periodic solutions. Our approach is analytical and the results are obtained using the averaging theory and some algebraic techniques. Keywords: extended Duffing–Van der Pol oscillator; periodic solution; non-autonomous systems; averaging theory 2010 AMS Subject Classifications: Primary: 34C07; 34C15; 34C25; 34C29; 37C60 1. Introduction 1.1 Setting the problem A large number of non-autonomous chaotic phenomena in physics, engineering, mechanics and biology, among others, are described by second-order differential systems of the form ẍ = g(x, ẋ, t) + γ (t), (1) where g(x, ẋ, t) is a continuous function and γ (t) is some external force. For instance, in biology, system (1) models the FHN neuron oscillator, and in engineering, system (1) is a model to the horizontal platform system. The specific topic addressed in this paper concerns another particular case of Equation (1), namely, an extension of the forced Van der Pol equation with external excitation. Van der Pol’s system also plays an important role in many applications in areas such as engineering, biology, physics and seismology (see [2] and references therein). *Corresponding author. Email: jllibre@mat.uab.cat © 2015 Taylor & Francis mailto:jllibre@mat.uab.cat International Journal of Computer Mathematics 1359 The system associated with the forced Van der Pol equation with external excitation is characterized by system (1) with functions g and γ in the form g(x, ẋ, t) = ρ0(1 − x2)ẋ − d dx V(x), γ (t) = δ0 cos(ωt), (2) where ρ0 is the damping parameter, V(x) is the potential function, and δ0 and ω are the amplitude and angular frequency of the driving force γ (t), respectively. We assume that ρ0 is non-negative and δ0 and ω are positive. The potential V(x) can be approximated by a finite Taylor expansion in series V2(x) = 1 2ω2 0x2, V4(x) = 1 2ω2 0x2 + 1 4α0x4, V6(x) = 1 2ω2 0x2 + 1 4α0x4 + 1 6λ0x6, where ω0 and λ0 are non-zero and α0 is a real number. Almost all papers on forced excited Van der Pol systems deal with the potential V2. How- ever, some papers concerning the potential V4 have shown a lot of interesting behaviours [7,12,25,27,28]. This case is usually referred to as Duffing–Van der Pol oscillator or φ4-Van der Pol oscillator. Nevertheless, more recently, some papers were published taking into account the potential V6 meanly addressing the problem of chaos control [13,17,22]. The dynamics considering the potential V6 is more complex and rich than the corresponding cases consid- ering the potentials V2 and V4 [24,26]. This is the case that we will regard. This case is quoted in the literature as the extended Duffing–Van der Pol oscillator or φ6-Van der Pol oscillator. In this paper, we will give an analytical treatment to system (1) in order to study its periodic solutions considering functions (2) and the potential V6. Indeed, calling y = ẋ we obtain ẋ = y, ẏ = −ω2 0x + ρ0y − α0x3 − ρ0x2y − λ0x5 + δ0 cos(ωt). (3) System (3) becomes simpler if we perform a rescaling s = ωt in the time t and another one y = ω0Y in the spatial variable y. In fact, calling again the new time s by t and the variable Y by y, after the rescaling, we have ẋ = y, ẏ = −x + ρy − αx3 − ρx2y − λx5 + δ cos t, (4) where the new parameters ρ, α, λ and δ are the respective old ones divided by ω2 0. The study of the periodic solutions of the non-linear non-autonomous 2π -periodic differential system (4) will be the objective of this paper. When α, λ and δ are zero, Equation (4) is referred to as unforced Van der Pol equation and has a unique stable periodic solution for ρ positive. Furthermore, if ρ is large, this periodic solution remains, and it describes a periodic oscillatory behaviour called relaxation oscillation. On the other hand, by considering δ non-zero, system (4) can present non-linear attractors from simple periodic solutions up to multi-periodicity and chaos. In [10], we can find a didactical discussion about these objects and some results using classical techniques are obtained for a special case of system (4) using particular values of ρ, δ, ω and initial conditions. There exists an exhaustive list of papers in the literature studying the properties of system (4) when δ is zero. For δ positive, many open questions remain mainly due to the difficulty of integrating the system. In particular, in [19], Ma et al. study some aspects of robust practical synchronization for system (1) and apply the results to a particular case of system (4). In the same direction, in [14], Leung investigates synchronization processes between chaotic attractors 1360 R.D. Euzébio and J. Llibre considering α = λ = 0 in system (4). In [2,11], many aspects of system (4) for the case where α and λ are zero and ρ is large are studied. Furthermore, in [8], Egami and Hirano provide suffi- cient conditions to the existence of one periodic solution under some analytical hypotheses for a general forced Van der Pol system considering ρ equal to zero. Periodic solutions for system (4) were found in [17], where Liu and Yamaura investigate chaos control. Besides, in [5,6], results on the existence of periodic solutions for an autonomous special case of the extended Duffing–Van der Pol oscillator were obtained. In [16,30,31], we can also find some results on systems similar to system (4). In this paper, we are concerned with periodic solutions of system (4). We will present sufficient conditions in order that this system possesses one or three periodic solutions, and we will provide conditions on the parameters ρ, α, λ and δ for which these bounds are realizable. In addition, we will prove that it is not possible to obtain different bounds of periodic solutions using the methodology presented in this paper. The phase space of our non-autonomous differential system is (x, y, t) and its Poincaré map is defined in the (x, y)-space. When possible the stability of the periodic solutions will be studied. As usual a periodic solution is stable if the eigenvalues of the fixed point associated with its Poincaré map have a negative real part; otherwise, the periodic solution is unstable. Inside the unstable periodic solutions there are two types: the unstable saddle periodic solutions having an eigenvalue with a negative real part and the other with a positive real part, and the repeller with both eigenvalues having a negative real part. We should note that the method used here for studying periodic solutions can be applied to any periodic non-autonomous differential system as done in [9,18]. In these papers, the authors applied the method used in this paper in order to guarantee the existence of periodic solutions in a periodic FitzHugh–Nagumo system and in the Vallis system, respectively. The paper is organized as follows. In Section 1.2, the main results are stated and compared with other results. Also, some important points on the results are clarified. In Section 2.1, we prove the results. In Section 2.2, some aspects of the results are pointed out. Lastly, Section 3 is devoted to give a brief summary of the results that we use from the averaging theory. 1.2 Statement of the main results In this section, we present our results. We will behave under an analytical approach and for this reason the results are valid for a large range of the parameters of system (4), different from the major part of the results dealing with numerical techniques. In addition, it is important to note that we are concerned with harmonic solutions in the sense that they do not bifurcate from periodic solutions with multiple periods. We have the following results. Theorem 1 Consider ε > 0 sufficiently small, (ρ, δ, α, λ) = (εr, εd, εn2 a, εn3�) with n2, n3 > 1 and 81d2 > 48r2 > 0. Then system (4) has a 2π -periodic solution (x(t, ε), y(t, ε)) such that (x(0, ε), y(0, ε)) → ( 0, 2(36r)1/3 3� + 1 3 ( 6 r )1/3 � ) , when ε → 0, where � = (9d + √ 81d2 − 48r2)1/3. Moreover, for r sufficiently small, this periodic solution is stable. International Journal of Computer Mathematics 1361 Theorem 2 Consider ε > 0 sufficiently small, (α, δ, ρ, λ) = (εa, εd, εn1 r, εn3�) with n1, n3 > 1 and α �= 0. Then system (4) has a 2π -periodic solution (x(t, ε), y(t, ε)) such that (x(0, ε), y(0, ε)) → (( 4d 3a )1/3 , 0 ) , when ε → 0. Theorem 3 Consider ε > 0 sufficiently small and (λ, δ, ρ, α) = (ε�, εd, εn1 r, εn2 a) with n1, n2 > 1. Then system (4) has a 2π -periodic solution (x(t, ε), y(t, ε)) such that (x(0, ε), y(0, ε)) → (( 8d 5� )1/5 , 0 ) , when ε → 0. We remark that the stability of the periodic solutions of Theorems 2 and 3 cannot be decided with the real part of the eigenvalues of their Poincaré map because these real parts are zero. Theorem 4 Consider ε > 0 sufficiently small and (x, y, ρ, δ, α, λ) = (εX , εY , εr, ε2d, εn2 a, εn3�) with ρ �= 0. Then system (4) has a 2π -periodic solution (x(t, ε), y(t, ε)) such that (x(0, ε), y(0, ε)) → ( O(ε), − δ ερ + O(ε) ) , when ε → 0. Moreover, this periodic solution is unstable. Note that in Theorem 4 we have a periodic solution that comes from infinity. As far as we know, this kind of behaviour also has not been observed in the papers concerned with system (4). Nevertheless, this behaviour is common when we perform a rescaling in the parameters and variables as we made previously and can be observed also in [9,18]. Theorem 5 Consider ε > 0 sufficiently small, (α, λ, δ, ρ) = (εa, ε�, εd, εn1 r) with n1 > 1, αλ < 0 and D = 53747712a5d2 78125�7 + 20480d4 �4 . Then system (4) has a 2π -periodic solution if D = 0 and three 2π -periodic solutions if D < 0. Moreover, there are values of α, λ and δ such that for all ρ the system realizes the number of these periodic solutions. Theorem 6 Consider ε > 0 sufficiently small, (ρ, α, δ, λ) = (εr, εa, εd, εn3�) with n3 > 1 and −324a4r2 − 9a2d2r2 + 36a2r4 − d2r4 �= 0. Consider also the values 1 = 324a4 + d2r2 + 9a2(d2 − 4r2), 1362 R.D. Euzébio and J. Llibre 2 = 2187a4d4 + 27d4r4 − 16d2r6 + 18a2(27d4r2 − 72d2r4 + 32r6). The system rx0(−4 + x2 0 + y2 0) − 3ay0(x 2 0 + y2 0) = 0, 4d − ry0(−4 + x2 0 + y2 0) − 3ax0(x 2 0 + y2 0) = 0 has one solution (x0 0, y0 0) if 1 2 = 0 or 2 > 0, and three solutions (xi 0, yi 0) if 2 < 0 for i = 1, 2, 3; we assume that all of them satisfy 27α2((xi 0) 2 + (yi 0) 2)2 + ρ2(−4 + (xi 0) 2 + (yi 0) 2)(−4 + 3(xi 0) 2 + 3(yi 0) 2) �= 0. Then if 1 2 = 0 or 2 > 0, system (4) has a 2π -periodic solution (x0(t, ε), y0(t, ε)) such that (x0(0, ε), y0(0, ε)) → (x0 0, y0 0) when ε → 0. Additionally if 2 < 0 system (4) has three 2π -periodic solutions (xi(t, ε), yi(t, ε)) such that (xi(0, ε), yi(0, ε)) → (xi 0, yi 0) when ε → 0 for i = 1, 2, 3. Moreover, there are values of ρ, α and δ such that for all λ the system realizes the number of these periodic solutions. Theorem 7 Consider ε > 0 sufficiently small, (α, ρ, λ, δ) = (εa, εr, ε�, εd), αρ �= 0, (r2 − 3a2)(r2 − 9a2) �= 0 and C = 4(3a + 10�)(540a3� − 9a2(d2 − 600�2) − 120a�(d2 − 150�2) + 50�2(−7d2 + 400�2))r6 + 6(a + 5�)(6a − d + 20�)(6a + d + 20�)r8 �= 0, D = 2066242608a14d6 − 3125d8r12 − 531441a12(3125d8 + 96d6r2 + 1536d4r4 − 1024d2r6) + 354294a10(3125d8r2 + 616d6r4 − 1600d4r6) + 18a2r10(9375d8 + 5000d6r2 − 88000d4r4 + 102400d2r6 − 32768r8) + 2916a6r6(15625d8 + 11700d6r2 − 45632d4r4 + 23552d2r6 − 2048r8) − 6561a8r4(46875d8 + 24832d6r2 − 66816d4r4 + 12288d2r6 + 4096r8) + 81a4r8(−46875d8 − 36000d6r2 + 275200d4r4 − 246784d2r6 + 61440r8). The system 2rx0(−4 + x2 0 + y2 0) − y0(x 2 0 + y2 0)(6a + 5�(x2 0 + y2 0)) = 0, 8d + 8ry0 − (x2 0 + y2 0)(6ax0 + 2ry0 + 5�x0(x 2 0 + y2 0)) = 0 has one solution (x0 0, y0 0) if D < 0, and three solutions (xi 0, yi 0) if D > 0 for i = 1, 2, 3; we assume that all of them satisfy 4ρ2(−4 + (xi 0) 2 + (yi 0) 2)(−4 + 3(xi 0) 2 + 3(yi 0) 2) + ((xi 0) 2 + (yi 0) 2)2(6α + ρλ((xi 0) 2 + (yi 0) 2))(18α + 25λ(xi 0) 2 + 3(yi 0) 2) �= 0. Then, if D < 0 system (4) has a 2π -periodic solution (x0(t, ε), y0(t, ε)) such that (x0(0, ε), y0(0, ε)) → (x0 0, y0 0) when ε → 0. Furthermore, if D > 0 system (4) has three 2π - periodic solutions (xi(t, ε), yi(t, ε)) such that (xi(0, ε), yi(0, ε)) → (xi 0, yi 0) when ε → 0 for i = 1, 2, 3. Moreover, there are values of α, ρ, λ and δ for which the system realizes the number of these periodic solutions. International Journal of Computer Mathematics 1363 Theorem 8 Consider ε > 0 sufficiently small, (ρ, λ, δ, α) = (εr, ε�, εd, εn2 a) with n2 > 1. Consider the numbers C = 4000d3�3r4 + 7200000d�5r4 + 60d3�r6 + 72000d�3r6 + (50d4�2r4 − 220000d2�4r4 − 8000000�6r4 − 2800d2�2r6 + 160000�4r6 − 6d2r8 + 2400�2r8), N2 = −4500�2r4 + 3r6 3125d2�4 , N3 = 5859375d4�6 + 18750000d2�6r2 + 87500�4(−7d2 + 1200�2)r4 + 25�2(−7d2 + 155600�2)r6 + 15600�2r8 + 9r10, N4 = 1171875d8�6(1250000�4 + 8500�2r2 + r4) − 2500d6�4(312500000000�8 + 5343750000�6r2 + 74625000�4r4 + 358375�2r6 + 69r8) − 3000d4�2r2(218750000000�10 − 21625000000�8r2 − 686625000�6r4 − 3715000�4r6 − 650�2r8 − 3r10) − 1600�2r6(400000000000�10 − 144900000000�8r2 − 2817000000�6r4 − 8880000�4r6 + 38700�2r8 + 27r10) + 12d2r4(−50000000000000�12 − 32250000000000�10r2 − 662400000000�8r4 − 2803000000�6r6 + 5410000�4r8 + 19200�2r10 + 9r12) N5 = 48828125d8�6 − 4687500d6�4r4 + 6400�2r10(1600�2 + r2) − 16d2r8 (2000000�4 + 2900�2r2 + r4) + d4r6(27500000�4 + 60000�2r2 + 27r4) M5 = −25d4�2 + d2(110000�4 + 1400�2r2 + 3r4) + 400(10000�6 − 200�4r2 − 3�2r4). We assume that CM5 �= 0. The system 2rx0(−4 + x2 0 + y2 0) − 5�y0(x 2 0 + y2 0) 2 = 0, 8d − 2ry0(−4 + x2 0 + y2 0) − 5�x0(x 2 0 + y2 0) 2 = 0 has one solution (x0 0, y0 0) if N2 ≤ 0, N3 ≥ 0 or N4 ≥ 0 and N5 > 0, and three solutions (xi 0, yi 0) if N5 < 0 for i = 1, 2, 3; we assume that all of them satisfy 125λ2((xi 0) 2 + (yi 0) 4) + 4ρ2(−4 + (xi 0) 2 + (yi 0) 2)(−4 + 3(xi 0) 2 + 3(yi 0) 2) �= 0. Hence, if N2 ≤ 0, N3 ≥ 0 or N4 ≥ 0 and N5 > 0, then system (4) has a 2π -periodic solution (x0(t, ε), y0(t, ε)) such that (x0(0, ε), y0(0, ε)) → (x0 0, y0 0) when ε → 0. If N5 < 0 then system (4) has three 2π -periodic solutions (xi(t, ε), yi(t, ε)) such that (xi(0, ε), yi(0, ε)) → (xi 0, yi 0) when ε → 0 for i = 1, 2, 3. Moreover, there are values of ρ, λ and δ such that for all α system (4) has one or three periodic solutions. We note that in Theorems 5, 6, 7 and 8 we do not say anything about the kind of stability of the periodic solutions because we do not have the explicit expressions of the real part of the eigenvalues of their Poincaré maps. We remark that the periodic solutions provided in Theorems 1 until 8 exist when we take small values for the parameters of system (4) obeying some relations among them. For instance, Theorem 1 states the existence of one periodic solution for system (4) if each parameter of this system is small and α and λ are much smaller than ρ and δ. This assertion becomes more clear 1364 R.D. Euzébio and J. Llibre if we observe the replacement done in each theorem and take into account that ε is sufficiently small. As far as we know, periodic solutions of system (4) whose parameters have this character- istic have not been observed in the literature. Besides, it seems that the simultaneous bifurcation of three harmonic periodic solutions in the extended Duffing–Van der Pol system is also new. Furthermore, we note that the periodic solutions presented in the present paper are different from those stated in [8]. Indeed, in the referred paper, the authors ask for V ′ 6(0) < 0, and in our case, we have V ′ 6(0) = 0. 2. Proof and discussion of the results 2.1 Proof of the results In order to apply the averaging theory described in Section 3 in systems (4), we start performing a rescaling of the variables x and y and of the parameters ρ, α, λ and δ as follows: x = εm1 X , y = εm2 Y , ρ = εn1 r, α = εn2 a, λ = εn3�, δ = εn4 d, (5) where ε is positive and sufficiently small and mi and nj are non-negative integers, for i = 1, 2 and j = 1, 2, 3, 4. We recall that since δ > 0, ρ ≥ 0 and λ �= 0 we have d > 0, r ≥ 0 and � �= 0. In the new variables (X , Y), system (4) is written as dX dt = ε−m1+m2 Y , dY dt = −εm1−m2 X + εn1 rY − ε2m1+n1 rX 2Y − ε3m1−m2+n2 aX 3 − ε5m1−m2+n3�X 5 + ε−m2+n4 d cos t. (6) In order to have non-negative powers of ε, we must impose the conditions m1 = m2 = m and n4 ≥ m, (7) where m is a non-negative integer. Hence, with conditions (7), system (6) becomes dX dt = Y , dY dt = −X + εn1 rY − ε2m+n1 rX 2Y − ε2m+n2 aX 3 − ε4m+n3�X 5 + ε−m+n4 d cos t. (8) In this paper, we will find periodic solutions of system (8) depending on the parameters r, a, � and d and the powers m, ni of ε, for i = 1, 2, 3, 4. Then, we will go back through rescaling (5) to ensure the existence of periodic solutions in system (4). From now on we assume that the values n1, m2 + n2 2, m2 + n2 3 and n4 − m are positive and observe that considering these conditions each power of ε in system (8) becomes positive. The reason for a such assumption will be explained later on in Section 2.2. Now we shall apply the averaging theory described in Section 3. Thus following the notation of the mentioned section and denoting again the variables (X , Y) by (x, y), we have x = (x, y)T, and system (13) International Journal of Computer Mathematics 1365 corresponding to system (8) can be written as ẋ = F0(t, x) = (y, −x)T. (9) So the solution x(t, z) = (x(t, z), y(t, z)) of system (9) such that x(0, z) = z = (x0, y0) is x(t, z) = x0 cos t + y0 sin t, y(t, z) = y0 cos t − x0 sin t. It is clear that the origin of coordinates of R 2 is a global isochronous centre for system (9) whose circular periodic solutions starting on z = (x0, y0) ∈ R 2 are parametrized by the above functions x(t, z) and y(t, z). Therefore, through every initial condition (x0, y0) in R 2 passes a 2π -periodic solution of system (9). Physically speaking, system (9) models a simple harmonic oscillator and its solutions (x(t, z), y(t, z)) describe the wave behaviour of this oscillator. In this direction, the problem of perturbation of the global centre (9) is equivalent to the problem of perturbation of a simple har- monic oscillator introducing a damping parameter as external force and considering a potential V by taking small values for ρ, α, λ and δ as stated in Theorems 1–8. We note that using the notation of Section 3, the fundamental matrix Y(t, z) of system (9) satisfying that Y(0, z) is the identity of R 2 is written as Y(t, z) = ( cos t sin t − sin t cos t ) . From the averaging theory we are interested in the simple zeros of the function f (z) = (f1(z), f2(z)) = ∫ 2π 0 Y−1(t, z)F1(t, x(t, z, 0)) dt, (10) where F1(t, x) are the terms of order 1 on ε of the vector field associated with system (8) (y, −x + εn1 ry − ε2m+n1 rx2y − ε2m+n2 ax3 − ε4m+n3�x5 + ε−m+n4 d cos t)T. Now we prove our results. Proof of Theorem 1 First we take n1 = n4 = 1, m = 0 and consider n2, n3 > 1. So rescaling (5) becomes (x, y, ρ, δ, α, λ) = (X , Y , εr, εd, εn2 a, εn3 d) and calling again (X , Y) by (x, y) we verify the hypotheses of Theorem 1. Now the vector field of system (4) becomes (y, −x + εry − εrx2y − εn2 ax3 − εn3�x5 + εd cos t)T. Hence, as we stated before, the terms of order 1 in this vector field are given by F1(t, x) = (0, ry − rx2y + d cos t). Then since z = (x0, y0), the function f (x0, y0) = (f1(x0, y0), f2(x0, y0)) given in Equation (10) turns into the form f1(x0, y0) = ∫ 2π 0 −(sin t)(d cos t + r(y0 cos t − x0 sin t) − r(y0 cos t − x0 sin t)(x0 cos t + y0 sin t)2) dt 1366 R.D. Euzébio and J. Llibre = −1 4 πrx0(−4 + x2 0 + y2 0), f2(x0, y0) = ∫ 2π 0 (cos t)(d cos t + r(y0 cos t − x0 sin t) − r(y0 cos t − x0 sin t)(x0 cos t + y0 sin t)2) dt = 1 4 (4d − ry0(−4 + x2 0 + y2 0)). The zeros of functions f1 and f2 are a pair of conjugate complex vectors and a real pair (x0 0, y0 0) satisfying x0 0 = 0 and y0 0 = 2 3 (36r)1/3 (9d + √ 81d2 − 48r2)1/3 + 1 3 ( 6 r )1/3 (9d + √ 81d2 − 48r2)1/3. Note that y0 0 is real because 81d2 > 48r2 > 0 by assumption. We observe that the expression of y0 0 = y0 0(r, d) does not change when we go back through rescaling (5) taking r = ε−1ρ and d = ε−1δ according to the hypotheses of Theorem 1. Moreover, we denote � = (9d + √ 81d2 − 48r2)1/3 and observe that the matrix M = (mij) = ∂(f1, f2)/∂(x0, y0) at (x0 0, y0 0) is diagonal and its elements are m11 = −πr�+/36 and m22 = πr�−, where �± = −12 ± 61/324r2/3 �2 ± 62/3�2 r2/3 . The determinant �0 of M written in a power series of r around r = 0 is �0 = 3 8 (2π2d(4d)1/3)r2/3 + (πr)2 3 + O(r7/3), which is positive for |r| sufficiently small because d > 0. So the averaging theory described in Section 3 guarantees the existence of a 2π -periodic solution (x(t, ε), y(t, ε)) such that (x(0, ε), y(0, ε)) tends to (x0 0, y0 0) when ε → 0. On the other hand, the trace �0 of M at (x0 0, y0 0) is �0 = −22/3π(2d)2/3r1/3 − 2πr 3 + O(r4/3). Thus for |r| sufficiently small, the trace �0 of M is negative. So the real part of the eigenvalues of the matrix M are both negative and then the periodic solution (x(t, ε), y(t, ε)) is stable for each ε sufficiently small. So Theorem 1 is proved. � Proof of Theorem 2 Now we take n2 = n4 = 1 and consider n1, n3 > 1. So we are under the hypotheses of Theorem 2. Now, taking m = 0 and doing rescaling (5), the vector field of system (8) is written as (y, −x + εn1 ry − εn1 rx2y − εax3 − εn3�x5 + εd cos t)T. International Journal of Computer Mathematics 1367 Thus, we have F1(t, x) = (0, −ax3 + d cos t), and f (x0, y0) = (f1(x0, y0), f2(x0, y0)) becomes f1(x0, y0) = ∫ 2π 0 −(sin t)(d cos t − a(x0 cos t + y0 sin t)3) dt = 3 4 aπy0(x 2 0 + y2 0), f2(x0, y0) = ∫ 2π 0 (cos t)(d cos t − a(x0 cos t + y0 sin t)3) dt = 1 4 π(4d − 3ax0(x 2 0 + y2 0)). The real zero (x0 0, y0 0) of f1 and f2 is (x0 0, y0 0) = (( 4d 3a )1/3 , 0 ) . Again the values (x0 0, y0 0) does not depend on ε when we go back to the original parameters α and δ through rescaling (5). Besides, the matrix M = ∂(f1, f2)/∂(x0, y0) at (x0 0, y0 0) now is written as M = 1 ε ⎛ ⎜⎜⎜⎝ 0 (( 3 4 ) αδ2 )1/3 π −3 (( 3 4 ) αδ2 )1/3 π 0 ⎞ ⎟⎟⎟⎠ , and then the determinant �0 of M is �0 = 3π2 2ε2 ( 9α2δ4 2 )1/3 > 0. Therefore, the averaging theory implies the existence of a 2π -periodic solution (x(t, ε), y(t, ε)) such that (x(0, ε), y(0, ε)) tends to (x0 0, y0 0) when ε → 0. � We observe that since the trace �0 of M is zero and the determinant �0 is positive, then the real part of both eigenvalues of M is zero. So we cannot decide about the stability of the periodic solution of Theorem 2. In what follows, we see that the same occurs in Theorem 3. Proof of Theorem 3 We start considering n3 = n4 = 1, n1, n2 > 1 and m = 0. Therefore, we are under the assumptions of Theorem 3. Doing rescaling (5) the vector field of system (8) becomes (y, −x + εn1 ry − εn1 rx2y − εn2 ax3 − ε�x5 + εd cos t)T. So we conclude that F1(t, x) = (0, −�x5 + d cos t) and then we get f1(x0, y0) = ∫ 2π 0 −(sin t)(d cos t − �(x0 cos t + y0 sin t)5) dt = 5 8 �πy0(x 2 0 + y2 0) 2, 1368 R.D. Euzébio and J. Llibre f2(x0, y0) = ∫ 2π 0 (cos t)(d cos t − �(x0 cos t + y0 sin t)5) dt = 1 8 π(8d − 5�x0(x 2 0 + y2 0) 2). Now the real zero (x0 0, y0 0) of f1 and f2 is (x0 0, y0 0) = (( 8d 5� )1/5 , 0 ) . We note that as in Theorems 1 and 2 the root (x0 0, y0 0) does not depend on ε after we perform rescaling (5). However, using Equation (5), the matrix M = ∂(f1, f2)/∂(x0, y0) at (x0 0, y0 0) is M = 1 ε ⎛ ⎜⎜⎜⎝ 0 (( 5 8 ) λδ4 )1/5 π −5 (( 5 8 ) λδ4 )1/5 π 0 ⎞ ⎟⎟⎟⎠ , whose determinant is �0 = 5π2 2ε2 ( 25λ2δ8 2 )1/5 > 0. Then the averaging method states the existence of a 2π -periodic solution (x(t, ε), y(t, ε)) such that (x(0, ε), y(0, ε)) tends to (x0 0, y0 0) when ε → 0. The stability of this periodic solution is unknown because the trace �0 of M is zero and the determinant �0 > 0 is positive just as in Theorem 2. � Proof of Theorem 4 In what follows, we assume that n1 = m = 1 and n4 = 2. Now the expression of the vector field of system (8) is (y, −x + εry − ε3rx2y − εn2+2ax3 − εn3+4�x5 + εd cos t)T, and then F1(t, x) = (0, ry + d cos t). In this case, we obtain f1(x0, y0) = ∫ 2π 0 −(sin t)(d cos t + r(y0 cos t − x0 sin t)) dt = πrx0, f2(x0, y0) = ∫ 2π 0 (cos t)(d cos t + r(y0 cos t − x0 sin t)) dt = π(d + ry0). It is immediate that the only zero (x0 0, y0 0) of f1 and f2 in this case is (x0 0, y0 0) = (0, d/r). So after going back through rescaling (5) taking r = ε−1ρ and d = ε−2δ, we have (x0 0, y0 0) = ( 0, − δ ερ ) . In addition, the matrix M = ∂(f1, f2)/∂(x0, y0) at (x0 0, y0 0) is diagonal and its elements are both given by ρπ/ε. Consequently, M has determinant �0 = (ρπ/ε)2 > 0 and trace �0 = 2ρπ/ε > International Journal of Computer Mathematics 1369 0. Then the averaging theory described in Section 3 assures the existence of a 2π -periodic solu- tion (x(t, ε), y(t, ε)) such that (x(0, ε), y(0, ε)) tends to (x0 0, y0 0) when ε → 0. Moreover, since ρ is positive this periodic solution is unstable. Then, Theorem 4 is proved. � Proof of Theorem 5 Now we take m = 0, n2 = n3 = n4 = 1 and n1 > 1. The expression of the vector field of system (8) in this case is (y, −x + εn1 ry − εn1 rx2y − εax3 − ε�x5 + εd cos t)T, and we get F1(t, x) = (0, −ax3 − �x5 + d cos t). Hence, function f = (f1, f2) is written as f1(x0, y0) = ∫ 2π 0 −(sin t)(d cos t − a(x0 cos t + y0 sin t)3 − �(x0 cos t + y0 sin t)5) dt = 1 8 πy0(x 2 0 + y2 0)(6a + 5�(x2 0 + y2 0)), f2(x0, y0) = ∫ 2π 0 (cos t)(d cos t − a(x0 cos t + y0 sin t)3 − �(x0 cos t + y0 sin t)5) dt = 1 8 π(8d − x0(x 2 0 + y2 0)(6a + 5�(x2 0 + y2 0))). In order to find the zeros of f (x0, y0), we compute a Gröbner basis {bk(x0, y0), k = 1, 2, 3} in the variables x0 and y0 for the set of polynomials {f̄1(x0, y0), f̄2(x0, y0)}, where f̄1,2 = (8/π)f1,2. Then, we look for the zeros of each bk , k = 1, 2, 3. It is a known fact that the zeros of a Gröbner basis of {f̄1(x0, y0)), f̄2(x0, y0)} are the zeros of f̄1 and f̄2, and consequently, zeros of f1 and f2 too. For more information about a Gröbner basis, see [1,15]. The Gröbner basis is b1(x0, y0) = dy0, b2(x0, y0) = (6ax2 0 + 5�x4 0)y0 + (6a + 10�x2 0)y 3 0 + 5�y5 0, b3(x0, y0) = −8d + 6ax3 0 + 5�x5 0 + (6ax0 + 10�x3 0)y 2 0 + 5�x0y4 0. The only zero of b1(x0, y0) is y0 0 = 0 since d = ε−1δ > 0. So replacing y0 0 into b3(x0, y0) and simplifying the new expressions, we reduce our problem to find a zero of the polynomial p(x0) = 1 5� b3(x0, 0) = −8d r� + 6a 5� x3 0 + x5 0. (11) As we know from algebra, there is no general formula to provide the roots of a quintic poly- nomial as polynomial (11). Nevertheless, some techniques may make this task easier. An elegant one can be found in [29]. In this paper, a way is provided to study the number of roots (complex and real, including multiplicity) of a polynomial of any order by only performing some calcula- tions on the coefficients of the considered polynomial. We will present a brief summary of the algorithm of this method in the appendix. By the Fundamental Theorem of Algebra, polynomial (11) has five roots taking into account their multiplicities and once we want to apply the averaging theory, we are only interested in the real simple roots of polynomial (11). Indeed, applying the method due to [29] and summarized 1370 R.D. Euzébio and J. Llibre in the appendix, we obtain D4 = 384ad2 �3 , D5 = 53747712 78125 a5d2 �7 + 20480 d4 �4 . Therefore as described in the appendix, since αλ = ε2a� < 0, by hypothesis we have that D4 is negative. So polynomial (11) has a unique simple real root if D5 = D = 0 according to statement (5) of the appendix. On the other hand, statement (3) says that polynomial (11) has three sim- ple real roots if D5 = D < 0. Moreover, D and D4 do not change when we replace (a, �, d) by (ε−1α, ε−1λ, ε−1δ) using rescaling (5). Now using the property on zeros of the Gröbner basis cited previously, the function f (x0, y0) has a zero (x0 0, 0) if D = 0 and three zeros (xi 0, 0) if D < 0, with i = 1, 2, 3. Besides the matrix M = ∂(f1, f2)/∂(x0, y0) evaluated at y0 0 = 0 after rescaling (5) is M = 1 ε ⎛ ⎜⎝ 0 π 8 (6αx2 0 + 5λx4 0) −π 8 (18αx2 0 + 25λx4 0) 0 ⎞ ⎟⎠ , whose determinant is �(x0) = 1 ε2 ( 27 16 (απx2 0) 2 + 15 4 αλ(πx3 0) 2 + 125 64 (λπx4 0) 2 ) . Thus by computing a Gröbner basis for polynomials {�(x0), f1(x0, 0), f2(x0, 0)} in the variables α, δ, λ and x0, we get the polynomials b̄1(α, δ, λ, x0) = δ and b̄2(α, δ, λ, x0) = 6α + 5λx2 0. So b̄1(α, δ, λ, x0) > 0 because δ > 0 and it means that we cannot have f1, f2 and � equal to zero simultaneously. Consequently, �(xi 0) �= 0 for all i = 0, 1, 2, 3. Then using the averaging theory we conclude the existence of a 2π -periodic solution if D = 0 and three 2π -periodic solutions if D < 0. Now we exhibit the values of α, λ and δ for which system (4) has one or three 2π -periodic solutions. On the one hand, if consider a = 25 18 , � = −1 and d = 5 12 , we obtain a� = − 25 18 < 0 and D = 0. On the other hand, taking d = 18a/31 and � = −42a/155, we get a� = −42a2/155 < 0 and D = − 1425339825408 823543 < 0. So Theorem 5 is proved. � We note that by assumption in Theorem 5 we have a� < 0. So considering the notation of the appendix, the coefficients D2 = −6a/5� and D3 = 384ad2/�3 corresponding to polynomial (11) are positive and D4 is negative. Thus the only possible configurations of roots of polynomial (11) are those one listed in statements (3) and (5) of the appendix. This fact means that we cannot have two or five periodic solutions when we consider m = 0, n2 = n3 = n4 = 1 and n1 > 1. Proof of Theorem 6 We start fixing the values m = 0, n1 = n2 = n4 = 1 and n3 > 1. Now the vector field of system (8) is written as (y, −x + εry − εrx2y − εax3 − εn3�x5 + εd cos t)T, International Journal of Computer Mathematics 1371 and then F1(t, x) = (0, ry − rx2y − ax3 + d cos t). With this expression of F1, function (10) is as follows: f1(x0, y0) = ∫ 2π 0 −(sin t)(d cos t + r(y0 cos t − x0 sin t) − r(y0 cos t − x0 sin t)(x0 cos t + y0 sin t)2 − a(x0 cos t + y0 sin t)3) dt = −1 4 π(rx0(−4 + x2 0 + y2 0) − 3ay0(x 2 0 + y2 0)), f2(x0, y0) = ∫ 2π 0 (cos t)(d cos t + r(y0 cos t − x0 sin t) − r(y0 cos t − x0 sin t)(x0 cos t + y0 sin t)2 − a(x0 cos t + y0 sin t)3) dt = 1 4 π(4d − ry0(−4 + x2 0 + y2 0) − 3ax0(x 2 0 + y2 0)). In order to find zeros (x∗ 0, y∗ 0) of f (x0, y0) as before, we compute a Gröbner basis {bk(x0, y0), k = 1, . . . , 14} in the variables x0 and y0 now for the set of polynomials {f̄1(x0, y0), f̄2(x0, y0)} with f̄1,2 = ∓(4/π)f1,2. We will look for zeros of two elements of the Gröbner basis for the polynomials {f̄1(x0, y0), f̄2(x0, y0)} in the variables x0 and y0. These polynomials are b1(x0, y0) = 144a2dr3 − 4d3r3 + (108a2d2r2 + 144a2r4 − 4d2r4)y0 + 144a2dr3y2 0 + (81a4d2 + 18a2d2r2 + d2r4)y3 0, b2(x0, y0) = 216a3dr2 + (−324a4r2 − 9a2d2r2 + 36a2r4 − d2r4)x0 + (27a3d2r + 216a3r3 + 3ad2r3)y0 + (243a5d + 54a3dr2 + 3adr4)y2 0. We must observe that b1(x0, y0) depends only on y0. Then for each zero y∗ 0 of b1(y0), the second polynomial b2(x0, y0) provides a zero x∗ 0 associated with y∗ 0 because the coefficient of x0 in b2(x0, y0) is not zero by the hypothesis. Now we will look for zeros of b1(y0). Indeed, the discriminant of the cubic polynomial b1(y0) is written as = −16d2r6(324a4 + d2r2 + 9a2(d2 − 4r2))2 × (2187a4d4 + 27d4r4 − 16d2r6 + 18a2(27d4r2 − 72d2r4 + 32r6)) = −16d2r6 2 1 2. So if 1 or 2 is zero, then is zero and consequently b1(y0) has a simple real root y0 0 because this polynomial has no root with multiplicity 3. Actually b′′′ 1 (y0) = 6(81a4d2 + 18a2d2r2 + d2r4) > 0. The same occurs if 2 > 0, since this condition implies < 0. On the other hand, if 2 < 0, then > 0 and consequently polynomial b1(y0) has three simple real roots yi 0, i = 1, 2, 3. Additionally, replacing each value yi 0 into b2(x0, y0) we obtain the respective val- ues xi 0, for = 0, 1, 2, 3. We note that coming back through rescaling (5), the signs of , 1 and 2 does not change because each monomial composing has the same degree. Now we will verify the condition M = det((∂f /∂z)(xi 0, yi 0)) �= 0 for i = 0, 1, 2, 3. In fact M is −π 4ε (−6αx0y0 + ρ(−4 + 3x2 0 + y2 0) 2ρx0y0 − 3α(x2 0 + 3y2 0) 2ρx0y0 + 3α(3x2 0 + y2 0) 6αx0y0 + ρ(−4 + x2 0 + 3y2 0) ) , whose determinant � now is �(x0, y0) = π2 16ε2 (27α2(x2 0 + y2 0) 2 + ρ2(−4 + x2 0 + y2 0)(−4 + 3x2 0 + 3y2 0)). 1372 R.D. Euzébio and J. Llibre However for each i = 0, 1, 2, 3 the determinant �(xi 0, yi 0) �= 0 by hypothesis. Therefore, using the averaging theory described in Section 3, system (4) has a 2π -periodic solution if 1 2 = 0 or 2 > 0 and three 2π -periodic solutions if 2 < 0. Now we present concrete values of α, ρ and δ for which we have one or three periodic solu- tions. Indeed taking a = r = 1 and d = 6, we have 2 = 3452544 > 0. Meantime if we consider d = 6a and r = − 9a√ −9 + 4 √ 6 , we obtain 2 = −34012224a8/(9 − 4 √ 6)2 < 0. We also observe that the values (r, a, d) = (( 555 4 ) √ 154073 9622 , 185, 22) and (r, a, d) = (585 √ 3, 195, 585 √ 2) make 1 and 2 equal to zero, respectively. This ends the proof of Theorem 6. � Proof of Theorem 7 For this case, we assume m = 0 and ni = 1, i = 1, 2, 3, 4. Then, the vector field of system (8) is written as (y, −x + εry − εrx2y − εax3 − ε�x5 + εd cos t)T, and we have F1(t, x) = (0, ry − rx2y − ax3 − �x5 + d cos t). So it follows that f1(x0, y0) = ∫ 2π 0 −(sin t)(d cos t + r(y0 cos t − x0 sin t) − r(y0 cos t − x0 sin t) (x0 cos t + y0 sin t)2 − a(x0 cos t + y0 sin t)3 − �(x0 cos t + y0 sin t)5) dt = −1 8 π(2rx0(−4 + x2 0 + y2 0) − y0(x 2 0 + y2 0)(6a + 5�(x2 0 + y2 0)), f2(x0, y0) = ∫ 2π 0 (cos t)(d cos t + r(y0 cos t − x0 sin t) − r(y0 cos t − x0 sin t) (x0 cos t + y0 sin t)2 − a(x0 cos t + y0 sin t)3 − �(x0 cos t + y0 sin t)5) dt = 1 8 π(8d + 8ry0 − (x2 0 + y2 0)(6ax0 + 2ry0 + 5�x0(x 2 0 + y2 0))). We will find zeros of f (x0, y0) through the zeros of a Gröbner basis {bk(x0, y0)} of it. As before, instead of f (x0, y0), the Gröbner basis will be related with the functions f̄1,2 = ∓(8/π)f1,2. One of the elements of the Gröbner basis in the variables x0 and y0 is b1(x0, y0) = 144a2dr5 − 4d3r5 + 960ad�r5 + 1600d�2r5 + (108a2d2r4 + 960ad2�r4 + 2000d2�2r4 + 144a2r6 − 4d2r6 + 960a�r6 + 1600�2r6)y0 + (120ad3�r3 + 500d3�2r3 + 144a2dr5 + 1200ad�r5 + 2400d�2r5)y2 0 + (81a4d2r2 + 540a3d2�r2 + 900a2d2�2r2 + 18a2d2r4 + 300ad2�r4 + 900d2�2r4 + d2r6)y3 0 + (−450a2d3�2r − 1500ad3�3r + 50d3�2r3)y4 0 + 625d4�4y5 0. We note that again b1(x0, y0) depends only on y0. Moreover, the Gröbner basis has another element b2(x0, y0) which is linear on x0 and the coefficients depend on y0, where the coefficient C of x0 is not zero by hypotheses. Its expression is too large and we will omit it here. It can be easily obtained by an algebraic manipulator. So for each zero y∗ 0 of b1(y0), the second polynomial b2(x0, y0) provides a zero x∗ 0 associated with y∗ 0. International Journal of Computer Mathematics 1373 Now we will look for zeros of b1(y0). Indeed using again the appendix, the conditions on a, r, d and � that provide three zeros, (xi 0, yi 0), where i = 1, 2, 3 for b1(y0), are D5 = −D1 5D2 5D < 0, where D is given in Theorem 7 and D1 5 and D2 5 are the values D1 5 = 110075314176a8r20 d20(−9a2 + r2)28 , D2 5 = (59049a10d4 + 16r14 + 36a2r10(11d2 − 4r2) − 81a4r6(d4 − 12d2r2 + 16r4) − 6561a8(3d4r2 − 4d2r4) + 729a6(3d4r4 − 28d2r6 + 16r8))2. By hypothesis ar = ε−2αρ �= 0 and D2 5 �= 0, then D1 5 and D2 5 are positive. Therefore, in order to have D5 negative, we must have D positive. On the other hand, if D5 is positive then according to statement (2) of the appendix, we have a unique zero (x0 0, y0 0) for b1(y0) because the value D2 = −1296a2r4(−3a2 + r2)2 d2(−9a2 + r2)4 is negative if (r2 − 3a2)(r2 − 9a2) �= 0. In addition by considering a = ε−1α, r = ε−1ρ and d = ε−1δ, we get M = −π 8ε ( 2ρ(−4 + 3x2 0 + y2 0) − �1 �2 − �4 − 5λ�3(x2 0 + 5y2 0) �2 + �4 + 5λ�3(5x2 0 + y2 0) 2ρ(−4 + x2 0 + 3y2 0) − �1 ) , where �1 = 2x0y0(3α + 5�(x2 0 + y2 0)), �2 = −4ρx0y0, �3 = x2 0 + y2 0 and �4 = 6α(x2 0 + 3y2 0). Moreover, determinant � of M is written as �(x0, y0) = 1 ε2 (4ρ2(−4 + x2 0 + y2 0)(−4 + 3x2 0 + 3y2 0) + (x2 0 + y2 0) 2(6α + ρλ(x2 0 + y2 0)(18α + 25λx2 0 + 3y2 0))). However, by hypothesis �(xi 0, yi 0) �= 0 for i = 0, 1, 2, 3 and hence we conclude the first part of Theorem 7. Now we provide values for α, ρ, λ and δ for which we have one or three 2π -periodic solutions. First, we observe that taking (a, r, d) = (1, 2, 8 3 ), we obtain D = −2770035802112 6561 < 0. However, if we consider (a, r, d, �) = ( 1√ 3 , 1, 1, − 1 (5 √ 3) ) , then D = 2752. So Theorem 7 is proved. � Now we will prove that under the conditions m = 0 and ni = 1, i = 1, 2, 3, 4, and the assump- tions of Theorem 7, we cannot have two or five periodic solutions by using the averaging theory and rescaling (5). Indeed, we observe that the only possibility to have five periodic solutions is that statement (1) of the appendix holds. However, this condition does not hold since D2 is negative for any a, r and d satisfying the hypotheses of Theorem 7. In addition, we note that the only condition that provides two periodic solutions according to the appendix is statement (7), and it needs D5 = D4 = E2 = 0. Nevertheless, we start con- sidering this condition and will find a Gröbner basis for the set of polynomials {D4, D5, E2} 1374 R.D. Euzébio and J. Llibre in the variables a, d and r. First, we will consider the factor D2 5 of D5 because D1 5 �= 0 by hypothesis. In this case, we obtain 55 polynomials in the Gröbner basis, where the first one is g(d, r) = r62ḡ(d, r) and ḡ(d, r) = 91125d12 − 438800d10r2 + 919360d8r4 − 1038464d6r6 + 660992d4r8 − 225280d2r10 + 32768r12. Since r �= 0 we must find a zero of ḡ(d, r). So solving ḡ(d, r) = 0 in the variable r we obtain six pairs of values ri ± = ±d √ qi 0, where i = 1, . . . , 6 and for each i the value qi 0 is the ith root of the polynomial q(x) = 91125 − 438800x + 919360x2 − 1038464x3 + 660992x4 − 225280x5 + 32768x6. This polynomial has only complex roots. Indeed, if we study the function q′(x) which is given by q′(x) = −438800 + 1838720x − 3115392x2 + 2643968x3 − 1126400x4 + 196608x5, we can apply the method described in the appendix to show that the respective values D2 and D5 related to q′(x) are negative and positive, respectively. Therefore, q′(x) has only one real root q′ 0 whose approximate value is q′ 0 = 106703 100000 . Additionally, the approximate value of q(x) evaluated in q′ 0 is q(q′ 0) = 16558687 10000 > 0. So if we observe that the coefficient of x6 is positive, it follows that the minimum value of q(x) is positive. Then each ri ±, i = 1, . . . , 6, is a complex value and the correspondent Gröbner basis has no zeros. Then we cannot have a common real zero of D4, D2 5 and E2 and so statement (7) in the appendix does not hold. Consequently, it is not possible to obtain two periodic solutions. The proof considering the factor D of D5 instead of D2 5 leads to the same polynomial g(d, r) and then we have the same conclusion. Proof of Theorem 8 In what follows we take m = 0, n1 = n3 = n4 = 1 and n2 > 1. So the vector field of system (8) becomes (y, −x + εry − εrx2y − εn2 ax3 − ε�x5 + εd cos t)T, and we obtain F1(t, x) = (0, ry − rx2y − �x5 + d cos t). So f (x0, y0) is written as f1(x0, y0) = ∫ 2π 0 −(sin t)(d cos t + r(y0 cos t − x0 sin t) − r(y0 cos t − x0 sin t) (x0 cos t + y0 sin t)2 − �(x0 cos t + y0 sin t)5) dt = −1 8 π(2rx0(−4 + x2 0 + y2 0) − 5�y0(x 2 0 + y2 0) 2), f2(x0, y0) = ∫ 2π 0 (cos t)(d cos t + r(y0 cos t − x0 sin t) − r(y0 cos t − x0 sin t) (x0 cos t + y0 sin t)2 − �(x0 cos t + y0 sin t)5) dt = 1 8 π(8d − 2ry0(−4 + x2 0 + y2 0) − 5�x0(x 2 0 + y2 0) 2). Again we will find zeros of f (x0, y0) through the roots of a Gröbner basis {bk(x0, y0)} of f (x0, y0). We will find a Gröbner basis of the functions f̄1,2 = ∓(8/π)f1,2. International Journal of Computer Mathematics 1375 An element of the Gröbner basis in the variables x0 and y0 is b1(x0, y0) = −4d3r5 + 1600dL2r5 + (2000d2�2r4 − 4d2r6 + 1600L2R6)y0 + (500d3�2r3 + 2400d�2r5)y2 0 + (900d2�2r4 + d2r6)y3 0 + 50d3�2r3y4 0 + 625d4�4y5 0. As before, b1(x0, y0) depends only on y0 and another element b2(x0, y0) of the Gröbner basis is linear on x0 with coefficients depending on y0, where the coefficient C of x0 in b2(x0, y0) is not zero by hypotheses. We will not present the expression of b2(x0, y0) in order to avoid large expressions. Then for each zero y∗ 0 of b1(y0) we have a second zero x∗ 0 through b2(x0, y0) related to y∗ 0. In order to apply the method described in [29] and the appendix for b1(y0), we will perform the translation y0 = ϕ − (2r3/d(125�)2). With this translation we obtain a new polynomial b∗ 1 in the form b∗ 1(ϕ) = ϕ5 + 4500�2r4 − 3r6 3125d2�4 ϕ3 + 2k1 k̄1d3�6 ϕ2 + 4r4k2 k̄2d4�8 ϕ − 4r5k3 k̄3d5�10 , where k1 = 156250d2�4r3 + 750000�4r5 − 13500�2r7 + r9, k̄1 = 390625, k2 = 39062500d2�6 − 390625�4(d2 − 80�2)r2 − 1500000�4r4 + 13500�2r6 + 3r8, k̄2 = 48828125, k3 = 48828125d4�6 − 781250d2(25000�8 − 500�6r2 + 3�4r4) + 2r4 (156250000�6 − 3750000�4r2 + 22500�2r4 + 9r6), k̄3 = 30517578125. We observe that the translation performed does not change the number or kind of the zeros of the original polynomial b1. Now we will apply the method of the appendix for b∗ 1. Indeed, we have D2 = N2, D3 = − 48r6 9765625d6�10 N3, D4 = − 16r12 3814697265625d12�18 N4, D5 = 256r20 1490116119384765625d20�26 N5M5. Therefore, if D5 is negative then b∗ 1 and consequently b1 have exactly three zeros, and hence function f = (f1, f2) has exactly three zeros (xi 0, yi 0), i = 1, 2, 3. On the other hand, if D5 is positive and one of the values D2, D3 or D4 is non-positive, then f has exactly one zero (x0 0, y0 0). Moreover, using rescaling (5) the matrix M is the same as the one of the proof of Theorem 7 taking a = 0. In addition, determinant � of M is written as �(x0, y0) = 1 ε2 (125λ2(x2 0 + y4 0) + 4ρ2(−4 + x2 0 + y2 0)(−4 + 3x2 0 + 3y2 0)). By hypothesis �(xi 0, yi 0) �= 0 for i = 0, 1, 2, 3, then we have the first part of Theorem 8 proved. 1376 R.D. Euzébio and J. Llibre Now we exhibit values of ρ, λ and δ for which we have either one or three periodic solutions. First, we consider (r, �, d) = (√ 31, 1, 1 125 √ 1547875 − 226851 √ 35 ) . With these values we obtain N3 = 1 625 (206837701817900 − 10419995302263 √ 35), N5 = 2(−649491478051728715458244 + 109639832554794027558525 √ 35) 48828125 , which are positive. On the other hand, considering (r, �, d) = ( √ 31/4, 1, 1) we get N5 = − 258261575 1048576 < 0. This ends the proof of Theorem 8. � Now we will show that if m = 0, n1 = n3 = n4 = 1 and n2 > 1, then we cannot have two periodic solutions by using the averaging theory described in Section 3 and rescaling (5). Indeed in order that function b∗ 1(ϕ), given in the proof of Theorem 8, has two periodic solutions, we need statement (7) of the appendix. We start considering the conditions N5 = N4 = E2 = 0 and will see that these equalities imply that N3 cannot be negative. We claim that N5 and M5 are factors of D5 as we saw in the last proof. Also we will not present the expression of E2 in order to avoid its large expression. It can be obtained by an algebraic manipulator such as Mathematica. Additionally, we note that N4 = N5 = 0 implies D4 = D5 = 0, which are necessary conditions for having exactly two real roots. We start replacing the condition N4 = N5 = E2 = 0 for other ones easier when r > 0. It is easy to see that if r = 0 we have N3 = 0 and consequently we cannot get two roots for b∗ 1. Now we obtain a Gröbner basis for the set of polynomials {N4, N5, E2} in the variables d, r and �. This basis has 12 elements, where 2 of them are g1(r, �) = 125 · 108�8 + 2 · 108�6r2 − 100000�4r4 − 10400�2r6 − 3r8, g2(d, r, �) = 25d4�2 − 110000d2�4 − 4000000�6 − 1400d2�2r2 + 80000�4r2 − 3d2r4 + 1200�2r4. In addition, the resultant of these two polynomials with respect to the variable r has the factor g3(d, �) = d8 − 1240000d6�2 − 269440000d4�4 + 6144 · 106d2�6 + 28672 · 108�8. Consequently, we will replace the problem of finding a zero of N4, N5 and E2 by the problem of finding a zero of g1, g2 and g3. Moreover, we suppose that we have a common zero between each gi and N3, i = 1, 2, 3. With this supposition and computing the resultant between g1 and N3, we obtain a new polynomial g4 in the variables r and d given by g4(r, d) = 1423828125d16 + 25974000000d14r2 − 47933388000000d12r4 − 478810245120000d10r6 + 1711832243200000d8r8 − 125591060860108800d6r10 − 1334616713986048000d4r12 + 5191865574606503936d2r14 − 4405603330689073152r16. However, the resultant between g3 and g4 is written as Kr128, where K is a positive constant. This polynomial cannot be zero since r ≥ 0. It means that when we consider each gi zero, N3 must be International Journal of Computer Mathematics 1377 positive or negative, exclusively. But taking the values (r, �, d) = (r0, �0, 1) we have gi ≡ 0, for each i = 1, 2, 3 and N3 is approximately 8738483 2500 > 0, where r0 and �0 are roots of the polynomials r̄(x) = −125 − 44400x2 + 184640x4 − 237568x6 + 86016x8, and l̄(x) = 1 − 1240000x2 − 269440000x4 + 6144000000x6 + 2867200000000x8, respectively. So N3 is positive when gi ≡ 0, for each i = 1, 2, 3, and then D3 is negative. Hence using the factor N5 of D5 we cannot obtain two periodic solutions. The proof that cannot exists two periodic solutions for the case that we consider the factor M5 of D5 instead of N5 is similar and we will omit it here. We remark that we cannot prove analytically the non-existence of five periodic solutions for system (8) considering m = 0, n1 = n3 = n4 = 1 and n2 > 1. Actually, using the algebraic manipulator Mathematica, we obtained evidences that this number of periodic solutions cannot happen, but an analytical treatment is not trivial. 2.2 Discussion of the results Now we discuss some aspects of the averaging method presented in Section 3 in order to find periodic solutions in system (4). In fact by using the averaging theory and rescaling (5), we shall see that we cannot prove that system (4) has periodic solutions except in the cases presented in Theorems 1–8. We start studying the possible values to the powers of ε in system (8). This is important because these powers play an important role in the averaging theory described in Section 3 because they determine the terms of order 1 that we are interested in. Indeed, we note that each one of the five different powers of ε in system (8) must be non-negative. Therefore, considering these powers as zero or positive, we have 32 possible combinations of these 5 different powers of ε. Actually some of them are not algebraically possible. Table 1 exhibits only the possible case. Observing Table 1 we see that in fact we have only 18 cases. Moreover, in order to apply the averaging theory for one of the 18 cases, we must integrate the equations of the non-perturbed part of the vector field of system (8). This is not a simple task because the major part of these equations are non-linear and non-autonomous. Indeed in each case from 1 to 8 and from 11 to 16 we could not integrate the equations even using the algebraic manipulator Mathematica or Maple. The non-integrable cases from Table 1 are listed as follows: Case 1: F0(t, x) = (y, −x + ry − ax3 − rx2y − �x5 + d cos t), Case 2: F0(t, x) = (y, −x + ry − ax3 − rx2y − �x5), Case 3: F0(t, x) = (y, −x + ry − ax3 − rx2y + d cos t), Case 4: F0(t, x) = (y, −x + ry − ax3 − rx2y), Case 5: F0(t, x) = (y, −x + ry − rx2y − �x5 + d cos t), Case 6: F0(t, x) = (y, −x + ry − rx2y − �x5), Case 7: F0(t, x) = (y, −x + ry − rx2y + d cos t), Case 8: F0(t, x) = (y, −x + ry − rx2y), Case 11: F0(t, x) = (y, −x − ax3 − �x5 + d cos t), Case 12: F0(t, x) = (y, −x − ax3 − �x5), Case 13: F0(t, x) = (y, −x − ax3 + d cos t), Case 14: F0(t, x) = (y, −x − ax3), Case 15: F0(t, x) = (y, −x − �x5 + d cos t), Case 16: F0(t, x) = (y, −x − �x5). 1378 R.D. Euzébio and J. Llibre Table 1. Here, when appear ∅, this means that the corresponding cases cannot occur. Possible combination of powers of ε for the system (9) n1 = 0 2m + n1 = 0 2m + n2 = 0 4m + n3 = 0 −m + n4 = 0 C1 −m + n4 > 0 C2 4m + n3 > 0 −m + n4 = 0 C3 −m + n4 > 0 C4 2m + n2 > 0 4m + n3 = 0 −m + n4 = 0 C5 −m + n4 > 0 C6 4m + n3 > 0 −m + n4 = 0 C7 −m + n4 > 0 C8 2m + n1 > 0 2m + n2 = 0 ∅ 2m + n2 > 0 4m + n3 = 0 ∅ 4m + n3 > 0 −m + n4 = 0 C9 −m + n4 > 0 C10 n1 > 0 2m + n1 = 0 ∅ 2m + n1 > 0 2m + n2 = 0 4m + n3 = 0 −m + n4 = 0 C11 −m + n4 > 0 C12 4m + n3 > 0 −m + n4 = 0 C13 −m + n4 > 0 C14 2m + n2 > 0 4m + n3 = 0 −m + n4 = 0 C15 −m + n4 > 0 C16 4m + n3 > 0 −m + n4 = 0 C17 −m + n4 > 0 C18 Note: The notation Ci, i = 1, . . . , 18 enumerate each possible case. In particular, we observe that cases 12, 14 and 16 turn the non-perturbed part of system (8) into a Hamiltonian system. In [3], Buica and Llibre provide a method to apply the averaging theory in planar systems when the system is Hamiltonian, but the expressions become also too much complicated and again we could not integrate the equations in cases 12, 14 and 16. In cases 9, 10 and 17, we can integrate the expressions but in these cases the hypotheses of the averaging method do not apply. In these cases, we have the following expressions for F0(t, x). Case 9: F0(t, x) = (y, −x + ry + d cos t), Case 10: F0(t, x) = (y, −x + ry), Case 17: F0(t, x) = (y, −x + d cos t). Solution x(t) in case 9 is C− e(1/2)t(r−√−4+r2) + C+ e(1/2)t(r+√−4+r2) − ( d r ) sin t, where C∓ = ∓2d + r((±r + √−4 + r2)x0 ∓ 2y0) 2r √−4 + r2 . In order that x(t) be periodic, we must have C− and C+ equal to zero. But these conditions are verified only if we take x0 = 0 and y0 = −d/r. So x0 and y0 are fixed and then the unperturbed system ẋ = F0(t, x) has no sub-manifold of periodic solutions. Consequently, we cannot apply the averaging theory in this case. On the other hand, in case 10 solution x(t) is written as C− e(1/2)t(r−√−4+r2) + C+ e(1/2)t(r+√−4+r2) − ( d r ) sin t, International Journal of Computer Mathematics 1379 where C∓ = (±r + √−4 + r2)x0 ∓ 2y0 2r √−4 + r2 . Again we need to choose C− and C+ equal to zero, but this happens only if (x0, y0) = (0, 0). So we cannot apply the averaging theory because (0, 0) is the equilibrium point of the non-perturbed part of system (8) in case 10. Finally, in case 17 solutions x(t) and y(t) have the form x(t) = x0 cos t + y0 sin t + d 2 t cos t, y(t) = y0 cos t − x0 sin t + d 2 (t cos t + sin t). It is immediate to note that these solutions are non-periodic because d is positive. Case 18 provides positive results and Theorems 1–8 are based on this case. In fact, we observe that in case 18 the vector field of system (8) is (y, −x + εn1 ry − ε2m+n1 rx2y − ε2m+n2 ax3 − ε4m+n3�x5 + ε−m+n4 d cos t)T. As we commented before, the periodic solutions of our problem correspond to the simple zeros of function (10). So in order to calculate the zeros of the function f (z) given in Equation (10), we must determine F1(t, x), where F1(t, x) is determined by the terms of order 1 on ε of the vector field of system (8). However, we see that the expression of F1(t, x) depends on the values m and ni, i = 1, 2, 3, 4. In fact, if m is positive and observing condition (7), the only terms of the vector field of order 1 on ε can be generated by powers n1 and −m + n4, where these two values can be 1 or greater than 1. This implies four possibilities. On the other hand, if m is zero the powers of ε in the vector field of system (8) depend on n1, n2, n3 and n4. Again, each one of these powers can be 1 greater than 1. So, if m is zero we have 16 possibilities for F1. Then case 18 has 20 subcases corresponding to the different possibilities of F1(t, x). These subcases are presented in Table 2. Table 2. Possible expressions for F1(t, x) when m2 + n2 2, m2 + n2 3, n4 − m and n1 are positives. Sub. Conditions Second coordinate of F1(t, x) 1 m = 0, ni = 1, i = 1, 2, 3, 4 ry − ax3 − rx2y − �x5 + d cos t 2 m = 0, ni = 1, i = 1, 2, 3, n4 > 1 ry − ax3 − rx2y − �x5 3 m = 0, ni = 1, i = 1, 2, 4, n3 > 1 ry − ax3 − rx2y + d cos t 4 m = 0, ni = 1, i = 1, 3, 4, n2 > 1 ry − rx2y − �x5 + d cos t 5 m = 0, ni = 1, i = 2, 3, 4, n1 > 1 −ax3 − �x5 + d cos t 6 m = 0, n1 = n2 = 1, n3, n4 > 1 ry − ax3 − rx2y 7 m = 0, n1 = n3 = 1, n2, n4 > 1 ry − rx2y − �x5 8 m = 0, n1 = n4 = 1, n2, n3 > 1 ry − rx2y + d cos t 9 m = 0, n2 = n3 = 1, n1, n4 > 1 −ax3 − �x5 10 m = 0, n2 = n4 = 1, n1, n3 > 1 −ax3 + d cos t 11 m = 0, n3 = n4 = 1, n1, n2 > 1 −�x5 + d cos t 12 m = 0, n1 = 1, ni > 1, i = 2, 3, 4 ry − rx2y 13 m = 0, n2 = 1, ni > 1, i = 1, 3, 4 −ax3 14 m = 0, n3 = 1, ni > 1, i = 1, 2, 4 −�x5 15 m > 0, n4 = 1, ni > 1, i = 1, 2, 3 d cos t 16 m > 0, ni > 1, i = 1, 2, 3, 4 0 17 m > 0, n1 = 1, −m + n4 = 1 ry + d cos t 18 m > 0, n1 = 1, −m + n4 > 1 ry 19 m > 0, n1 > 1, −m + n4 = 1 d cos t 20 m > 0, n1 > 1, −m + n4 > 1 0 Notes: We exhibit only the second coordinate of F1(t, x) because the first one has no terms depending on ε. 1380 R.D. Euzébio and J. Llibre Theorems 1–8 correspond to subcases 8, 10, 11, 17, 5, 3, 1 and 4, respectively. These are the only subcases of case 18 where the averaging method provides positive results. Indeed, in cases 15, 16, 19 and 20, function (10) does not have zeros. On the other hand, in cases 2, 6, 7, 13, 14, 18, the only zero of function (10) is (x0, y0) = (0, 0) that corresponds to the equilibrium point of the system, and consequently, in these subcases the system does not have periodic solutions. In cases 9 and 12, function (10) has real zeros different from (0, 0), but they are non-isolated and then the Jacobian of function (10) at the zero is zero, and consequently, the averaging theory cannot be applied. 3. The averaging theory for periodic solutions Now we present the basic results on the averaging theory of first order that we need to prove our results. Consider the problem of bifurcation of T-periodic solutions from differential systems of the form ẋ = F0(t, x) + εF1(t, x) + ε2R(t, x, ε), (12) with ε = 0 to ε �= 0 sufficiently small. Here the functions F0, F1 : R × � → R n and R : R × � × (−εf , εf ) → R n are C2, T-periodics in the first variable and � is an open subset of R n. One of the main assumptions is that the unperturbed system ẋ = F0(t, x) (13) has a manifold of periodic solutions. A solution to this problem is given using the averaging theory. Indeed, assume that there is an open set V with V̄ ⊂ D ⊂ � and such that for each z ∈ V̄ , x(·, z, 0) is T-periodic, where x(·, z, 0) is the solution of the unperturbed system (13) with x(0) = z. Answer to the problem of bifurcation of T-periodic solutions from x(·, z, 0) is given in the following theorem. Theorem 9 We assume that there exists an open set V with V̄ ⊂ D and such that for each z ∈ V̄ , x(·, z, 0) is T-periodic and consider the function f : V̄ → R n given by f (z) = ∫ T 0 Y−1(t, z)F1(t, x(t, z, 0)) dt. Then the following statements hold. (a) If there exists a ∈ V with f (a) = 0 and det((∂f /∂z)(a)) �= 0, then there exists a T-periodic solution ϕ(·, ε) of system (12) such that ϕ(0, ε) → a as ε → 0. (b) The type of stability of the periodic solution ϕ(·, ε) is given by the eigenvalues of the Jacobian matrix M = ((∂f /∂z)(a)). For a proof of Theorem 9(a), see [4, Corollary 1]. In fact, the result of Theorem 9 is a classical result due to Malkin [20] and Roseau [21]. For a shorter proof of Theorem 9(a), see [4]. For additional information on the averaging theory, see [23]. Disclosure statement No potential conflict of interest was reported by the authors. International Journal of Computer Mathematics 1381 Funding The first author is supported by the FAPESP-BRAZIL grants 2010/18015-6, 2012/05635-1 and 2013/25828-1. The second author is partially supported by MINECO/FEDER grants MTM2008-03437 and MTM2013-40998-P, an AGAUR grant number 2014SGR568, an ICREA Academia, grants FP7-PEOPLE-2012-IRSES 318999 and 316338, FEDER-UNAB-10-4E-378 and a CAPES grant 88881. 030454/2013-01 do Programa CSF-PVE. References [1] W.W. Adams and P. Loustaunau, An Introduction to Gröbner Bases, Graduate Studies in Mathematics, Vol. 3, American Mathematical Society, Providence, RI, 1994. [2] K. Bod, C. Edwards, J. Guckenheimer, S. Guharay, K. Hoffman, J. Hubbard, R. Oliva, and W. Weckesser, The forced Van der Pol equation II: Canards in the reduced system, SIAM J. Appl. Dyn. Syst. 2 (2003), pp. 570–608. [3] A. Buica and J. Llibre, Averaging methods for finding periodic orbits via Brouwer degree, Bull. Sci. Math. 128 (2004), pp. 7–22. [4] A. Buica, J.P. Françoise, and J. Llibre, Periodic solutions of nonlinear periodic differential systems with a small parameter, Comm. Pure Appl. Anal. 6 (2007), pp. 103–111. [5] A. Chen and G. Jiang, Periodic solution of the Duffing-Van der Pol oscillator by homotopy perturbation method, Int. J. Comput. Math. 87 (2010), pp. 2688–2696. [6] Y.M. Chen and J.K. Liu, Uniformly valid solution of limit cycle of the Duffing-van der Pol equation, Mech. Res. Comm. 36 (2009), pp. 845–850. [7] W.O. Criminale, T.L. Jackson, and P.W. Nelson, Limit cycle-strange attractor competition, Stud. Appl. Math. 112 (2004), pp. 133–160. [8] C. Egami and N. Hirano, Periodic solutions for forced van der Pol type equations, Math. Econ. 1264 (2001), pp. 159–172. [9] R.D. Euzébio and J. Llibre, Periodic solutions of El Ni no model through the Vallis differential system, Discrete Contin. Dyn. Syst. Ser. A 34 (2014), pp. 3455–3469. [10] T.H. Fay, The forced Van der Pol equation, Int. J. Math. Educ. Sci. Techol. 40 (2009), pp. 669–677. [11] J. Guchkenheimer, K. Hoffman, and W. Weckesser, The forced Van der Pol equation I: The slow flow and its bifurcations, SIAM J. Appl. Dyn. Syst. 2 (2003), pp. 1–35. [12] F.M.M. Kakmeni, S. Bowong, C. Tochawoua, and E. Kaptouom, Strange attractors and chaos control in a Duffing– Van der Pol oscillator with two external periodic forces, J. Sound Vib. 227 (2004), pp. 783–799. [13] F.M.M. Kakmeni, S. Bowong, C. Tchawoua, and E. Kaptouom, Chaos control and synchronization of a φ6-Van der Pol oscillator, Phys. Lett. A 322 (2004), pp. 305–323. [14] H.K. Leung, Synchronization dynamics of coupled van der Pol systems, Phys. A 321 (2003), pp. 248–255. [15] H. Li, Gröbner Bases in Ring Theory, World Scientific Publishing, Hackensack, NJ, 2011. [16] Z.-G. Li, W. Xu, and X.-Y. Zhang, Analysis of chaotic behavior in the extended Duffing-Van der Pol system subject to additive non-symmetry biharmonical excitation, Appl. Math. Comput. 183 (2006), pp. 858–871. [17] D. Liu and H. Yamaura, Chaos control of a φ6-Van der Pol oscillator driven by external excitation, Nonlinear Dyn. 68 (2012), pp. 95–105. [18] J. Llibre and C. Vidal, Periodic solutions of a periodic FitzHugh-Nagumo differential system, preprint (2014). [19] M. Ma, J. Zhou, and J. Cai, Practical synchronization of second-order nonautonomous systems with parameter mismatch and its applications, Nonlinear Dyn. 69 (2012), pp. 1285–1292. [20] I.G. Malkin, Some Problems of the Theory of Nonlinear Oscillations (Russian), Gosudarstv. Izdat. Tehn.-Teor. Lit., Moscow, 1956. [21] M. Roseau, Vibrations non linéaires et théorie de la stabilité (French), Springer Tracts in Natural Philosophy, Vol. 8, Springer-Verlag, Berlin-New York, 1966. [22] H. Salarieh and A. Alasty, Control of stochastic chaos using sliding mode method, J. Comput. Appl. Math. 225 (2009), pp. 135–145. [23] J.A. Sanders, F. VerhlstT, and J. Murdock, Averaging Method in Nonlinear Dynamical Systems, Applied Mathe- matical Science, Vol. 59, Springer, New York, 2007. [24] M.S. Siewe, F.M.M. Kakmeni, and C. Tchawoua, Resonant oscillation and homoclinic bifurcation in a φ6-Van der Pol oscillator, Chaos Solitons and Fractals 21 (2004), pp. 841–853. [25] W. Szemplinska-Stupnicka and J. Rudowski, Neimark bifurcation, almost-periodicity and chaos in the forced van der Pol-Duffing system in the neighbourhood of the principal resonance, Phys. Lett. A 192 (1994), pp. 201–206. [26] R. Tchoukuegno, B.R.N. Nbendjo, and P. Woafo, Resonant oscillations and fractal basin boundaries of a particle in a φ6-potential, Phys. A 304 (2002), pp. 362–378. [27] Y. Ueda and N. Akamatsu, Chaotically transitional phenomena in the forced negative-resistance oscillator, IEEE Trans. Circuits Syst. 28 (1981), pp. 217–224. [28] A. Venkatesan and M. Lakshmanan, Bifurcation and chaos in the double-well Duffing–Van der Pol oscillator: Numerical and analytical studies, Phys. Rev. E 56 (1997), pp. 6321–6330. [29] L. Yang, Recent advances on determining the number of real roots of parametric polynomials, J. Symbolic Comput. 28 (1999), pp. 225–242. 1382 R.D. Euzébio and J. Llibre [30] X. Yang, W. Xu, and Z. Sun, Effect of Gaussian white noise on the dynamical behaviours of an extended Duffing-Van der Pol oscillator, Int. J. Bifur. Chaos Appl. Sci. Eng. 16 (2006), pp. 2587–2600. [31] X. Yang, W. Xu, and Z. Sun, Effect of bounded noise on the chaotic motion of a Duffing Van der pol oscillator in a φ6 potential, Chaos Solitons Fractals 27 (2006), pp. 778–788. Appendix. Root classification for a quintic polynomial In this section, we present a brief summary of the results about the number and multiplicities of the real/complex roots for a quintic polynomial with arbitrary coefficients presented in [29]. Indeed, consider the polynomial P(x) = x5 + px3 + qx2 + ux + v. So, the following table gives the number of real and complex roots and multiplicities of roots of P(x) in all cases: (1) D5 > 0 ∧ D4 > 0 ∧ D3 > 0 ∧ D2 > 0 {1, 1, 1, 1, 1} (2) D5 > 0 ∧ (D4 ≤ 0 ∨ D3 ≤ 0 ∨ D2 ≤ 0) {1} (3) D5 < 0 {1, 1, 1} (4) D5 = 0 ∧ D4 > 0 {2, 1, 1, 1} (5) D5 = 0 ∧ D4 < 0 {2, 1} (6) D5 = 0 ∧ D4 = 0 ∧ D3 > 0 ∧ E2 �= 0 {2, 2, 1} (7) D5 = 0 ∧ D4 = 0 ∧ D3 > 0 ∧ E2 = 0 {3, 1, 1} (8) D5 = 0 ∧ D4 = 0 ∧ D3 < 0 ∧ E2 �= 0 {1} (9) D5 = 0 ∧ D4 = 0 ∧ D3 < 0 ∧ E2 = 0 {3} (10) D5 = 0 ∧ D4 = 0 ∧ D3 = 0 ∧ D2 �= 0 ∧ F2 �= 0 {3, 2} (11) D5 = 0 ∧ D4 = 0 ∧ D3 = 0 ∧ D2 �= 0 ∧ F2 = 0 {4, 1} (12) D5 = 0 ∧ D4 = 0 ∧ D3 = 0 ∧ D2 = 0 {5} where D2 = −p, D3 = −12p3 − 45q2 + 40pu, D4 = −4p3q2 − 27q4 + 12p4u + 117pq2u − 88p2u2 + 160u3 − 40p2qv − 300quv + 125pv2, D5 = −4p3q2u2 − 27q4u2 + 16p4u3 + 144pq2u3 − 128p2u4 + 256u5 + 16p3q3v + 108q5v − 72p4quv − 630pq3uv + 560p2qu2v − 1600qu3v + 108p5v2 + 825p2q2v2 − 900p3uv2 + 2250q2uv2 + 2000pu2v2 − 3750pqv3 + 3125v4, E2 = 16p4q2 − 48p5u + 60p2q2u + 160p3u2 + 900q2u2 − 1100p3qv − 3375q3v + 1500pquv + 625p2v2, F2 = 3q2 − 8pu. The polynomials Di, i = 2, 3, 4, 5, E2 and F2 form a discriminant system which is sufficient for the classification of roots of the polynomial P(x), which is described by the right column of the table. For instance, {1, 1, 1} means three real simple roots and a pair of complex roots, and {3, 1, 1} means a real triple root plus two real simple roots. 1. Introduction 1.1. Setting the problem 1.2. Statement of the main results 2. Proof and discussion of the results 2.1. Proof of the results 2.2. Discussion of the results 3. The averaging theory for periodic solutions Disclosure statement Funding References Appendix