One Hat Cyber Team

View File Name : simsek.tex

\documentclass[12pt]{article}
\usepackage{times}
\newcommand{\dis}{\displaystyle}
\newcommand{\lsim}{\stackrel{<}{_\sim}}
\newcommand{\gsim}{\stackrel{>}{_\sim}}
\usepackage{epsf}
\setlength{\topmargin}{-1.5cm} \setlength{\textheight}{23.5cm}
\setlength{\oddsidemargin}{0.cm} \setlength{\textwidth}{17.cm}
\def\beq{\begin{equation}}
\def\eeq{\end{equation}}
\def\bea{\begin{eqnarray}}
\def\eea{\end{eqnarray}}
\def\ve{\vert}
\def\ve{\vert}
\def\ver{\right|}
\def\vel{\left|}
\def\ver{\right|}
\def\nnb{\nonumber}
\def\ga{\left(}
\def\dr{\right)}
\def\aga{\left\{}
\def\adr{\right\}}
\def\rar{\rightarrow}
\def\nnb{\nonumber}
\def\lla{\langle}
\def\rra{\rangle}
\def\ba{\begin{array}}
\def\ea{\end{array}}
\def\bea{\begin{eqnarray}}
\def\eea{\end{eqnarray}}
\def\Bgll{$B_{d}\rightarrow \,\gamma\, \ell^+ \ell^-$}
\def\tep{$b \rar d \ell^+ \ell^-$}
\def\ds{\displaystyle}
\title{Stability Boundaries of Disturbed Van der Pol Oscillator
}
\author{\vspace{1cm}\\
         {\bf \"{O}zlem Ye\c{s}ilta\c{s}}
         \thanks{E-mail address:
        yesiltas@gazi.edu.tr} \, \, and \, \,
        {\bf Mehmet \c{S}im\c{s}ek }
        \thanks{E-mail address:
        msimsek @gazi.edu.tr}\,\,
     %}
\\
\\
{\small \sl ${}$Gazi University, Faculty of Arts and Science,
Department of Physics, 06500, }
\\  {\small \sl  $\ $Teknikokullar, Ankara, Turkey}}



\begin{document}
%\setlength{\baselineskip}{24pt}
\maketitle
%\setlength{\baselineskip}{7mm}

\begin{abstract}
In this study, the approximate solutions are presented for
non-linear multiple-degree of freedom  vibrating systems. As an
example disturbed Van der Pol oscillator with periodic
coefficients is solved by using multiple time scale analysis
without secular terms grow and periodic stability conditions are
obtained. Then, using Von Zeipel Method stability considerations
are presented. It is shown that the transition curves are same
with those obtained by both methods.
\end{abstract}

PACS No: 03.65.-w,63.50.+x,46.30.-i,31.15.+q,02.90.+p

\newpage
%
\section{Introduction}
Although there has been a considerable amount of research study on
non-linear vibrating systems, in general, exact analytic solutions
of non-linear differential equations are possible only a limited
classes of systems. Then, we have to resort to approximate method
to solve many of the real systems. Multiple-Scale Perturbation
Theory (MSPT) is one of the effective techniques of the
approximate method that can be applied to the many problems in
physics and natural sciences[1]. Reformulating the perturbation
series to get through secular growth, MSPT is an extremely useful
method for solving such perturbation physical problems with a
small parameter $\varepsilon $. MSPT is applicable to both linear
and non-linear differential equations through the classical and
quantum mechanical problems [2,3]. The technique has been
developed by Sturrock ,Frieman and Nayfeh [4,5,6]. The main goal
underlying MSPT is that dynamical systems have distinct
characteristic physical behavior at different length or slow/ fast
time scales[7]. Eliminating  the secular terms in the fast time
variable,$T_{0}=t$,it leads to the well-known solvability
conditions. For two scale problems, the solutions oscillate on a
time scale of order $t$ with an amplitude and phase which drifts
on a time scale of order $\epsilon t$. Starting from the first
order, free terms which are solutions of the unperturbed equations
emerge in every order in the expansion of the approximate
solution. If the constraints are not obeyed, then the analysis may
lead to wrong results, or allow only trivial solutions.

The study of non-linear oscillators has been important in the
development of the theory of dynamical systems. Van der Pol and
Van der Mark [8] studied a simple non-linear electronic circuit(a
neon tube was the non-linear element) and experimentally found,
but were not much interested in, noisy behavior that can be
identified as chaos. Recently, there are many studies about Van
der Pol equation [9,10,11]. The mechanism for the action of time
delay in a non-autonomous system such as VdP-Duffing oscillator
with excitation is investigated by Xu and Chung [12]. Dynamics of
relaxation oscillations of VdP equation is investigated using an
analytical method requiring only the connection at a point of the
two dynamical fast and slow regions [13]. We study the form of
disturbed motion,

$  $

$
%\begin{eqnarray}
\ddot{y}(t)-\varepsilon \left(
1-x^{2}(t)\right)\dot{y}(t)+(1+2\varepsilon x(t)\dot{x}(t) )
y(t)=0.
%\end{eqnarray}
$

$ $


 Straightforward application of perturbation theory to non-linear
equations of motion in classical mechanics gives rise to secular
terms that increase unboundedly with time even for periodic motion
and unphysical terms also appear in the application of time
dependent perturbation theory to quantum mechanical systems[14].

This paper is organized as follows:

In section 2, as discussing first variational equation of VdP and
eliminating the first order derivative in equations(4-5),
disturbed VdP oscillator, which is a second order non-linear
differential equation with periodic coefficients is underlined.
Disturbed VdP equation is solved in section 3, without secular
terms grow is solved to first order by using MSPT as one of the
effective singular perturbation methods. Transition curves and
steady states are obtained. In addition, for more accurate results
are obtained in order to find stability boundary to second order
and steady states. Using Von Zeipel's method, new
momenta-coordinates for reduced VdP equation are obtained in
section 4, and then, it is shown that the stability boundaries are
same with those obtained both method. The study end
conclusions(Section 5).


\section{Disturbed Van der Pol Oscillator}
In this section, before introducing MSPT, Van der Pol (VdP)
equation is discussed in Saaty's way [15]. Consider well-known VdP
equation, which is studied for many years:

\begin{eqnarray}
\ddot{x}(t)-\epsilon \left( 1-x^{2}(t)\right)\dot{x}(t)+x(t)=0
\end{eqnarray}

This equation has  one limit cycle for $~\varepsilon >~0$ [15].
One can find singular points of VdP (1) with

\begin{eqnarray}
\frac{dy}{dx}=-\frac{x-\varepsilon(1-x^{2})y}{y}.
\end{eqnarray}

