Copyright © 2009 J. L. López et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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.
1. Introduction
A dynamical system whose time evolution 
is determined
by the differential equation,
(1)with 
a real
parameter, and the dot denoting the time derivative is called the van der Pol
oscillator [1]. For 
, and due to the nonlinear term 
, the system accumulates energy in the region 
and dissipates
this energy in the region 
. This constraint implies the existence of a stable
periodic motion (limit cycle [2]) when 
. If the nonlinearity is increased, the dynamics in
the time domain runs from near-harmonic oscillations when 
to relaxation
oscillations when 
, making it a good model for many practical situations
[3, 4]. The closed curve
representing this oscillation in the plane 
is quasi
circular when 
and a sharp
figure when 
. For 
, the dynamics is dissipative in the region 
and
amplificative for 
. Under these conditions, the periodic motion is still
possible but unstable. In this case, the limit cycle can be derived from that
one with 
taking into
account the symmetry 
. Therefore, it is enough to study the case 
to obtain also
the behavior of the system for 
.
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 (
) 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, 
, which is the generalization of the van der Pol
system when 
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 
. 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 can be easily computed with a classical Runge-Kutta method. The
results of this computational calculation 
are shown in
Table 1. Let us observe that 
for 
, with the asymptotic values: 
. An upper bound rigorously established in [22] for the amplitude 
is 
. However, as it was also signaled in that work, the
maximum of 
is 
, and it is
obtained for 
. In view of this result, Odani [22] conjectured that the amplitude of the limit cycle of the van
der Pol equation is estimated by 
for every 
.
Table 1: The value
represents the
amplitude
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
.
A closed formula for the amplitude 
as function of 
is unknown. By
inspecting Table 1, one could propose as a good solution the constant
amplitude
(2) for the whole
range of the parameter 
. Knowing that the experimental upper bound for 
is 2.02342,
that is, 
for every 
, the error made with this approximation, 
, would be about 
.
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 
. This is, up to order 
,
(3)This expansion agrees for small 
with the
computational calculation 
presented in
Table 1. For 
, the error is
less than 
, for 
, the error is
bigger than 
, and for 
, the formula has
lost its validity, and the error is bigger than 
. In [25], the asymptotic dependence of the amplitude on 
is given for
sufficiently large 
,
(4)Compared with 
, this formula generates the amplitude with an error
bigger than 
for 
. The formula starts to have validity for 
with an error
around 
, that passes to be less than 
when 
. In Figures 1(a) and 1(b), formulas (3) and (4) are,
respectively, plotted versus 
in their
regions of validity.
Figure 1: Comparison of the “experimental” amplitude
(dotted curve)
with (a) amplitude
given by
formula (
3) (solid line) and (b) amplitude
given by
formula (
4) (solid line).
We propose here two formulas, 
and 
, that approximate the amplitude of the van der Pol
limit cycle for every 
: the first one with an error less than 
and the second one with an error less than 
. 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)that derives from expression
(38). The error obtained, 
, with this formula is less than 
for every 
. Figure 2(a)
shows 
in comparison
with the “experimental” amplitude 
. The maximum of 
is 
, and it is taken
for 
.
Figure 2: (a) Amplitude
given by
formula (
5) (solid curve) and experimental amplitude
. (b) Amplitude
given by
formula (
6) (solid curve) and experimental amplitude
(dotted curve),
and
given in
logarithmic scale.
The second formula is
(6)that derives from expression
(41). The error obtained, 
, with this formula is less than 
for every 
. The maximum of 
is 
, and it is taken
for 
. The plot of 
in comparison
with the amplitude 
obtained
computationally can be seen in Figure 2(b). Let us observe that it can be guessed
from Figure 2(b) that 
for every 
.
The expansion of expression (6) for small 
gives
(7)For large 
, we obtain
(8)Let us note the different
scaling of this last expression (with behavior 
) with respect to
expansion (4) (with behavior 
). Taking into
account that these approximations start to be valid when 
, it can be easily seen that this difference is
negligible, in fact less than 
. Nevertheless, we have tried to modify formula (6)
to obtain the correct scaling 
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” points 
by the
amplitudes 
and 
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
amplitude 
obtained with
formula (6) is an upper bound for the amplitude 
of the van der
Pol limit cycle, that is,
(9)
3. The HAM and the van der Pol Equation
The generalization of the van der Pol equation, 
, is called the Liénard equation. In coordinates 
, it reads
(10)where we consider that 
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 variable 
, and suppose that 
is a function
of 
, the acceleration 
can be
rewritten as 
, where the prime denotes the derivative respect to 
. The result is a new equation with the time variable
eliminated:
(11)When 
is even
[26, 27], the limit cycles of this
equation are closed trajectories of amplitude 
in the 
plane, with 
and 
.
After the rescaling 
to the new
variables 
, the limit cycles are transformed in solutions of the
equation
(12)where the prime is now denoting
the derivative with respect to 
.
From [26, 27], we know that an approximate solution of the positive
branch solution of (12), for any value of 
, is
(13)where 
is a negative
root of 
with 
, and 
is a positive
root of 
with 
. The number of possible stable limit cycles in the
strongly nonlinear regime (for large 
) is obtained
from the number of positive roots of 
when 
runs over the
set of all the negative roots of 
. We take the positive roots 
selected by the
algorithm explained in [26, 27] as the corresponding approximate amplitudes. Let us
observe that the initial approximate limit cycle 
is justified
because it recovers the exact circular form when 
, and also, when 
, it becomes exactly the two-piecewise limit cycle
proposed and plotted in [26]. Moreover, it must be noticed that 
makes sense
only when 
, and then the results obtained from this initial
guess only will have validity for 
.
After eliminating the time variable, the nonlinear
differential equation 
is of order
one. Then, to construct the homotopy [8, 9], we build a linear differential equation of order one
for the whole range 
of what is
called the homotopy parameter 
.
We define an auxiliary linear operator 
by
(14)with the
property
(15)where 
is a constant
function (the kernel of operator 
), and 
is a (homotopy)
parameter explained below.
From (12) we define a nonlinear
operator
(16)and then construct the
homotopy
(17)where 
is a nonzero
auxiliary parameter. Setting 
, we have the zero-order deformation
equation
(18)subject to the boundary
conditions
(19)where 
is an embedding
parameter. When the parameter 
increases from
0 to 1, the solution 
varies from 
to 
, and 
varies from 
to 
. Assume that 
and 
are analytic in 
and can be
expanded in the Maclaurin series of 
as
follows:
(20)where
(21)Notice that series (20)
contain the auxiliary parameter 
, which has influence on their convergence regions.
Assume that 
is properly
chosen such that all of these Maclaurin series are convergent at 
. Hence at 
we
have
(22)At the 
th-order
approximation, we have the analytic solution of (12),
namely,
(23) The parameter 
is free and can
be chosen arbitrarily; in particular, it can be a function of 
. Nevertheless, the function 
is the exact
limit cycle solution of the Liénard equation (12) in the limits 
and 
. This means that the general solution 
of (18)
should tend to 
for any 
in those limits
of 
, loosing its dependence on 
. Then, a reasonable property for 
would be to
vanish in those limits of 
. Hence, the solution of (18) would be the exact
solution 
for any value
of the parameter 
in the limits 
and 
. A simple function 
satisfying
these properties is
(24) Differentiating (18) and (19) 
times with
respect to 
, then setting 
, and finally dividing by 
, we obtain the 
th-order
deformation equation
(25)subject to the boundary
conditions
(26)where 
is defined
by
(27)At each iteration, we have two
unknowns: 
, in (15), and 
. These unknowns are obtained by considering the
boundary conditions (26) as follows.
At zero order, we obtain
(28)At 1st-order, we
obtain
(29)where we have imposed the
condition 
choosing the
integration constant, 
, appropriately. At this moment 
is free but we
can fix the value of 
by imposing 
:
(30)This is a nonlinear equation
for 
. When 
is an even
polynomial of degree 
, it is a polynomial equation of degree 
in the variable 
. Therefore, at this order, the number of limit cycles
is the number of positive roots of this equation, at most 
(see [27] for a longer discussion of
this point).
Then, for every one of the 
solutions of
the above equation, that is, for every one of the limit cycles of the system, we
proceed order by order in 
to obtain
higher order corrections to the amplitude 
and to the
shape of the limit cycle 
.
For example, at 2th-order, we obtain
(31)
(32) Therefore, a first approximation to the shape of the
limit cycle is 
, and a first-order approximation to the amplitude of
the limit cycle is 
.
At 3th-order, we have
(33)
(34)Therefore, a second-order
approximation to the shape of the limit cycle is 
, and a second-order approximation to the amplitude of the limit cycle is 
.
Moreover, at every order in 
the appropriate
function 
is indicated by
the experiment, and in some sense, 
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 have 
and
(35) The expressions for 
, and 
can be
calculated as explained in Section 3.1. From (30) we get
(36) From (32) we obtain
(37) Therefore, a first approximation to the amplitude of
the limit cycle is
(38)As it was advanced at the
beginning of this section, note that this expression is only valid for positive 
. Choosing 
and 
in
(24),
(39)the formula 
given in (5)
is obtained. The values for 
and 
are selected by
guess and check.
From (34) we get
(40) Therefore, a second approximation to the amplitude of
the limit cycle is
(41)Choosing 
and 
in
(24),
(42)the formula 
shown in (6)
is finally obtained.
As an example, and to conclude this section, we plot
in Figure 3 the form of the van der Pol limit cycle when 
, and 
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
and
given in
(
42). Black line: the experimental limit cycle. Brown, blue, and red lines:
the curves
and
, 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 
of the van der
Pol limit cycle. In addition to the well-known constant approximation 
that generates
an error less than 
, we establish a family of recursive formulas that are
valid for the whole range of the parameter 
. Two of them, 
and 
, have been explicitly given. The first one, 
, produces an error less than 
, and the second
one, 
, reduces the error to less than 
. Moreover, 
is conjectured
to be an upper bound of 
.
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.
Acknowledgment
The Gobierno
of Navarra, Spain, Res. 07/05/2008 is acknowledged by its financial
support.
- B. van der Pol, “On relaxation oscillations,” Philosophical Magazine, vol. 2, pp. 978–992, 1926.
- A. A. Andronov, A. A. Vitt, and S. E. Khaikin, Theory of Oscillators, Dover, New York, NY, USA, 1st edition, 1989.
- B. van der Pol and M. van der Mark, “The beating of the heart considered as relaxation oscillation and an electric model of the heart,” L'Onde électrique, vol. 7, pp. 365–392, 1928.
- R. López-Ruiz and Y. Pomeau, “Transition between two oscillation modes,” Physical Review E, vol. 55, no. 4, pp. R3820–R3823, 1997.
- K. Odani, “The limit cycle of the van der Pol equation is not algebraic,” Journal of Differential Equations, vol. 115, no. 1, pp. 146–152, 1995.
- A. W. Jordan and P. Smith, Nonlinear Ordinary Differential Equations, Oxford University Press, Oxford, UK, 1st edition, 1977.
- Y. Ye, Theory of Limit Cycles, vol. 66 of Translations of Mathematical Monographs, American Mathematical Society, Boston, Mass, USA, 1986.
- S. J. Liao, The proposed homotopy analysis technique for the solutions of non-linear problems, Ph.D. thesis, Shanghai Jiao Tong University, Shanghai, China, 1992.
- S. J. Liao, Beyond Perturbation: Introduction to the Homotopy Analysis Method, Chapman & Hall/CRC Press, Boca Raton, Fla, USA, 2003.
- S.-J. Liao, “An analytic approximate approach for free oscillations of self-excited systems,” International Journal of Non-Linear Mechanics, vol. 39, no. 2, pp. 271–280, 2004.
- W. Wu and S.-J. Liao, “Solving solitary waves with discontinuity by means of the homotopy analysis method,” Chaos, Solitons & Fractals, vol. 26, no. 1, pp. 177–185, 2005.
- T. Hayat and M. Khan, “Homotopy solutions for a generalized second-grade fluid past a porous plate,” Nonlinear Dynamics, vol. 42, no. 4, pp. 395–405, 2005.
- T. Hayat, M. Khan, and M. Ayub, “On non-linear flows with slip boundary condition,” Zeitschrift für Angewandte Mathematik und Physik, vol. 56, no. 6, pp. 1012–1029, 2005.
- T. Hayat, M. Khan, M. Sajid, and M. Ayub, “Steady flow of an Oldroyd 8-constant fluid between coaxial cylinders in a porous medium,” Journal of Porous Media, vol. 9, no. 8, pp. 709–722, 2006.
- T. Hayat, M. Khan, and M. Ayub, “The effect of the slip condition on flows of an Oldroyd 6-constant fluid,” Journal of Computational and Applied Mathematics, vol. 202, no. 2, pp. 402–413, 2007.
- M. Khan, Z. Abbas, and T. Hayat, “Analytic solution for flow of Sisko fluid through a porous medium,” Transport in Porous Media, vol. 71, no. 1, pp. 23–37, 2008.
- S. Abbasbandy, “The application of homotopy analysis method to solve a generalized Hirota-Satsuma coupled KdV equation,” Physics Letters A, vol. 361, no. 6, pp. 478–483, 2007.
- S. Abbasbandy, “Homotopy analysis method for heat radiation equations,” International Communications in Heat and Mass Transfer, vol. 34, no. 3, pp. 380–387, 2007.
- S. Abbasbandy, Y. Tan, and S. J. Liao, “Newton-homotopy analysis method for nonlinear equations,” Applied Mathematics and Computation, vol. 188, no. 2, pp. 1794–1800, 2007.
- Y. Tan and S. Abbasbandy, “Homotopy analysis method for quadratic Riccati differential equation,” Communications in Nonlinear Science and Numerical Simulation, vol. 13, no. 3, pp. 539–546, 2008.
- S. Abbasbandy, J. L. López, and R. López-Ruiz, “The homotopy analysis method and the Liénard equation,” http://arxiv.org/abs/0805.3916.
- K. Odani, “On the limit cycle of the Liénard equation,” Archivum Mathematicum, vol. 36, no. 1, pp. 25–31, 2000.
- J. L. López and R. López-Ruiz, “Approximating the amplitude and form of limit cycles in the weakly nonlinear regime of Liénard systems,” Chaos, Solitons & Fractals, vol. 34, no. 4, pp. 1307–1317, 2007.
- J. L. López and R. López-Ruiz, “The limit cycles of Liénard equations in the weakly nonlinear regime,” to appear in Far East Journal of Dynamical Systems (2009); also available at http://arxiv.org/abs/nlin/0605025.
- M. L. Cartwright, “Van der Pol's equation for relaxation oscillation,” in Contributions to the Theory of Non-Linear Oscillations II, S. Lefschetz, Ed., vol. 29 of Annals of Mathematics Studies, pp. 3–18, Princeton University Press, Princeton, NJ, USA, 1952.
- J. L. López and R. López-Ruiz, “The limit cycles of Liénard equations in the strongly nonlinear regime,” Chaos, Solitons & Fractals, vol. 11, no. 5, pp. 747–756, 2000.
- R. López-Ruiz and J. L. López, “Bifurcation curves of limit cycles in some Liénard systems,” International Journal of Bifurcation and Chaos, vol. 10, no. 5, pp. 971–980, 2000.
