3.2. Analysis of Delayed Differential Equations

In this tutorial, you will learn how to implement and solve delayed differential equation (DDE) systems via PyRates. DDEs take the general form

\[\dot{\mathbf{y}} = \mathbf{f}(\mathbf{y}, t, \mathbf{y_{ au}})\]

with state vector \(\mathbf{y}\), vector field \(f\), time point \(t\) and state history \(\mathbf{y_{ au}} = {\mathbf{y}(\tau) : \tau \leq t}\). Here, we will address the case of implementing and solving a DDE system with discrete delays of the form

\[\dot{\mathbf{y}} = \mathbf{f}(t, \mathbf{y}(t), \mathbf{y}(t-\tau)),\]

with discrete delay \(\tau\). As a specific example, we will use a modified version of the Van der Pol oscillator with additional delayed feedback:

\[ \begin{align}\begin{aligned}\dot x(t) &= z(t),\\\dot z(t) &= \mu (1-x(t)^2) z - x(t) + k x(t-\tau),\end{aligned}\end{align} \]

with delayed feedback strength k. For a detailed introduction of the Van der Pol oscillator model, see the gallery example in the model introduction section.

Below, we will show how such a DDE model can be implemented and how it can be solved via the third-party DDE solver ddeint.

import numpy as np

from pyrates import CircuitTemplate, OperatorTemplate, clear
import matplotlib.pyplot as plt

3.2.1. Step 1: DDE definition

DDEs can be defined in two different ways in PyRates. The first way would be the usage of the past(y, tau) function call that indicates to PyRates that \(y(t-tau)\) should be evaluated instead of \(y(t)\). An example of how to implement an OperatorTemplate representing the delayed Van der Pol equations is provided below.

# define parameters
k = -2.0
tau = 5.0

# define operator
eqs = [
    "x' = z",
    "z' = mu*(1-x**2)*z - x + k*past(x,tau)"
]
variables = {
    'x': 'output(0.0)',
    'z': 'variable(1.0)',
    'mu': 1.0,
    'k': k,
    'tau': tau
}
op = OperatorTemplate(name='vdp_delayed', equations=eqs, variables=variables)

This OperatorTemplate could then be used to define a PyRates model (see the examples in the model introduction section of the gallery for a tutorial on how to do that).

Alternatively, it is also possible to add an edge to a CircuitTemplate, including a discrete delay. This edge will then automatically translated into a corresponding past call by PyRates. Note that this only works if the edge source variable is a variable defined by a differential equation. So in the case of the Van der Pol model, both \(x\) and \(z\) would be valid variables for such an edge definition. Below, we show how to do that using the Van der Pol model that is already implemented in PyRates:

vdp = CircuitTemplate.from_yaml("model_templates.oscillators.vanderpol.vdp")
vdp.update_template(edges=[('p/vdp_op/x', 'p/vdp_op/inp', None, {'weight': k, 'delay': tau})], in_place=True)

3.2.2. Step 2: DDE simulation

We can perform a simple simulation using a forward Euler algorithm to solve the DDE as follows:

# define simulation time
T = 100.0
step_size = 1e-2

# solve DDE via forward Euler
results = vdp.run(step_size=step_size, simulation_time=T, outputs={'x': 'p/vdp_op/x'}, solver='euler',
                  in_place=False, clear=True)

# plot resulting signal
plt.plot(results)
plt.show()

To use a different solver, for example a Runge-Kutta algorithm with automatic step-size adaptation, you would either have to use backend='julia' and solver='julia_dde' as options to the run method, or you could generate the run function instead and use any DDE solver available in the backend of your choice. Below, we demonstrate how to solve the DDE via the ddeint Python package. The latter is not part of the PyRates prerequisites and thus has to be manually installed via pip install ddeint.

from ddeint import ddeint

# generate DDE function
func, args, _, _ = vdp.get_run_func(func_name='vdp_dde_run', file_name='vdp_dde', solver='scipy', step_size=step_size,
                                    vectorize=False, backend='default', in_place=False)

# define ddeint wrapper function
t, y, hist, *fargs = args
def dde_run(y, t):
    return func(t, y(t), y, *fargs)

# solve DDE via ddeint
eval_time = np.linspace(0, T, num=int(T/step_size))
res2 = ddeint(func=dde_run, g=hist, tt=eval_time)

# compare results
fig, ax = plt.subplots(figsize=(10, 5))
ax.plot(eval_time, res2[:, 0])
ax.plot(results)
ax.set_xlabel('time')
ax.set_ylabel('x')
plt.legend(['euler', 'ddeint'])
plt.show()

As can be seen, the model dynamics agree qualitatively, with quantitative differences being caused by the different solving algorithms. The most important feature of the DDE function generated by PyRates is that it contains a call to a history function in order to calculate the delayed state variable versions. These history function calls vary in their syntax between backends in order to allow to interface the most common DDE solvers in each programming language. In the case of the default backend, the generated DDE function looks as follows:

f = open('vdp_dde.py', 'r')
print('')
print(f.read())
f.close()

In order to interface any given DDE solver, a wrapper function as defined above may have to be provided to the solver to makes sure that the correct variable is passed to the hist argument.

clear(vdp)

Gallery generated by Sphinx-Gallery