On the Qualitative Behavior of a Class of Generalized Liénard Planar Systems

We study the problem of existence/nonexistence of limit cycles for a class of Liénard generalized differential systems in which, differently from the most investigated case, the function F depends not only on x but also on the y-variable. In this framework, some new results are presented, starting from a case study which, actually, already exhibits the most significant properties. In particular, the so-called “superlinear case” presents some new phenomena of escaping orbits which will be discussed in detail.


Introduction
In this paper we study the qualitative behavior of the trajectories of a planar system of the form ẋ = y − F(x, y) y = −g(x).  [26] in 1928 for the study of the equationẍ Although there is an enormous literature concerning (1.1) and its generalizations as (see the monographs [14,24,34,44,47], as well as the articles [8,15,16,25,27,29,[36][37][38][39][40] and the quotations therein) or ẋ = h(y) − F(x) with the additional assumption h(0) = 0, h(y)y > 0 for y = 0, (see the monographs [44,47] as well as the articles [1,2,13,[19][20][21][22]42,43,45,46] and the references therein), as far as we know, no systematic study has been performed for (S). This because the classical assumptions on F, g, h guarantee that the origin is the unique equilibrium point and the trajectories move clockwise around it in the plane. For equation (S), additional conditions on F(x, y) would be required in order to obtain a similar behavior. Moreover, the problem of possible blow-up in finite forward time of the solutions should be considered. On the other hand, the presence of a term of the form F(x, y) in the first equation gives some more flexibility, if one looks for examples of rich dynamics, for instance, exhibiting a prescribed number of limit cycles. Previous contributions on the study of system (S) can be found in the work of Cartwright and Swinnerton-Dyer [6,7] concerning the boundedness of the solutions of non-autonomous equations of the forṁ and, for the autonomous case, in [3,32], concerning the uniqueness of limit cycles (see also [30], where a special case of (S), with peculiar assumptions, was considered) and in [31], dealing with a φ-Laplacian Liénard equation which, studied in the phase-plane, leads to a system of the form (S) with x and y inverted. The plan of the paper is the following. In Sect. 2 we discuss some basic facts of equation (S) related to the use of the energy of the associated Duffing equations as a Lyapunov function. Moreover, we present a first result of multiple limit cycles which shows the above mentioned flexibility in the use of the function F(x, y).
The other sections are focused on the study of the case

F(x, y) = λB(y)A(x)
where A(x) satisfies the standard assumptions on F(x) in the classical case. In particular, In Sect. 3, a case study is analyzed for the choice which already exhibits the mean features of this problem that, as far as we know, have not been investigated before. Section 4 presents the main results, namely a theorem of existence/nonexistence of limit cycles for system (S) (see Theorem 6) as well as results about existence and uniqueness (see Theorem 4 and Theorem 5).

Preliminary Results
Throughout the paper the following assumptions are made on system ẋ = y − F(x, y) y = −g(x). (S) We suppose that F : R × R → R and g : R → R are locally Lipschitz continuous functions, in order to guarantee the uniqueness of the solutions for the associated initial value problems. We also assume g(0) = 0, g(x)x > 0 for x = 0 and F 0 (y) := y − F(0, y), vanishes only at y = 0, therefore the origin is the only singular point of (S). When F ≡ 0, our system reduces to the classical Duffing one, namely It is well known (cf. [14,24]) that system (D) has the Hamiltonian structurė and the origin is a global center if and only if G(x) → +∞ as |x| → +∞. Taking the energy H as a Lyapunov function for system (S), we obtain, for its time-derivative along the trajectories, the expressionḢ (x, y) = −g(x)F(x, y), (2.1) which generalizes the classical one for system (1.2). Therefore, the orbits of (S) enter the ones of (D) which are the level lines of H , when g(x)F(x, y) > 0, while they exit when this product is negative. In view of this property and in the light of the classical Poincaré example (see [29, §3.3]), a first result giving examples of a prescribed number of periodic solutions is the following.
system (S) has exactly n limit cycles which become arbitrarily large as n grows. Conversely, if G(x) → +∞ as |x| → +∞, given any energy level E with 0 < E < min{G(−∞), G(+∞)}, then, for system (S) has exactly n limit cycles which become arbitrarily close to the origin as n grows.
Proof The proof is straightforward using the above mentioned properties of the orbits following from (2.1).
A similar construction has been recently proposed in [5] for the relativistic and the curvature generalized Liénard equations Clearly, it is also possible to produce an arbitrarily large number of small limit cycles approaching the origin even when G(x) → +∞ as |x| → +∞. In this case we can start from any positive energy level E. As a side remark, we observe that this result shows the flexibility of working with systems like (S), because analogous results of non-perturbative nature (like those referring to small-amplitude limit cycles) are very hard to obtain for the classical Liénard systems (see [44, §7]). At this point, following again (2.1) we assume at first the following conditions. There exist α < 0 < β such that g(x)F(x, y) < 0, for x ∈ ]α, β[ and x = 0 and for all y = 0.
These assumptions fit with with the classical conditions considered in the Liénard systems. Moreover, we also assume G(x) → +∞ as |x| → +∞.
As a consequence, the origin is a source and outside the strip [α, β] × R, trajectories of system (S) are guided by those of system (D). Hence, as in the classical case, the problem is that of controlling the orbits in the strip. However, the phase-portrait in the strip is now more complicated due to the possible presence of two phenomena, namely: (a) the existence of several branches of the vertical isocline, that is the set y = F(x, y), (b) the existence of blow-up solutions, given by possible superlinear terms in y.
In order to control the above defined situations, we focus our attention to the case of with A such that There exist α, β with α < 0 < β such that A(x) is strictly increasing for x ≤ α and x ≥ β and, moreover, For the function B assume that We observe again that such choice seems natural, because, when B ≡ 1, the obtained Liénard system presents a nonlinearity with a typical cubic-shaped curve, which is common in most papers and is crucial in order to prove uniqueness of the limit cycle (see [14,17,18,24,32,[39][40][41] and the references quoted therein). On the other hand, the assumptions on B(y) allow to consider power-like nonlinearities as B(y) = |y| p . This will be our starting point and will be treated in the next section.

A Case Study
In this section we focus our attention on the equation with λ > 0 and p > 0.
Here we are mainly interested in the case p ≥ 1. This, because for 0 < p < 1, |y| p is sublinear and therefore the geometry is not very far from the classical one. Moreover, the equation of the vertical isocline may be split in a simple way, namely the x-axis and a graph of the function y =F(x), so that the situation is similar to the one of the classical case and the problem can be treated in a similar way.
We shall discuss separately the two cases study p = 1 and p > 1.

The case p = 1.
We start with an analysis of the system for p = 1. In this case, for simplicity, we write (H λ ) = (H 1,λ ). System (H λ ) can be split as two systems of the form and Both systems have an Hamiltonian structure, with Hamiltonians Using the fact that max{( , we find that for λ > 0, the equation λ(x 2 − 1)x − 1 = 0 has a unique zero if and only if 0 < λ < λ * := 3 √ 3/2. In this case, we denote by x = a λ such unique solution, with a λ > 1. Hence we have that (a) (b) Fig. 1 The functions Φ ± λ are the primitives of The graphs are produced using [48] Φ + λ is well defined on (−∞, a λ [ and, without loss of generality, we can take, among all the primitives of Figure 1a shows the typical where −a λ < −1 is the unique real root of 1 + λ(x 2 − 1)x. Figure 1b shows the above symmetry.
With reference to the graph of Φ + λ we can describe the dynamics of system (H + λ ) for 0 < λ < λ * . Let c λ := Φ + λ (−∞). For every energy level c with 0 < c < c λ , the energy level line is a closed curve in the phase-plane (x, y) which is also a periodic orbit of the system (H + λ ), considered as a differential system defined on the whole plane. This orbit intersects the x-axis on two points u ± (λ, c) where On the other hand, for c ≥ c λ , the energy level line E + λ (x, y) = c is an unbounded curve, symmetric with respect to the x-axis which, restricted to y ≥ 0, is the graph of a function with an asymptote at y = √ 2(c − c λ ), increasing on (−∞, 0] and decreasing on [0, u + (λ, c)], where 0 < u + (λ, c) < a λ , with Φ + λ (u + (λ, c)) = c. Analogously, with reference to the graph of Φ − λ we can describe the dynamics of system (H − λ ) for 0 < λ < λ + . Let d λ := Φ − λ (+∞). By symmetry, we know that d λ = c λ . For every energy level d with 0 < d < d λ , the energy level line is a closed curve in the phase-plane (x, y) which is also a periodic orbit of the system (H − λ ), considered as a differential system defined on the whole plane. This orbit intersects the x-axis On the other hand, for d ≥ d λ , the energy level line E − λ (x, y) = d is an unbounded curve, symmetric with respect to the x-axis which, restricted to y ≤ 0, is the graph of a function with an asymptote at y = Now, to get the global dynamics of (H λ ), we can glue the orbits of (H + λ ) on the upper half-plane with those of (H − λ ) on the lower half-plane. Let us consider system (H λ ) for 0 < λ < λ * . Without loss of generality, we can assume to start from an initial point P 0 := (x 0 , y 0 ) with y 0 > 0. First of all, we observe that if x 0 = a λ thenẋ = 0 andẏ = constant = −x 0 < 0, so that the trajectory hits the x-axis in finite time. For convenience, we suppose now that y 0 > 0 is sufficiently large and x 0 = a λ . It is easy to check that in this situation, the trajectory moves fast close to the vertical line x = a λ and then goes down to the x-axis almost parallel to the line x = a λ . Thus, we can suppose that we enter in the lower half-plane from a point of the x-axis P 1 := (x 1 , 0) with x 1 close to a λ . Let now d 1 := Φ − λ (x 1 ). We have 0 < d 1 < d λ and now we follow in the lower half-plane the lower part of the closed orbit of energy d 1 we get a new point P 2 := (x 2 , 0) on the orbit of (H λ ) where the x-axis is hit again. From P 2 we continue with the orbit path of (H + λ ) which is the upper-half of a closed orbit which hits again the x-axis at a point P 3 := (x 3 , 0) with x 3 = u + (λ, c) for c := Φ + λ (x 2 ). Continuing in this manner, by induction, we obtain a sequence of points P i = (x i , 0) on the x-axis with x 1 > x 3 > x 5 , · · · > 0 and x 2 < x 4 < · · · < 0 which shows that the semi-orbit γ + (P 0 ) is winding. Since we already know that the orbits unwind from the origin, we conclude, in virtue of the Poincaré-Bendixson theorem, that the orbit γ + (P 0 ) tends to a limit cycle Γ . The intersection of the limit cycle Γ with the x-axis is made by two symmetric points P − = (−x, 0) and P + := (x, 0) with P + being the limit of the subsequence (P k ) k for k odd and P − the limit of the subsequence (P k ) k for k even. The uniqueness of the limit cycle can be deduced by the fact that the graphs of Φ + λ and Φ − λ intersect at a unique point for x > 0 (and a unique symmetric one for x < 0).
We observe that the above analytical argument works only in the case p = 1 because only in this case we can exploit the symmetry and the Hamiltonian structure for the two phase-portraits which are glued together. In order to overcome this difficulty, we present now a different and more general argument to prove the existence and the uniqueness of a limit cycle which can be extended to more general situations. This approach was actually already present in the above mentioned work of Liénard [26] and then by Levinson and Smith [25], and Sansone [33]. Such approach was also treated in [3,17,18,32,40,41]. The idea is to observe that the integral of the time-derivativeḢ (cf. (2.1)) along a closed trajectory must be zero and then to show that this can happen only once. A detailed proof for F(x, y) = F(x) may be found in [41].
Here we give a sketch of the proof, because the argument is similar. At first we observe that, from (2.1), we haveḢ = −λ|y|(x 2 − 1)x 2 and therefore the energy function H is increasing inside the strip ] − 1, 1[ ×R and hence, due to the symmetry property of system (H λ ), all the limit cycles must intersects the vertical lines x = −1 and x = 1. Assume (by contradiction) that there are two limit cycles Γ 1 and Γ 2 (with Γ 1 Fig. 2 From the phase-portrait of (H λ ), the symmetry between (H + λ ) for y ≥ 0 and (H − λ ) for y ≤ 0 appears evident. The simulation corresponds to the case λ = 1. The graph is produced using [49] in the internal region of Γ 2 ) and compute the line integrals Γ i −g(x(t))F(x(t), y(t)) dt for i = 1, 2. Splitting these integrals along the orbit paths inside and outside the strip ]−1, 1[ ×R, we observe the following: (1) Inside the strip, the integral takes the form Notice that the denominators in the above integrals do not vanish for 0 < λ < λ * = 3 √ 3/2. As y does not appear in the integrals, the variation of energy in the strip is the same for both limit cycles.
(2) outside the strip, the integral takes the form Now one must take into account the intervals of integration, in the y-variable, which are different for the two limit cycles. However, we observe that, since the function x 3 − x is monotone increasing for |x| ≥ 1 (positive for x > 1 and negative for x < 1), we can follows step by step the proof in [41, Section 2] and, in virtue of the clockwise orientation of the limit cycles, we find that the contribution to the energy along the external limit cycle is less than the contribution in the internal one.
As a final balance, we get that which is a contradiction and shows also that the unique limit cycle is stable (Figs. 2 and 3). Summarizing this discussion, we have obtained the following result.

Theorem 1 The system
has a unique limit cycle for 0 < λ < λ * = 3 √ 3/2, while, for λ > λ * there are no limit cycles. The figures represent the dynamics for system (H 1,λ ), for λ less than and λ greater than the critical value λ * . The graphs are produced using [49] Remark 1 Following the dynamics for λ λ * we can say that, at this value, the limit cycle collapses at ±∞ in a vertical strip. In this light we produce simulations for λ = 2 < λ * < λ = 3 (see Fig. 4 below). Therefore, this result can be viewed as a bifurcation from infinity in the strip ] − 1, 1[ and this clearly does not appear in the classical case of the Van der Pol equation.

The Case p > 1.
We consider now system (H p,λ ) for p > 1. In order to understand the dynamical behavior of such class of systems, we deal, at first, with the particular case The infinite-isocline and the zero-isocline for system (H 2,1 ), with the direction of the vector field. The graph is produced using [48] since the other systems exhibit the same dynamical features. Now, the vertical isocline will play a crucial role. It is given by the union of the x-axis with the graph of .
Standard phase-plane arguments show that the trajectories, as before, are symmetric with respect to the origin and the directions of the vector field are the following (see Fig. 5).
The symmetric regions are positively invariant and play a crucial role because all the trajectories entering such regions become unbounded in the y-component. Notice that, however, there is no blow-up in forward time and the solutions entering V 1 or V 2 are defined on R + . Indeed, if (x(t), y(t)) ∈ V i (for i = 1, 2), for t ≥ 0, then |x(t)| ≤ 1 and hence, from the second equation in system (H p,λ ) we have |y (t)| ≤ 1.
Unboundedness in backward time occur also in the negatively invariant regions Fig. 5), but this latter fact does not involve extra difficulties (and sometimes it may help) in the search of the limit cycle.
An inspection of the vector field in a neighborhood of the origin and observing also thaṫ H (x, y) = −λy 2 x 2 (x 2 −1) > 0 for 0 < |x| < 1 and y = 0, shows that the origin is a source and therefore small orbits are unwinding. However, for λ > 0 large, all these orbits will enter the positively invariant regions V 1 or V 2 and ultimately will be unbounded, so that there are no limit cycles (the precise details and estimates justifying our assertion on the nonexistence of limit cycles are given in Step 1 of the proof of Theorem 2, below). On the other hand, for λ small enough, the are trajectories exiting the negatively invariant region W 1 that reach the negative y-axis, without entering the "forbidden region" V 1 . Such trajectories eventually intersect the x-axis at some x < 0. Due to the symmetry of the system, any such trajectory is bounded above in the upper plane by a symmetric one and therefore is winding. As usual, the Poincaré-Bendixson theorem guarantees the existence of at least a limit cycle. The fact Fig. 6 The movement of the infinite-isocline when λ varies. The graph represent (from the external to the internal ones) three different isoclines for λ = 1, 2, 3 and show the monotonicity. The graphs are produced using [48] that this behavior of the orbits is possible for λ is small, as well the fact that the trajectories become unbounded for λ large, is granted by the fact that the infinite isoclines in the strip ] − 1, 1[×R monotonically move with respect to the parameter λ (see Fig. 6).
In view of the above discussion, now we are in position to establish the following result.
Proof In the introductory part preceding the statement of the theorem, we have already described the qualitative analysis leading to the strategy of the proof. We provide now the technical estimates which justify our previous assertions. To this end, we split the proof in three steps.
Step 1. Nonexistence of limit cycles for λ large. Considering the previously introduced energy function H (x, y) = 1 2 y 2 + x 2 ) associated to the (linear) autonomous Duffing equatioṅ x = y,ẏ = −x, we find thatḢ HenceḢ (x, y) ≥ 0, with ∇ H (x, y) = 0 on all the points of the circumference H (x, y) = 1/2. This implies that the disc D := {(x, y) : H (x, y) ≤ 1/2} (the closed disc of center the origin and radius one) is a negatively invariant region for the dynamical system associated with (H 2,λ ) (for each λ > 0). Hence the Lyapunov-LaSalle theory (cf. [23]) guarantees that for each initial point P 0 ∈ D the negative semi-orbit γ − (P 0 ) tends to the origin (for t → −∞). In other words the disc D is contained in the region of repulsivity of the origin. Now, it is sufficient to take λ > 0 sufficiently large so that at least one of the regions V 1 or V 2 (in our case both the regions, by the symmetry of the system) intersects the disc D (see Fig. 7 for an illustration of this fact). Clearly, if we take any initial point , we have that the orbit γ (P 0 ) tends to the origin in the negative time, while its positive semi-orbit enters V and is unbounded in the future. By construction, the orbit γ (P 0 ) prevents the possibility of the existence of a closed orbit (which should surround the origin). At this point the Boincaré-Bendixson theory implies that all the nontrivial positive trajectories must be unbounded. In this manner we have proved that there existsλ > 0 such that for each λ >λ all the nontrivial trajectories are ultimately unbounded.
Step 2. Existence of limit cycles for λ small. Let 7 Taking λ > 0 sufficiently large we move down the "forbidden" region V 1 and up the symmetric one V 2 (in dark) until they cross the disc D. The graph is produced using [48] be the coordinates of the minimum point of 1/(x 3 − x) on the interval ] − 1, 0[ . Consider also the line which is below the region V 1 and has y m /2λ as a minimum value. By construction, we have thatẋ .
Consider the initial point Clearly y 0 > 0 for λ small, namely, 0 < λ < y m /4. The slope of the trajectory departing from P 0 satisfies the relation hence, an elementary and standard estimate shows that the positive trajectory departing from P 0 remains in the region R and crosses the y-axis at some point P 1 = (y 1 , 0) with y 0 < y 1 < y 0 + 2. Consider now the point P 2 = (x 1 , y 1 ) of intersection of the horizontal line y = y 1 with the boundary of the region W 1 , given by the equation y = 1/(λ(x 3 − x)) for x > 1 and observe that the segment P 1 P 2 is crossed by the vector field from the exterior to the interior, asẋ ≥ 0 andẏ < 0 (for x > 0). Consider now the negative semi-orbit departing from P 0 and let C(P 0 ) be the circumference of center the origin and passing through P 0 of equation SinceḢ (x, y) ≤ 0 for all x ≤ −1, we have that the region outside the circumference C(P 0 ) in the half-plane x ≤ −1 is negatively invariant is for the dynamical system associated with (H 2,λ ) (for each λ > 0), with respect to the open set x < −1. Now we are in position to consider the following Claim For λ > 0 sufficiently small the circumference C(P 0 ) crosses the curve y = We observe that we look for an intersection point in the third quadrant, therefore we study the system ⎧ ⎨ looking for a solution x in the domain Clearly, from system (3.3) we get the equation Such equation has always solutions in the desired range for λ sufficiently small. It is enough to observe that and then take a suitable (large in absolute value)x < −1 such that η := (x 3 −x)y m /2 < −1 and observe that and therefore the claim is proved. The claim we have just proved guarantees that the negative semi-trajectory departing from P 0 crosses the line y = 1/(λ(x 3 − x)) at some point with P −1 = (x −1 , y −1 ) with x −1 < −1 and y −1 < 0.
To conclude the proof, we just repeat symmetrically the same construction starting from the point Q 0 = (1, −y 0 ) and obtain the corresponding (and symmetric) points Q 1 , Q 2 and Q −1 . In this manner, we have determine a piecewise smooth Jordan curve Υ defined as follows: from P −1 to P 0 and from P 0 to P 1 along an orbit-path of the trajectory through P 0 , next, the segment P 1 P 2 , then the arc of curve y = 1/(λ(x 3 − 1)) for x > 1 from P 2 to Q −1 and, finally, the symmetric orbit-path from Q −1 to Q 1 (passing through Q 0 ), the segment Q 1 Q 2 and the arc of curve y = 1/(λ(x 3 − 1)) for x < −1 from Q 2 to P −1 . The curve Υ is the boundary of a positively invariant region containing the origin as the unique equilibrium point, which is also a source. Therefore, the Poincaré-Bendixson theorem guarantees the existence of at least a limit cycle bounded by Υ . In this manner we have proved that there existsλ > 0 such that, for each 0 < λ <λ, system (H 2,λ ) has at least a limit cycle. Figure 8 illustrates the geometry involved in the proof of Step 2.
Step 3. The existence of λ * 2 . This last assertion follows by defining λ * 2 := sup{k > 0 : (H 2,λ ) has a bounded nontrivial orbit ∀ λ : 0 < λ < k} and observing that if for a certain λ 1 all the solutions of (H 2,λ ) are unbounded, then the same holds for any λ 2 > λ 1 . This because the unbounded solutions of (H 2,λ ) either eventually enter in the regions V 1 and V 2 or they correspond to the separatrices contained, respectively in ] − 1, 0[ and in ]0, 1[ and very close (below/above) to these regions. Now, increasing λ has the effect of moving down the region V 1 and up the region V 2 (and the corresponding separatrices). This in turn implies that all the obits which were escaping to infinity for λ 1 are also escaping to infinity for λ 2 . This completes the proof of our result. Fig. 8 Taking λ > 0 sufficiently small we move up the "forbidden" region V 1 and down the symmetric one V 2 (in dark) and an invariant region bounded by W 1 and W 1 can be constructed. The graph is produced using [48] As a comment to this lengthy proof, we observe that we choose to describe in detail all the steps, because the same argument will work in the general case and this will avoid tedious computations in the proof of Theorem 6, where the geometry is exactly the same except for the symmetry (Figs. 9 and 11.).

Remark 3
The same argument, with minor changes, allows to extend the result of Theorem 2 to system (H p,λ ) for every p > 1. In this case, we find a constant λ * p such that the system has a limit cycle if and only if 0 < λ < λ * p . The same result applies also to more general functions A(x) (see also Theorem 6). Example 2, considered at the end of the article, shows the applicability of our approach to other situations.

The problem of the uniqueness of the limit cycle
We conclude our analysis by discussing the problem of the uniqueness of the limit cycle for p > 1. Although there is a numerical evidence of the uniqueness, unfortunately, the previous proof does not work. It is interesting to point out the issues even because this will give a proof for the case 0 < p < 1.
Let p > 1. As before, a computation ofḢ = −λ|y| p (x 2 − 1)x 2 and the symmetry property of system (H p,λ ), show that the origin is a source and all the limit cycles must intersects the vertical lines x = −1 and x = 1. Again, assume (by contradiction) that there are two limit cycles Γ 1 and Γ 2 (with Γ 1 in the internal region of Γ 2 ) and compute the line integrals Γ i −g(x(t))F(x(t)) dt for i = 1, 2. Splitting these integrals along the orbit paths inside and outside the strip ] − 1, 1[ ×R, we observe, in the same way, the following: (a) (b) Fig. 11 The figures represent the behavior of the components x(t) (left panel) and y(t) (right panel) for system (H 2,λ ) with λ = 1.475 as in the right panel of Fig. 10. The graphs are produced using [48] Outside the strip, the integral takes the form Arguing as before, we find that the contribution to the energy along the external limit cycle is less than the contribution of the internal one. This is true not only for p > 1, but for every p > 0. Inside the strip, the integral takes the form Clearly, a limit cycle cannot intersect the vertical isocline and enter in the "forbidden" regions that before were indicated by V 1 and V 2 for p = 2. Hence, for 0 < λ < λ * p , the denominators in the above integrals do not vanish. Let us consider the first of the above integrals and write it as (for the second integral the situation is symmetric and so we do not repeat the discussion). However, now a problem appears, because, as one can see, for p > 1 the integrand is increasing with respect to y. This means that in the strip ] − 1, 1[×R + , the part of the outer limit cycle gains more energy than the part of the inner limit cycle. This does not allow to repeat the above proof and show that the complete integral along the outer limit cycle is less than the integral along the inner one. This is also evident from the numerical simulations, where we see that in the strip the limit cycles tend to have an arbitrary large amplitude and the function H (x, y) is actually the square of the distance from the origin.

The case 0 < p < 1.
At first we notice that even if in this case B(y), and therefore F(x, y), is not locally Lipschitz continuous at y = 0, nevertheless there is still uniqueness of the trajectories. Indeed, near the origin (which could be the only point in which the uniqueness might fail) standard phase inspection shows that trajectories are clockwise and the origin is a source. Concerning the existence of limit cycles, we stress that the proof for the case p > 1 does not work because a crucial point in that case was the fact there were trajectories coming from the negatively invariant regions W 1 and W 2 and such regions are not present in this case. To overcome this question, we use a different (still geometrical) approach.
From (H p,λ ) we consider the system It is worth to note that this change of variables is not smooth at w = 0. However, being the origin a source, this is not a problem and the uniqueness of the solutions is always granted (the uniqueness also follows from more general results in [12] regarding differential equations with even a discontinuous right-hand side, provided that the vector field presents a surface of discontinuity which is crossed transversally, as is in our case). Finally, from (3.10), we consider the associated system It is well known (see, for instance, [2]) that such a system has at least a limit cycle. As the slopes for (3.10) and (3.11) are the same, due to the fact that the factor L(w) simplifies, the two systems have the same phase-portrait. However, this argument fails at the x-line, where L(w) = 0. Indeed, one may observe that (3.10) as w = 0 as a whole line of critical points. But this is not a problem if we compare the two system for w > 0 (that is y > 0) and w < 0 (that is y < 0) respectively. More in detail, let Γ be a limit cycle of system (3.11). We split Γ in two arcs, Γ + for w > 0 and Γ − for w < 0, ignoring the two points where Γ intersects the x-axis. Clearly, the same arcs Γ ± are orbits of system (3.10). At this point we apply on these arcs the transformation (x, w) → (x, y) considered in (3.9), obtaining two new arcs Γ ± 1 which are orbits of system (3.8). Cancelling |y| p in the second equation of (3.8), we are back to system (H p,λ ). Finally, we glue the closures of Γ + 1 and Γ − 1 (in a smooth way due to the fact that the vector field on the x-axis is vertical), obtaining the desired limit cycle for (H p,λ ).
Let us consider now the problem of the uniqueness of the limit cycle. We just observe that the previous computation, which was failing or the case p > 1, now works for 0 < p < 1. Indeed, the integral (3.7) in decreasing in the strip [−1, 1] (and a symmetric argument works for y < 0). Hence, the computations produced for (3.4)-(3.5)-(3.6) can be repeated and we have the following result.
A comment to Theorem 3 is given by the following Fig. 12 concerning the case p = 1/2 and λ = 2.

The General Case
As in the previous section, we consider the case and A, B, g satisfying regularity conditions in order that the uniqueness of the solutions for the initial value problems is guaranteed.
In the light of the previous section, we consider two cases for the function B(y), namely, the sublinear case and the superlinear one.

The sublinear case: B(y)/y → 0 for y → ±∞
As already discussed, if B(y) is not locally Lipschitz continuous at y = 0, this is not a problem because the origin is a source in virtue of the sign of A(x) in a neighborhood of x = 0. As before, we replace system (S1) with and transform system (4.1) to The orbit-paths of (4.2) for w = 0 (that is y = 0) are the same of those of Observe that the sublinear condition implies that h(w) → ±∞ for w → ±∞. Hence, in order to enter in the setting of the existence result in [19,42] (see also [46] for a survey on this topic) for system (4.3), we just need to assume that is strictly monotone increasing with lim y→0 y B(y) = 0.
At this point we can replicate the argument for (H p,λ ) with 0 < p < 1 and get the existence of at least a limit cycle for the original system. Thus we have the following result.

Example 1
To show an example of applicability of Theorem 4 we consider the system ẋ = y − λ 2 π arctan(|y| 1/2 + y 2 )(x 2 − x − 1)ẋ y = −x 3 . Figure 13 illustrates the associated dynamics for λ = 1. Fig. 13 The figure represents the dynamics for system (S1), in the phase-plane (panel left) and for component x(t) in the time variable (panel right). This example concerns the sublinear case with B(y) = 2 π arctan(|y| 1/2 + y 2 ), A(x) = (x 2 − x − 1)x and g(x) = x 3 . The simulation suggests the uniqueness of the limit cycle, but in order to rigorously prove this fact, one should first check the validity of the property (S ) introduced below. The graphs are produced using [48] As discussed in [9,20,22,46] the topic concerning the uniqueness of the limit cycle is more delicate. We show how our argument developed in Sect. 3 for system (H p,λ ) can be affectively applied to our situation as well. We recall that the problem of the uniqueness of the limit cycle for the classical Liénard equation (1.2) has been studied in [40,41], by considering the following classical property:

(a) (b)
All the possible limit cycles cross both the lines x = α and x = β, which we referred to R. Conti [10] in [41]. It is worth to note, and it is a pleasure, that Professor Roberto Conti was the thesis advisor of the first author, and was a mentor and a guide for both authors when they were moving their first steps in this fascinating mathematical field area. As remarked in [46, p.1192], this crossing condition (implicitly used in the proof of the result in [20], further corrected in [22] and [46]) is crucial also for the proof of the uniqueness of the limit cycles for (4.3) or, equivalently, (1.3), otherwise some counterexamples can be constructed. In this setting we present the following result.
Proof By property (S ) all the possible limit cycles cross the strip [α, β] × R.
Assume (by contradiction) that there are two limit cycles Γ 1 and Γ 2 (with Γ 1 in the internal region of Γ 2 ) and compute the line integrals Γ i −g(x(t))F(x(t), y(t)) dt for i = 1, 2. Splitting these integrals along the orbit paths inside and outside the strip [α, β] × R, we obtain: Outside the strip, the integrals (3.4) take the form Since the function A(x) is monotone increasing and is positive for x > β and negative for x < α, we find that the contribution to the energy along the external limit cycle is less than the contribution in the internal one. Inside the strip the integrals (3.5) and (3.6), become, respectively dx, for y < 0. Now, using the assumption (B2) we find again that the contribution to the energy along the external limit cycle is less than the contribution in the internal one. In conclusion we get that which is a contradiction.
In our case, in order to fulfill the implicit assumption (S ), we need some symmetry properties which, for example are satisfied if B(y) is even and A(x) and g(x) are odd. If these additional assumptions hold, combining Theorem 4 and Theorem 5 we get that: System (S1) has exactly one limit cycle which is asymptotically stable.

Remark 4
One may argue that assuming the above mentioned symmetry conditions one can get directly the uniqueness of the limit cycle from the auxiliary system (4.3). This because the results on the number of limit cycles for (4.3) transfer to our case via the change of variables considered in the proof of Theorem 4. This argument is correct. However, we observe that our uniqueness result can be applied to situations where the above quoted existence results for (4.3) are not applicable. This occurs, for instance, when the growth conditions for F(x) or G(x) at infinity are not satisfied. On the other hand, limit cycles may in principle exist and Theorem 5 gives the uniqueness. A result related to Theorem 4 was previously obtained in [32,Theorem 1], for the study of a system of the form (S1) in the sublinear case. This also shows the effectiveness of the method of energy integration introduced in the pioneering work of Liénard and developed by several authors since then.

The superlinear case: y/B(y) → 0 for y → ±∞
This case is more intriguing, due to the presence of new branches of the infinite isocline as it was shown in Figs. 5 and 6 and now the parameter λ > 0 will play a crucial role. As far as we know, this interesting phenomenon has never been investigated in this context before.
From now on, besides (B0) and (g0), we suppose that (A0) holds with the function A(x) having exactly three zeros α < 0 < β, so that we assume the following Observe again that this assumption is common for the classical Liénard equation (1.2). Moreover, assume that B(y) y is strictly monotone increasing and B(y) y → 0 for y → 0.
Under the above assumption, the map η(y) := B(y)/y is an increasing bijection of the real line with η(0) = 0 and therefore the vertical isocline is given by the union of the x-axis with the graph of This graph, in view of the hypotheses on A(x), exhibits a behavior similar to the one of the case study with p > 1, depicted in Fig. 5. As in Sect. 3 (cf. Fig. 5), we can define the positively invariant regions
Our next and final result is an extension of Theorem 2 which was proved for the "case study" (H 2,λ ) to general systems of the form (S1). For the proof of our previous theorem we used the fact that B(y) = |y| p , for p > 1 (even if only the case p = 2 was discussed in detail). In Theorem 2 we considered a special form of odd A(x), namely A(x) = x 3 − x and g(x) = x. In the general superlinear case studied in this section we try to avoid as possible any condition of symmetry. However, checking the proof of Step 2 in Theorem 2 one realizes that the crucial fact in the proof of the intersection of the circumference C(P 0 ) with the curve which is the boundary of the region W 2 involves a comparison between the 1/A(x) on ] − 1, 0[ with A(x) on (−∞, −1[ . Now, there two terms are transformed through the function η −1 , hence a condition relating the behavior of η −1 at −∞ with η −1 at +∞ will be required. Accordingly, we assume the following condition. Condition (RV ) is strongly related to the usual assumption considered in the theory of Karamata's Regularly Varying Functions. Classical regularly varying function usually involve a limit of the form f (σ s)/ f (s) and appear in different areas of interest such as real analysis, probability theory and asymptotic theory for ODEs (see the books [28,35] and the references therein). In this framework the following main existence result holds.
Proof The proof follows the same argument already described in detail along the proof of the "case study" Theorem 2. Therefore, we will concentrate only on the small modifications which are needed in the present case.
Step 1. Nonexistence of limit cycles for λ large. We consider the energy function H (x, y) = .

Consider now the region
The set D is negatively invariant and for each P 0 ∈ D it holds that γ − (P 0 ) tends to the origin. Using the fact that μ λ , ν λ → 0 as λ → +∞, we can determine a value of λ such that for all the larger values a region V (with V = V 1 or V = V 2 ) intersects D (a situation completely similar to that depicted in Fig. 7). Taking P 0 ∈ V ∩ D we find that the orbit γ (P 0 ) tends to the origin in negative time and is unbounded. This prevents the possibility of having limit cycles. In this manner we have proved that there existsλ > 0 such that for each λ >λ all the nontrivial trajectories are ultimately unbounded.
Step 2. Existence of limit cycles for λ small. First of all, we introduce the line , for α < x < 0 which is below the region V 1 and has v M,λ : as a minimum value. By construction, we have thaṫ which is the analogue of (3.2).

Consider the initial point
with K > √ 2G(α) a fixed constant. Clearly y 0 > 0 for λ small, as v M,λ +∞ as λ → 0 + . The slope of the trajectory departing from P 0 satisfies the relation as long as (x, y(x)) remains in the rectangle. Hence, the above a priori estimate shows that the positive trajectory departing from P 0 remains in the region R and crosses the y-axis at some point P 1 = (y 1 , 0) with y 0 < y 1 < y 0 + K . Consider now the point P 2 = (x 1 , y 1 ) of intersection of the horizontal line y = y 1 with the boundary of the region W 1 , given by the equation y = κ λ (x) for x > β and observe that the segment P 1 P 2 is crossed by the vector field from the exterior to the interior, asẋ ≥ 0 andẏ < 0 (for x > 0). Considering now the negative semi-orbit γ − (P 0 ) departing from P 0 we have that such a trajectory remains in the region outside the level line as long as x < α. We claim that γ − (P 0 ) crosses the negatively invariant region W 2 for λ > 0 sufficiently small. To this end, it is sufficient to show that the system (which is the analogue of (3.3) to our situation) has solution x in the domain where G −1 l is the inverse of G restricted on the negative real numbers. From the above system we get the equation To prove that this equation has solutions in the desired range for λ sufficiently small, we proceed as follows.
On the other hand, for any fixedx < α we have that Using condition (RV ) with s := 1/(λ|A(x)|) and σ := |A(x)|/2 A(x M ), we obtain that 2 and observe that we can have η < −1 if we takex < α sufficiently negative (using also the fact that |A(x)| → +∞ for x → ±∞). Therefore the existence of an intersection point P −1 of γ − (P 0 ) with the boundary of W 2 is guaranteed. From now on, we proceed in a complete symmetrical manner to determine points Q −1 , Q 0 , Q 1 and Q 2 starting from a point Q 0 = (β, z 0 ) above the region V 2 and, like in the proof of Theorem 2, we arrive at the construction of a Jordan curve Υ which is the boundary of a positively invariant region surrounding the origin and bounding a limit cycle. In this manner we have proved that there existsλ > 0 such that, for each 0 < λ <λ, system (S2) has a limit cycle.
We skip the discussion concerning Step 3 because the geometry of (S2) in the strip α < x < β is the same as that of (H 2,λ ) in −1 < x < 1. This concludes the proof of our theorem.

Example 2
To show an example of applicability of Theorem 6 we consider the system ẋ = y − λ(tanh(y 2 + y + 1))y 2 (x 2 − 2x − 3)ẋ y = −x. Figures 14, 15, 16 illustrate, respectively, the associated dynamics for λ = 0.0522 (existence of a limit cycle) and for λ = 0.0523 (escape of all the nontrivial orbits). Clearly, this is just a numerical simulation which shows the above phenomenon which is proved analytically. Precise estimates from above and below on Λ * are in any case obtained within the limits of the numerical approximations.
The "escape effect" exhibited by a single orbit in Fig. 15 is even more dramatically illustrated considering the whole set of trajectories in Fig. 16.
Another approach to prove the existence of limit cycles for λ ∈]0, Λ * [ is based on a classical bifurcation technique, dating back to Poincaré (according to Lefschetz [24,) and Liénard to prove the existence of limit cycles bifurcating form a periodic orbit of a center when a parameter multiplying the nonlinear terms is small. This method (a) (b) Fig. 14 The figure represents the dynamics for system (S1), in the phase-plane (panel left) and for component x(t) in the time variable (panel right). This example concerns the superlinear case with B(y) = (tanh(y 2 + y + 1))y 2 , A(x) = (x 2 − 2x − 3)x and g(x) = x. The simulation suggests that for these parameters the limit cycle is unique. In the phase-plane we also have plotted part of the infinite-isocline near the trajectories. The limit cycle is possible because "fortunately" the trajectories avoid (it is narrow escape !) the regions V 1 and V 2 . The graphs are produced using [48] (a) (b) Fig. 15 The figure represents the dynamics for system (S1), in the phase-plane (panel left) and for component x(t) in the time variable (panel right). This example concerns the superlinear case with B(y) = (tanh(y 2 + y + 1))y 2 , A(x) = (x 2 − 2x − 3)x and g(x) = x. In this case there is no limit cycle because all the nontrivial solutions enter eventually the regions the regions V 1 and V 2 . In the phase-plane we also have plotted part of the infinite-isocline near the trajectories. The graphs are produced using [48] Fig. 16 Planar dynamics for system (S1). It is evident that the trajectories are unwinding from the origin and eventually become unbounded entering the regions V 1 or V 2 . Our simulation shows that a large region around the origin (excluding (0, 0)) is attracted to V 2 . The graph is produced using [49] for the same value of the parameter λ = 0.0523 as in Fig. 15 and the same functions B(y), A(x) and g (x) was successfully applied by Duff and Levinson [11] (see also [34]) to produce multiple limit cycles bifurcating from circular orbits of the harmonic oscillator and represents a very useful technique to construct specific examples of multiplicity results for Liénard or Rayleigh equations. We refer also to [4,Remark 2.4] for another application of this method. Due to space limitations, we do not enter into these details, however we think that it is interesting to mention the fact that this averaging approach provides an independent proof for the existence of limit cycles.
Funding Open access funding provided by Università degli Studi di Udine within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.