Formulas for the amplitude of the van der Pol limit cycle

The limit cycle of the van der Pol oscillator, $\ddot{x}+ \epsilon (x^2-1) \dot{x} + x =0$, is studied in the plane $(x,\dot{x})$ by applying the homotopy analysis method. A recursive set of formulas that approximate the amplitude and form of this limit cycle for the whole range of the parameter $\epsilon$ is obtained. These formulas generate the amplitude with an error less than 0.1%. To our knowledge, this is the first time where an analytical approximation of the amplitude of the van der Pol limit cycle, with validity from the weakly up to the strongly nonlinear regime, is given.


Introduction
The differential equation, x(t) + (x 2 − 1)ẋ(t) + x(t) = 0, t ≥ 0, (1.1) with a real parameter and the dot denoting the derivative with respect to time t, is called the van der Pol oscillator [1]. For > 0, and due to the nonlinear term (x 2 − 1)ẋ, the system accumulates energy in the region | x |< 1 and dissipates this energy in the region | x |> 1. This constraint implies the existence of an stable periodic motion (limit cycle [2]) when > 0. If the nonlinearity is increased, the dynamics in the time domain runs from near-harmonic oscillations when → 0 to relaxation oscillations when → ∞, making it a good model for many practical situations [3,4]. The closed curve representing this oscillation in the plane (x,ẋ) is quasi circular when → 0 and a sharp figure when → ∞. For < 0, the dynamics is dissipative in the region | x |< 1 and amplificative for | x |> 1. Under these conditions, the periodic motion is still possible but unstable. In this case, the limit cycle can be derived from that one with > 0 taking into account the symmetry ( , x(t)) → (− , −x(−t)). Therefore, it is enough to study the case > 0 to obtain also the behavior of the system for < 0.
Different standard methods (perturbative, non-perturbative, geometrical) [1,2,3,4,5,6,7] have permitted to study extensively the limit cycle of the van der Pol equation, in the weakly ( → 0) and in the strongly ( → ∞) nonlinear regimes. However, investigations giving analytical information of this object in the intermediate regime (0 ∞) are lacking in the literature. In this paper, it is our aim to fill in this gap by applying to the Eq. (1.1) the homotopy analysis method (HAM) introduced by Liao [8,9] in the nineties. This method has been shown to be powerful to solve different nonlinear problems [10,11,12,13,14,15,16]. In particular, it has been applied in [17] to Liénard equation,ẍ + f (x)ẋ + x = 0, which is the generalization of the van der Pol system when f (x) is an arbitrary function. As the interest in that work [17] was the amplitude and the frequency of the periodic motions, the calculations were performed with the time variable being explicit. As here we are only interested in the amplitude and form of the limit cycles, the time dependence of the solutions can be omitted, and we can work directly in the phase space (x,ẋ). Our results for the van der Pol limit cycle are presented in Section 2. In Section 3 it is explained how these results have been obtained from our specific application of the HAM to the van der Pol equation. Some conclusions are given in the last Section.
2 The amplitude of the van der Pol limit cycle Nowadays, the amplitude of the van der Pol limit cycle a can be easily approximated with a classical Runge-Kutta method. The results of this computational calculation a E are shown in Table 1. Let us observe that a E ( ) > 2 for > 0, with the asymptotic values: a E ( → 0) = a E ( → ∞) = 2. An upper bound rigorously established in [18] for the amplitude a is 2.3233. However, as it was signalled there [18], the maximum of a E is 2.02342 and it is obtained for = 3.3. In view of this result, Odani [18] conjectured that the amplitude of the limit cycle of the van der Pol equation is estimated by 2 < a( ) < 2.0235 for every > 0. If more precision is needed, the different analytical expansions in that have been found for the amplitude in the weakly and in the strongly nonlinear regimes can be considered. Evidently, the error becomes very large in the regions where these approximations are not valid. In [17,19,20], a recursive perturbation approximation is used to find the formula for the amplitude when → 0. This is, up to order O( 8 ), This expansion agrees for small with the computational calculation a E presented in Table 1. For 0 < < 2 the error is less than 1%, for ≈ 4 the error is bigger than 10%, and for ≈ 6 the formula has lost its validity and the error is bigger than 50%. In [21], the asymptotic dependence of the amplitude on is given for sufficiently large > 0, Compared with a E , this formula generates the amplitude with an error bigger than 15% for < 2. The formula starts to have validity for ≈ 10 with an error around 1%, that passes to be less than 0.1% when > 50. In Figs. 1.a and 1.b, formulas (2.2) and (2.3) are respectively plotted versus in their regions of validity.
We propose here two formulas, a R1 ( ) and a R2 ( ), that approximate the amplitude of the van der Pol limit cycle for every > 0. The first one with an error less than 0.1%, and the second one with an error less than 0.05%. These formulas have been obtained by applying the HAM to the van der Pol equation. The details of the derivation of these formulas are given in Section 3.
The first formula is , that derives from expression (3.30). The error obtained, |a( ) − a R1 ( )|/a( ), with this formula is less than 0.1% for every > 0. Fig. 2.1 shows a R1 in comparison with the 'experimental' amplitude a E . The maximum of a R1 is 2.02317 and it is taken for = 3.3.
The expansion of expression (2.5) for small gives For large , we obtain Let us note the different scaling of this last expression (with behavior −1 ) respect to expansion (2.3) (with behavior −1.33 ). Taking into account that these approxi- mations start to be valid when > 100, it can be easily seen that this difference is negligible, in fact less than 0.01%. Nevertheless, we have tried to modify formula (2.5) to obtain the correct scaling −1.33 of expression (2.3), but the increase of the error in other regions of the parameter does not recommend this possibility.
In summary, let us remark the exceptional fit of the 'experimental' points a E by the amplitudes a R1 and a R2 generated with formulas (2.4) and (2.5), respectively (see Figs. 2a and 2b). Moreover, by inspection of Fig. 2b, let us finish this section by posing the following Conjecture: the amplitude a R2 ( ) obtained with formula (2.5) is an upper bound for the amplitude a( ) of the van der Pol limit cycle, that is The HAM and the van der Pol equation The generalization of the van der Pol equation,ü + (u 2 − 1)u + u = 0, is called the Liénard equation. In coordinates (u, v), it reads where we consider that f (u) is an even function. The HAM has been applied to this equation in [17] working explicitly in the time domain. We proceed now to apply the HAM to this system by following a different strategy, namely, by omitting the time variable of Eq. (3.1).

