Maximum Amplitude of Limit Cycles in Liénard Systems

We establish sufficient criteria for the existence of a limit cycle in the Liénard system ˙ x = y − εF (x), ˙ y = −x, where F (x) is odd. In their simplest form the criteria lead to the result that, for all finite nonzero ε, the amplitude of the limit cycle is less than ρ and 0 a ρ u, where F (a) = 0 and u 0 F (x)dx = 0. We take the van der Pol oscillator as a specific example and establish that for all finite, nonzero ε, the amplitude of its limit cycle is less than 2.0672, a value whose precision is limited by the capacity of our symbolic computation software package. We show how the criterion for the upper bound can be extended to establish a bound on the amplitude of a limit cycle in systems where F (x) contains both odd and even components. We also show how the criteria can be used to establish bounds for bifurcation sets.


I. INTRODUCTION
Oscillatory phenomena occur widely in nature [1] and they also arise in engineered systems [2].Some of the simplest mathematical descriptions of nonlinear sustained oscillations are based on limit cycles in second-order systems [3].An area of early research in relation to specific nonlinear systems was the identification of whether limit cycles could occur and, if so, the determination of their amplitude and frequency [4][5][6].The behavior of coupled systems of nonlinear oscillators has also attracted attention, particularly since the development of the Kuramoto model for synchronisation in oscillator ensembles [7].In some analyses of coupled oscillators it has been sufficient to consider only phase aspects of the dynamical behavior and ignore amplitude aspects [8,9], but in other systems of coupled oscillators it has been observed that the dynamical behavior also depends significantly on oscillator amplitude [10].A recent development in relation to the analysis of oscillatory dynamics which are stable under continuous perturbations includes a generalization of the limit-cycle concept, the "chronotaxic limit cycle" [11].
In this paper we return to the issue of limit-cycle amplitude.However, rather than seek an approximation to the amplitude of oscillations for a specific set of parameters, we seek to establish general bounds on the amplitude of oscillations.Such general bounds can serve as a useful cross-check on numerical results.They may be useful in some engineering contexts, because the maximum amplitude of oscillations can have a significant impact on critical aspects of system behavior.For example, it appears that ferrous alloys have distinct fatigue limits (i.e., for specific alloys there appears to be a maximum amplitude of cyclic stress that can be applied to the material without causing fatigue failure) [12].Other examples where the maximum amplitude of motion is important are provided by earthquakes and animal limb movement in locomotion [13].Furthermore, the determination of bounds on oscillation amplitudes is related to a fundamental mathematical problem [14]: What is the maximum number of limit cycles for coupled pairs of first-order differential equations?
There has been much interest [15] in the control of nonlinear dynamical systems where bifurcation can occur, i.e., the occurrence of a significant change in behavior as a system parameter is smoothly changed.A successful controller for such systems would typically be required to delay the onset of a bifurcation or moderate its effects in order to achieve some target dynamical performance.Our method to establish bounds on limit-cycle amplitudes generates sufficient conditions for the existence of limit cycles.A by-product is that these conditions provide a nonperturbative way of establishing bounds on bifurcation sets in relation to limit cycles, and this may assist in the development of systems for real-time bifurcation control.
A popular choice of model for simulating the dynamics of nonlinear physical phenomena is the van der Pol oscillator [16].Its behavior is governed by the equation where ε is a positive scaling parameter.This equation was introduced in the early 1920s by van der Pol to model discharges in triode circuits.It has subsequently been applied to a vast range of phenomena, including vortex shedding [17], human circadian rhythms [18], and tunnel diode oscillators [2].
The van der Pol oscillator can exhibit a wide range of behavior.For example, with ε 1 the time variation of x for a single oscillator is nearly sinusoidal [5], but with ε 1 the time variation of x resembles a square wave [19].For the latter case van der Pol introduced the term relaxation oscillation [20].The popular appeal of the van der Pol oscillator as a modeling tool may relate to the wide range of behavior it can exhibit, the simplicity of its formulation, and the large canon of literature reporting its successful use.
The van der Pol oscillator is a simple example from a class of nonlinear oscillators known as Liénard oscillators.In the late 1920s Liénard investigated the equation in relation to electrical circuits [4].Here f (x) is even, ε is a positive scaling parameter, and ω is a constant.In recognition of this, Eq. ( 2) bears his name.Because Liénard systems have received much attention from mathematicians, various authors have generated transformations to convert other two-dimensional systems into Liénard form.See, for example, Ref. [21].However, we note that often the definition of the Liénard form is extended to where g(x) is not limited to ω 2 x.In the case where f (x) = x 2 − 1 and ω = 1, Eq. ( 2) becomes the equation for the van der Pol oscillator, Eq. (1).Equation ( 2) may also be expressed as two coupled firstorder equations, Other pairs of coupled first-order equations have been developed to generate the same system dynamics as Eq. ( 4) (see Ref. [22] and the references therein).In this paper we use the version due to Liénard [4], where F (x) = x 0 f (s)ds and, for consistency with Refs.[23] and [19], we set ω = 1.εF (x) is known as the damping term.In the particular case of the van der Pol oscillator F (x) = x 3  3 − x in Eq. ( 5).The main focus of this paper is the determination of bounds for the amplitude (bounds for the maximum x coordinate) of the limit cycle of Eq. ( 5), subject to reasonable physical constraints on εF (x).
In general it is not possible to find an analytic solution to an arbitrary differential equation [19].Instead, qualitative techniques have been developed to enable important characteristics of solutions to be deduced in the absence of a formal solution.
A qualitative device we use to establish our results is the phase plane [19].For every pair of values x,y the planar vector field E = ( dx dt , dy dt ) is tangent to the solution of dx dt = P (x,y), The plane corresponding to the planar vector field E = ( dx dt , dy dt ) is known as the "phase plane" and, in the case that the coupled equations are given by Eqs.(5), the phase plane is known as the Liénard plane.To permit the use of some standard vector results we utilize the Liénard plane but acknowledge the existence of a z axis perpendicular to the Liénard plane.This z axis together with the x and y axes of the Liénard plane forms a right-handed coordinate system.At all times we impose the condition dz dt = 0.For Eq. ( 6) equilibrium points in the phase plane occur where P (x,y) and Q(x,y) are simultaneously zero.A trajectory is the locus of temporally successive points in the phase plane and trajectories may only cross at equilibrium points.A limit cycle is an isolated, closed trajectory.Trajectories adjacent to the limit cycle either spiral in towards the limit cycle or spiral away from the limit cycle.Every limit cycle must enclose at least one equilibrium point [19].The sustained oscillations of Liénard oscillators correspond to limit cycles in the Liénard plane.All limit cycles in Liénard systems must be concentric, because there is only one equilibrium point [i.e., (0,0)] and trajectories cannot cross.Furthermore, TABLE I. Amplitudes of limit cycles in the van der Pol equation (1), calculated by numerical integration for different values of ε, according to Zonneveld [27] and Odani [26]. in representations such as Eq. ( 5) they possess the symmetry that the phase plane is unaffected by the substitutions (x,y) → (−x,−y).The Poincaré-Bendixson theorem (see, for example, Ref. [19]) states that if a trajectory remains in a closed bounded region D of the phase plane, where D consists of ordinary points of a system (P ,Q) = ( dx dt , dy dt ), then either the trajectory is closed (i.e., a limit cycle), approaches a closed trajectory, or terminates at an equilibrium point.If an annular region D, which contains no equilibrium points, can be found such that all trajectories which cross the boundaries of the D region do so inwards, then the existence of at least one limit cycle in the region D can be inferred and thus bounds on its amplitude can be established.Practically, it can be very difficult to determine regions D which satisfy these requirements.
It is known [24] that Eq. (1) has a single limit cycle for all finite, nonzero ε and Odani [25] has shown that this limit cycle is not algebraic.Odani [26] has commented that it is a very difficult problem to determine the maximal amplitude of the limit cycle of Eq. ( 1) and, by implication, it is a difficult problem to determine the maximal amplitude of the limit cycles of Liénard equations generally.Over 50 years ago, in the context of numerical analysis, Lefschetz [6] pointed out the desirability of establishing the amplitude of the van der Pol limit cycle.The early numerical issues have been resolved but we have not found reports of any simple, general method to identify upper bounds on the amplitudes of limit cycles in Liénard systems.
A wide variety of methods has been used to approximate the amplitude of the van der Pol limit cycle.They include numerical integration [27] (see Table I), harmonic balance [28,29], variational analysis [30], the Melnikov method [31], empirical analyis [23], piecewise approximation [32], and bifurcation analysis [33].Most of the methods (Melnikov, piece-wise approximation, bifurcation, and harmonic balance) are valid only for specific ranges of the strength parameter ε.The enhanced harmonic balance method proposed by Buonomo [29] generates series in ε which can, in principle, be developed up to any order.Lopez et al. [34] have applied the homotopy analysis method to obtain a recursive set of formulae that approximate the amplitude and form of the van der Pol limit cycle for the whole range of the strength parameter ε.
Some of the methods used to approximate the amplitude of the limit cycle of the van der Pol oscillator Eq. ( 1) have also been employed with success to approximate the amplitude of limit cycles in other Liénard systems where F (x) is an odd polynomial of degree greater than 3.Among these methods are the variational method of Depassier and Mura [30] and the empirical, nonperturbative method of Giacomini and Neukirch [23].
In this paper we establish sufficient criteria for the existence of a limit cycle in the system given by Eq. ( 5), subject to reasonable physical constraints.In their simplest form the criteria lead to an upper bound on the amplitude of the limit cycle, valid for all finite, nonzero, ε.We illustrate how the bounds developed with F (x) odd can be applied to systems where the damping term contains both odd and even components and how the criteria can be used to provide lower bounds for bifurcation sets in Liénard systems.

