# Formulas for the Amplitude of the van der Pol Limit Cycle through the Homotopy Analysis Method

Authors: J. L. López; S. Abbasbandy; R. López-Ruiz
Journal: Scholarly Research Exchange
 (2009)

---

## Abstract

The limit cycle of the van der Pol oscillator,x¨+ε(x2−1)x˙+x=0, is studied in the plane (x,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 ε 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.


---

## Body

## 1. Introduction

A dynamical  system whose time evolutionx(t) is determined
by the differential equation,(1)x¨(t)+ϵ(x2−1)x˙(t)+x(t)=0,t≥0,with ϵ a real
parameter, and the dot denoting the time derivative is called the van der Pol
oscillator [1]. For ϵ>0, and due to the nonlinear term ϵ(x2−1)x˙, the system accumulates energy in the region |x|<1 and dissipates
this energy in the region |x|>1. This constraint implies the existence of a 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,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,
nonperturbative, geometrical) [1–7] have been used to study 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 of ϵ are lacking in
literature. In this paper, it is our aim to fill in this gap by applying to
(1) the homotopy analysis method (HAM) introduced by Liao [8, 9] in the nineties. This method
has been shown to be very useful to solve different nonlinear problems
[10–20]. In particular, it has been
applied in [21] to
Liénard equation, x¨+ϵf(x)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 [21] 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,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 amplitudea of the van der
Pol limit cycle can be easily computed with a classical Runge-Kutta method. The
results of this computational calculation aE are shown in
Table 1. Let us observe that aE(ϵ)>2 for ϵ>0, with the asymptotic values: aE(ϵ→0)=aE(ϵ→∞)=2. An upper bound rigorously established in [22] for the amplitude a is 2.3233. However, as it was also signaled in that work, the
maximum of aE is 2.02342, and it is
obtained for ϵ=3.3. In view of this result, Odani [22] conjectured that the amplitude of the limit cycle of the van
der Pol equation is estimated by2<a(ϵ)<2.0235for everyϵ>0.Table 1
The valueaE represents the
amplitude a of the van der
Pol limit cycle obtained with five decimal digits by integrating directly (1) with a fourth order Runge-Kutta method for the indicated values of ϵ.
ϵ0.10.51.01.52.02.5aE2.000102.002492.008622.015222.019892.02235ϵ3.03.54.05.01050aE2.023302.023372.022962.021512.014282.00295A closed formula for the amplitudea as function of ϵ is unknown. By
inspecting Table 1, one could propose as a good solution the constant
amplitude(2)a¯=2  for the whole
range of the parameter ϵ. Knowing that the experimental upper bound for a is 2.02342,
that is, 2<a(ϵ)<2.02342 for every ϵ>0, the error made with this approximation, (a(ϵ)−a¯)/a(ϵ), would be about 1%.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 [21, 23, 24], a recursive perturbation
approximation is used to find the formula for the amplitude when ϵ→0. This is, up to order 𝒪(ϵ8),(3)a(ϵ)ϵ→0=2+196ϵ2−1033552960ϵ4+101968955738368000ϵ6.This expansion agrees for small ϵ with the
computational calculation aE 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 [25], the asymptotic dependence of the amplitude on ϵ is given for
sufficiently large ϵ>0,(4)a(ϵ)ϵ→∞=2+0.7793ϵ−4/3.Compared with aE, 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 Figures 1(a) and 1(b), formulas (3) and (4) are,
respectively, plotted versus ϵ in their
regions of validity.Comparison of the “experimental” amplitudeaE(ϵ) (dotted curve)
with (a) amplitude a(ϵ)ϵ→0 given by
formula (3) (solid line) and (b) amplitude a(ϵ)ϵ→∞ given by
formula (4) (solid line).
(a)(b)We propose here two formulas,aR1(ϵ) and aR2(ϵ), 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(5)aR1(ϵ)=2+1.737ϵ2(8π+9ϵ)(4+ϵ2),that derives from  expression
(38). The error obtained, |a(ϵ)−aR1(ϵ)|/a(ϵ), with this formula is less than 0.1% for every ϵ>0. Figure 2(a)
shows aR1 in comparison
with the “experimental” amplitude aE. The maximum of aR1 is 2.02317, and it is taken
for ϵ=3.3.(a) AmplitudeaR1(ϵ) given by
formula (5) (solid curve) and experimental amplitude aE(ϵ). (b) Amplitude aR2(ϵ) given by
formula (6) (solid curve) and experimental amplitude aE(ϵ) (dotted curve),
and ϵ given in
logarithmic scale.
(a)(b)The second formula is(6)aR2(ϵ)=2+0.74958ϵ2(8π+9ϵ)(4+ϵ2)+ϵ2(75.3562+43.0023ϵ+28.1589ϵ2+8.3479ϵ3)(8π+9ϵ)2(4+ϵ2)2,that derives from expression
(41). The error obtained, |a(ϵ)−aR2(ϵ)|/a(ϵ), with this formula is less than 0.05% for every ϵ>0. The maximum of aR2 is 2.02346, and it is taken
for ϵ=3.29482. The plot of aR2 in comparison
with the amplitude aE obtained
computationally can be seen in Figure 2(b). Let us observe that it can be guessed
from Figure 2(b) that aR2(ϵ)>aE(ϵ) for every ϵ>0.The expansion of expression (6) for small ϵ gives(7)aR2(ϵ)ϵ→0=2+0.01491ϵ2+𝒪(ϵ3).For large ϵ, we obtain(8)aR2(ϵ)ϵ→∞=2+0.18648ϵ−1+𝒪(ϵ−2).Let us note the different
scaling of this last expression (with behavior ϵ−1) with respect to
expansion (4) (with behavior ϵ−1.33). Taking into
account that these approximations 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 (6)
to obtain the correct scaling ϵ−1.33 of expression
(4) 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” pointsaE by the
amplitudes aR1 and aR2 generated with
formulas (5) and (6), respectively, (see Figures 2(a) and 2(b)). Moreover, by
inspection of Figure 2(b), let us finish this section by posing the following.Conjecture
The
amplitudeaR2(ϵ)
obtained with
formula (6) is an upper bound for the amplitude  a(ϵ)
of the van der
Pol limit cycle, that is,(9)aR2(ϵ)>a(ϵ)>2for everyϵ>0.


