van der Pol Oscillator

from Essential Mathematica for Students of Science
© James J. Kelly, 1998
Last revised: January, 2006

Introduction

An interesting feature of many physical or biological systems is the appearance of limit-cycle oscillations.  The classic example is the van der Pol equation

(((∂^2) x[t])/(∂(t^2)))==-x[t]+ε (1 - x[t]^2) ((∂x[t])/(∂t))

where ε is a "damping parameter".  When ε->0, the equation reduces to a simple-harmonic oscillator with unit frequency.  [Note that we can always choose the time unit to normalize the frequency to unity.]  For ε ≠ 0 there is an interesting nonlinear velocity-dependent force which has a damping effect for large displacements and an anti-damping effect for small displacements; thus, this nonconservative force supplies energy for small displacements or drains energy for large displacements.  The net effect is to drive the system toward a stable limit cycle for which the damping and anti-damping effects are balanced over the course of each cycle.  This equation was originally used by Balthasar van der Pol in the 1920s to study nonlinearity in electrical circuits containing vacuum tubes and iron-core inductors, but more recently has found application to lasers and biological systems as well.

It is useful to re-express the second-order differential equation as a pair of coupled first-order equations of the form

{x'[t]==y[t],y'[t]==-x[t] + ε (1 - x[t]^2) y[
                  t]}

where a point in phase space is represented by the "position" vector r = {x, y} and the corresponding phase velocity is then v = (∂_t) r = {(x^′), 
            (y^′)}.  The equations of motion now take the form

(∂_t) x[t] == (Subscript[v, x])[x, y]
(∂_t) y[t] == 
      (Subscript[v, y])[x, y]

where it is important to recognize that there is no explicit time dependence in the phase velocity.   Thus, the driving terms can be described as a time-independent velocity field in phase space.  Systems of this type are often described as autonomous because, in the absence of explicit time dependence, the phase point at any time uniquely determines the subsequent trajectory.  Thus, each point in phase space belongs to a unique trajectory and trajectories do not intersect.

The behavior of autonomous systems can be analyzed in several steps.  First, the fixed points (Subscript[r, fp]) for the system are determined by the solutions to the system of equations v == 0.  Second, near any fixed-point we can linearize the equations of motion by performing a series expansion of the velocity field.  The resulting linearized equations of motion can now be expressed in the form

(∂_t) r[t]  ==  
      A·(r[t] - (Subscript[r, fp]))

where A = (Subscript[((∂((Subscript[v, i])/∂)(Subscript[r, j]))), fp]) is the stability matrix composed of derivatives of velocity with respect to position evaluated at the fixed point.  Next, we construct displacement eigenvectors (Subscript[δ, i]) satisfying