If~$|x|<<1$,~ approximately $\frac{dy}{dx}\simeq
-\frac{x-\varepsilon y}{y}$ and, hence the origin~ $x=y=0$~is a
singular and unstable spiral point.
If~$|x|>>1$,$\frac{dy}{dx}\simeq -\frac{x-\varepsilon x^2 y}{y}$
 and the origin is also a singular point. If $%
\varepsilon =0$, ~$y=A~\cos t$~ is a periodic solution with
arbitrary A. After eliminating secular terms, one can find
approximate solution up to first order terms in ~$\varepsilon $ as
[15]

\begin{eqnarray}
 y(t) \simeq
2cost+\frac{\varepsilon}{4}(3sint-sin3t)+O(\varepsilon^2)
\end{eqnarray}

In order to test the stability \ of the approximate periodic
solution $a\sin t$ ,\  applying the variational method to VdP
equation with $\delta x=y$ , then (1) turns into the form of the
disturbed motion [15],

\begin{eqnarray}
\ddot{y}-\varepsilon \left( 1-x^{2}\right)\dot{y}+(1+2\varepsilon
x\dot{x} ) y=0
\end{eqnarray}

Write $y=e^{v(t)}u(t)$ in (4) and then put
$v(t)=-\frac{\varepsilon }{2}\left[ \left(
\frac{a^{2}}{2}-1\right) t-\frac{a^{2}}{4}\sin 2t\right] $, the
term $\dot{y} $ disappears [15],

\begin{eqnarray}
\ddot{u}(t)+\left( 1+\frac{1}{2}\varepsilon a^{2}\sin (2t)-\frac{3}{4}%
\varepsilon ^{2}\left( \frac{1-a^{2}\sin ^{2}t}{2}\right)
^{2}\right) u(t)=0
\end{eqnarray}

If the terms  in $\varepsilon ^{2}\,\,$ are neglected and $t$ is
transformed into  $t=t+\frac{\pi }{4}\,,$ (5) reduces to the
well-known Mathieu equation:

\begin{eqnarray}
\ddot{u}(t)+\left( 1+\frac{\varepsilon a^{2}}{2}\cos 2t\right)
u(t)=0
\end{eqnarray}

The Mathieu equation is an example of a  differential equation
with periodic coefficients[1,15,16,17].


\section{Multiple Scale analysis of Reducable Van der Pol Equation  }

It is possible to solve (5), disturbed VdP oscillator equation, by
using MSPT, so we seek the perturbative solution to (5) having two
variables the short-time scale $t$ and the long-time scale $\tau
=\varepsilon \,t$ when $\varepsilon $ is sufficiently small. We
would like to investigate the solutions of (5) by disturbing this
equation. Because to terminate \ more\ secular terms and unstable
solutions near $\frac{1}{4}$[1] , time variable is changed as $\
t\rightarrow \frac{t}{2}+\frac{\pi }{4}$ \ and the (5) turns into

\begin{eqnarray}
\ddot{u}(t)+\left(\frac{1}{4}+\frac{1}{8}\varepsilon a^{2}\cos
t-\frac{3}{16}\varepsilon ^{2}\left( 1-\frac{a^{2}}{2}(1+\sin
t)\right) ^{2}\right)u(t)=0.
\end{eqnarray}

Applying the multiple-scale perturbation theory we find the
boundaries between the regions in the $\left( \Delta ,\varepsilon
\right) $
plane for which all solutions to general equation are stable. So we put $%
\Delta $ parameter\ in (7)

\begin{eqnarray}
\ddot{u}(t)+\left( \frac{\Delta}{4} +\frac{1}{8}\epsilon a^{2}\cos
t- \frac{3}{16}\varepsilon ^{2}\left( 1-\frac{a^{2}}{2}(1+\sin
t)\right) ^{2}\right) u(t)=0,
\end{eqnarray}
and expand it\ as a power series in [16] $\epsilon $. By writing \
$\Delta
=1+\epsilon \delta _{1}+\epsilon ^{2}\delta _{2}+...$  determine the $%
\Delta $ so that the resulting expansion is periodic.
Consequently, the resulting expressions for $\Delta $ define the
transition curves separating stability from instability

\begin{eqnarray}
\ddot{u}(t)+\left( \frac{1}{4}+\epsilon (\frac{\delta_{1}}{4}+\frac{1}{8}%
a^{2}\cos t)+\epsilon ^{2}(\frac{\delta_{2}}{4}-\frac{3}{16}\left( 1-\frac{a^{2}}{2}%
(1+\sin t)\right) ^{2})\right) u(t)=0 ,
\end{eqnarray}

is the general equation. An uniform expansion is in the form :

\begin{eqnarray}
u(t)=U_{0}(t,\tau )+\varepsilon U_{1}(t,\tau )+\varepsilon
^{2}U_{2}(t,\tau )+O(\varepsilon ^{3})
\end{eqnarray}

The chain rule for the partial differentiation to compute
derivatives of $u(t)$ is, with using $\frac{d\tau
}{dt}=\varepsilon $ $\quad $

\begin{eqnarray}
\ddot{u}(t)&=&\frac{\partial ^{2}U_{0}(t,\tau )}{\partial t^{2}%
}+\varepsilon \,\left( 2\frac{\partial ^{2}U_{0}(t,\tau
)}{\partial \tau\partial t}+\frac{\partial ^{2}U_{1}(t,\tau
)}{\partial t^{2}}\right) \nnb \\
&&+\qquad \qquad \qquad \varepsilon ^{2}\left( \frac{\partial ^{2}U_{2}}{%
\partial t^{2}}+2\frac{\partial ^{2}U_{1}}{\partial \tau \partial t}+\frac{%
\partial ^{2}U_{0}}{\partial \tau ^{2}}\right) +O(\epsilon ^{3}).
\end{eqnarray}

Substituting (10) and (11) in (9) , partial differential equations
are obtained for the dependent variables $U_{0}(t,\tau ),
U_{1}(t,\tau ),..$ So,in this solution there are less secular
terms to all order in the perturbation series. The first three
terms of the series are given by

\begin{eqnarray}
\frac{\partial ^{2}U_{0}(t,\tau )}{\partial
t^{2}}+\frac{1}{4}U_{0}(t,\tau)=0,
\end{eqnarray}

\begin{eqnarray}
\frac{\partial ^{2}U_{1}(t,\tau )}{\partial t^{2}%
}+\frac{1}{4}U_{1}(t,\tau )=-2\frac{\partial ^{2}U_{0}(t,\tau
)}{\partial \tau \partial
t}-(\frac{\delta_{1}}{4}+\frac{1}{8}a^{2}\cos t)U_{0}(t,\tau ),
\end{eqnarray}

