Back to the main Index

Van der Pol oscillator

Matthias König

The Van der Pol oscillator is a non-conservative oscillator with non-linear damping. It is one of the standard models analysed in dynamics. In this tutorial the Van der Pol oscillator is analysed using SBML.


The Van der Pol oscillator was originally proposed by the Dutch electrical engineer and physicist Balthasar van der Pol while he was working at Philips. Van der Pol found stable oscillations, which he called relaxation-oscillations[2] and are now known as limit cycles, in electrical circuits employing vacuum tubes. When these circuits were driven near the limit cycle they become entrained, i.e. the driving signal pulls the current along with it. Van der Pol and his colleague, van der Mark, reported in the September 1927 issue of Nature that at certain drive frequencies an irregular noise was heard. This irregular noise was always heard near the natural entrainment frequencies. This was one of the first discovered instances of deterministic chaos.

The Van der Pol equation has a long history of being used in both the physical and biological sciences. For instance, in biology, Fitzhugh and Nagumo extended the equation in a planar field as a model for action potentials of neurons.


The Van der Pol oscillator evolves in time according to the second-order differential equation

$$\frac{d^2x}{dt^2} - µ(1-x^2)\frac{dx}{dt} + x = 0$$

with initial conditions $x=−2$ and $\frac{dx}{dt}=0$. $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.

To create a SBML file, we need to convert the second-order equation to two first-order equations by defining the velocity $\frac{dx}{dt}$ as a new variable $y$:


The initial conditions are now $x=−2$ and $y=0$.


  • antimony
  • roadrunner
  • matplotlib

Notebook settings

In a first step general settings for the notebook are defined.

In [1]:
from __future__ import print_function, division
# print settings for notebook 
%matplotlib inline

import matplotlib
matplotlib.rcParams;  # available global parameters

matplotlib.rcParams['figure.figsize'] = (9.0, 6.0)
matplotlib.rcParams['axes.labelsize'] = 'medium'
font = {'family' : 'sans-serif',
        'weight' : 'normal', # bold
        'size'   : 14}
matplotlib.rc('font', **font)

SBML Model

Now the model description of the Van der Pol oscillator in a standard format for computational models, the Systems Biology Markup Language (SBML) is created [Hucka2003]. This allows to analyse the model in a multitude of tools supporting SBML, to name a few

  • cy3sbml

The SBML for the above system of ordinary differential equations (ODEs) is created with Antimony: A modular human-readable, human-writeable model definition language [Smith?].

We set the initial conditions $x=−2$ and $y=0$ and the damping parameter $mu=0$, i.e. no dampening.

In [2]:
# create SBML
import antimony