\(TraditionalForm\`\(\(\ \)\(A·δ\_i == \(λ\_i\) δ\_i\)\)\)

where (Subscript[λ, i]) is an eigenvalue that governs the time dependence of displacements in the direction of (Subscript[δ, i]).  Finally, the displacement of the phase point can be represented by the expansion

r[t] - 
        (Subscript[r, fp]) == (∑_i) 
              (Subscript[a, i]) (e^((Subscript[λ, i]) 
                t)) (Subscript[δ, i])

where the constant coefficients (Subscript[a, i]) represent the initial displacement

r[0] - 
        (Subscript[r, fp]) == (∑_i) (Subscript[a, i]) (Subscript[δ, i])

The dynamical behavior of the system near its fixed points can now be classified according to the properties of the eigenvalues of the stability matrix.  For simplicity, consider a two-dimensional system such that the eigenvalues are either both real or are a pair of complex conjugates.  Although one can classify all possibilities rather completely, for the present purposes it is sufficient to consider the case of complex eigenvalues.  If the real part of an eigenvalue is negative, the corresponding fixed point is considered stable because small initial displacements are damped exponentially, whereas if the real part is positive the fixed point is unstable because small initial displacements grow exponentially.  However, if it should happen that for larger displacements there exists a closed curve in phase space on which the normal component of the phase velocity is everywhere inward, the system will eventually find itself trapped on a limit cycle.  This occurs because points inside the limit cycle are driven outward while points outside are driven inward, such that the trajectory relaxes exponentially toward the limit cycle.  This phenomenon is explored in this notebook.  Furthermore, this system provides an example of a relaxation oscillator for which stress accumulates slowly and then is relieved rapidly.  This behavior is illustrated by strobe plots in phase space.

Initialization

ClearAll["Global`*"];
Off[General::spell, General::spell1];
$TextStyle = {FontFamily -> "Times", 
          FontSize -> 12};

$DefaultFont={"Times", 12}; 
$TextStyle={FontFamily->"Times", FontSize->12,FontSlant->"Italic"};

Needs["Utilities`Notation`"];
Needs["Graphics`PlotField`"];
Needs["Graphics`Graphics`"]

Off[Symbolize::"boxSymbolExists"];

Symbolize[(Subscript[x, 0])]; Symbolize[(Subscript[y, 0])]; Symbolize[(Subscript[t, max])];
Symbolize[(Subscript[x, max])]; Symbolize[(Subscript[y, max])]

Linear stability analysis

Define phase-space position and velocity

position = {x, y};
velocity = {y, ε (1 - (x^2)) y  -  x}

Locate fixed point

Although it is obvious by inspection that the only fixed point lies at the origin, perhaps it is worth using a procedure that can be generalized to other systems.

fixedPoints  =  
      Solve[MapThread[Equal, {velocity,  Table[0, {Length[velocity]}]}], 
        position]

Construct stability matrix

stabilityMatrix  =  
      Outer[D, velocity, position] /. 
          fixedPoints[[1]] // 
        Simplify

Examine eigenvalues

evalues  =  Eigenvalues[stabilityMatrix]

Plot[Evaluate[Re[evalues]], {ε, -4, 4}, 
        PlotStyle -> {RGBColor[1, 0, 0], 
            RGBColor[0, 0, 1]}];

Thus, we can identify four distinct regimes.  1) For -2 < ε < 0 the fixed-point is stable because the complex eigenvalues have negative real parts, such that nearby initial points will spiral inward toward the origin.   2) For 0  <  ε < 2 the fixed-point is unstable because the complex eigenvalues have positive real parts, such that nearby initial points will spiral outward.  3) For ε < -2 both eigenvalues are real and negative, producing exponential damping of initial displacements.  4) For ε > 2 both eigenvalues are real and positive, producing exponential growth of initial displacements.

Limit-cycle oscillations

The most interesting behavior occurs near an unstable fixed point.  Although the initial motion tends to carry the system away from the fixed point, if restoring forces are encountered for larger displacements, the system can find itself trapped on a closed path through phase space known as a limit cycle.  This occurs because points inside the limit cycle are driven outward, while points outside are driven inward.  This behavior can be illustrated by solving the equations of motion numerically.

Define display functions

Here we define several simple display functions that can be used to examine various aspects of limit-cycle behavior.  Although it is a relatively easy to generalize these functions to handle any pair of autonomous first-order equations, we leave that exercise to the student.

r[t]  =  position /. {x -> x[t], y -> y[t]};
v[t]  =  
        velocity  /. {x -> x[t],  y -> y[t]};

eq1 = MapThread[ 
        Equal, {D[r[t], t],  v[t]} ]

trajectory[ε_, x0_, y0_, tmax_, opts___Rule]:=NDSolve[{x'[t]==y[t],y'[t]==-x[
                        t] + ε (1 - x[t]^2) y[t],x[0] == x0,y[0] == y0},{x, y},{t, 0, tmax},opts,MaxSteps -> (10^4)][[1]]

velocityField[εε_, xmax_, ymax_, 
          opts___Rule]  :=  
        PlotVectorField[
          velocity /. ε -> εε, {x, -xmax, xmax}, {y, -ymax, ymax}, opts];

phaseCurve[x0_,  y0_, ε_, tmax_, opts___Rule]  :=  
        ParametricPlot[
          Evaluate[r[t] /. trajectory[ε, x0, y0, tmax]], {t, 0, 
            tmax}, opts];

StrobePlot[x0_,  y0_, ε_, tmin_, tmax_, tinc_, 
        opts___]  :=  
      Module[{sol}, sol = r[t] /. trajectory[ε, x0, y0, tmax]; 
        ListPlot[Table[sol, {t, tmin, tmax, tinc}], opts];]

xyCurves[x0_,  y0_, ε_, tmax_]  := 
        Plot[Evaluate[r[t] /. trajectory[ε, x0, y0, tmax]], {t, 
            0, tmax}, 
          PlotStyle -> {RGBColor[1, 0, 0], 
              RGBColor[0, 0, 1]}];

phaseAndField[x0_,  y0_, ε_, tmax_, xmax_, ymax_]  := 
      Show[{velocityField[ε, xmax, ymax, 
            DisplayFunction -> Identity], 
          phaseCurve[x0, y0, ε, tmax, 
            DisplayFunction -> Identity]}, 
        DisplayFunction -> $DisplayFunction]

