Welcome to pyoddtool’s documentation!¶
The Python package pyoddtool provides an Python interface for the C++ programm oddtool (which is included in this package).
You can find the oddtool project page here: https://github.com/leoluecken/oddtool
See the oddtool documentation for information about the used numerical method.
Requirements¶
This package was tested under Ubuntu and openSUSE. It requires a C++ compiler (e.g. gcc 4.9 or higher) and the C++ boost library. You will also need numpy and scipy.
Installing pyoddtool¶
Download and extract the file pyoddtool_v0.0.2.tar.gz.
You will find the following directory structure:
pyoddtool_v0.0.2
├─ setup.py
└─ pyoddtool
├─ __init__.py
└─ oddtool
Open a terminal in the pyoddtool_v0.0.2
directory and run
python setup.py install
That’s it.
If you have installed both Python 2.7 and Python 3 on your computer, then the command for Python 3 is probably python3
. That is you need to run
python3 setup.py install
if you want to be able to use pyoddtool with Python 3.
Installing pyoddtool on a TU Berlin computer¶
If you try to install the package on a computer at the Institute of Mathematics of the TU Berlin, there probably occurs an error (permission denied). In order to use the package in this case, you need to do the following:
Go to the .local
directory in your home directory. (Enable “show hidden files” if you cannot see the .local
directory.) Here you create the directories
lib
└─ python3.4
└─ sitepackages
if they do not exist yet. (Check the version number of your Python and use this number instead of 3.4 if you have a different version.) Now, copy the pyoddtool directory into the folder sitepackages
:
sitepackages
└─ pyoddtool
├─ __init__.py
└─ oddtool
Now you should be able to use pyoddtool with Python 3.
The class DDE¶

class
pyoddtool.
DDE
(time_grid, dimension, rhs, rhs_parameters, delays, history, additional_lines=None, t_start=None, t_end=None, h_max=1, h_min=1e05, terminal_output=False)¶ The class DDE uses the C++ software “oddtool” to solve delay differential equations.
Parameters:  time_grid (numpy.ndarray or None) – An array which contains the time points for which the states
of the system should be stored. If
time_grid
isNone
, the solution states at the time points chosen by oddtool’s adaptive stepper will be stored.  dimension (int) – The dimension of the system.
 rhs (list of strings or single string) – Contains the right hand sides
of the system (i.e. the dynamics) written in C++ code. The list
must have the length
dimension
. Ifdimension
is1
,rhs
can also be a string. For a system of two delaycoupled MackeyGlass oscillatorsrhs
would be['p[0] * xd[0][1] / (1 + pow(xd[0][1], p[1]))  p[2] * x[0]', 'p[0] * xd[0][0] / (1 + pow(xd[0][0], p[1]))  p[2] * x[1]']
. The C++ variable namesp
,x
andxd
are reservered for the parameter array, the system state array and the delayed state array. The valuesp[0]
,p[1]
, … must be provides by the below argument rhs_parameters. The variablesx[0]
,x[1]
, … are the first, second, … component of the current systems state. The variablexd[j][i]
denotes the \(i\)th component of the system state at time \(t  \tau _j\).  rhs_parameters (numpy.ndarray of numpy floats or numpy ints) – An array which contains the parameters for the right hand sides.
For the example with two delaycoupled MackeyGlass oscillators this array could
look as follows:
array([2.2, 10.0, 1.0])
.  delays (numpy.ndarray of numpy floats or numpy ints or single float or int) – An array which contains the constant delay times.
If there is only one delay time, this parameter can also be a float.
For the example with two coupled MackeyGlass oscillators this parameter could
look as foolows:
array([20])
. Since there is only one delay time for this system, the parameter can also be a simple float or int:20
.  history (list of strings or single string) – Contains history terms written in C++ code. The history term are functions
for the system states at and before the initial time point.
The list must have the length
dimension
. For the example with two coupled MackeyGlass oscillators this list could look as follows:['abs(cos(t)) + 1', '1 + sin(t)']
.  additional_lines (None or string) – If a string is provided, it will be inserted in the rhs function in oddtool’s C++ main file. (See example 2 below.)
 t_start (float, int or None) – If