\begin{eqnarray}
\frac{\partial ^{2}U_{2}(t,\tau )}{\partial t^{2}
}+\frac{1}{4}U_{2}(t,\tau )=&-2\frac{\partial ^{2}U_{1}(t,\tau
)}{\partial \tau \partial
t}-(\frac{\delta_{1}}{4}+\frac{1}{8}a^{2}\cos t)U_{1}(t,\tau
)-\frac{\partial ^{2}U_{0}}{
\partial \tau ^{2}}  \nnb \\
&-\left( \frac{\delta_{2}}{4}-\frac{3}{16}\left(
1-\frac{a^{2}}{2}(1+\sin t)\right) ^{2}\right) U_{0}(t,\tau )
\end{eqnarray}

The form of general solution to (12) can be chosen as

\begin{eqnarray}
U_{0}(t,\tau )=A_{0}(\tau )e^{it/2}+A_{0}^{\star }(\tau )e^{-it/2}
\end{eqnarray}

In order to determine $~\tau ~~(\tau =\varepsilon t)$ dependent
coefficients$~A_{0}(\tau )$,~ $A_{0}^{\star }(\tau )$  substitute $%
U_{0}(t,\tau )\;$into the right hand side of (13). One can see
that $e^{it/2}$ and$\ \ e^{-it/2}$ are the solutions of the left
-hand side equation i.e.
homogenous equation $\frac{\partial ^{2}U_{1}(t,\tau )}{\partial t^{2}}+%
\frac{1}{4}U_{1}(t,\tau )=0$ . Therefore, the right-hand side
contains terms that produce secular terms in $U_{1}$. For a
uniform expansion, these terms must be eliminated, this is
proceeding by setting coefficients of \ $e^{it/2}$ and$\ \
e^{-it/2}$ equal to zero. So, \ $A_{0}(\tau )$ satisfy:

\begin{eqnarray}
-i \frac{\partial A_{0}}{\partial \tau }-
\frac{a^2}{16}A_{0}^{\star }-\frac{\delta_{1}}{4} A_{0} =0,
\end{eqnarray}

\begin{eqnarray}
i \frac{\partial A_{0}^{\star }}{\partial \tau }-
\frac{a^2}{16}A_{0}-\frac{\delta_{1}}{4} A_{0}^{\star } =0.
\end{eqnarray}

Setting $A_{0}(\tau )=w_{0}(\tau )+iv_{0}(\tau )$ we find:

\begin{eqnarray}
-w_{0}^{^{\prime }}(\tau )-\frac{\delta_{1}}{4}v_{0}(\tau
)+\frac{a^{2}}{16} v_{0}(\tau )=0,
\end{eqnarray}

\begin{eqnarray}
v_{0}^{^{\prime }}(\tau )-\frac{\delta_{1}}{4}w_{0}(\tau
)-\frac{a^{2}}{16}w_{0}(\tau )=0,
\end{eqnarray}

from (18) and (19)

\begin{eqnarray}
v_{0}^{^{\prime \prime }}(\tau
)+\frac{1}{4}(\delta^{2}_{1}-\frac{a^{4}}{16})v_{0}=0,
\end{eqnarray}

and the solutions of \ $w_{0}(\tau )$,$v_{0}(\tau )$ are:

\begin{eqnarray}
v_{0}(\tau )=C_{1}\cos \lambda \tau +C_{2}\sin \lambda \tau,
\end{eqnarray}

\begin{eqnarray}
w_{0}(\tau )=\frac{\lambda }{\delta _{1}+\frac{a^{2}}{4}}(
-C_{1}\sin \lambda \tau +C_{2}\cos \lambda \tau )
\end{eqnarray}
where
\begin{eqnarray}
\lambda=\frac{1}{2}\sqrt{\delta _{1}^{2}-\frac{a^{4}}{16}}.
\end{eqnarray}


 Here,instability occurs if $\ \delta _{1}^{2}-\frac{a^{4}}{16}$ is
negative. Thus, $\left\vert \delta _{1}\right\vert
>\frac{a^{2}}{4}$ gives stable solutions and $\left\vert \delta
_{1}\right\vert <\frac{a^{2}}{4}$ gives unstable solutions.Near
$\epsilon =0$ , the $\Delta $ is given as:

\begin{eqnarray}
\Delta =1\pm \frac{a^{2}}{4}\epsilon +O(\epsilon ^{2}),\qquad
\epsilon \longrightarrow 0.
\end{eqnarray}

If the initial conditions are specified as $U(0,0)=1$ and $\dot{U}%
(0,0)=0$ and $U_{0}(0,0)=1,\dot{U}_{0}(0,0)=0$ , then, $C_{1}=0,C_{2}=%
\frac{\delta _{1}+\frac{a^{2}}{4}}{2\lambda }$ and $U_{0}(t,\tau
):$

\begin{eqnarray}
U_{0}(t,\tau )=\left( \cos \lambda \tau \cos
\frac{t}{2}-\frac{\delta _{1}+ \frac{a^{2}}{4}}{\lambda }\sin
\lambda \tau \sin \frac{t}{2}\right)
\end{eqnarray}

Equation (25) is a solution of Mathieu equation also and it is
obvious that with choosing appropriate parameters in this
equation,the result is same obtained by Bender [1]\ .The
differential equation in the present case is the simple harmonic
oscillator with an external periodic force which is its amplitude
have damping factor for~$\varepsilon >0$. Furthermore, one can
find the periodic solutions of partial differential (13) as,
\begin{eqnarray}
U_{1}(t,\tau ) &=&A_{1}(\tau )e^{\frac{it}{2}}+A_{1}^{\ast }(\tau
)e^{-\frac{it}{2}} \nnb \\
 &&-\frac{a^{2}}{512\lambda }\left[ \left(
a^{2}+4\delta _{1}\right) \sin \lambda \tau \sin
\frac{3t}{2}+4\lambda \cos \lambda \tau \cos \frac{3t}{2} \right]
\end{eqnarray}
where$~A_{1}(\tau )~$ and $~A_{1}^{\ast }(\tau )$ are unknown
integral coefficients yet. To find these coefficients,used (14)is
used. However, inserting $U_{0}(t,\tau )$ , and~ $U_{1}(t,\tau
)$~into (14), an another resonance case as would occur in this
equation. At this step, the coefficients of ~$cos\frac{t}{2}~$ and
~$sin\frac{t}{2}~$ be zero. It is usually interested in the
solution of a differential equation if the initial conditions or
right side of the equation is changed. The existence theorem
guarantees a unique solution of the differential equation for each
choice of initial conditions. Substitute the $U_{1}(t,\tau )$ and
$U_{0}(t,\tau )$ and trigonometric formulas into the right
hand-side of (14) and $A_{1}=w_{1}+iv_{1}$,

\begin{eqnarray}
 -16384 \frac{\partial v_{1}}{\partial \tau
}-\eta_{1} \sin \lambda \tau-\eta_{2}\cos\lambda \tau=0,
\end{eqnarray}

