where x is the position coordinate---which is a function of the time t, and μ is a scalar parameter indicating the nonlinearity and the strength of the damping. This unforced differential equation was first introduced by Lord Rayleigh in 1883 to model the oscillations of a clarinet reed. However, the Dutch electrical engineer and physicist Balthasar van der Pol (1889--1959), who
was known as Balth, was the first to initiate modern experimental and
theoretical dynamics in the laboratory during the 1920s and
1930s. Balth was the first who used the above model in his
research. He investigated electrical circuits by employing vacuum tubes and found that they have stable oscillations, now called limit cycles. He was awarded the Institute of Radio Engineers (now the IEEE) Medal of Honor in 1935. The asteroid 10443 van der Pol was named after him. Van der Pol built a number of electronic circuit models of the human heart to study the range of stability of heart dynamics. His investigations with adding an external driving signal were analogous to the situation in which a real heart is driven by a pacemaker.
In the September 1927 issue of the British journal Nature, Van der Pol and his colleague J. van der Mark published an article entitled "Frequency Demultiplication" in which they reported an "irregular noise" before transition from one subharmonical regime to another one. By reconstructing his electronic tube circuit, we now know that they had discovered deterministic chaos. Their paper is probably one of the first experimental reports of chaos. In particular, they considered the forced oscillations of a triode, but in the field of relaxation oscillations:
\[
\ddot{v} - \alpha \left( 1- v^2 \right) + \omega_0^2 v = \omega_1^2 E\,\sin \left( n \omega_1 t \right) \qquad\mbox{with} \quad \varepsilon = \frac{\alpha}{\omega_0} \gg 1 .
\]
They showed that for n = 3, the above forced van der Pol equation admits an exact solution \( v(t) = 2\,\cos \left( \omega_0 t \right) \) for E = 2ε. It therefore appears that the frequency of this solution is three times the forcing one.
During the 1930s, Russian mathematicians Nikolai M. Krylov (1879--1955) and Nikolai N. Bogolyubov (1909--1992), inscribed the slowly varying amplitudes method unveiled a few years earlier by van der Pol in the framework of a theory that now bares their names. The Krylov--Bogolyubov averaging method is based on the averaging principle when the exact differential equation of the motion is replaced by its averaged version. Their work attracted quite a bit of attention with full support of the famous French mathematician Jacques Hadamard (1865--1963). They put mathematical foundation into van der Pol experimental research.
In 1945, M.L. Cartwright and J.E. Litlewood published the article where they showed that the van der Pol equation admits singular solutions for large values of parameter μ. Later in 1949,
N. Levinson analytically analyzing the van der Pol equation confirms
that it has singular solutions. On the basis of numerical
simulations, it was observed that the transition from periodic
oscillations to chaotic ones is possible. Further history of
analyzing van der Pol oscillator could be found in the article by
Ginoux and Letellier.
To analyze stability of solutions to van der Pol equation, we rewrite his
second order differential equation
If μ is positive, then the derivative of V with respect to
the van der Pol equation is also positive for |x| < 1, so
the origin is unstable. This happens because energy is added to the
system when the displacement is small, i.e. |x| < 1. The
Lyapunov
theorem guarantees that every trajectory starting inside the
unit circle tends to a closed loop.
On the other hand, when x is large, the cubic
term x3 becomes dominant and the damping becomes
positive. In this case, \( \dot{V}(x,y) \le 0
. \) Energy is dissipated for large displacement
corresponding to |x| > 1. Therefore, the system is expected
to be restricted in some area around the fixed point and the system
approaches a limit cycle.
Using the Liénard's transformation\( y = x - x^3 / - \dot{x} /\epsilon , \)
the van der Pol equation can be rewritten as
Actually, the van der Pol equation is a particular case of more general Liénard equation:
\[
\ddot{x} + f(x) \,\dot{x} + g(x) =0 ,
\]
named after the French physicist Alfred-Marie Liénard (1869--1958). The Van der Pol equation has no exact, analytic solution, but it has a limit cycle.
Example:
In order to obtain information regarding the approach to the limit cycle,
we use a powerful perturbation method called the method of averaging. We
begin with a system of the more general form:
\[
\begin{split}
x &= a(t)\,\cos \left( t + \psi (t) \right) , \\
\dot{x} &= -a(t) \,\sin \left( t + \psi (t) \right) .
\end{split}
\]
Our motivation for this ansatz is that when µ is zero, the van der Pol
equation has its solution
of the above form with a and ψ constants. For small values of µ we expect
the same form of the solution to be approximately valid, but now a and ψ
are expected to be slowly varying functions
of t. Differentiating equations for x and its
derivative, we obtain the relation for unknown parameters a and ψ:
So far our treatment has been exact and is essentially the procedure of variation
of parameters which is used to obtain particular solutions to nonhomogenous
linear differential equations.
Now we introduce the following approximation: since µ is small,
\( \dot{a}
\quad\mbox{and}\quad \dot{\psi} , \)
are also small. Hence a(t) and ψ(t) are slowly varying functions of t.
Thus over one cycle of oscillations the quantities a(t) and ψ(t) on the right
hand sides of the latter can be treated as nearly constant, and thus
these right hand sides may be replaced by their averages:
The right hand side of the latter is zero. The averaging in the
previous equation for \( \dot{\psi} \) can be
done using the following trigonometric identities:
Our next step in analyzing the trajectories of the van der Pol equation is to linearize the corresponding system of first order equations. Although Mathematica has a dedicated command JacobianMatrix, we calculate it directly. First, we build a matrix of its slope functions.
vdot[t_] = mu (1 - x[t]^2)*x'[t] - x[t]
Out[1]= -x[t] + mu (1 - x[t]^2) Derivative[1][x][t]
M = {{0, 1}, {D[vdot[t], x[t]], D[vdot[t], x'[t]]}}
Out[2]= {{0, 1}, {-1 - 2 mu x[t] Derivative[1][x][t], mu (1 - x[t]^2)}}
The following example demonstrates a series of Mathematica codes for solving numerically the van der Pol equation. Note for numerical procedure NDSolve it does not matter whether the equation is homogeneous or driven.
We start with a particular numerical value of parameter μ = 0.5. the initial conditions are chosen very close to zero.
Then we take another initial condition far away from the origin.
Example:
Chaos synchronization has received increasing attention as an interesting phe-
nomenon of practical importance in coupled chaotic oscillators, due to its theoretical challenge and its great potential applications in secure communication, chemical reaction, biological systems and so on. Various types of synchronization have been studied.
One particular characteristic in the van der Pol model is that its phase depends on initial conditions. Therefore, if two van der Pol
oscillators run with different initial conditions, their trajectory will finally circulate on the same limit cycle with different phases φ1 and φ2. The chaos synchronization problem can be formulated as follows: Given a chaotic system, which is considered as the master (or driving) system, and another identical system, which is considered as the slave (or response)
system, the aim is to force the response of the slave system to synchronize
the master system in such a way that the dynamical behaviors of these two
systems be identical after a transient time.
■
Liénard's theorem can be used to prove that the system has a limit
cycle. Applying the Liénard transformation \( y = x - x^3 / 3 -
\dot{x}/\mu , \) where the dot indicates the time derivative, the
Van der Pol oscillator can be written in its two-dimensional form:
The forced, or driven, Van der Pol oscillator takes the 'original' function and adds a driving function Asin(ωt) to give a differential equation of the form:
One can also write a time-independent Hamiltonian formalism for the
Van der Pol oscillator by augmenting it to a four-dimensional
autonomous dynamical system using an auxiliary second-order nonlinear
differential equation as follows:
Note that the dynamics of the original Van der Pol oscillator is not
affected due to the one-way coupling between the time-evolutions of x
and y variables. A Hamiltonian H for this system of equations can be
shown to be
where \( p_x = \dot{y} + \mu \left( 1 - x^2
\right) y \) and \( p_y = \dot{x}
\) are the conjugate momenta corresponding to x and y, respectively.
Van der Pol, B. and Van der Mark, J., Frequency demultiplication, Nature, 120, Number 3019, 1927, pp. 363--364.
M.L. Cartwright and J.E. Litlewood, "On Non‐Linear Differential
Equations of the Second Order: I. the Equation \(
y'' - k (1 - y^2) y' +y = b\lambda k\,\cos (\lambda t + \alpha
), \) k large", Journal of the London
Mathematical Society, 20, Issue 3, pp. 180--189, 1945; https://doi.org/10.1112/jlms/s1-20.3.180
Norman Levinson,
A Second Order Differential Equation with Singular Solutions
Norman Levinson, Annals of Mathematics,
Second Series, Vol. 50, No. 1 (Jan., 1949), pp. 127-153 (27 pages).
W. A. Albarakati, N. G. Lloyd, &
J. M. Pearson, Transformation
to Liénard form, Electronic Journal of Differential
Equations, Vol. 2000
(2000), No. 76, pp. 1--11.