Examine behavior for various initial conditions and relaxation parameters

Case 1: initial point inside limit cycle

phaseAndField[1, 1, 1, 20, 3, 3];

phaseAndField[0.1, 0, 1, 40, 2.5, 2.5];

xyCurves[0.1, 0, 1, 40];

Even systems which start very near the fixed point evolve toward the limit cycle quite quickly.  Notice that when the relaxation parameter ε is relatively small, the oscillation period is near the 2π period of the driving term.

Case 2: initial point outside limit cycle

phaseAndField[3, 1, 1, 20, 3.5, 3];

xyCurves[3, 1, 1, 20];

Case 3: large ε

The shape of the limit cycle depends strongly upon ε.

phaseCurve[1, 1, 3, 20];

phaseCurve[0.1, 0, 3, 20];

phaseAndField[0.1, 0, 3, 20, 3, 5];

phaseCurve[-3, -3, 3, 20];

It is also important to recognize that limit cycles for large ε are traversed at rather uneven rates and that their periods are much longer than the underlying 2π period for free oscillations (ε = 0).

xyCurves[-3, -3, 3, 30];

StrobePlot[-3, -3, 3, 15, 24, 0.1, 
        PlotRange -> All];

These figures illustrate the typical behavior of a relaxation oscillator for which stress accumulates slowly and is then released rapidly.  Thus, intervals of relatively slow variation alternate with brief intervals of rather rapid variation, with the ratio between these intervals being governed by a relaxation parameter.  The characteristics of relaxation oscillators are investigated in the next section.

(((∂^2) x[t])/(∂(t^2)))==ε (1 - x[t]^2) ((∂x[t])/(∂t))

(((∂^2) x[t])/(∂(t^2)))+ (x[t]^2 - β) ((∂x[t])/(∂t))+x[t] == 0

Relaxation oscillators

Introduction

There are an infinite variety of transformations which represent a second-order differential equation as a pair of first-order differential equations.  Further insight into the behavior of the van der Pol equation can be obtained by representing it in the form of a relaxation oscillator described by a pair of first-order equations of the form

{α ((∂x[t])/(∂t))==y[t] - f[x],   ((∂y[t])/(∂t))==-x[t]}

where the relaxation parameter 0 < α ≪ 1 is small and where f[x] is a continuous function.  An important feature of this representation is that the null clines are readily apparent.  A null cline is a curve in phase space upon which one of the phase variables is constant.  For this system it is clear that