HAM applied to Liénard equation
If we define the variable v =u, and suppose that v is a function of u, the acceleration u can be rewritten as v v, where the prime denotes the derivative respect to u. The result is a new equation with the time variable eliminated: When f (u) is even [22,23], the limit cycles of this equation are closed trajectories of amplitude a in the (u, v) plane, with −a ≤ u ≤ a and v(a) = v(−a) = 0.
After the re-scaling (u, v) = (ax, ay) to the new variables (x, y), the limit cycles are transformed in solutions of the equation where the prime is now denoting the derivative respect to x.
From [22,23], we know that an approximate solution of the positive branch solution of (3.3), for any value of > 0, is withã > 0. The number of possible stable limit cycles in the strongly non-linear regime (for large ) is obtained from the number of positive roots of F (x 0 ) − F (x) when x 0 runs over the set of all the negative roots of f (x 0 ) = 0. We take the positive rootsã selected by the algorithm explained in [22,23] as the corresponding approximate amplitudes. Let us observe that the initial approximate limit cycleỹ(x) recovers the exact circular form when → 0, and also, when → ∞, it becomes exactly the two-piecewise limit cycle proposed in [22]. Moreover, it must be noticed thatỹ(x) has sense only when > 0, and then the results obtained from this initial guess only will have validity for > 0.
After eliminating the time variable, the non-linear differential equation yy + yf (ax)+ x = 0 is of order one. Then, to construct the homotopy [8,9], we build a linear differential equation of order one for the whole range [0, 1] of what is called the homotopy parameter p.
We define an auxiliary linear operator L by with the property where C 1 is a constant, and p is a (homotopy) parameter explained below.