A. Overview
We will develop sufficient criteria for the existence of a limit cycle in a system given by Eq. ( 5) with the addition of simplifying criteria for ε and F (x) [Eq.(7), below].Our method shares some features with the geometric path-integral approach presented by Jordan and Smith [19] but involves the additional process of aggregating the path integrals.In order to aggregate the integrals we employ a demichord technique which was introduced by Liénard [4].The use of the demichord technique permits us to convert path integrals to simple integrals and thereby change the order of integration of a double integral.This in turn allows us to isolate, to an integral with respect to x, the sign change related to the occurrence of a limit cycle.[The main result that Liénard derived using the demichord technique was that in the Liénard plane for every trajectory of Eq. ( 5) the y amplitude exceeds the x amplitude].
Based on the above we will present two general equations [Eqs.(26) and (27)] in relation to the amplitude of the corresponding limit cycle.A root of one of the equations gives an upper bound for the limit-cycle amplitude and a root of the other gives a lower bound for the limit-cycle amplitude.We then utilize techniques from the theory of rotated vector fields [3] to develop upper bounds for limit-cycle amplitudes in systems where F (x) has both odd and even components.
The analysis is undertaken with respect to the Liénard plane, but for some vector results we make reference to a right-handed coordinate system with unit vectors ( î , ĵ, k) along the x,y,z directions respectively (i.e., î and ĵ in the Liénard plane and k perpendicular to the Liénard plane).In the generic upperbound analysis we consider Liénard equations of the form with ε > 0 and finite, E = ( dx dt , dy dt ), F (x) = x 0 f (s)ds, and subject to the following conditions on f (x) and F (x).
Conditions II A: (i) f (x) is continuous and even, (ii) u 0 F (x)dx = 0 for some u > 0, (iii) for |x| < u,F (x) = 0 only at x = 0 and x = ±a, for some a such that 0 < a < u, and (iv) for 0 < x < a,F (x) < 0 and for a < x < u,F (x) > 0.
The vector field E = ( dx dt , dy dt ) for this system is unaltered by the substitution (x,y) → (−x,−y) so if (x,y) is a point on the limit cycle in the Liénard plane, then so also is (−x,−y).In particular, if (0,y) is a point on the limit cycle, then so is (0,−y).There is only one equilibrium point, i.e., (0,0).The extremum of each trajectory with respect to variations in y occurs on the line y = εF (x), which divides the region −u < x < u of the phase plane into an upper portion and a lower portion.The extrema of each trajectory with respect to variations in x all occur on the y axis and trajectories rotate clockwise about the origin.The direction of rotation can be confirmed by considering the vector cross product of a vector along the positive x axis with E (i.e., î × E).A resultant vector parallel to the vector k in a right-handed coordinate system corresponds to a counterclockwise rotation from the x axis to the vector field E. However, for the system given by ( 7), ( î × E) = −x k so the trajectories rotate clockwise.

