pydelay v 0.1.1

Introduction

pydelay is a program which translates a system of delay differential equations (DDEs) into simulation C-code and compiles and runs the code (using scipy weave). This way it is easy to quickly implement a system of DDEs but you still have the speed of C. The Homepage can be found here:

http://pydelay.sourceforge.net/

It is largely inspired by PyDSTool.

The algorithm used is based on the Bogacki-Shampine method [1] which is also implemented in Matlab’s dde23 [2].

We also want to mention PyDDE – a different python program for solving DDEs.

License

pydelay is licensed under the MIT License.

Citing

You can cite pydelay through the arXiv.:

@misc{FLU09a,
    title = {pydelay -- a python tool for solving delay differential equations},
    author = {Flunkert, V. and Sch{\"o}ll, E. },
    year = {2009},
    howpublished = {{\tt arXiv:0911.1633~[nlin.CD]}},
}

Installation and requirements

Unix:

You need python and python headers files (in debian/ubuntu these are in the package python-dev), numpy and scipy and the gcc-compiler.

To plot the solutions and run the examples you also need matplotlib.

To install pydelay download the latest tar.gz from the website and install the package in the usual way:

cd pydelay-$version
python setup.py install

When the package is installed, you can get some info about the functions and the usage with:

pydoc pydelay

TODO easy_install

For Arch linux there is a PKGBUILD .

Windows:

The solver has not been tested on a windows machine. It could perhaps work under cygwin.

An example

The following example shows the basic usage. It solves the Mackey-Glass equations [3] for initial conditions which lead to a periodic orbit (see [4] for this example).

# import pydelay and numpy and pylab
import numpy as np
import pylab as pl
from pydelay import dde23

# define the equations
eqns = {
    'x' : '0.25 * x(t-tau) / (1.0 + pow(x(t-tau),p)) -0.1*x'
    }

#define the parameters
params = {
    'tau': 15,
    'p'  : 10
    }

# Initialise the solver
dde = dde23(eqns=eqns, params=params)

# set the simulation parameters
# (solve from t=0 to t=1000 and limit the maximum step size to 1.0)
dde.set_sim_params(tfinal=1000, dtmax=1.0)

# set the history of to the constant function 0.5 (using a python lambda function)
histfunc = {
    'x': lambda t: 0.5
    }
dde.hist_from_funcs(histfunc, 51)

# run the simulator
dde.run()

# Make a plot of x(t) vs x(t-tau):
# Sample the solution twice with a stepsize of dt=0.1:

# once in the interval [515, 1000]
sol1 = dde.sample(515, 1000, 0.1)
x1 = sol1['x']

# and once between [500, 1000-15]
sol2 = dde.sample(500, 1000-15, 0.1)
x2 = sol2['x']

pl.plot(x1, x2)
pl.xlabel('$x(t)$')
pl.ylabel('$x(t - 15)$')
pl.show()

[source code, hires.png, pdf]

Usage

Defining the equations, delays and parameters

Equations are defined using a python dictionary. The keys are the variable names and the entry is the right hand side of the differential equation. The string defining the equation has to be a valid C expression, i.e., use pow(a,b) instead of a**b etc.

Delays are written as (t-delay), where delay can be some expression involving parameters and numbers but not (yet) involving the time t or the dynamic variables:

eqns = {
    'y1': '- y1 * y2(t-tau) + y2(t-1.0)',
    'y2': 'a * y1 * y2(t-2*tau) - y2',
    'y3': 'y2 - y2(t-(tau+1))'
  }

Complex variables can be defined by adding ':c' or ':C' in the eqn-dictionary. The imaginary unit can be used through 'ii' in the equations:

eqns = {
    'z:c': '(la + ii*w0 + g*pow(abs(z),2) )*z + b*(z(t-tau) - z(t))',
}

Parameters are defined in a separate dictionary where the keys are the parameter names, i.e.,:

params = {
    'a'  : 0.2,
    'tau': 1.0
}

Setting the history

The history of the variables is stored in the dictionary dde23.hist. The keys are the variable names and there is an additional key 't' for the time array of the history.

There is a second dictionary dde23.Vhist where the time derivatives of the history is stored (this is needed for the solver). When the solver is initialized, i.e.,:

dde = dde23(eqns, params)

the history of all variables (defined in eqns) is initialized to an array of length nn=101 filled with zeros. The time array is evenly spaced in the interval [-maxdelay, 0].

It is possible to manipulate these arrays directly, however this is not recommended since one easily ends up with an ill-defined history resulting for example in segfaults or false results.

Instead use the following methods to set the history.

dde23.hist_from_funcs(dic, nn=101)

Initialise the histories with the functions stored in the dictionary dic. The keys are the variable names. The function will be called as f(t) for t in [-maxdelay, 0] on nn samples in the interval.

This function provides the simplest way to set the history. It is often convenient to use python lambda functions for f. This way you can define the history function in place.

If any variable names are missing in the dictionaries, the history of these variables is set to zero and a warning is printed. If the dictionary contains keys not matching any variables these entries are ignored and a warning is printed.

Example: Initialise the history of the variables x and y with cos and sin functions using a finer sampling resolution:

from math import sin, cos

histdic = {
    'x': lambda t: cos(0.2*t),
    'y': lambda t: sin(0.2*t)
}

dde.hist_from_funcs(histdic, 500)

dde23.hist_from_arrays(dic, useend=True)

Initialise the history using a dictionary of arrays with variable names as keys. Additionally a time array can be given corresponding to the key t. All arrays in dic have to have the same lengths.

If an array for t is given the history is interpreted as points (t,var). Otherwise the arrays will be evenly spaced out over the interval [-maxdelay, 0].

If useend is True the time array is shifted such that the end time is zero. This is useful if you want to use the result of a previous simulation as the history.

If any variable names are missing in the dictionaries, the history of these variables is set to zero and a warning is printed. If the dictionary contains keys not matching any variables (or 't') these entries are ignored and a warning is printed.

Example::

t = numpy.linspace(0, 1, 500)
x = numpy.cos(0.2*t)
y = numpy.sin(0.2*t)

histdic = {
    't': t,
    'x': x,
    'y': y
}
dde.hist_from_arrays(histdic)

Note that the previously used methods hist_from_dict, hist_from_array and hist_from_func (the last two without s) have been removed, since it was too easy to make mistakes with them.

The solution

After the solver has run, the solution (including the history) is stored in the dictionary dde23.sol. The keys are again the variable names and the time 't'. Since the solver uses an adaptive step size method, the solution is not sampled at regular times.

To sample the solutions at regular (or other custom spaced) times there are two functions.

dde23.sample(tstart=None, tfinal=None, dt=None)

Sample the solution with dt steps between tstart and tfinal.

tstart, tfinal

Start and end value of the interval to sample. If nothing is specified tstart is set to zero and tfinal is set to the simulation end time.

dt

Sampling size used. If nothing is specified a reasonable value is calculated.

Returns a dictionary with the sampled arrays. The keys are the variable names. The key 't' corresponds to the sampling times.

dde23.sol_spl(t)

Sample the solutions at times t.

t

Array of time points on which to sample the solution.

Returns a dictionary with the sampled arrays. The keys are the variable names. The key 't' corresponds to the sampling times.

These functions use a cubic spline interpolation of the solution data.

Noise

Noise can be included in the simulations. Note however, that the method used is quite crude (an Euler method will be added which is better suited for noise dominated dynamics). The deterministic terms are calculated with the usual Runge-Kutta method and then the noise term is added with the proper scaling of $\sqrt{dt}$ at the final step. To get accurate results one should use small time steps, i.e., dtmax should be set small enough.

The noise is defined in a separate dictionary. The function gwn() can be accessed in the noise string and is a Gaussian white noise term of unit variance. The following code specifies an Ornstein-Uhlenbeck process.:

eqns = { 'x': '-x' }
noise = { 'x': 'D * gwn()'}
params = { 'D': 0.00001 }

dde = dde23(eqns=eqns, params=params, noise=noise)

You can also use noise terms of other forms by specifying an appropriate C-function (see the section on custom C-code).

Custom C-code

You can access custom C-functions in your equations by adding the definition as supportcode for the solver. In the following example a function f(w,t) is defined through C-code and accessed in the eqn string.:

# define the eqn f is the C-function defined below
eqns = { 'x': '- x + k*x(t-tau) + A*f(w,t)' }
params = {
    'k'  : 0.1,
    'w'  : 2.0,
    'A'  : 0.5,
    'tau': 10.0
}

mycode = """
double f(double t, double w) {
    return sin(w * t);
}
"""

dde = dde23(eqns=eqns, params=params, supportcode=mycode)

When defining custom code you have to be careful with the types. The type of complex variables in the C-code is Complex. Note in the above example that w has to be given as an input to the function, because the parameters can only be accessed from the eqns string and not inside the supportcode. (Should this be changed?)

Using custom C-code is often useful for switching terms on and off. For example the Heaviside function may be defined and used as follows.:

# define the eqn f is the C-function defined below
eqns = { 'z:c': '(la+ii*w)*z - Heavi(t-t0)* K*(z-z(t-tau))' }
params = {
    'K'  : 0.1 ,
    'w'  : 1.0,
    'la' : 0.1,
    'tau': pi,
    't0' : 2*pi
}

mycode = """
double Heavi(double t) {
    if(t>=0)
        return 1.0;
    else
        return 0.0;
}
"""
dde = dde23(eqns=eqns, params=params, supportcode=mycode)

This code would switch a control term on when t>t0. Note that Heavi(t-t0) does not get translated to a delay term, because Heavi is not a system variable.

Since this scenario occurs so frequent the Heaviside function (as defined above) is included by default in the source code.

Use and modify generated code

The compilation of the generated code is done with scipy.weave. Instead of using weave to run the code you can directly access the generated code via the function dde23.output_ccode(). This function returns the generated code as a string which you can then store in a source file.

To run the generated code manually you have to set the precompiler flag\ #define MANUAL (uncomment the line in the source file) to exclude the python related parts and include some other parts making the code a valid stand alone source file. After this the code should compile with g++ -lm -o prog source.cpp and you can run the program manually.

You can specify the history of all variables in the source file by setting the for loops after the comment\ /* set the history here ...  */.

Running the code manually can help you debug, if some problem occurs and also allows you to extend the code easily.

Another example

The following example shows some of the things discussed above. The code simulates the Lang-Kobayashi laser equations [5]

import numpy as np
import pylab as pl
from pydelay import dde23

tfinal = 10000
tau = 1000

#the laser equations
eqns = {
    'E:c': '0.5*(1.0+ii*a)*E*n + K*E(t-tau)',
    'n'  : '(p - n - (1.0 +n) * pow(abs(E),2))/T'
}

params = {
    'a'  : 4.0,
    'p'  : 1.0,
    'T'  : 200.0,
    'K'  : 0.1,
    'tau': tau,
    'nu' : 10**-5,
    'n0' : 10.0
}

noise = { 'E': 'sqrt(0.5*nu*(n+n0)) * (gwn() + ii*gwn())' }

dde = dde23(eqns=eqns, params=params, noise=noise)
dde.set_sim_params(tfinal=tfinal)

# use a dictionary to set the history
thist = np.linspace(0, tau, tfinal)
Ehist = np.zeros(len(thist))+1.0
nhist = np.zeros(len(thist))-0.2
dic = {'t' : thist, 'E': Ehist, 'n': nhist}

# 'useend' is True by default in hist_from_dict and thus the
# time array is shifted correctly
dde.hist_from_arrays(dic)

dde.run()

t = dde.sol['t']
E = dde.sol['E']
n = dde.sol['n']

spl = dde.sample(-tau, tfinal, 0.1)

pl.plot(t[:-1], t[1:] - t[:-1], '0.8', label='step size')
pl.plot(spl['t'], abs(spl['E']), 'g', label='sampled solution')
pl.plot(t, abs(E), '.', label='calculated points')
pl.legend()

pl.xlabel('$t$')
pl.ylabel('$|E|$')

pl.xlim((0.95*tfinal, tfinal))
pl.ylim((0,3))
pl.show()

[source code, hires.png, pdf]

Module Reference

dde23.__init__(eqns, params=None, noise=None, supportcode='', debug=False)

Initialise the solver.

eqns

Dictionary defining for each variable the derivative. Delays are written as as (t-...) example:

eqns = {
    'y1': '- y1 * y2(t-tau1) + y2(t-tau2)',
    'y2': 'a * y1 * y2(t-tau1) - y2',
    'y3': 'y2 - y2(t-tau2)'
    }

You can also directly use numbers or combination of parameters as delays:

eqns = {
    'x1': '-a*x1 + x1(t - 1.0)',
    'x2': 'x2-b*x1(t-2.0*a+b)
    }

At the moment only constant delays are supported.

The string defining the equation has to be a valid C expression, i.e., use pow(a,b) instead of a**b etc. (this might change in the future):

eqns = {'y': '-2.0 * sin(t) * pow(y(t-tau), 2)'}

Complex variable can be defined using :C or :c in the variable name. The imaginary unit can be used through ii in the equations:

eqns = {'z:C': '(-la + ii * w0) * z' }

params

Dictionary defining the parameters (including delays) used in eqns. example:

params = {
   'a'   : 1.0, 
   'tau1': 1.0, 
   'tau2': 10.0
   }

noise

Dictionary for noise terms. The function gwn() can be accessed in the noise string and provides a Gaussian white noise term of unit variance. example:

noise = {'x': '0.01*gwn()'}

debug

If set to True the solver gives verbose output to stdout while running.

dde23.set_sim_params(tfinal=100, AbsTol=9.9999999999999995e-07, RelTol=0.001, dtmin=9.9999999999999995e-07, dtmax=None, dt0=None, MaxIter=1000000000.0)

tfinal

End time of the simulation (the simulation always starts at t=0).

AbsTol, RelTol

The relative and absolute error tolerance. If the estimated error e for a variable y obeys e <= AbsTol + RelTol*|y| then the step is accepted. Otherwise the step will be repeated with a smaller step size.

dtmin, dtmax

Minimum and maximum step size used.

dt0

initial step size

MaxIter

maximum number of steps. The simulation stops if this is reached.

dde23.hist_from_arrays(dic, useend=True)

Initialise the history using a dictionary of arrays with variable names as keys. Additionally a time array can be given corresponding to the key t. All arrays in dic have to have the same lengths.

If an array for t is given the history is interpreted as points (t,var). Otherwise the arrays will be evenly spaced out over the interval [-maxdelay, 0].

If useend is True the time array is shifted such that the end time is zero. This is useful if you want to use the result of a previous simulation as the history.

If any variable names are missing in the dictionaries, the history of these variables is set to zero and a warning is printed. If the dictionary contains keys not matching any variables (or 't') these entries are ignored and a warning is printed.

Example::

t = numpy.linspace(0, 1, 500)
x = numpy.cos(0.2*t)
y = numpy.sin(0.2*t)

histdic = {
    't': t,
    'x': x,
    'y': y
}
dde.hist_from_arrays(histdic)

dde23.hist_from_funcs(dic, nn=101)

Initialise the histories with the functions stored in the dictionary dic. The keys are the variable names. The function will be called as f(t) for t in [-maxdelay, 0] on nn samples in the interval.

This function provides the simplest way to set the history. It is often convenient to use python lambda functions for f. This way you can define the history function in place.

If any variable names are missing in the dictionaries, the history of these variables is set to zero and a warning is printed. If the dictionary contains keys not matching any variables these entries are ignored and a warning is printed.

Example: Initialise the history of the variables x and y with cos and sin functions using a finer sampling resolution:

from math import sin, cos

histdic = {
    'x': lambda t: cos(0.2*t),
    'y': lambda t: sin(0.2*t)
}

dde.hist_from_funcs(histdic, 500)

dde23.output_ccode()

dde23.run()

run the simulation

class pydelay.dde23(eqns, params=None, noise=None, supportcode='', debug=False)

This class translates a DDE to C and solves it using the Bogacki-Shampine method.

Attributes of class instances:

For user relevant attributes:

self.sol

Dictionary storing the solution (when the simulation has finished). The keys are the variable names and 't' corresponding to the sampled times.

self.discont

List of discontinuity times. This is generated from the occurring delays by propagating the discontinuity at t=0. The solver will step on these discontinuities. If you want the solver to step onto certain times they can be inserted here.

self.rseed

Can be set to initialise the random number generator with a specific seed. If not set it is initialised with the time.

self.hist

Dictionary with the history. Don’t manipulate the history arrays directly! Use the provided functions to set the history.

self.Vhist

Dictionary with the time derivatives of the history.

For user less relevant attributes:

self.delays

List of the delays occurring in the equations.

self.chunk

When arrays become to small they are grown by this number.

self.spline_tck

Dictionary which stores the tck spline representation of the solutions. (see scipy.interpolate)

self.eqns

Stores the eqn dictionary.

self.params

Stores the parameter dictionary.

self.simul

Dictionary of the simulation parameters.

self.noise

Stores the noise dictionary.

self.debug

Stores the debug flag.

self.delayhashs

List of hashs for each delay (this is used in the generated C-code).

self.vars

List of variables extracted from the eqn dictionary keys.

self.types

Dictionary of C-type names of each variable.

self.nptypes

Dictionary of numpy-type names of each variable.

[1]	Bogacki, P. and Shampine, L. F., A 3(2) pair of Runge - Kutta formulas, Applied Mathematics Letters 2, 4, 321 ISSN 0893-9659, (1989).

[2]	Shampine, L. F. and Thompson, S., Solving DDEs in Matlab, Appl. Num. Math. 37, 4, 441 ( 2001)

[3]	Mackey, M. C. and Glass, L. (1977). Pathological physiological conditions resulting from instabilities in physiological control system. Science, 197(4300):287-289.

[4]	http://www.scholarpedia.org/article/Mackey-Glass_equation

[5]	Lang, R. and Kobayashi, K. , External optical feedback effects on semiconductor injection laser properties, IEEE J. Quantum Electron. 16, 347 (1980)