\begin{eqnarray}
1024\lambda \left( 4\delta _{1}+a^{2}\right) w_{1}+\eta _{3}\sin
\lambda \tau +\eta _{4}\cos \lambda \tau =0.
\end{eqnarray}

Here $A_{1}(\tau )=-A_{1}^{\ast }(\tau )$ .$\ \eta _{1},\eta _{2}$
are constants and given as respectively:

\begin{eqnarray}
\eta _{1} &=&\left( \lambda -36\right) a^{6}+\left( 96+16\lambda
\delta _{1}-576\delta _{1}\right) a^{4}-\left( 6656-8192\lambda
^{2}\right) \delta
_{1}  \nnb \\
&&+\left( 1536\delta _{1}+512\delta _{2}-512\lambda ^{2}-96\right)
a^{2},
\end{eqnarray}

\begin{eqnarray}
\eta _{2}=384a^{2}\lambda \left( 1-2a^{2}\right),
\end{eqnarray}

\begin{eqnarray}
\eta _{3}=24a^{2}\left( a^{4}-2a^{2}+16a^{2}\delta _{1}-32\delta
_{1}\right),
\end{eqnarray}

\begin{eqnarray}
\eta _{4}=64\lambda \left( -24-128\lambda ^{2}+128\delta
_{2}+24a^{2}-9a^{4}\right).
\end{eqnarray}

\begin{eqnarray}
A_{1}(\tau ) &=&iC_{1}+\sin \lambda \tau \left( -\frac{\eta _{3}}{
1024\lambda \left( a^{2}+4\delta _{1}\right) }-\frac{i\eta _{2}}{
16384\lambda ^{2}}\right) +  \nnb \\
&&\cos \lambda \tau \left( -\frac{\eta 4}{1024\lambda \left(
a^{2}+4\delta _{1}\right) }+\frac{i\eta _{1}}{16384\lambda
^{2}}\right)
\end{eqnarray}

\begin{eqnarray}
A_{1}^{\ast }(\tau ) &=&-iC_{1}+\sin \lambda \tau \left(
-\frac{\eta _{3}}{ 1024\lambda \left( a^{2}+4\delta _{1}\right)
}+\frac{i\eta _{2}}{
16384\lambda ^{2}}\right)   \nnb \\
&&+\cos \lambda \tau \left( -\frac{\eta 4}{1024\lambda \left(
a^{2}+4\delta _{1}\right) }-\frac{i\eta _{1}}{16384\lambda
^{2}}\right)
\end{eqnarray}

Here is the $C_{1}$ is an arbitrary constant and we say $C_{1}=0$
.$U_{1}(t,\tau )$ can be written again:

\begin{eqnarray}
U_{1}(t,\tau ) &=&\alpha (\eta _{3}\sin \lambda \tau \cos
\frac{t}{2}+\eta _{4}\cos \lambda \tau \cos \frac{t}{2})+\beta
(\eta _{2}\sin \lambda \tau
\sin \frac{t}{2}-\eta _{1}\cos \lambda \tau \sin \frac{t}{2})  \nnb \\
&&-\frac{a^{2}}{512\lambda }\left[ \left( a^{2}+4\delta
_{1}\right) \sin \lambda \tau \sin \frac{3t}{2}+4\lambda \cos
\lambda \tau \cos \frac{3t}{2} \right]
\end{eqnarray}
where $\alpha =-\frac{1}{512\lambda \left( a^{2}+4\delta _{1}\right) }%
, ~~~ \beta =\frac{1}{8192\lambda ^{2}}.$\ Thus, we illustrated up
to first-order calculation. However, using similar algorithm, it
is easy to calculate in a power series in~$\varepsilon $~ up to a
high order and its coefficients accurately determined.

%\section{Higher Order Corrections To The Stability Boundary}

In addition  to calculations  it is needed to examine the higher
order terms to determine the location of the stability boundary
more precisely than (24). Moreover one can write the $U_{0}$ and
$U_{1}$ states more accurate and simple than in (25) and (34). To
do this, turn back to (20) :

\begin{eqnarray}
v_{0}^{^{\prime \prime }}(\tau )+\frac{1}{4}\left( \delta
_{1}^{2}-\frac{a^{4}}{16} \right) v_{0}(\tau )=0,
\end{eqnarray}

 Write $v_{0}(\tau )$ again but this time in this form:

\begin{eqnarray}
v_{0}(\tau )=A\exp \left( \pm i\sqrt{\delta
_{1}^{2}-\frac{a^{4}}{16}}\tau \right)
\end{eqnarray}

where $A$ is a constant. Thus , a new time scale for the problem
occurs: assume that $\delta _{1}=\frac{a^{2}}{4}+\delta
_{2}\epsilon $ in (37). Then $v_{0}(\tau )$ becomes approximately $%
K\exp \left( \pm i\sqrt{\frac{a^{2}}{2}\epsilon \delta _{2}}\tau
\right)
=K\exp \left( \pm i\sqrt{\frac{a^{2}}{2}\delta _{2}}\epsilon ^{3/2}t\right) $%
, which advances that a new time scale $\sigma =\epsilon ^{3/2}t$
must be introduced. Therefore substitute $\Delta$ into (8), and it
is noted that $\sigma =\epsilon ^{3/2}t $

\begin{eqnarray}
\ddot{U}(t)+\left( \frac{1}{4}+\epsilon
(\frac{a^{2}}{4}+\frac{1}{8} a^{2}\cos (t))+\epsilon ^{2}(2\delta
_{2}-\frac{3}{16}\left( 1-\frac{a^{2}}{2 }(1+\sin t)\right)
^{2})\right) U(t)=0.
\end{eqnarray}

It means expanding\ new $u(t)$ series and second order
differential operator as:

\begin{eqnarray}
U(t)\simeq U_{0}(t,\sigma )+\epsilon ^{\frac{1}{2}}U_{1}(t,\sigma
)+\epsilon U_{2}(t,\sigma )+\epsilon ^{\frac{3}{2}}U_{3}(t,\sigma
)+\epsilon ^{2}U_{4}(t,\sigma )+...
\end{eqnarray}

\begin{eqnarray}
\ddot{U}(t)=&\frac{\partial ^{2}U_{0}}{\partial t^{2}}+\epsilon
^{\frac{1}{2}}\frac{\partial ^{2}U_{1}}{\partial t^{2}}+\epsilon
\frac{
\partial ^{2}U_{2}}{\partial t^{2}}+\epsilon ^{\frac{3}{2}}(\frac{\partial ^{2}U_{3}}{\partial
t^{2}}+2\frac{
\partial ^{2}U_{0}}{\partial t\partial \sigma })+\epsilon ^{2}(\frac{
\partial ^{2}U_{4}}{\partial t^{2}}+2\frac{\partial ^{2}U_{1}}{\partial t\partial \sigma })+...
\end{eqnarray}

Putting these equations into (38) and equate powers of $%
\epsilon ^{\frac{1}{2}}$:

\begin{eqnarray}
\epsilon ^{0}~~~:~~~\frac{\partial ^{2}U_{0}}{\partial t
^{2}}+\frac{U_{0}}{4}=0,
\end{eqnarray}

\begin{eqnarray}
\epsilon ^{\frac{1}{2}}~~~:~~~\frac{\partial ^{2}U_{1}}{\partial
t^{2}}+\frac{U_{1}}{4}=0,
\end{eqnarray}

\begin{eqnarray}
\epsilon ^{1}~~~:~~~\frac{\partial ^{2}U_{2}}{\partial t^{2}}+\frac{U_{2}}{4}=-%
\frac{a^{2}}{4}\left( 1+e^{it}+e^{-it}\right) U_{0},
\end{eqnarray}

\begin{eqnarray}
\epsilon ^{\frac{3}{2}}~~~:~~~\frac{\partial ^{2}U_{3}}{\partial
t^{2}}+ \frac{U_{3}}{4}=-2\frac{\partial ^{2}U_{0}}{\partial
t\partial \sigma }- \frac{a^{2}}{4}\left( 1+e^{it}+e^{-it}\right)
U_{1},
\end{eqnarray}

\begin{eqnarray}
\epsilon ^{2} &~~~:~~~&\frac{\partial ^{2}U_{4}}{\partial t^{2}}+\frac{U_{4}}{4}=-2%
\frac{\partial ^{2}U_{1}}{\partial t\partial \sigma }-\frac{a^{2}}{4}%
\left( 1+e^{it}+e^{-it}\right) U_{2}  \nnb \\
&&-(2\delta{2}-\frac{3}{16}\left( 1-\frac{a^{2}}{2}(1+\sin
t)\right) ^{2})U_{0}.
\end{eqnarray}

The solutions of (41) and (42) are respectively:

\begin{eqnarray}
U_{0}(t,\sigma )=A_{0}(\sigma )e^{\frac{it}{2} }+A_{0}^{\ast
}(\sigma )\!e^{\frac{-it}{2}},
\end{eqnarray}

\begin{eqnarray}
U_{1}(t,\sigma )=A_{1}(\sigma )e^{\frac{it}{2} }+A_{1}^{\ast
}(\sigma )\!e^{\frac{-it}{2}}.
\end{eqnarray}

As following the same procedure, from the right hand side of (43)
,secularity
removes and one can find $A_{0}(\sigma )=-A_{0}^{\ast }(\sigma )$ and say $%
A_{0}(\sigma )=iB_{0}(\sigma ).$ Now  $U_{2}$ can be solved as:

\begin{eqnarray}
U_{2}(t,\sigma )=A_{2}(\sigma )e^{\frac{it}{2}%
}+A_{2}^{\ast }(\sigma
)\!e^{\frac{-it}{2}}+\frac{a^{2}}{8}A_{0}e^{\frac{3it}{2}}-\frac{a^{2}}{8}A_{0}e^{\frac{-3it}{2}}
\end{eqnarray}

Vanishing the coefficients of $\ e^{\pm \frac{it}{2}}~$in (44)\
gives:

\begin{eqnarray}
 \frac{\partial B_{0}}{\partial \sigma
}= \frac{a^2}{4} \left( A_{1}(\sigma )+A_{1}^{\ast }(\sigma
)\right),
\end{eqnarray}

From (42-46)

\begin{eqnarray}
 -i\frac{d}{d\sigma}(A_{1}+A_{1})=(4\delta_{2}+\frac{3a^{2}}{8})A_{0}
\end{eqnarray}

from (49)

\begin{eqnarray}
 B_{0}^{^{\prime \prime }}(\tau )+\frac{a^{2}}{4}\left(4 \delta _{2}+\frac{3a^{2}}{8}%
\right) B_{0}=0,
\end{eqnarray}

the solution of  the $B_{0}$ is in form

\begin{eqnarray}
 B_{0}(\sigma )=C_{4}\cos \mu \sigma +C_{5}\sin \mu \sigma,
\end{eqnarray}

where $\mu=\sqrt{\delta_{2}+\frac{3a^{2}}{32}}$. We can choose
$a=\sqrt{2}$ as an arbitrary parameter. If we look at the
stability boundaries ,
stability occurs when $\delta _{2}<\frac{3}{16}$ and instability occurs when $\delta _{2}>%
\frac{3}{16}$. Thus, the higher order stability boundary is given
by:

\begin{eqnarray}
 \Delta =1\pm\frac{1}{2}\epsilon -\frac{3}{16}\epsilon
^{2}+O(\epsilon ^{3}),\qquad \epsilon \longrightarrow 0.
\end{eqnarray}

Then $A_{1}(\sigma )$ is given as

\begin{eqnarray}
 A_{1}(\sigma )=4\mu \left( -C_{4}\sin \mu \sigma +C_{5}\cos \mu
\sigma \right).
\end{eqnarray}

Using the initial conditions as $U_{0}(0)=0, \dot{U}_{0}
(0)=1,\quad U_{1}(0,0)=0,\dot{U}_{1}(0)=0$ , $U_{0}$ and $U_{1}:$

\begin{eqnarray}
 U_{0}(t,\sigma )=2\cos \mu \sigma \sin \frac{t}{2},
\end{eqnarray}

\begin{eqnarray}
U_{1}(t,\sigma )=-8\mu \sin \mu \sigma \cos \frac{t}{2}.
\end{eqnarray}

Using similar algorithm in here, it is easy to calculate in a
power series in~$\varepsilon $~ up to a high order again.

\section{Stability Boundaries by Using Von Zeipel Method}

The classical method of generating the canonical transformations
is called Von Zeipel's method. The problem with this method is the
awkward mixture of odd and new variables that has to be
unscrambled. To find higher order approximations ,Von Zeipel [18]
evaluated a technique . The main idea of the technique is to
expand the generating function S in powers of a small parameter
$\epsilon$ . First of all by using generalized momentum vector p
and coordinate vector q, one can write following canonical
equations of motion [17]

\begin{eqnarray}
\dot{q}=\frac{\partial H}{\partial p}
\end{eqnarray}

\begin{eqnarray}
\dot{P_{i}}=-\frac{\partial H}{\partial q}
\end{eqnarray}

Under a transformation from $\textbf{q}$ and $\textbf{p}$ to
$\textbf{Q}(\textbf{q},\textbf{p},t)$ and
$\textbf{P}(\textbf{q},\textbf{p},t)$ are transformed into :

\begin{eqnarray}
\dot{Q_{i}}=f_{i}(\textbf{P},\textbf{Q},t)
\end{eqnarray}

\begin{eqnarray}
\dot{P_{i}}=g_{i}(\textbf{P},\textbf{Q},t)
\end{eqnarray}

Here is the function K(\textbf{P},\textbf{Q},t)

\begin{eqnarray}
f_{i}=\frac{\partial K}{\partial P_{i}}, g_{i}=-\frac{\partial
K}{\partial Q_{i}}
\end{eqnarray}

\begin{eqnarray}
Q_{i}=\frac{\partial K}{dP_{i}}, P_{i}=-\frac{\partial K}{dQ_{i}}
\end{eqnarray}

Q and P are called canonical variables and the canonical
transformations can be generated using S(\textbf{p},\textbf{q},t)
which is generating function[19] :

\begin{eqnarray}
P_{i}=\frac{\partial S}{\partial q_{i}}, Q_{i}=-\frac{\partial
S}{\partial p_{i}}
\end{eqnarray}

When the equations above are solved for

\begin{eqnarray}
\textbf{q}=\textbf{q(\textbf{P},\textbf{Q},t)}
\end{eqnarray}

and

\begin{eqnarray}
\textbf{p}=\textbf{p(\textbf{P},\textbf{Q},t)}
\end{eqnarray}

K is given as:

\begin{eqnarray}
K(P,Q,t)=H(p(P,Q,t),q(P,Q,t),t)+\frac{\partial S}{\partial t}
\end{eqnarray}

Also S must satisfy the Hamilton-Jacobi equation:

\begin{eqnarray}
H(\frac{\partial S}{dq_{1}},..,\frac{\partial S}{\partial
q_{N}},q_{1},q_{2},...,q_{N},t)+\frac{\partial S}{\partial t}=0.
\end{eqnarray}

Complete solutions of the equation are not available for general
H. If $H=H_{0}+H_{1}$ where $H_{1}$ is small compared to $H_{0}$
and a complete solution $
S_{0}(P_{1},...,P_{N},q_{1},...,q_{N},t)$ is available for

\begin{eqnarray}
H_{0}(\frac{\partial S_{0}}{\partial q_{1}},..,\frac{\partial
S_{0}}{\partial q_{N}},q_{1},q_{2},...,q_{N},t)+\frac{\partial
S_{0}}{\partial t}=0.
\end{eqnarray}

We can use a generating function $
S=S_{0}(P_{1},...,P_{N},q_{1},...,q_{N},t)$ where
$\dot{P_{i}}=-\frac{\partial K}{\partial Q_{i}}$,
$\dot{Q_{i}}=-\frac{\partial K}{\partial P_{i}}$ and

\begin{eqnarray}
K=H_{0}+\tilde{H}+\frac{\partial S_{0}}{\partial t}=\tilde{H}
\end{eqnarray}

The method is the same as the generalized method of averaging[17].
Stern showed [20] for Hamiltonian systems, Kruskal's technique is
[21] equivalent to Von Zeipel's technique. The system under
discussion is described by the Hamiltonian

\begin{eqnarray}
H(\textbf{p},\textbf{q},t)=\sum(\varepsilon^{n}H_{n}(\textbf{p},\textbf{q},t)),
\varepsilon\ll 1
\end{eqnarray}

If $S_{0}=S_{0}(\textbf{P},\textbf{q},t)$ is a complete solution
of the Hamilton-Jacobi equation, $\varepsilon\ll1$

\begin{eqnarray}
H_{0}[\frac{\partial S_{0}}{\partial
\textbf{q}},\textbf{q},t]+\frac{\partial S_{0}}{\partial t}=0.
\end{eqnarray}

and the equations $\textbf{p}=\textbf{p}(\textbf{P},\textbf{Q},t)$
, $\textbf{q}=\textbf{q}(\textbf{P},\textbf{Q},t)$ are the
solutions of

\begin{eqnarray}
p_{i}=\frac{\partial S_{0}}{\partial q_{i}}
\end{eqnarray}

\begin{eqnarray}
Q_{i}=\frac{\partial S_{0}}{\partial P_{i}}
\end{eqnarray}

Assume $\textbf{P}$ and $\textbf{Q}$ to be time varying and
generating function $S=S_{0}(\textbf{P},\textbf{q},t)$ to
transform from the canonical system \textbf{p} and \textbf{q} to
the canonical system \textbf{P} and \textbf{Q}, the Hamiltonian is
transformed into

\begin{eqnarray}
\tilde{H}=\sum(\varepsilon^{n}H_{n}[\textbf{p},\textbf{q},t])+\frac{\partial
S_{0}}{\partial t} \sum(\varepsilon^{n}\tilde{H}_{n})
\end{eqnarray}

Hence $\textbf{P}$ and $\textbf{Q}$  are stated by the variational
equations

\begin{eqnarray}
\dot{\textbf{P}}=-\sum(\varepsilon^{n}\frac{\partial
\tilde{H}_{n}}{\partial \textbf{Q}})
\end{eqnarray}

\begin{eqnarray}
\dot{\textbf{Q}}=\sum(\varepsilon^{n}\frac{\partial
\tilde{H}_{n}}{\partial \textbf{P}})
\end{eqnarray}

Using the generating function S, to determine an approximate
solution to (75) and (76) to any order, introduce a transformation
from the canonical system $\textbf{P}$ and $\textbf{Q}$ to the new
canonical system

\begin{eqnarray}
S=\sum(P_{i}^{*}Q_{i})+\sum(\varepsilon_{n}S_{n}(P^{*},\textbf{Q},t))
\end{eqnarray}

then

\begin{eqnarray}
P_{i}=P^{*}_{i}+\sum(\varepsilon^{n}\frac{\partial S_{n}}{\partial
Q_{i}})
\end{eqnarray}

so $\tilde{H}$ is transformed into

\begin{eqnarray}
K\equiv\sum\varepsilon^{n}K_{n}(P^{*},Q,t)=\sum\varepsilon^{n}\tilde{H_{n}}+\sum\varepsilon^{n}\frac{\partial
S_{n}}{\partial t}
\end{eqnarray}

If the right-hand side of (79) is expanded for small $\varepsilon$
and then equate the coefficients depend on $\varepsilon$ on both
sides to obtain

\begin{eqnarray}
K_{1}=\tilde{H_{1}}+\frac{\partial S_{1}}{\partial t}
\end{eqnarray}

\begin{eqnarray}
K_{2}=\tilde{H_{2}}+\sum\frac{\partial S_{1}}{\partial
Q_{i}}\frac{\partial \tilde{H}_{1}}{\partial P_{i}}+\frac{\partial
S_{2}}{\partial t}
\end{eqnarray}

If this procedure is applied to the (8) with using $a=\sqrt{2}$

\begin{eqnarray}
\ddot{U}+(\frac{\omega^{2}}{4}+\frac{1}{4}\varepsilon
cost-\frac{3}{16}\varepsilon^{2}sin^{2}t)U=0
\end{eqnarray}

If we write $H_{0}$ and $\tilde{H}$ for (82)

\begin{eqnarray}
H_{0}=\frac{1}{2}(p^{2}+\frac{\omega^{2}}{4}q^{2})
\end{eqnarray}

\begin{eqnarray}
\tilde{H}=\frac{1}{2}(\frac{1}{4}\varepsilon
cost-\frac{3}{16}\varepsilon^{2}sin^{2}t)q^{2}
\end{eqnarray}

(83) and (84) are the hamiltonians for disturbed VdP oscillator
equation. The Hamilton-Jacobi equation corresponding to the case
$\varepsilon=0$ is

\begin{eqnarray}
\frac{1}{2}[(S'(q))^{2}+\frac{\omega^{2}}{4}q^{2}]+\frac{\partial
S}{\partial t }=0
\end{eqnarray}

The equation above can be solved by separation of variables ;

\begin{eqnarray}
S=S_{1}(q)+\sigma(t)
\end{eqnarray}

Then (86) separates into $\dot{\sigma}=-\alpha$ and
$\sigma=-\alpha t$

We find $S_{1}$ generator and $\beta$ new coordinate and $q$ as:

\begin{eqnarray}
S_{1}=-\alpha t+\int\sqrt{2\alpha-\frac{\omega^{2}}{4}q^{2}dq}
\end{eqnarray}

\begin{eqnarray}
\beta=-t+\frac{2}{\omega}arcsin(\frac{q\omega}{2\sqrt{2\alpha}})
\end{eqnarray}

\begin{eqnarray}
q=\frac{2\sqrt{2\alpha}}{\omega}cos(\frac{\omega(t+\beta)}{2})
\end{eqnarray}

Hence $\alpha$ and $\beta$ are canonical variables with respect to

\begin{eqnarray}
\tilde{H}=\frac{4\alpha}{\omega^{2}}cos^{2}(\frac{\omega(t+\beta)}{2})[\frac{1}{4}\varepsilon
cost-\frac{3}{16}\varepsilon^{2}sin^{2}t]
\end{eqnarray}

and $\dot{\alpha}=-\frac{\partial \tilde{H}}{\partial \beta }$ ,
$\dot{\alpha}=\frac{\partial \tilde{H}}{\partial \alpha }$ to
determine an approximate solution to these equations, introduce a
near-identity transformation from $\alpha$ and $\beta$ to
$\alpha^{*}$ and $\beta^{*}$ using generating function

\begin{eqnarray}
S=\alpha^{*}\beta+\varepsilon
S_{1}(\alpha^{*},\beta,t)+\varepsilon^{2}S_{2}(\alpha^{*},\beta,t)+...
\end{eqnarray}

Hence

\begin{eqnarray}
\alpha=\alpha^{*}+\varepsilon\frac{\partial {S_{1}}}{\partial
\beta }+...
\end{eqnarray}

Using (91) and (92), $K$ can be written as

\begin{eqnarray}
K=\varepsilon K_{1}+\varepsilon^{2}K_{2}+...
\end{eqnarray}

Hence $K_{1}$ and $K_{2}$ are

\begin{eqnarray}
K_{1}=\frac{\partial S_{1}}{\partial
t}+\frac{4\alpha^{*}}{\omega^{2}}
cos^{2}(\frac{\omega(t+\beta)}{2})cost
\end{eqnarray}

\begin{eqnarray}
K_{2}=\frac{\partial S_{2}}{\partial
t}+\frac{4}{\omega^{2}}\frac{\partial S_{1}}{\partial
\beta}cos^{2}(\frac{\omega(t+\beta)}{2})
cost-\frac{3\alpha^{*}}{4\omega^{2}}cos^{2}(\frac{\omega(t+\omega\beta)}{2})
sin^{2}t
\end{eqnarray}

Using some trigonometric relations, the equations(94) above can be
written in these forms:

\begin{eqnarray}
K_{1}=\frac{\partial S_{1}}{\partial
t}+\frac{\alpha^{*}}{\omega^{2}}(2cost+cos((\omega+1)t+\omega\beta)+cos((\omega-1)t)+\omega\beta)
\end{eqnarray}

in the same way one can write $K_{2}$ .In the case of
$\omega\neq1$ all the terms on the right-hand side of (96) are
fast varying. Hence $K_{1}=0$ and $S_{1}$ is

\begin{eqnarray}
S_{1}=-\frac{\alpha^{*}}{\omega^{2}}(2sint+\frac{1}{\omega+1}sin((\omega+1)t+\omega\beta)
+\frac{1}{\omega-1}sin((\omega-1)t+\omega\beta))
\end{eqnarray}


In the case of $\omega\approx 1$ , $cos((\omega-1)t+\omega\beta)$
is slowly varying because  $S_{1}$ is singular at $\omega\approx1$
as it is seen from (100). By  equating $K_{1}$ to long-period in
(100),we have

\begin{eqnarray}
K_{1}=\frac{\alpha^{*}}{\omega^{2}} cos((\omega-1)t+\omega\beta)
\end{eqnarray}

Substituting $S_{1}$ into (95) one can easily get  $K_{2}$ and
equating the $K_{2}$  to the long terms in this equation ,we have

\begin{eqnarray}
K_{2}=-(\frac{1}{\omega^{3}(\omega+1)}+\frac{3}{16\omega^{2}})\alpha^{*}
\end{eqnarray}

Therefore $K$ can be written to second order:

\begin{eqnarray}
K=\frac{\alpha^{*}\varepsilon}{\omega^{2}}cos((\omega-1)t+\omega\beta)
-\frac{\alpha^{*}\varepsilon^{2}}{\omega^{2}}(\frac{3}{16}+\frac{1}{\omega(\omega+1)})
\end{eqnarray}

It is obvious that  $\alpha$ and  $\beta$ in terms of $\alpha^{*}$
$\beta^{*}$ :

\begin{eqnarray}
\alpha=\alpha^{*}-\frac{\varepsilon\alpha^{*}}{\omega(\omega+1)}
cos ((\omega+1)t+\omega\beta)
\end{eqnarray}

\begin{eqnarray}
\beta=\beta^{*}-\frac{2\varepsilon}{\omega^{2}}(sin
t+\frac{1}{2(\omega+1)}sin((\omega+1)t+\omega\beta))
\end{eqnarray}

Removing the dependence of $K$ on $t$ by changing from
$\alpha^{*}$ and $\beta^{*}$ to $\alpha^{'}$ and $\beta^{'}$ by
using

\begin{eqnarray}
S^{'}=\alpha^{'}((\omega-1)t+\omega\beta^{*})
\end{eqnarray}

and

\begin{eqnarray}
\alpha^{*}=\frac{\partial S^{'}}{\partial\beta^{*}}
\end{eqnarray}

\begin{eqnarray}
\beta^{'}=\frac{\partial S^{'}}{\partial \alpha^{'}}
\end{eqnarray}

Now $K$ can be written as  $K^{'}=K+\frac{\partial S^{'}}{\partial
t}$ and:

\begin{eqnarray}
K=\frac{\varepsilon\alpha^{'}}{\omega}cos\beta^{'}-
\frac{\varepsilon^{2}\alpha^{'}}{\omega}(\frac{3}{16}+\frac{1}{\omega}(\omega+1))
+(\omega-1)\alpha^{'}
\end{eqnarray}

Therefore

\begin{eqnarray}
\dot{\alpha}^{'}=\frac{\varepsilon\alpha^{'}}{\omega}sin\beta^{'}
\end{eqnarray}

\begin{eqnarray}
\dot{\beta}^{'}=-\frac{\varepsilon}{\omega}cos\beta^{'}-\frac{3\varepsilon^{2}}{16\omega}
-\frac{\varepsilon^{2}}{\omega^{2}(\omega+1)}+\omega-1
\end{eqnarray}

The solution of (107-108) is

\begin{eqnarray}
ln \alpha^{'}=ln
[\frac{\varepsilon}{\omega}cos\beta^{'}+(\frac{3}{16\omega}+
\frac{1}{\omega^{2}(\omega+1)})\varepsilon^{2}+\omega-1]
\end{eqnarray}

\begin{eqnarray}
\frac{\varepsilon}{\omega}=|\omega-1+\varepsilon^{2}(\frac{3}{16\omega}+\frac{1}{\omega^{2}(\omega+1)})|
\end{eqnarray}

so the transition curves in agreement with those obtained in
sections(3-4) by multiple-scale are:

\begin{eqnarray}
\frac{\omega^{2}}{4}\rightarrow
 1+\frac{1}{2}\varepsilon-\frac{3}{16}\varepsilon^{2}
\end{eqnarray}

\section{Conclusions}

In this study, we examined the disturbed Van der Pol oscillator
equation and applied MSPT and Von Zeipel method to the  equation
which is turned into a second order \ differential equation with
periodic\ coefficients. For convenient solutions, necessary
transforms are applied to this variational equation. We used
$\left( \Delta ,\epsilon \right) $ parameter space in the equation
in order to examine the stability conditions through boundaries.
As it has already been shown[2,3], the use of MSPT allows for the
elimination of all harmonic and subharmonic solution related
secular
producing terms appearing in the evolution of the $U_{0}(t,\tau )$ and $%
U_{1}(t,\tau )$ . So the method worked effectively for terminating
the secular terms grow and more accurate results have been
obtained for sub-harmonics by using different time scales. \ By
comparing the results obtained for two\ different time scales,
stability and instability conditions depend on constants and one
can see if \ the motion is strictly periodic or not. Our results
show that solutions of terminated nonlinear structure of Van der
Pol equation are periodic and stable and are similar to solutions
of the original VdP equation[17]. To determine a first
approximation to our Hamiltonian system for disturbed VdP equation
, Von Zeipel method which is a perturbation method for classical
Hamiltonian systems using an averaging procedure in phase space is
used. Transition curves are in agreement with those obtained by
MSPT methods. As a result, we can say that this work contains
useful tools for extension of these applications in the case of
disturbed motions.

\newpage
\section{References}

1.  $C.M.Bender, S.A. Orszag, Advanced Mathematical Methods For
   Scientists And Engineers,McGraw-Hill Book Company,New York
   (1978).$

2. $Carl M. Bender and Luis M.A.Bettencourt,Multiple-scale
    analysis of quantum systems,Physical Review D,54,
    7710-7723,1996.$

3.  $C. M. Bender and L. M. A. Bettencourt ,Multiple-scale
    Analysis of the Quantum Anharmonic Oscillator,Physical Review
    Letters  77, 4114 (1996).$

4.  $P. A. Sturrock, Nonlinear effects in electron plasmas, Proc.
    Roy. Soc. (London),A242,277-299, 1957.$

5.  $E. A. Frieman, On a new method in the theory of irreversible
    processes,J.Math. Phys. ,4,410-418, 1963.$

6.  $A. H. Nayfeh, Nonlinear oscillations in a hot electron
    plasma,   Phys. Fluids,8, 1896-1898, 1965.$

7.  $Peter B. Kahn and Yair Zarmi, arxiv:nlin/0206045$

8.  $Van der Pol and Van der Mark, Nature 120,363(1927).$

9.  $R.E. Mickens, Fractional Van der Pol equations,Journal of
   Sound and Vibration(2003),259(2),457-460.$

10. $A.Y.T. Leung and Q.C. Zhang, Complex Normal Form For Strongly
    Non-Linear Vibration Systems  Exemplified by Duffing-Van der Pol
    Equation,Journal of Sound and Vibration(1998),213(5),907-914.$

11. $Je-Chiang Tsai,Tai-Yih Tso, the Critical Phase Curve of Van
    der Pol Equation, Journal of Taiwan Normal Univ.: Mathematics,
    Science and Technology,2001,46,pp 1-12.$

12.  $J. Xu, K.W. Chung, Effects of Time Delayed Position Feedback
    on a Van der Pol Oscillator,Physica D3096(2003)1-23.$

13.  $P.E. Phillipson, P. Schuster, Dynamics of Relaxation
    Oscillations, International Journal of
    Bifurcation\ \ and Chaos, vol.11,No.5(2001) 1471-1482.$

14. $F. M. Fernandez and A. Pathak, Frequency operator for
    anharmonic oscillators by perturbation theory, J. Phys. A: Math. Gen.
    36(2003),5061-5066.$

15. $T. L. Saaty,J. Bram, Nonlinear Mathematics,Mc-Graw-Hill Book
    Company,New York (1964).$

16. $P. Glenndining, Stability,Instability and Chaos: an
    introduction to the theory of non-linear  differential
    equations,Cambridge University Press(1994).$

17. $A. H. Nayfeh, Perturbation Methods, John Wiley and
    Sons,NY(1973)$

18. $H. Von Zeipel , Movements of minor planets, Ark. Mat. Astron.
    Fysik,Stockholm,11,No.1,1-58,No7,1-62, 1916.$

19. $H. Goldstein, Classical Mechanics, Addison-Wesley, 1965.$

20. $D. P. Kruskal's perturbation method, J. Math.
    Phys.,11,2771-2775,1970.$

21. $M. Kruskal, in Dynamical systems,Theory and Applications,
    Lecture Notes in Physics 38,ed. by J. Moser
    (Springer-Verlag,Berlin,1975).$
\end{document}