time_grid
isNone
,t_start
will be used as the first time of the solution. Iftime_grid
is notNone
,t_start
does not matter.  t_end (float, int or None) – If
time_grid
isNone
,t_end
will be used as the last time of the solution. Iftime_grid
is notNone
,t_end
does not matter.  h_max (float or int) – Maximal time step for the integration. It should be
smaller than the distance of time points in the parameter
time_grid
.  h_min (float or int) – Minimal time step for the integration. It should be
smaller than or equal to
h_max
.  terminal_output (bool) – If
True
, oddtool’s terminal output will be printed.
Attributes:  solution (numpy.ndarray) – The solution values and time points are stored here
after running the method
solve()
. The time points are stored in the 0th columns, and the ith component of system states are stored in the ith column of the array.

solve
()¶ This method integrates the DDE. The solution will be stored in the attribute
DDE.solution
.
 time_grid (numpy.ndarray or None) – An array which contains the time points for which the states
of the system should be stored. If
Example¶
You may use pyoddtool to simulate two delaycoupled MackeyGlass oscillators
In this example we use \(\sigma = 2.2, n = 10\) and \(\gamma = 1\) as parameters.
Open a Python console and import pyoddtool:
>>> import pyoddtool
We want to get the solution at each time point \(0, 0.1, 0.2, \ldots , 99.9, 100\). Thus we create a numpy array with these time points:
>>> import numpy as np
>>> t = np.arange(0, 100.1, 0.1)
The right hand side of the system equation must be expressed as a list of strings containing C++ code:
>>> cpp_formulas = ['p[0] * xd[0][1] / (1 + pow(xd[0][1], p[1]))  p[2] * x[0]',
... 'p[0] * xd[0][0] / (1 + pow(xd[0][0], p[1]))  p[2] * x[1]']
In the above code snippet x[0]
and x[1]
denote the system states \(x_1(t)\) and \(x_2(t)\), and xd[0][0]
and xd[0][1]
are they delayed states \(x_1(t\tau)\) and \(x_2(t\tau)\). The parameters \(\sigma, n\) and \(\gamma\) are represented by p[0]
, p[1]
and p[2]
. The parameter values must be provided in a numpy.array:
>>> parameters = np.array([2.2, 10.0, 1.0])
Moreover, we need to provide a function which describes the system states at the initial time point and before. This function must be written in C++ too. Again, we need a list of strings:
>>> history_functions = ['abs(cos(t)) + 1', '1 + sin(t)']
Now, we can initiate an instance of pyoddtool.DDE
. Note that we also need to provide the system’s dimension and the delay:
>>> dde = pyoddtool.DDE(time_grid=t, dimension=2, rhs=cpp_formulas,
... rhs_parameters=parameters, delays=20, history=history_functions)
We call the method solve()
to do the computation:
>>> dde.solve()
Now, we can plot the solution:
>>> import matplotlib.pyplot as plt
>>> x1 = dde.solution[:, 1]
>>> x2 = dde.solution[:, 2]
>>> plt.plot(t, x1, title='Solution x_1')
>>> plt.show()
>>> plt.plot(t, x2, title='Solution x_2')
>>> plt.show()
>>> plt.plot(x1, x2, title='Phase portrait')
>>> plt.show()
In this case dde.solution
contains the values at the given time points that we obtained from a cubic spline interpolation of the “original” solution which was recieved from oddtool.
The parameters time_grid
, rhs_parameters
, delays
, history
, t_start
, t_end
, h_max
and h_min
can be accessed directly. That is we can for example change the time parameters and the maximal step size and compute another solution using the new parameters:
>>> dde.time_grid = None
>>> dde.t_start = 0
>>> dde.t_end = 100
>>> dde.h_max = 0.5
>>> dde.solve()
Since dde.time_grid
is None
in this case, dde.solution
contains the system states at the time points chosen by the oddtool stepper. That is dde.solution
contains the “original” values which were returned by oddtool instead of interpolated values.
Complete Python script of the example¶
# import matplotlib
# matplotlib.rcParams["backend"] = "qt4Agg"
# uncomment the above lines if and only if you run
# the example script on a TUBerlin computer
import matplotlib.pyplot as plt
import numpy as np
import pyoddtool
t = np.arange(0,100.1,0.1)
N = 2
rhs = ['p[0] * xd[0][1] / (1 + pow(xd[0][1], p[1]))  p[2] * x[0]',
'p[0] * xd[0][0] / (1 + pow(xd[0][0], p[1]))  p[2] * x[1]']
rhs_params = np.array([2.2, 10.0, 1.0])
delays = 20
hist = ['abs(cos(t)) + 1', '1 + sin(t)']
dde = pyoddtool.DDE(t, N, rhs, rhs_params, delays, hist)
dde.solve()
solution_interpolated = dde.solution
dde.time_grid = None
dde.t_start = 0
dde.t_end = 100
dde.h_max = 0.5
dde.solve()
solution_pure = dde.solution
plt.plot(solution_interpolated[:, 0], solution_interpolated[:, 1])
plt.plot(solution_interpolated[:, 0], solution_interpolated[:, 2])
plt.title("Two coupled MackeyGlass oscillators: interpolated solution")
plt.show()
plt.plot(solution_interpolated[:, 1], solution_interpolated[:, 2])
plt.title("Two coupled MackeyGlass oscillators: phase portrait of interpolated solution")
plt.show()
plt.plot(solution_pure[:, 0], solution_pure[:, 1])
plt.plot(solution_pure[:, 0], solution_pure[:, 2])
plt.title("Two coupled MackeyGlass oscillators: of pure solution")
plt.show()
plt.plot(solution_pure[:, 1], solution_pure[:, 2])
plt.title("Two coupled MackeyGlass oscillators: phase portrait pure solution")
plt.show()
Example 2¶
Sometimes it is not enough to provide just one simple mathematical term for the right hand side. If we want to run a simulation of two MackeyGlass oscillators which are delaycoupled for \(0 \leq t < 75\) but not coupled for \(t \geq 75\), we need to use the parameter additional_lines
. We set
add_lines = "int coupling;\n" + \
"if(t < 75){\n" + \
" coupling = 1;\n" + \
"}\n" + \
"else{\n" + \
" coupling = 0;\n" + \
"}"
dde.additional_lines = add_lines
and
dde.rhs = ['p[0] * (coupling * xd[0][1] + (1  coupling) * xd[0][0]) / (1 + pow(xd[0][1], p[1]))  p[2] * x[0]',
'p[0] * (coupling * xd[0][0] + (1  coupling) * xd[0][1]) / (1 + pow(xd[0][0], p[1]))  p[2] * x[1]']
Note that C++ variables (such as coupling
in our example) have to be declared, which we did in the first line of the C++ code string add_lines
above. One can also use p[k]
, x[i]
or x[j][i]
in additional_lines
. For each time step, the additional lines will be executed before the evaluation of rhs
.
With these modifications (and t_end = 150
) the above complete python skript changes to:
# import matplotlib
# matplotlib.rcParams["backend"] = "qt4Agg"
# uncomment the above lines if and only if you run
# the example script on a TUBerlin computer
import matplotlib.pyplot as plt
import numpy as np
import pyoddtool
t = np.arange(0,150.1,0.1)
N = 2
rhs = ['p[0] * (coupling * xd[0][1] + (1  coupling) * xd[0][0]) / (1 + pow(xd[0][1], p[1]))  p[2] * x[0]',
'p[0] * (coupling * xd[0][0] + (1  coupling) * xd[0][1]) / (1 + pow(xd[0][0], p[1]))  p[2] * x[1]']
rhs_params = np.array([2.2, 10.0, 1.0])
add_lines = "int coupling;\n" + \
"if(t < 75){\n" + \
" coupling = 1;\n" + \
"}\n" + \
"else{\n" + \
" coupling = 0;\n" + \
"}"
delays = 20
hist = ['abs(cos(t)) + 1', '1 + sin(t)']
dde = pyoddtool.DDE(t, N, rhs, rhs_params, delays, hist, add_lines)
dde.solve()
solution_interpolated = dde.solution
dde.time_grid = None
dde.t_start = 0
dde.t_end = 150
dde.h_max = 0.5
dde.solve()
solution_pure = dde.solution
plt.plot(solution_interpolated[:, 0], solution_interpolated[:, 1])
plt.plot(solution_interpolated[:, 0], solution_interpolated[:, 2])
plt.title("Two coupled MackeyGlass oscillators: interpolated solution")
plt.show()
plt.plot(solution_interpolated[:, 1], solution_interpolated[:, 2])
plt.title("Two coupled MackeyGlass oscillators: phase portrait of interpolated solution")
plt.show()
plt.plot(solution_pure[:, 0], solution_pure[:, 1])
plt.plot(solution_pure[:, 0], solution_pure[:, 2])
plt.title("Two coupled MackeyGlass oscillators: of pure solution")
plt.show()
plt.plot(solution_pure[:, 1], solution_pure[:, 2])
plt.title("Two coupled MackeyGlass oscillators: phase portrait pure solution")
plt.show()