From (3.3) we define a nonlinear operator
and then construct the homotopy where h is a nonzero auxiliary parameter. Setting H[φ(x; p), A(p)] = 0, we have the zero-order deformation equation where y n (x) = 1 n! ∂ n φ(x; p) ∂p n p=0 , a n = 1 n! ∂ n A(p) ∂p n p=0 .
Notice that series (3.11) contain the auxiliary parameter h, which has influence on their convergence regions. Assume that h is properly chosen such that all of these Maclaurin series are convergent at p = 1. Hence at p = 1 we have At the N th-order approximation, we have the analytic solution of Eq. The parameter h is free and can be chosen arbitrarily, in particular, it can be a function of . Nevertheless, the functionỹ(x) is the exact limit cycle solution of the Liénard equation (3.3) in the limits → 0 and → ∞. This means that the general solution y(x, h) of Eq. (3.9) should tend toỹ(x) for any p in those limits of , loosing its dependence on h. Then, a reasonable property for h would be to vanish in those limits of . Hence, the solution of (3.9) would be the exact solutioñ y(x) for any value of the parameter p in the limits → 0 and → ∞. A simple function h( ) satisfying these properties is where R n (x) is defined by and At each iteration, we have two unknowns, C 1 , in Eq. (3.6), and a n . These unknowns are obtained by considering the boundary conditions (3.15) as follows.
(3.17) At 1th-order, we obtain where we have imposed the condition y 1 (−1) = 0 choosing the integration constant, C 1 , appropriately. At this moment a 0 is free, but we can fix the value of a 0 by imposing y 1 (1) = 0: This is a non linear equation for a 0 . When f (x) is an even polynomial of degree 2n, it is a polynomial equation of degree n in the variable a 2 0 . Therefore, at this order, the number of limit cycles is the number of positive roots of this equation, at most n (see [23] for a longer discussion of this point).
Then, for every one of the a 0 solutions of the above equation, i.e. for every one of the limit cycles of the system, we proceed order by order in p to obtain higher order corrections to the amplitude a and to the shape of the limit cycle y(x). At any nth-order, with n ≥ 2, we obtain the nth-order correction y n to the shape of the limit cycle: where B n,m (a 1 , ..., a n ) are the partial ordinary Bell polynomials [ [24], p. 190]. We have B 0,0 = 1 and, for n = 1, 2, 3, ..., B n+1,0 = 0 and they satisfy the recurrence B n,m (a 1 , ..., a n ) = n−1 k=m−1 a n−k B k,m−1 (a 1 , ..., a n ). (3.21) We recall that B n,m (a 1 , ..., a n ) is linear in a n .
At nth-order, we can obtain also the nth-order correction a n−1 to the amplitude of the limit cycle by imposing the condition y n (1) = 0. The result is a linear equation for a n−1 with solution: 22) where the star in the sum symbol means that the term (k, m) = (n − 1, 1) is excluded.
For example, at 2th-order we obtain and Therefore, a first approximation to the shape of the limit cycle is y 0 (x) + y 1 (x), and a first order approximation to the amplitude of the limit cycle is a = a 0 + a 1 .
At 3th-order, we have and (3.26) Therefore, a second order approximation to the shape of the limit cycle is y 0 (x) + y 1 (x) + y 2 (x) and a second order approximation to the amplitude of the limit cycle is a = a 0 + a 1 + a 2 .
Moreover, at every order in p the appropriate function h( ) is indicated by the experiment.

Results for the van der Pol equation
For the van der Pol system, we have .
Therefore, a first approximation to the amplitude of the limit cycle is .    the formula a R2 ( ) shown in (2.5) is finally obtained.
As an example, and to conclude this section, we plot in Fig. 3 the form of the van der Pol limit cycle when y 0 (x), y 0 (x) + y 1 (x) and y 0 (x) + y 1 (x) + y 2 (x) are used as approximated limit cycles. It is not banal to recall here that the reconstruction of the shape of limit cycles is not a problem less difficult than that of the calculation of their amplitudes.
In this work, the conjecture (2.8) has been posed. This is a consequence of the different formulas here presented for approximating the amplitude a of the van der Pol limit cycle. In addition to the well known constant approximationā = 2 that generates an error less than 1%, we establish a family of recursive formulas that are valid for the whole range of the parameter . Two of them, a R1 ( ) and a R2 ( ), have been explicitly given. The first one, a R1 , produces an error less than 0.1% and the second one, a R2 , reduces the error to less than 0.05%. Moreover, a R2 is conjectured to be an upper bound of a.
As far as we know, this is the first time where an analytical approximation of the amplitude of the van der Pol limit cycle, with validity from the weakly up to the strongly nonlinear regime, is proposed.