solve an ODE function of a previous time step (delay differential equation)

Question

I have this set of differential equations:

dy/dt = a*y   - b*x*y  
dx/dt = b*x*y - c*y(t - t_0)

t_0 is a constant time and the term is neglected when t<t_0. How do I solve this in python using numpy/scipy, given initial conditions and all coefficients?

EDIT: the y(t-t_0) is the value of y at the time t-t_0, NOT y times t-t_0

Check out https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html — PilouPili, Sep 09 '18 at 12:03
Then your system is a delay-differential equation and you need to use a DDE solver. While these exist in matlab, the standard packages for python do not cover these equations. — Lutz Lehmann, Sep 09 '18 at 13:32

score 3 · Answer 1 · answered Sep 10 '18 at 15:13

It seems that performing the interpolation on the global variable sol_y works also:

import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt

def dudt(t, u, params):
    x, y = u
    a, b, c, t0 = params

    dydt = a*y   - b*x*y  
    if t <= t0:
        dxdt = b*x*y
    else:
        dxdt = b*x*y - c*get_old_y(t-t0)

    return [dxdt, dydt]

def get_old_y(old_t):
    return np.interp(old_t, sol_t, sol_y)

def jac_dudt(t, u, params):
    x, y = u
    a, b, c, t0 = params
    jac = [[ b*y, b*x-c],
           [-b*y, a-b*y]]
    return jac

# parameters
t0 = 1
params = 1, 1, 2, t0

u0 = [1, 2]

t_end = 3*t0
dt = 0.05

# integration
r = ode(dudt, jac_dudt).set_integrator("vode",
                            method="adams",
                            with_jacobian=True)

r.set_initial_value(u0, 0).set_f_params(params).set_jac_params(params)

sol_t, sol_x, sol_y = [], [], []                             
while r.successful() and r.t < t_end:
    r.integrate(r.t + dt)
    sol_x.append(r.y[0])
    sol_y.append(r.y[1])
    sol_t.append(r.t)

# graph
plt.plot(sol_t, sol_x, '-|', label='x(t)')
plt.plot(sol_t, sol_y, '-|', label='y(t)')
plt.legend(); plt.xlabel('time'); plt.ylabel('solution');

The output graph, with example parameters, is:

If I have, say, another equation with a variable z, that depended on t, and whose derivative depended on x, y and t, then how would I go about finding the Jacobian for more than 2 equations and plotting a curve? — theduckgoesquark, Sep 10 '18 at 15:18
@theduckgoesquark The jacobian is not necessary, it will be estimated numerically if not provided, you can try without it first. (The jacobian will be a 3x3 matrix with the corresponding derivatives) — xdze2, Sep 10 '18 at 15:27

Merlin1896 · Answer 2 · 2018-09-09T13:49:08.047

In an earlier version of the question, the problem stated only a simple system of ODE. It was then changed to a delay differential equation and the answer below is no longer valid. I leave it for future reference.

To solve a system with delay, additional python packages have to be used. For example the package JiTCDDE allows to solve this kind of equations. A related question was asked here: Solve ODE in Python with a time-delay

Old Answer

The scipy function ode might be what you are looking for:

Lets first define the two equation systems. One for t<t0 and one for t>t0. We call these functions f and f2. Additionally, we also calculate the Jacobian matrix, which can later be used by the integrator.

def f(t,y,a,b,c,t_0):                                  
    return [b*y[0]*y[1]-c*(t-t_0), a*y[1]-b*y[0]*y[1]]

def f2(t,y,a,b,c):                           
    return [b*y[0]*y[1], a*y[1]-b*y[0]*y[1]] 
def jac_f(t,y,a,b):                            
    return [[b*y[1],b*y[0]],[-b*y[1],a-b*y[1]]]

Then we import ode and call the integrator two times. The first time, we integrate from our start value (I set it to t=0) until we reach t0 and then start a second integration with the equation system valid for t>t0. We pass the last calculated values as initial conditions to the integrator and continue our integration until we reach t=4 (arbitrarily chosen).

from scipy.integrate import ode 
y_res = []                             
t_list = []
a,b,c,t0=1,1,1,1
y0=[1,2]
t_start=0
t_fin=t0
dt=0.01
r=ode(f2,jac_f).set_integrator("vode", method="adams", with_jacobian=True)
r.set_initial_value(y0, t_start).set_f_params(a,b).set_jac_params(a,b)
while r.successful() and r.t < t_fin:
    r.integrate(r.t+dt)
    y_res.append(r.y)
    t_list.append(r.t)
y0=y_res[-1]
t_start=t0
t_fin=4
dt=0.01
r=ode(f,jac_f).set_integrator("vode", method="adams", with_jacobian=True)
r.set_initial_value(y0, t_start).set_f_params(a,b,c,t0).set_jac_params(a,b)
while r.successful() and r.t < t_fin:
    r.integrate(r.t+dt)
    y_res.append(r.y)
    t_list.append(r.t)

We can now plot the resulting curves:

import matplotlib.pyplot as plt
yy=np.stack(y_res)
plt.plot(t_list, yy[:,0], label="x(t)")
plt.plot(t_list, yy[:,1], label="y(t)")
plt.legend()
plt.show()

We get this nice graph:

could you please explain the "vode" and the "adams" and the integrate code in brief? — theduckgoesquark, Sep 09 '18 at 12:58
also, what if i have more equations that do now have the t-t_0 dependency? — theduckgoesquark, Sep 09 '18 at 13:02
What you can do to have only one function is to use `c*max(0,t-t0)`. As that is continuous the solver should have no problem integrating over that point, even if it has to regulate the step size down around `t0`. — Lutz Lehmann, Sep 09 '18 at 13:05
@LutzL, that is actually a nicer way to solve this than my attempt. @ theduckgoesquark, you can read up on different integration methods in various books and publications. If you have a specific question on how these methods work, please ask a new question or maybe head over to math.SE. A starting point for reading about vode is the publication of Brown et al.: https://computation.llnl.gov/casc/nsde/pubs/207532.pdf — Merlin1896, Sep 09 '18 at 13:14
but i want a function of t-t_0, not the value of t-t_0 itself. How would I go about that? — theduckgoesquark, Sep 09 '18 at 13:19

solve an ODE function of a previous time step (delay differential equation)

2 Answers2

Old Answer