# 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¶

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
└─ site-packages


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 site-packages:

site-packages
└─ 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=1e-05, 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 is None, 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. If dimension is 1, rhs can also be a string. For a system of two delay-coupled Mackey-Glass oscillators rhs 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 names p, x and xd are reservered for the parameter array, the system state array and the delayed state array. The values p[0], p[1], … must be provides by the below argument rhs_parameters. The variables x[0], x[1], … are the first, second, … component of the current systems state. The variable xd[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 delay-coupled Mackey-Glass 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 Mackey-Glass 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 Mackey-Glass 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 is None, t_start will be used as the first time of the solution. If time_grid is not None, t_start does not matter. t_end (float, int or None) – If time_grid is None, t_end will be used as the last time of the solution. If time_grid is not None, 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. solution (numpy.ndarray) – The solution values and time points are stored here after running the method solve(). The time points are stored in the 0-th columns, and the i-th component of system states are stored in the i-th column of the array.
solve()

This method integrates the DDE. The solution will be stored in the attribute DDE.solution.

## Example¶

You may use pyoddtool to simulate two delay-coupled Mackey-Glass oscillators

\begin{align}\begin{aligned}\dot{x}_1 &= \frac{\sigma x_2(t-\tau )}{1 + x_2(t-\tau )^n} - \gamma x_1(t),\\\dot{x}_2 &= \frac{\sigma x_1(t-\tau )}{1 + x_1(t-\tau )^n} - \gamma x_2(t).\end{aligned}\end{align}

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 TU-Berlin 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 Mackey-Glass oscillators: interpolated solution")
plt.show()

plt.plot(solution_interpolated[:, 1], solution_interpolated[:, 2])
plt.title("Two coupled Mackey-Glass 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 Mackey-Glass oscillators: of pure solution")
plt.show()

plt.plot(solution_pure[:, 1], solution_pure[:, 2])
plt.title("Two coupled Mackey-Glass 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 Mackey-Glass oscillators which are delay-coupled 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" + \
"}"


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 TU-Berlin 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 Mackey-Glass oscillators: interpolated solution")
plt.show()

plt.plot(solution_interpolated[:, 1], solution_interpolated[:, 2])
plt.title("Two coupled Mackey-Glass 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 Mackey-Glass oscillators: of pure solution")
plt.show()

plt.plot(solution_pure[:, 1], solution_pure[:, 2])
plt.title("Two coupled Mackey-Glass oscillators: phase portrait pure solution")
plt.show()