B. Line integrals in the Liénard plane
Consider a polynomial H (x,y) which satisfies the following conditions.
Conditions II B: Clearly H (x,y) = x 2 + y 2 satisfies these conditions.As shown in Sec.II C, there is a systematic method to generate other polynomials which possess these properties.The polynomials H (x,y) can be used to generate the vectors where k is a unit vector along the positive z axis.Using the vector identity ∇ Eq. ( 8) it can be readily shown that for all Consider integrals along paths in the Liénard plane as illustrated in Fig. 1.Each path consists of a trajectory γ of the system (7) in the right half-plane from a point A (0,y A ) on the positive y axis through the extremum of γ at B (x,εF (x)) back to a point on the negative y axis at C (0,y C ).The path is then closed back along the y axis to A. Using the divergence theorem and Eq. ( 9) the integral along the closed path ABCA of the flux due to the vector field H(x,y) is 0. We designate the flux through the portion of the path along the trajectory γ as I ABC and the flux through the portion of the path along the y axis as I CA .Thus The outward normal to the exterior of the path ABCA along CA is − î so the aggregate flux of the vector field H through path segment CA along the y axis is given by Using Eq. ( 10) and the property that − î • ( k × ∇H ) x=0 is a monotonic odd function of y only, we see that the net flux, I CA , through the portion of the path CA along the y axis is dependent only on the y coordinates of C and A. If |y A | = |y C |, then the net flux through the path CA along the y axis will be 0; the net flux, I ABC , through the path ABC will also be 0 and the trajectory in the right-hand plane will be part of a limit cycle.If |y A | > |y C |, then I CA < 0,I ABC > 0 and on the path ABC the trajectory is spiraling inwards so it is not part of a limit cycle.If |y A | < |y C |, then I CA > 0,I ABC < 0 and on the path ABC the trajectory is spiraling outwards and so it is not part of a limit cycle.
Since the trajectory is rotating clockwise about the z axis the outward normal along the portion of the path ABC in the right half-plane is given by ( k × E).The aggregate flux, I ABC , through the path ABC is given by Consider the line through x = X,(0 X < x B ) and two trajectories γ i and γ k as illustrated in Fig. 2. For y > εF (x) all trajectories whose extremum (with respect to variations in y) occurs for x > X cross the line x = X in the direction of increasing x.Since trajectories cannot cross each other [except at the equilibrium point (0,0)], we deduce that if OA i < OA k , then the y coordinates of the points where γ i and γ k cross x = X satisfy M i < M k .Similarly, for y < εF (x), all trajectories whose extremum (with respect to variations in y) occurs for x > X cross the line x = X in the direction of decreasing x and the y coordinates of the points where γ i and and for all x = X,(0 X < x B ) the length of the chord of the trajectory γ i at X is less than the length of the chord of the trajectory γ k at X.
We now use a result due to Liénard [4] to convert the line integral I ABC to a simple integral.For a general trajectory γ let M(x) be the y coordinate of the portion of the trajectory in the right half-plane for y > εF (x) and let L(x) be the y coordinate of the portion of the trajectory for y < εF (x) and let Thus η(x) is the length of the demichord at x of the portion of trajectory γ in the right half-plane.It can readily be shown (see Ref. [4]) for the trajectory γ in the right half-plane with extremum (x B ,εF (x B )) that (i) η(x) 0 for 0 < x x B .η(x) = 0 when x = x B and y = εF (x B ), i.e., η(x) = 0 at the extremum of the trajectory with respect to variations in y, (ii This property of dη dx is derived in Ref. [4].We include it here for completeness but do not use it in our analysis, and (iii) It can also easily be shown (again see Ref. [4]) that for an arbitrary function (x) which is a function of x only that the path integral along a path ABC whose start point A and end point C have the same x coordinate, By conditions II B, 1 x dH dt is a function of x only, so applying Eq. ( 14) to Eq. ( 11) with (x) = 1 x dH dt , where η 0 is the length of the demichord at x = 0.For every trajectory in the right-hand plane there is a unique extremum with respect to variations in y [i.e., (x B ,εF (x B ))], and for every positive x B there is a trajectory and associated path integral We now consider the aggregate J (x β ) of all integrals I ABC (i.e., for all trajectories) in the right half-plane which have an extremum, (x B ,εF (x B )) with x B x β in the right half-plane.If J (x β ) changes sign as x β gets progressively larger, then this indicates that, at some value of x B < x β , the value of I ABC has changed sign.This in turn implies that the trajectory through (x β ,εF (x β )) lies outside the limit cycle, Integration is a linear operation and, with ε finite, we assume that for all physically reasonable trajectories the integrals 0 η 0 1 x dH dt dη are finite.Under these conditions inversion of the order of integration is permitted.We also invert the limits of integration for η where the integrand, 1 x dH dt of the inner integral, is now a function of x only.Since η 0 for all the integrals, the sign of J (x β ) is determined by the sign of − We assemble these results for line integrals along trajectories in the Liénard right half-plane to establish sufficient conditions for the existence a limit cycle.
From this we deduce that all trajectories wholly contained in the region −a * < x < a * are spiralling outwards.This implies there can be no limit cycle wholly contained in the region −a * < x < a * .If, further, x β 0 1 x dH dt dx = 0 for some x β ,x β > a * , then J (x β ) = 0 and so for the trajectory through (x β ,ε(F (x β )),|y A | > |y C | and this trajectory is spiralling inwards.This implies there is a limit cycle through some point (x B ,εF (x B )), where a * < x B < x β .In particular, the maximum amplitude of the limit cycle satisfies a * < x B < x β .

Thus sequential sign changes in 1
x dH dt and x β 0 1 s dH dt ds provide a sufficient condition for the existence of a limit cycle.

C. Giacomini and Neukirch polynomials
Giacomini and Neukirch [23] introduced a sequence of polynomials h N (x,y) with the following form: where g n,N (x),0 n N − 1 are polynomials in x only and N is even.By enforcing the constraint that dh N (x,y) dt is a function of x only and observing the behavior of dh N (x,y) dt as x is varied, Giacomini and Neukirch were able to make empirical claims about the number and approximate location of limit cycles in Liénard systems [23].
Our subscript notation for polynomials h N (x,y) differs from that introduced in Ref. [23] because, for clarity in recurrence relations, we reserve N to designate the highest degree of y in the polynomial.The time derivative of h N is given by Substituting for dx dt and dy dt gives where and ġn,N (x) with 0 n N − 1 represents the derivative of polynomial g n,N (x) with respect to x. Combining Eqs. ( 20)- (22) gives and, if The condition that for 0 < n N − 1 the coefficient of y n should be zero in Eq. ( 23) gives the following recurrence relation: with initial conditions Following Ref. [31] we use the designation "GN polynomials" for polynomials h N (x,y) = Z which satisfy Eqs. ( 18), (24), and (25).We note that h 2 (x,y) (i.e., x 2 + y 2 ) is a polynomial irrespective of the form of F (x).However, in general, for N > 2,h N (x,y) is not a polynomial unless F (x) is a polynomial.The ordinary differential equations in the recurrence relation given by Eq. ( 24) are easily integrated with respect to x when F (x) is a polynomial.Symmetry conditions require the constants of integration for ġn,N to be zero for n odd.Following Ref. [23] we set all constants of integration to be 0 but note that symmetry does not require the constants of integration for ġn,N to be 0 for n even.We consider the choice of constants of integration further in Sec.IV.
Examples of GN polynomials are given in Table II and examples of their time derivatives are given in Table III.By construction each h N (x,y) satisfies conditions II B. Thus dt is finite for all finite x, and − î • ( k × ∇h N ) x=0 is a monotonic odd function of y.Accordingly, the results derived in Sec.II B for the polynomial H (x,y) also apply to the GN polynomials h N (x,y).

D. Bounds on limit-cycle amplitude using h 2 and h
For paths using h 2 = x 2 y 2 we consider the integral I ABC [see Eq. ( 11)] for the trajectory through point (x B ,εF (x B )) in the right half-plane, For all trajectories through points (x B ,εF (x B )) in the righthand plane such that |x B | < a,F (x B ) < 0 (see conditions II A), so |y A | < |y C | and the trajectories are spiraling outwards.This gives the result of Ref. [19], applicable for all finite nonzero ε, that no limit cycle can be wholly contained in the region −a < x < a.
We then consider the sum J (x β ) [see Eq. ( 17)].J (x β ) changes sign when From conditions II A J (x β ) = 0 when x β = u so, for all finite, nonzero ε, there is a trajectory through some point (x B ,εF (x B )), with a < x B < u, for which I ABC = 0. Thus for all finite, nonzero ε the system (7) subject to conditions II A has a limit cycle, and the amplitude of the limit cycle lies between a and u, where For paths using h 4 (x,y) (see Table II) we again consider the integral I ABC [see Eq. ( 11)] for the trajectory through point (x B ,εF (x B )) in the right half-plane.From these we deduce that the system (7) subject to conditions II A has a limit cycle, and the amplitude of the limit cycle lies between a 4 and u 4 , where 4x 2 εF (a 4 ) + 4ε 3 F 3 (a 4 ) + 4ε For ε 1 the term x 0 4ε 3 F 3 (s)ds will be dominant in Eq. ( 29) and for ε 1 the term ε[ x 0 4s 2 1 F (s 1 ) + 4 s 1 0 sF (s)ds)ds 1 ] will be dominant.The result for the lower bounds [i.e., Eqs. ( 26) and (28)] is effectively the same as that given in Ref. [23], except our result explicitly retains ε.The results given by Eqs. ( 27) and ( 29) are new.Similar results can be given for h N (x,y) for N > 4 but the number of terms in the expressions for the corresponding a N and u N becomes larger as N increases.

E. Rotated vectors
We now utilize elements from the theory of rotated vector fields [3] to extend our method.The extension allows us to establish an upper bound for limit-cycle amplitude in systems where the damping term has both odd and even components.If the minimum x value on the limit cycle is x − (in the left half-plane) and the maximum x value on the limit cycle is x + (in the right half-plane), then we seek a bound ρ on the amplitude such that ρ > x + and ρ > |x − |.
Let E = ( dx dt , dy dt ) = (P (x,y; α),Q(x,y; α)) with F (x) satisfies the conditions II A and P * (x) contains both odd and even components.We impose the following further conditions.
(iii) P * (0) = 0; this ensures that for every α the system (30) has only one equilibrium point in the Liénard plane for all α.
(iv) xP * (x) < 0 for 0 < |x| < u, where about the origin wholly contained within 0 |x| < u.Thus in , E cuts the circles of h 2 in an interior to exterior direction.Here we have used the notation of Sec.II D (i.e., h 2 (x,y) = x 2 + y 2 ).System (30) is just system (7) with the addition of an extra term αP * (x) to dx dt .Hence with α = 0 system (30) reduces to the system considered in Sec.II D and for all finite, nonzero is the region about the origin for which condition II E (v) applies.D is the annular region between the boundary of and the limit cycle 0 of Eq. ( 30) with α = 0.The arrows with straight tails illustrate the vector field ( E) α at the boundaries of D. The arrowheads without tails illustrate the vector field ( E) 0 on the limit cycle 0 .ε, there is a limit cycle with amplitude <u.From Sec.II B we know that the corresponding limit-cycle trajectory, which we designate 0 , is rotating clockwise.We now show that for α = 0 + α ( α > 0) the system given by Eq. ( 30) has a limit cycle α bounded by 0 .Consider the vector cross product The right-hand side of Eq. ( 31) corresponds to a vector perpendicular to the Liénard plane with magnitude determined by the sine of the angle, measured counterclockwise, between the vectors ( E) 0 and ( E) α , From condition II E (iv), in the region 0 |x| < u,xP * (x) 0, so the trajectories γ α of the system ( E) α are either tangent to 0 or are rotated clockwise with respect to 0 .Thus any trajectory γ α after entering the region about the origin bounded by 0 cannot run out again.This is a local result, similar to the global result given in Ref. [3] in relation to a complete family of rotated vector fields.
Although we cannot express 0 algebraically we select it as the outer boundary for an annular region D to which we will apply the Poincaré-Bendixson theorem.For the inner boundary of the annular region we select any level set of h 2 (x,y) wholly contained in .Regions D and are illustrated in Fig. 3 where we use arrows without tails to illustrate the vector field ( E) 0 and vectors with straight tails to illustrate the vector field ( E) α .By condition II E (v), all trajectories of E cross level sets of h 2 (x,y) in an interior to exterior direction (i.e., from the interior to the exterior of in Fig. 3).This in turn means that they cross the inner boundary of region D in an exterior-to-interior direction.Thus no trajectory which enters into the region D can ever leave the region.By the Poincaré-Bendixson theorem a limit cycle ( α ) must exist within the annular region D or on its boundary.In particular, α is bounded by 0 and thus the amplitude of α is less than u.
Using induction and the above result that a limit cycle α exists, these arguments can be extended to show that a limit cycle wholly bounded by 0 exists for any value of α > 0 provided that condition II E (v) remains valid.

F. Multiple limit cycles
In Sec.II A we do not impose conditions on the behavior of F (x) for |x| > u.We now extend the conditions of II A on system (7) for |x| > u by conditions II F.
Conditions II F: (v) (vii) For u < x < a 2 ,F (x) > 0 and for With these extra conditions the arguments of Sec.II B can be extended to show the existence of a second limit cycle, concentric with the first.The amplitude of the second limit cycle lies between a 2 and u 2 .This process can be used iteratively to provide sufficient conditions for further limit cycles.

A. Van der Pol and damping term ε F(x) = ε(x 3 /3 − x)
In this section we take the van der Pol oscillator Eq. (1) as an example to illustrate our results.In Sec.II B we showed that for a damping term satisfying conditions II A system (7) has a limit cycle whose amplitude lies between the value a such that | dH dt = 0| x=a and the value of u such that dt ds (here s is a dummy variable to represent x under integration) for even N up to 18.For h N ,N = 2,4,10,12,14,16, and 18, Table IV illustrates the x values at the zeros of the resultant terms.For ease of presentation, the terms are aggregated in their powers of ε.
With respect to the van der Pol oscillator we observe empirically that for every combination of even N (N 18) and odd power of ε in − 1 x dh N dt , the corresponding polynomial has a form that contains the key features of W A(x) and W B(x) as illustrated in Fig. 4. Empirically, we observe that for even N the term for each odd power of ε in −  dt ds = 0. Terms for h 6 and h 8 are omitted for ease of presentation.As N increases, for each degree of ε the x values for the zero of the corresponding polynomial associated with the upper bound progressively decrease (i.e., improve).Also as N increases, the largest x value for the zero of a polynomial for the upper bound progressively decreases, i.e., the zero which determines the maximum possible limit-cycle amplitude progressively decreases.(iii) W A(x) and W B(x) are both negative as x increases to their respective root and are positive thereafter. Let , where σ 1 and σ 2 are finite constants such that σ 1 σ 2 > 0. Applying the intermediate value theorem between x = A and x = B we deduce that W (x) has a zero for x > 0 at x = w, where A < w < B. The values A and B as illustrated in Fig. 4 correspond to the x values associated with a zero of a polynomial with specific N and specific degree of ε in Table IV.For example, consider the terms in Table III for N = 4.We take the term with degree 1 in ε, substitute F (x) = (x 3 /3 − x), and integrate, 4ε We take the corresponding term with degree 3 in ε and integrate, We observe empirically that, as the highest degree of y in the corresponding GN polynomial (i.e., N ) increases, for each degree of ε the x value for the zero of the corresponding polynomial associated with the upper bound progressively decreases, (i.e., improves).Thus, for example, for N = 4 to N = 18 we observe that the x value of the zero for the ε 3 polynomial progressively decreases from 2.1950 to 2.0617.We also observe that, as N increases, the largest value of x corresponding to a zero in the polynomials for the upper bound progressively decreases (i.e., the x value of the zero which determines the maximum possible limit-cycle amplitude progressively decreases).For the highest degree GN polynomial (N = 18) shown in Table IV, the largest value of x for the zero of a polynomial is 2.0672 and this occurs for ε 7 .(For comparison see the numerical approximations to the amplitude given in Table I.) From this we deduce that for all finite, nonzero ε, the amplitude of the limit cycle of the van der Pol oscillator is less than 2.0672.

B. Damping term with odd and even components
In this section we establish an upper bound for the amplitude of the system E = ( dx dt , dy dt ) with where satisfies the conditions II A and in particular u 0 F (x)dx = 0 for u = √ 6. P * (x) contains a factor x so P * (0) = 0. Also xP * (x) 0 for −u < x < u + u so xP * (x) satisfies conditions II E (iii) and II E (iv).
To show the presence of a limit cycle for this system and to determine an upper bound for its amplitude, we seek an annular region D such that all the trajectories of E which cut the inner boundary of D do so in an exterior to interior direction, and all the trajectories of E which cut the outer boundary of D do so from exterior to interior.For the case α = 0 we have shown in Sec.II D that there is a limit cycle, which we now designate as 0 , with amplitude less than u.We then select 0 for the outer boundary of our region D.
We now seek to establish the inner boundary of D. Using the notation for h 2 from Sec. II D and substituting u 2 = 6, For x small we neglect x 3 and higher-order terms in x and find E • ∇h 2 > 0 for ε > α(u 2 + u u).Accordingly, it is possible to construct a circular (i.e., h 2 ) inner boundary for D, using a sufficiently small circle that does not intersect 0 , such that all trajectories of system (36) cut the inner boundary of D from exterior to interior.Applying the Poincaré-Bendixson theorem to region D we find that, for ε > α(u 2 + u u), the system (36) has a limit cycle and its amplitude is less than √ 6.Our bound for the amplitude of 0 can be improved from √ 6 to 2.0617, using the result of Sec.III A. Likewise, our bound for the amplitude of the limit cycle for system (36) can be improved to 2.0617.

C. Bifurcation and damping term
In Ref. [35] Giacomini and Neukirch considered the Liénard system with ε(x 5 − μx 3 + x) as a damping function.
They determined that for ε = 0.1 and μ = 41 9 the system has two limit cycles but for ε = 8 there are no limit cycles.Giacomini and Neukirch reported the example to illustrate an erroneous prediction using Melnikov theory that there are two limit cycles for this system independent of ε.
Our condition [i.e., Eq. ( 27)] for an upper bound for the limit-cycle amplitude for all finite, nonzero ε requires u = 0 such that u 0 F (x)dx = 0.This corresponds to the condition u 0 x 5 − μx 3 + xdx = 0, which requires u 6  6 − μ u 4 4 + u 2 2 = 0. Real roots for u exist only for μ 48 9 (i.e., μ > 2.309) so for μ = 41 9 (i.e., μ = 2.134) we cannot establish a real value of u such that Eq. ( 27) is satisfied.Accordingly, and correctly, using h 2 (x,y) we are unable to establish sufficient conditions for the existence of a limit cycle for all nonzero ε.
We now use h 4 (x,y) and introduce a double subscript notation for the parameters a and u.The first subscript designates the degree N of the corresponding GN polynomial, and the second subscript indicates the zero being considered (e.g., the smallest positive value of a satisfying 1 x dh 4 dt = 0 is designated a 4,1 ).We find that for x > 0 each term in Eq. ( 28) has two zeros and for all ε there are values x = a 4,1 and x = a 4,2 that satisfy 1 x dh 4 dt = 0.In Eq. ( 29) with μ > 100 21 (i.e., μ > 2.182) there are two values u 4,i=1 and u 4,i=2 such that the term u 4,i 0 [4x 2 εF (x) + 4ε x 0 sF (s)ds]dx = 0, but there is no u 4 such that u 4 0 ε 3 F 3 (s)ds = 0 with μ < 2.2654.If ε 1 the former terms will dominate u 4 0 ε 3 F 3 (s)ds = 0 and for μ > 2.182 we can establish two distinct values of a 4,i and two distinct values of u 4,i such that Eqs. ( 28) and ( 29) are both satisfied.For ε 1 the term 4 x 0 ε 3 F 3 (s)ds = 0 will dominate and we cannot yet establish sufficient conditions for limit cycles to be present for μ < 2.2654.The largest intermediate value of ε such that there are two distinct values of u 4,i satisfying Eq. ( 29) provides a lower bound on the value of ε for which the system has two limit cycles for 2.182 < μ < 2.2654.Using h N (x,y) with N > 4 the ranges of ε or μ for which two limit cycles exist can be further refined.In this way sets of bounding parameters for the (saddle node) bifurcation of two limit cycles to no limit cycles can be found by considering the range of parameters (e.g., ε or μ) for which two positive values u N,i=1 ,u N,i=2 , exist such that u N,i 0 1 s dh N dt ds = 0.

A. Constants of integration for GN polynomials
In the recurrence relation Eq. (24) we selected the constants of integration for the integral of dg N,n  dx ,n even, to be 0 although this was not required on symmetry grounds.Let H 0 N (x,y) designate the GN polynomial of degree N in y developed by taking constants of integration to be 0. Let H * N (x,y) designate the GN polynomial of degree N in y generated with a nonzero constant of integration k for the integral of dg K,N dx .Then, applying the recurrence relation Eq. ( 24), we find H * N (x,y) = H 0 N (x,y) + kH 0 K (x,y).This result for H N (x,y) is just a linear combination of GN polynomials H 0 N (x,y) and H 0 K (x,y) derived with the rule that constants of integration should be 0.
In Sec.II C we dismissed the possibility of nonzero constants of integration for terms dg N,n  dx ,n odd, as being incompatible with the symmetry of the Liénard equation.However, the question arises, if the damping term in a system is not odd and thus the limit cycle of the system (if it exists) is altered by the change (x,y) → (−x,−y), could we admit nonzero constants of integration for n odd?For example, the choice of a constant of integration k,k = 0 for the integral of dg 2,1 dx generates polynomials H 2 (x,y) = y 2 + ky + x 2 = Z.These polynomials consist of circles with center (0,−k) in the Liénard plane.However, we note that there would be some choices of Z for which the corresponding circle would not enclose the equilibrium point (0,0).If the closed curve does not enclose the equilibrium point, then our arguments of Sec.II B about trajectories spiralling inwards or spiralling outwards do not apply.As a consequence, we would be unable to infer the existence of limit cycles from our criteria.Accordingly, our method does not admit nonzero constants of integration for n odd.In any case, in Sec.II E, we present a simple alternative approach to the problem of establishing bounds on the amplitude of the limit cycle for cases of damping terms with both odd and even components.

B. Retention of ω in the original Liénard equation
If the factor ω 2 from Eq. ( 2) is retained, then the system given by (5) becomes The recurrence relation Eq. ( 24) can be readily generalized to allow for this.The simplest corresponding GN polynomial is the ellipse It is easily verified that h 2 (x,y) of Eq. ( 39) conditions II B and There can be no limit cycle wholly contained in the region where the integrand of I ABC is of one sign (see Sec. II B).A change of sign for the integrand of I ABC occurs when 2εxF (x) changes sign [i.e., at x = a such that F (a) = 0].A change of sign for J (x β ) [see Eq. ( 17)] occurs when ω 2 dx changes sign [i.e., at x = u such that u 0 F (x)dx = 0].Thus there is a limit cycle whose amplitude ρ satisfies a < ρ < u, and the simplest bounds on the amplitude which we established for the corresponding limit cycle with ω = 1 remain valid for ω = 1.

C. Convergence failures
The homotopy analysis method is an analytic approximation method and over the past 20 years there have been many reports of its use in relation to nonlinear differential equations.For example, Lopez et al. [34] have used the technique to approximate the limit cycle of the van der Pol equation.Despite its widespread use, circumstances have been reported where application of the homotopy analysis method has failed to converge to correct values [36].Where homotopy analysis is applied to Liénard systems simple bounds for the amplitude of limit cycles, such as those presented in this paper, can provide a plausibility check on numerical results and reduce the likelihood of convergence problems propagating unnoticed.

D. Application to a Burridge-Knopoff model for elastic excitable media
Our result ( 27) can be used to place a bound on the amplitude of propagating wavefronts in a Burridge-Knopoff model of an excitable medium with elastic coupling [37], irrespective of the propagation speed of the front.These simple models were developed to describe the interactions between tectonic plates [38] in relation to earthquakes.Cartwright et al. [37] replaced the "stick slip" friction behavior typically used in such models with a creep-slip friction, which they designated as "asymptotically velocity weakening friction," leading to the equation Here χ (x,t) is the local longitudinal deformation of the surface of the upper plate in the reference frame of the lower plate, γ is the magnitude of the friction, c is the speed of sound, and ν represents the slip rate.They applied Eq. (42) to laboratory friction experiments and postulated that it might be applicable to electronic transmission lines and active optical waveguides.With ψ = ∂χ ∂t , equation (42) can be reduced to two firstorder partial differential equations, For some slip rates the system generates global oscillations (i.e., one surface executes temporally periodic motion with respect to the other surface).At slip rates close to, but above, the slip threshold, global oscillations become unstable.Small perturbations from the global oscillations grow to become synchronized pacemakers that emit fronts in both directions.Cartwright et al. investigated the equations qualitatively and considered solutions of the type ψ(x,t) = f (z), where z = x v + t, and v is the velocity of the front.As a tool to understand the behavior, they manipulated Eq. ( 43) so φ(ψ) = ψ 3 3 − ψ.
As a result, the model becomes a set of elastically coupled van der Pol oscillators.By assuming solutions of the form ψ(x,t) = f (z), where z = x v + t further manipulation leads to This is similar to the van der Pol equation ( 1), but with the nonlinearity scaled by μ = γ (1 − c 2 v 2 ) −1 and the x axis in the Liénard plane displaced by ν.The propagating fronts are then periodic solutions of the van der Pol equation but the value of the parameter μ is not known until the value of the front velocity v is chosen.Our results in Eqs. ( 26) and ( 27) provide bounds for the amplitude of such wave fronts irrespective of the chosen value of v. Our results could also be applied for different forms φ(ψ).

E. Application to bifurcation control
In Sec.III C we discussed how the behavior of zeros for the polynomials x 0 1 s dh N dt ds, as parameters are varied, can provide a nonperturbative method to determine lower bounds in relation to bifurcation sets for limit cycles.In applications involving real-time control, methods which involve simple polynomial root finding have potential speed benefits over both perturbative techniques and also over techniques involving integration forward in time.

F. Comparison with results of Giacomini and Neukirch
We acknowledge similarities between the technique described by Giacomini and Neukirch in a series of papers in the late 1990s [23,35,39,40], and our method.Both approaches are nonperturbative, both utilize GN polynomials, and both bring insights into the location of limit cycles and parameters associated with bifurcations.We highlight some significant similarities and differences between the two approaches: (i) The lower bound a N [i.e., the zero of 1 x dh N (x,y) dt ] developed by our method coincides with the approximation to the limit cycle developed by Giacomini and Neukirch.They acknowledge that their use of this value as an approximation to the amplitude of the limit cycle is based on empirical observation (i.e., they do not analytically demonstrate the value coincides with the existence of a limit cycle).By contrast our method delivers a pair of bounds (a N ,u N ) for the limit-cycle amplitude as part of sufficient conditions for the existence of a limit cycle.In the absence of an upper bound u N we are unable to definitely associate the existence of a zero, a N , of 1 x dh N (x,y) dt with the existence of a limit cycle.
(ii) Giacomini and Neukirch acknowledge that they are unable to explain analytically the empirical improvement they observe in the approximations to the limit-cycle amplitude as N increases.We have found that in particular cases the improvements in a N as N increases may be extremely slow.For example, with ε 1 the term with highest degree in ε dominates the location (or existence) of the value(s) of x for which 1 x dh N dt = 0.It can be seen from the recursion formula (24) that the term of highest degree of ε in 1 , where K is a constant, and that the value for which this term is zero is simply the value of x for which F (x) = 0.This is clearly illustrated in Table IV.Thus if ε 1, the value of x such that 1 x dh N dt = 0 is "trapped" near the value of x for which F (x) = 0 and does not rapidly converge to the limit-cycle amplitude as N increases.Empirically, we observe improvements in our upper bound, u N (see Table IV).In a few cases we are able to show analytically some convergence properties for a N and u N , but generally we, too, rely on empirical results.The particular cases where we can show improvements are all in the system (7) with conditions II A and F (x) convex for a < x < u.For the transition N = 2 → N = 4 we can show that that the bounds provided by both a N and u N improve (or remain the same) for all ε.To do this we use a combination of integration by parts, Jensen's inequality [41], and the intermediate value theorem.For system (7) with conditions II A and F (x) convex for a < x < u, we can also show that for N 2, u N tends to the recognized limit-cycle amplitude.
(iii) Both methods can generate bounds for the parameters associated with bifurcations.In many cases the method in Ref. [39] may converge more rapidly to the bifurcation set than do the zeros of the polynomials discussed in Sec.III C.However, our bounds for the bifurcation set are conservative, in the sense that they are based on sufficient conditions for the existence of a corresponding limit cycle.
(iv) In Sec.II E we show how our techniques can be applied to systems where the damping term has both odd and even components.By contrast, Neukirch [40] acknowledges that the technique used by Giacomini and Neukirch to approximate limit-cycle amplitudes (based on zeros for 1 x dh N dt ) does not work if the damping term in Eq. ( 5) is not an odd function of x.

V. CONCLUSION
In summary, we have introduced simple and sufficient conditions for the existence of limit cycles in Liénard systems.Our method is sufficiently simple to permit a valid, reasonably accurate upper bound for the amplitude of the limit cycle to be generated by inspection of the damping function of the Liénard system.Our method starts from geometric observations but primarily uses the techniques of vector algebra.The method does not rely on perturbation techniques and is applicable for all nonzero, finite values of the strength parameter ε.
We have shown that all systems of the form (7) which satisfy conditions II A possess a limit cycle.The maximal amplitude ρ of the limit cycle satisfies a < ρ < u, where F (a) = 0 and u 0 F (x)dx = 0. We have shown that, by utilizing a technique developed by Giacomini and Neukirch in a different context, it is possible to generate improved upper bounds for the amplitude of the limit cycle for the van der Pol oscillator Eq. (1).The criteria for establishing the existence of a limit cycle can also be used to develop bounds on bifurcation sets relating to the existence of limit cycles.
In many practical situations it is unlikely that the system being investigated will have the high symmetry of a Liénard equation, such as the van der Pol equation.Nevertheless, as shown in Sec.II E, there are many lower-symmetry systems for which bounds on the amplitude can be readily derived from those of a constituent Liénard system.The parameters of differential equations used to model physical systems are rarely known exactly and to assess systematically the consequences of parameter variations numerically could require extensive numerical simulations.By contrast, the results we present here can be used to provide simple bounds for the limit-cycle amplitude without the need for extensive numerical work.

FIG. 1 .
FIG. 1. (Color online) Right half-plane of the Liénard plane illustrating the integration path ABCA.The portion of the path ABC is a trajectory γ of the system; the portion CA lies along the y axis.The trajectory has its maximum x coordinate at the point B = (x,εF (x)).

FIG. 2 .
FIG. 2. (Color online) Right half-plane of the Liénard plane illustrating two trajectories γ i ,γ k crossing the line x = X.M k indicates the y coordinate of the upper portion of trajectory γ k at x = X,M i indicates the y coordinate of the upper portion of trajectory γ i at x = X,L k indicates the y coordinate of the lower portion of trajectory γ k at x = X, and L i indicates the y coordinate of the lower portion of trajectory γ i at x = X.

FIG. 3 .
FIG. 3. (Color online) Liénard plane illustrating the regions D and .is the region about the origin for which condition II E (v) applies.D is the annular region between the boundary of and the limit cycle 0 of Eq. (30) with α = 0.The arrows with straight tails illustrate the vector field ( E) α at the boundaries of D. The arrowheads without tails illustrate the vector field ( E) 0 on the limit cycle 0 .
= 0, where H satisfies conditions II B. In Sec.II C we showed that the GN polynomials h N (x,y) satisfy conditions II B. For the van der Pol equation we have used MAPLE10 symbolic computation software to generate the terms for ds also has a form that contains the key features of W A(x) and W B(x).The key features are (i) W A(0) = 0 and W B(0) = 0 (i.e., a root at x = 0), (ii) W A(x) and W B(x) have a single root (respectively at x = A and x = B) for x > 0, and

FIG. 4 .
FIG. 4. (Color online) W A(x) and W B(x) illustrate the form of each aggregate term in − x 0 1 s dh N dt ds for specific degree in ε and N = 2,4,6,8,10,12,14,16,18.In particular, each W A (0) = 0,W B(0) = 0,W A(x) < 0 for x 1,W B(x) < 0 for x 1 and W A(A) = 0,W B(B) = 0.The values A and B as illustrated correspond to the x values associated with a zero of a polynomial with specific N and specific degree in ε in the columns of Table IV marked "Upper."

4 = 4 4 4
0 gives roots x = 0 and x = ±2.1950.The positive root corresponds to the value in Table IV on the row ε 3 in the column "Upper h 4 ."Adding the terms with ε degree 1 and ε degree 3 we obtain − dt ds has the same form as W (x) with 4ε = σ 1 and 4ε 3 = σ 2 so applying the intermediate value theorem to dt ds between u = 2.1950 and u = 2.2361 we find for all finite, nonzero ε that dt ds = 0 for some u, where 2.1950 u 2.2361.All the polynomials for N = 2,4,10,12,14,16, and 18 and specific degree in ε have the same form as W A(x) and W B(x), so the value of x such that ds = 0 cannot lie outside the range of x values associated with zeros in Table IV for the corresponding value of N .

TABLE II .
GN Polynomials for N = 2 and N = 4.

TABLE III .
Derivatives of GN Polynomials for N = 2 and N = 4.

TABLE IV .
Locations of zeros derived from the corresponding GN polynomial.Values designated "Lower" correspond to the value of x for which 1 = 0 and values designated "Upper" correspond to the values of x for which dt