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]}},
}
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.
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]
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
}
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.
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)
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.
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.
Sample the solution with dt steps between tstart and tfinal.
Returns a dictionary with the sampled arrays. The keys are the variable names. The key 't' corresponds to the sampling times.
Sample the solutions at times t.
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 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).
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.
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.
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]
Initialise the solver.
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' }
Dictionary defining the parameters (including delays) used in eqns. example:
params = {
'a' : 1.0,
'tau1': 1.0,
'tau2': 10.0
}
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()'}
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)
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)
This class translates a DDE to C and solves it using the Bogacki-Shampine method.
Attributes of class instances:
For user relevant attributes:
For user less relevant attributes:
[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) |