\(TraditionalForm\`x == 
      0\ \ ⟹\ \ \(∂y\/∂t\) == 0\)
\(TraditionalForm\`y == 
      f[x]\ ⟹\ \ \(∂x\/∂t\) == 0\)

represent the null clines.  Fixed points occur where null clines intersect.  For this type of system, one of the null clines is simply the y-axis, and is trivial, while the nontrivial function f[x] specifies the dynamics.

Suppose that y is near the null cline f[x] , where nearness is gauged in terms of (α^(-1)); in other words, we consider the phase point to be near the null cline when |y - 
            f[x]| ≲ (α^(-1)).  Under those circumstances both components of the phase velocity will be of order unity.  Conversely, when |y - 
            f[x]| ≫ (α^(-1)) the velocity scaled by α becomes very large while the other velocity remains of order unity.  Thus, the relaxation parameter α regulates the rate at which various portions of the limit cycle are traversed.

Relaxation representation of van der Pol equation

First, we show that the van der Pol equation can be represented as a relaxation oscillator.  It is useful to define a function which establishes the equations for any choice of null cline.

RelaxationEquations[x_, y_, f_, 
        t_, α_] :=  {α D[x[t], t] == y[t] - f[x[t]], 
        D[y[t], t] == -x[t]}

Next, we show that a judicious choice of f[x] gives an equation

\(testEq = 
      RelaxationEquations[x, y, \((\(1\/3\) #\^3 - #)\) &, 
        t, α]\)

which can be transformed into the van der Pol equation by means of the variable changes \(TraditionalForm\`{t -> α\ t, α -> ε\^\(-2\)}\):

D[testEq[[1]], 
            t] /. (testEq[[2]] /. 
              Equal -> Rule) /. {(x^′)[
              t] -> ε (x^′)[
                t], (x^″)[
              t] -> (ε^2) (x^″)[
                t], α -> (ε^(-2))} // 
      Simplify

Thus, we have identified the characteristic null cline and have related the relaxation parameter to the damping parameter found in the original form van der Pol equation.  It is important to recognize that the time scale must be scaled by α to represent the van der Pol equation as a relaxation oscillator.

Characteristic behavior of relaxation oscillators

The following functions are designed to illustrate the characteristic behavior of relaxation oscillators by plotting the velocity field, the null clines, and phase trajectories together.  The first shows the phase trajectory as a continuous curve, while the second shows the phase velocities stroboscopically using points at equal time intervals.

\(\(\(Clear[RelaxationPlot];\)\(
\)
    \)\)
RelaxationPlot[f_, α_, xmax_, ymax_, x0_, y0_, tmax_, 
        opts___]  := 
      Show[{ParametricPlot[
            Evaluate[{x[t], y[t]} /. 
                NDSolve[
                  Join[RelaxationEquations[x, y, f, 
                      t, α], {x[0] == x0, y[0] == y0}], {x[
                      t], y[t]}, {t, 0, tmax}, MaxSteps -> (10^4)]], {t, 
              0, tmax}, DisplayFunction -> Identity], 
          PlotVectorField[{((y - f[x])/α), -x}, {x, -xmax, 
              xmax}, {y, -ymax, ymax}, DisplayFunction -> Identity], 
          Plot[f[x], {x, -xmax, xmax}, 
            PlotStyle -> RGBColor[1, 0, 0], 
            DisplayFunction -> Identity]}, opts, 
        DisplayFunction -> $DisplayFunction]

\(\(\(Clear[RelaxationStrobePlot];\)\(
\)
    \)\)
\(RelaxationStrobePlot[f_, α_, xmax_, ymax_, x0_, y0_, tmin_, 
        tmax_, tinc_, opts___]\  := 
      Module[{sol}, 
        sol = \(NDSolve[
              Join[RelaxationEquations[x, y, f, 
                  t, α], {x[0] == x0, y[0] == y0}], {x[t], 
                y[t]}, {t, 0, tmax}, 
              MaxSteps -> 
                10\^4]\)[[1]]; 
	
        Show[{Plot[f[x], {x, \(-xmax\), xmax}, 
              PlotStyle -> RGBColor[1, 0, 0], 
              DisplayFunction -> Identity], \ 
            ListPlot[Table[{x[t], y[t]} /. sol, {t, tmin, tmax, tinc}], 
              DisplayFunction -> Identity]}, opts, 
          DisplayFunction -> $DisplayFunction];]\)

We illustrate this function using the van der Pol oscillator for a fairly small relaxation parameter.  The null cline is a cubic equation with two turning points — the presence of these turning points is responsible for the appearance of limit cycle oscillations.

\(\(RelaxationPlot[\((\(1\/3\) #\^3 - #)\) &, 0.1, 4, 2, 0.1, 0. , 20, 
        PlotRange -> {{\(-3\), 3}, {\(-1.5\), 1.5}}];\)\)

\(\(RelaxationStrobePlot[\((\(1\/3\) #\^3 - #)\) &, 0.01, 4, 2, 0.1, 0. , 
        18, 20, 0.01, 
        PlotRange -> {{\(-3\), 3}, {\(-1.5\), 1.5}}];\)\)

The limit cycle now resembles a slightly curved parallelogram.  On the left side, x being negative requires the y velocity to be positive so that the phase point moves slowly upward and remains near the null cline.  However, once the phase point passes the peak of the null cline the positive velocity brings it into a region where the small relaxation parameter strongly amplifies the horizontal velocity.  The system traverses the top segment very rapidly until it passes through the right side of the null cline where the large positive horizontal velocity abruptly changes to a small negative value.  Now that the phase point is below the null cline, the vertical velocity is negative and it remains small as long as the system can remain near the null cline.  However, once the phase point overshoots the bottom of the null cline, it once again enters a regime with large horizontal velocity.  The system races to the left until crossing the null cline again where it slowly retraces the left segment of the limit cycle.  The smaller the value of α the more tightly the trajectory adheres to the null cline in its vertical segments and the flatter are its horizontal segments.

ε x'[t]==(1 - x[t]) x[t] + q y[t] - x[t] y[t]
	δ y'[t]==-q y[t] - x[t] y[t] + 2 f z[t]
	z'[t]==x[t] - z[t]

(∂_t) x[t] == 
      (Subscript[v, x])[x, y, ε]
(∂_t) y[t] == 
      (Subscript[v, y])[x, y, ε]

(((∂^2) x[t])/(∂(t^2)))==-x[t]+ε (1-(((((∂x[t])/(∂t))))^2)) ((∂x[t])/(∂t))