## 3. The HAM and the van der Pol Equation

The generalization of the van der Pol equation,u¨+ϵ(u2−1)u˙+u=0, is called the Liénard equation. In coordinates (u,v), it reads(10)u¨(t)+ϵf(u)u˙(t)+u(t)=0,t≥0,where we consider that f(u) is an even
function. The HAM has been applied to this equation in [21] 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 (10).


### 3.1. HAM Applied to Liénard Equation

If we define the variablev=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:(11)vv′+ϵf(u)v+u=0.When f(u) is even
[26, 27], 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 rescaling(u,v)=(ax,ay) to the new
variables (x,y), the limit cycles are transformed in solutions of the
equation(12)yy′+ϵyf(ax)+x=0,−1≤x≤1,y(−1)=y(1)=0,where the prime is now denoting
the derivative with respect to x.From [26, 27], we know that an approximate solution of the positive
branch solution of (12), for any value of ϵ>0, is(13)y˜(x)=1−x2+ϵa˜{0if−1≤x<x0a˜,F(a˜)−F(a˜x)ifx0a˜≤x≤1,where F′(x)=f(x), x0 is a negative
root of f(x): f(x0)=0 with x0<0, and a˜ is a positive
root of F(x0)−F(x): F(x0)=F(a˜) with a˜>0. The number of possible stable limit cycles in the
strongly nonlinear regime (for large ϵ) is obtained
from the number of positive roots of F(x0)−F(x) when x0 runs over the
set of all the negative roots of f(x0)=0. We take the positive roots a˜ selected by the
algorithm explained in [26, 27] as the corresponding approximate amplitudes. Let us
observe that the initial approximate limit cycle y˜(x) is justified
because it recovers the exact circular form when ϵ→0, and also, when ϵ→∞, it becomes exactly the two-piecewise limit cycle
proposed and plotted in [26]. Moreover, it must be noticed that y˜(x) makes 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 nonlinear
differential equationyy′+ϵ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ℒ by(14)ℒ[ϕ(x;p)]=∂∂xϕ(x;p),with the
property(15)ℒ[C1]=0,where C1 is a constant
function (the kernel of operator ℒ), and p is a (homotopy)
parameter explained below.From (12) we define a nonlinear
operator(16)𝒩[ϕ(x;p),A(p)]=ϕ(x;p)∂ϕ(x;p)∂x+ϵϕ(x;p)f(A(p)x)+x,and then construct the
homotopy(17)ℋ[ϕ(x;p),A(p)]=(1−p)ℒ[ϕ(x;p)−y˜(x)]+hp𝒩[ϕ(x;p),A(p)],where h is a nonzero
auxiliary parameter. Setting ℋ[ϕ(x;p),A(p)]=0, we have the zero-order deformation
equation(18)(1−p)ℒ[ϕ(x;p)−y˜(x)]+hp𝒩[ϕ(x;p),A(p)]=0,subject to the boundary
conditions(19)ϕ(−1;p)=0,ϕ(1;p)=0,where p∈[0,1] is an embedding
parameter. When the parameter p increases from
0 to 1, the solution ϕ(x;p) varies from y0(x)=y˜(x) to y(x), and A(p) varies from a0 to a. Assume that ϕ(x;p) and A(p) are analytic in p∈[0,1] and can be
expanded in the Maclaurin series of p as
follows:(20)ϕ(x;p)=∑n=0+∞yn(x)pn,A(p)=∑n=0+∞anpn,where(21)yn(x)=1n!∂nϕ(x;p)∂pn|p=0,an=1n!∂nA(p)∂pn|p=0.Notice that series (20)
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(22)y(x)=y0(x)+∑n=1+∞yn(x),a=a0+∑n=1+∞an.At the Nth-order
approximation, we have the analytic solution of (12),
namely,(23)y(x)≈YN(x):=∑n=0Nyn(x),a≈AN:=∑n=0Nan. The parameter h is free and can
be chosen arbitrarily; in particular, it can be a function of ϵ. Nevertheless, the function y˜(x) is the exact
limit cycle solution of the Liénard equation (12) in the limits ϵ→0 and ϵ→∞. This means that the general solution y(x,h) of (18)
should tend to y˜(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 (18) would be the exact
solution y˜(x) for any value
of the parameter p in the limits ϵ→0 and ϵ→∞. A simple function h(ϵ) satisfying
these properties is(24)h(ϵ)=bϵc+ϵ2,withb,c∈R. Differentiating (18) and (19) n times with
respect to p, then setting p=0, and finally dividing by n!, we obtain the nth-order
deformation equation(25)ℒ[yn(x)−χnyn−1(x)]+hRn(x)=0,(n=1,2,3,…),subject to the boundary
conditions(26)yn(−1)=0,yn(1)=0,where Rn(x) is defined
by(27)Rn(x)=1(n−1)!∂n−1𝒩[ϕ(x;p),A(p)]∂pn−1|p=0,χn={0,n≤1,1,n>1.At each iteration, we have two
unknowns: C1, in (15), and an. These unknowns are obtained by considering the
boundary conditions (26) as follows.At zero order, we obtain(28)y0(x)=y˜(x).At 1st-order, we
obtain(29)y1(x)=h2[1−x2−y02(x)]−hϵ∫−1xy0(t)f(a0t)dt,where we have imposed the
condition y1(−1)=0 choosing the
integration constant, C1, appropriately. At this moment a0 is free but we
can fix the value of a0 by imposing y1(1)=0:(30)∫−11y0(t)f(a0t)dt=0.This is a nonlinear equation
for a0. When f(x) is an even
polynomial of degree 2n, it is a polynomial equation of degree n in the variable a02. Therefore, at this order, the number of limit cycles
is the number of positive roots of this equation, at most n (see [27] for a longer discussion of
this point).Then, for every one of thea0 solutions of
the above equation, that is, 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).For example, at 2th-order, we obtain(31)y2(x)=y1(x)−hy0(x)y1(x)−hϵ∫−1x[y1(t)f(a0t)+ty0(t)f′(a0t)a1]dt,(32)a1=−[∫−11y1(x)f(a0x)dx][∫−11xy0(x)f′(a0x)dx]. Therefore, a first approximation to the shape of the
limit cycle is y0(x)+y1(x), and a first-order approximation to the amplitude of
the limit cycle is a=a0+a1.At 3th-order, we have(33)y3(x)=y2(x)−hy0(x)y2(x)−h2y12(x)−hϵ∫−1x[y2(t)f(a0t)+a1ty1(t)f′(a0t)+12a12t2y0(t)f′′(a0t)+ty0(t)f′(a0t)a2]dt,(34)a2=−[∫−11[y2(x)f(a0x)+a1xy1(x)f′(a0x)+12a12x2y0(x)f′′(a0x)]dx]/[∫−11xy0(x)f′(a0x)dx].Therefore, a second-order
approximation to the shape of the limit cycle is y0(x)+y1(x)+y2(x), and a second-order approximation to the amplitude of the limit cycle is a=a0+a1+a2.Moreover, at every order inp the appropriate
function h(ϵ) is indicated by
the experiment, and in some sense, h(ϵ) should be seen
as somewhat provisional if a better guess could be inferred.


### 3.2. Results for the van der Pol Equation

For the van der Pol system, we havef(x)=x2−1, F(x)=x3/3−x, x0=−1, a˜=2 and(35)y˜(x)=1−x2+ϵ2{0if−1≤x<−12,F(2)−F(2x)if−12≤x≤1. The expressions for y0(x), y1(x), and y2(x) can be
calculated as explained in Section 3.1. From (30) we get(36)∫−11y0(t)f(a0t)dt=9ϵ+8π64(a02−4)=0⇒a0=2. From (32) we obtain(37)a1(h)=−273ϵh35(8π+9ϵ). Therefore, a first approximation to the amplitude of
the limit cycle is(38)a(h)≃2−273ϵh35(8π+9ϵ).As it was advanced at the
beginning of this section, note that this expression is only valid for positive ϵ. Choosing b=−1.3 and c=4 in
(24),(39)h(ϵ)=−1.3ϵ4+ϵ2,the formula aR1(ϵ) given in (5)
is obtained. The  values for b and c are selected by
guess and check.From (34) we get(40)a2(h)=(3ϵh[191973ϵ2h−5040(83−21h)π+2ϵ(−226803+h(54837+77003π))])/19600(9ϵ+8π)2. Therefore, a second approximation to the amplitude of
the limit cycle is(41)a(h)≃2−273ϵh35(8π+9ϵ)+(3ϵh[191973ϵ2h−5040(83−21h)π+2ϵ(−226803+h(54837+77003π))])/19600(9ϵ+8π)2.Choosing b=−0.561 and c=4 in
(24),(42)h(ϵ)=−0.561ϵ4+ϵ2,the formula aR2(ϵ) shown in (6)
is finally obtained.As an example, and  to conclude this section, we plot
in Figure3 the form of the van der Pol limit cycle when y0(x), y0(x)+y1(x), and y0(x)+y1(x)+y2(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 less difficult problem than that of the
calculation of their amplitudes.Figure 3
Approximated shape of the van der
Pol limit cycle forϵ=5 and h(ϵ) given in
(42). Black line: the experimental limit cycle. Brown, blue, and red lines:
the curves y0(x), y0+y1(x) and y0+y1(x)+y2(x), respectively.


## 3.1. HAM Applied to Liénard Equation

If we define the variablev=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:(11)vv′+ϵf(u)v+u=0.When f(u) is even
[26, 27], 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 rescaling(u,v)=(ax,ay) to the new
variables (x,y), the limit cycles are transformed in solutions of the
equation(12)yy′+ϵyf(ax)+x=0,−1≤x≤1,y(−1)=y(1)=0,where the prime is now denoting
the derivative with respect to x.From [26, 27], we know that an approximate solution of the positive
branch solution of (12), for any value of ϵ>0, is(13)y˜(x)=1−x2+ϵa˜{0if−1≤x<x0a˜,F(a˜)−F(a˜x)ifx0a˜≤x≤1,where F′(x)=f(x), x0 is a negative
root of f(x): f(x0)=0 with x0<0, and a˜ is a positive
root of F(x0)−F(x): F(x0)=F(a˜) with a˜>0. The number of possible stable limit cycles in the
strongly nonlinear regime (for large ϵ) is obtained
from the number of positive roots of F(x0)−F(x) when x0 runs over the
set of all the negative roots of f(x0)=0. We take the positive roots a˜ selected by the
algorithm explained in [26, 27] as the corresponding approximate amplitudes. Let us
observe that the initial approximate limit cycle y˜(x) is justified
because it recovers the exact circular form when ϵ→0, and also, when ϵ→∞, it becomes exactly the two-piecewise limit cycle
proposed and plotted in [26]. Moreover, it must be noticed that y˜(x) makes 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 nonlinear
differential equationyy′+ϵ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ℒ by(14)ℒ[ϕ(x;p)]=∂∂xϕ(x;p),with the
property(15)ℒ[C1]=0,where C1 is a constant
function (the kernel of operator ℒ), and p is a (homotopy)
parameter explained below.From (12) we define a nonlinear
operator(16)𝒩[ϕ(x;p),A(p)]=ϕ(x;p)∂ϕ(x;p)∂x+ϵϕ(x;p)f(A(p)x)+x,and then construct the
homotopy(17)ℋ[ϕ(x;p),A(p)]=(1−p)ℒ[ϕ(x;p)−y˜(x)]+hp𝒩[ϕ(x;p),A(p)],where h is a nonzero
auxiliary parameter. Setting ℋ[ϕ(x;p),A(p)]=0, we have the zero-order deformation
equation(18)(1−p)ℒ[ϕ(x;p)−y˜(x)]+hp𝒩[ϕ(x;p),A(p)]=0,subject to the boundary
conditions(19)ϕ(−1;p)=0,ϕ(1;p)=0,where p∈[0,1] is an embedding
parameter. When the parameter p increases from
0 to 1, the solution ϕ(x;p) varies from y0(x)=y˜(x) to y(x), and A(p) varies from a0 to a. Assume that ϕ(x;p) and A(p) are analytic in p∈[0,1] and can be
expanded in the Maclaurin series of p as
follows:(20)ϕ(x;p)=∑n=0+∞yn(x)pn,A(p)=∑n=0+∞anpn,where(21)yn(x)=1n!∂nϕ(x;p)∂pn|p=0,an=1n!∂nA(p)∂pn|p=0.Notice that series (20)
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(22)y(x)=y0(x)+∑n=1+∞yn(x),a=a0+∑n=1+∞an.At the Nth-order
approximation, we have the analytic solution of (12),
namely,(23)y(x)≈YN(x):=∑n=0Nyn(x),a≈AN:=∑n=0Nan. The parameter h is free and can
be chosen arbitrarily; in particular, it can be a function of ϵ. Nevertheless, the function y˜(x) is the exact
limit cycle solution of the Liénard equation (12) in the limits ϵ→0 and ϵ→∞. This means that the general solution y(x,h) of (18)
should tend to y˜(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 (18) would be the exact
solution y˜(x) for any value
of the parameter p in the limits ϵ→0 and ϵ→∞. A simple function h(ϵ) satisfying
these properties is(24)h(ϵ)=bϵc+ϵ2,withb,c∈R. Differentiating (18) and (19) n times with
respect to p, then setting p=0, and finally dividing by n!, we obtain the nth-order
deformation equation(25)ℒ[yn(x)−χnyn−1(x)]+hRn(x)=0,(n=1,2,3,…),subject to the boundary
conditions(26)yn(−1)=0,yn(1)=0,where Rn(x) is defined
by(27)Rn(x)=1(n−1)!∂n−1𝒩[ϕ(x;p),A(p)]∂pn−1|p=0,χn={0,n≤1,1,n>1.At each iteration, we have two
unknowns: C1, in (15), and an. These unknowns are obtained by considering the
boundary conditions (26) as follows.At zero order, we obtain(28)y0(x)=y˜(x).At 1st-order, we
obtain(29)y1(x)=h2[1−x2−y02(x)]−hϵ∫−1xy0(t)f(a0t)dt,where we have imposed the
condition y1(−1)=0 choosing the
integration constant, C1, appropriately. At this moment a0 is free but we
can fix the value of a0 by imposing y1(1)=0:(30)∫−11y0(t)f(a0t)dt=0.This is a nonlinear equation
for a0. When f(x) is an even
polynomial of degree 2n, it is a polynomial equation of degree n in the variable a02. Therefore, at this order, the number of limit cycles
is the number of positive roots of this equation, at most n (see [27] for a longer discussion of
this point).Then, for every one of thea0 solutions of
the above equation, that is, 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).For example, at 2th-order, we obtain(31)y2(x)=y1(x)−hy0(x)y1(x)−hϵ∫−1x[y1(t)f(a0t)+ty0(t)f′(a0t)a1]dt,(32)a1=−[∫−11y1(x)f(a0x)dx][∫−11xy0(x)f′(a0x)dx]. Therefore, a first approximation to the shape of the
limit cycle is y0(x)+y1(x), and a first-order approximation to the amplitude of
the limit cycle is a=a0+a1.At 3th-order, we have(33)y3(x)=y2(x)−hy0(x)y2(x)−h2y12(x)−hϵ∫−1x[y2(t)f(a0t)+a1ty1(t)f′(a0t)+12a12t2y0(t)f′′(a0t)+ty0(t)f′(a0t)a2]dt,(34)a2=−[∫−11[y2(x)f(a0x)+a1xy1(x)f′(a0x)+12a12x2y0(x)f′′(a0x)]dx]/[∫−11xy0(x)f′(a0x)dx].Therefore, a second-order
approximation to the shape of the limit cycle is y0(x)+y1(x)+y2(x), and a second-order approximation to the amplitude of the limit cycle is a=a0+a1+a2.Moreover, at every order inp the appropriate
function h(ϵ) is indicated by
the experiment, and in some sense, h(ϵ) should be seen
as somewhat provisional if a better guess could be inferred.


## 3.2. Results for the van der Pol Equation

For the van der Pol system, we havef(x)=x2−1, F(x)=x3/3−x, x0=−1, a˜=2 and(35)y˜(x)=1−x2+ϵ2{0if−1≤x<−12,F(2)−F(2x)if−12≤x≤1. The expressions for y0(x), y1(x), and y2(x) can be
calculated as explained in Section 3.1. From (30) we get(36)∫−11y0(t)f(a0t)dt=9ϵ+8π64(a02−4)=0⇒a0=2. From (32) we obtain(37)a1(h)=−273ϵh35(8π+9ϵ). Therefore, a first approximation to the amplitude of
the limit cycle is(38)a(h)≃2−273ϵh35(8π+9ϵ).As it was advanced at the
beginning of this section, note that this expression is only valid for positive ϵ. Choosing b=−1.3 and c=4 in
(24),(39)h(ϵ)=−1.3ϵ4+ϵ2,the formula aR1(ϵ) given in (5)
is obtained. The  values for b and c are selected by
guess and check.From (34) we get(40)a2(h)=(3ϵh[191973ϵ2h−5040(83−21h)π+2ϵ(−226803+h(54837+77003π))])/19600(9ϵ+8π)2. Therefore, a second approximation to the amplitude of
the limit cycle is(41)a(h)≃2−273ϵh35(8π+9ϵ)+(3ϵh[191973ϵ2h−5040(83−21h)π+2ϵ(−226803+h(54837+77003π))])/19600(9ϵ+8π)2.Choosing b=−0.561 and c=4 in
(24),(42)h(ϵ)=−0.561ϵ4+ϵ2,the formula aR2(ϵ) shown in (6)
is finally obtained.As an example, and  to conclude this section, we plot
in Figure3 the form of the van der Pol limit cycle when y0(x), y0(x)+y1(x), and y0(x)+y1(x)+y2(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 less difficult problem than that of the
calculation of their amplitudes.Figure 3
Approximated shape of the van der
Pol limit cycle forϵ=5 and h(ϵ) given in
(42). Black line: the experimental limit cycle. Brown, blue, and red lines:
the curves y0(x), y0+y1(x) and y0+y1(x)+y2(x), respectively.


## 4. Conclusions

In this work, the conjecture (9) 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 a¯=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, aR1(ϵ) and aR2(ϵ), have been explicitly given. The first one, aR1, produces an error less than 0.1%, and the second
one, aR2, reduces the error to less than 0.05%. Moreover, aR2 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.



---

*Source: 854060-2009-04-22.xml*