model_id = 'van_der_pol'
# ----------------------------
model_str = '''
model {}
var species x = 2; 
var species y = 0;
const mu = 0;

J1: -> x; y
J2: -> y; mu *(1-x^2)*y - x

# ----------------------------

model_file = '../models/{}.xml'.format(model_id)
antimony.writeSBMLFile(model_file, model_id)


In [3]:
# show SBML

<?xml version="1.0" encoding="UTF-8"?>
<!-- Created by libAntimony version v2.8.1 on 2016-01-15 10:14 with libSBML version 5.12.1. -->
<sbml xmlns="" level="3" version="1">
  <model id="van_der_pol" name="van_der_pol">
      <unitDefinition id="dimensionless_litre">
          <unit kind="dimensionless" exponent="1" scale="0" multiplier="1"/>
          <unit kind="litre" exponent="1" scale="0" multiplier="1"/>
      <compartment sboTerm="SBO:0000410" id="default_compartment" spatialDimensions="3" size="1" constant="true"/>
      <species id="x" compartment="default_compartment" initialConcentration="2" substanceUnits="dimensionless_litre" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
      <species id="y" compartment="default_compartment" initialConcentration="0" substanceUnits="dimensionless_litre" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
      <parameter id="mu" value="0" units="dimensionless" constant="true"/>
      <reaction id="J1" reversible="true" fast="false">
          <speciesReference species="x" stoichiometry="1" constant="true"/>
          <modifierSpeciesReference species="y"/>
          <math xmlns="">
            <ci> y </ci>
      <reaction id="J2" reversible="true" fast="false">
          <speciesReference species="y" stoichiometry="1" constant="true"/>
          <modifierSpeciesReference species="x"/>
          <math xmlns="" xmlns:sbml="">
                <ci> mu </ci>
                  <cn sbml:units="dimensionless" type="integer"> 1 </cn>
                    <ci> x </ci>
                    <cn sbml:units="dimensionless" type="integer"> 2 </cn>
                <ci> y </ci>
              <ci> x </ci>


Validation of the SBML model and display of problems with the file.

In [4]:
import multiscale.sbmlutils.validation as validation
validation.validate_sbml(model_file, ucheck=False); # no unit checks

 filename : ../models/van_der_pol.xml
 file size (byte) : 2856
 read time (ms) : 0.938177
 c-check time (ms) : 0.416994
 validation error(s) : 0
 consistency error(s): 0
 validation warning(s) : 0
 consistency warning(s): 0 


After the model definition the model can be be simulated. We can now use the defined model for some simulations. We use roadrunner for the simulations.

In [5]:
import roadrunner
import multiscale.odesim.simulate.roadrunner_tools as rt
# load model in roadrunner
r = roadrunner.RoadRunner(model_file)
# what are the current selections
# make the integrator settings

['time', '[x]', '[y]']
<roadrunner.RoadRunner() { 
'this' : 0x426d8b0
'modelLoaded' : true
'modelName' : van_der_pol
'libSBMLVersion' : LibSBML Version: 5.11.0
'jacobianStepSize' : 1e-05
'conservedMoietyAnalysis' : false
'simulateOptions' : 
< roadrunner.SimulateOptions() 
'this' : 0x42dfe10, 
'reset' : 0,
'structuredResult' : 0,
'copyResult' : 1,
'steps' : 50,
'start' : 0,
'duration' : 5
'integrator' : 
< roadrunner.Integrator() >
      relative_tolerance: 0.00000001
      absolute_tolerance: 0.00000001
                   stiff: true
       maximum_bdf_order: 5
     maximum_adams_order: 12
       maximum_num_steps: 20000
       maximum_time_step: 0
       minimum_time_step: 0
       initial_time_step: 0
          multiple_steps: false
      variable_step_size: true


In [6]:
# Simulation 
print('mu = {}'.format(r['mu']))
tend = 50
s = r.simulate(0, tend)

# create plot
import matplotlib.pylab as plt
plt.plot(s['time'], s['[x]'], color='blue', label='x')
plt.plot(s['time'], s['[y]'], color='black', label='y')
plt.ylabel('x , y');

plt.plot(s['[x]'], s['[y]'], color="black")

mu = 0.0

Model behavior

Evolution of the limit cycle in the phase plane. The limit cycle begins as circle and, with varying $\mu$, become increasingly sharp. An example of a Relaxation oscillator.

The Van der Pol oscillator shows an interesting behavior depending on the dampin parameter $\mu$.

In [7]:
# add the additional values of interest to the selection
r.selections = ['time'] + ['[{}]'.format(sid) for sid in r.model.getFloatingSpeciesIds()] \
                            + r.model.getReactionIds()

['time', '[x]', '[y]']
['time', '[x]', '[y]', 'J1', 'J2']

In [8]:
# change the control parameter mu in the model
import numpy as np
results = []
mu_values = [0, 0.01, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0]
for mu in mu_values:
    s, gp = rt.simulate(r, t_start=0, t_stop=100, parameters={'mu': mu})

Integration time: 0.0222599506378
Integration time: 0.015722990036
Integration time: 0.0341429710388
Integration time: 0.0447070598602
Integration time: 0.048525094986
Integration time: 0.0532228946686
Integration time: 0.0546460151672
Integration time: 0.0552170276642
Integration time: 0.0539009571075
Integration time: 0.0552830696106

In [9]:
for k, mu in enumerate(mu_values):
    res = results[k]
    plt.plot(res['[x]'], res['[y]'], color='black')
plt.title('Phase plane of limit cycle')

(-7, 7)

In [10]:
# calculating the derivatives for given concentrations
def dxdt(rr, X, Y): 
    DX = np.zeros_like(X)
    DY = np.zeros_like(Y)
    for k, _ in np.ndenumerate(X):
        # print('X[k], Y[k]', X[k], Y[k])
        rr['[x]'], rr['[y]'] = X[k], Y[k] 
        DX[k], DY[k] = rr['J1'], rr['J2']
    return DX, DY

In [11]:
def phase_portrait(r, x=np.linspace(-6, 6, 20), y=np.linspace(-8, 8, 20), figsize = (5,8)):
    fig2 = plt.figure(figsize=figsize)
    ax2 = fig2.add_subplot(1,1,1)

    # quiverplot
    # define a grid and compute direction at each point
    X1 , Y1  = np.meshgrid(x, y)                    # create a grid
    DX1, DY1 = dxdt(r, X1, Y1)                   # compute J1 and J2 (use roadrunner as calculator)

    M = (np.hypot(DX1, DY1))                        # norm the rate 
    M[ M == 0] = 1.                                 # avoid zero division errors 
    DX1 /= M                                        # normalize each arrows
    DY1 /= M

    ax2.quiver(X1, Y1, DX1, DY1, M, pivot='mid')
    # ax2.xaxis.label = 'x'
    # ax2.yaxis.label = 'y'
    # ax2.legend()

In [12]:
# single trajectories
def trajectory(mu, x0=2.0, y0=0.0, color="black"):
    s, _ = rt.simulate(r, t_start=0, t_stop=100, 
                       parameters = {'mu' : mu}, init_concentrations={'x':x0, 'y':y0})
    plt.plot(s['[x]'], s['[y]'], color=color)

In [13]:
for mu in [0.3, 0.0, 0.3, 4.0]:
    r['mu'] = mu
    phase_portrait(r, figsize=(8,8))
    trajectory(mu, x0=2.3, y0=4.0)

Integration time: 0.0303640365601
Integration time: 0.0159130096436
Integration time: 0.0311679840088
Integration time: 0.0522689819336

Phase plane and trajectories

We can now analyse the phase plane and the trajectories for different damping values $mu$.

In [14]:
mu = 0.3  # 0.0, 0.3, 4.0

# create phase portrait
r['mu'] = mu

trajectory(mu, x0=2.3, y0=4.0)
trajectory(mu, x0=0.2, y0=0.4)
trajectory(mu, x0=-2.0, y0=6.0)
trajectory(mu, x0=-2.0, y0=-6.0)
# limit cycle
trajectory(mu, color="blue")

Integration time: 0.0304691791534
Integration time: 0.0283119678497
Integration time: 0.0305590629578
Integration time: 0.030464887619
Integration time: 0.0297479629517

In [ ]: