︠db934717-fcba-4795-96b4-33c223398b9dis︠ %hide %html

Function definition: dde_solve(dde, statevar, delayedvars, history, tmax, timestep)

Description: Numerically integrate a delay differential equation

A single-variable delay differential equation is an equation of the form X'(t) = f(X, Xτ1, Xτ2, ... , Xτn) where on the right hand side X denotes X(t) and each Xτ denotes X(t - τ) for some constant delay τ.

Inputs:

dde: The right hand side of the delay differential equation

statevar: The state variable to use (X in the example above)

delayedvars: A dictionary whose keys are the symbolic variables that appear in the DDE expression, which represent the value of the state variable some time ago, and whose values are the actual delays. For example, if the DDE is X'(t) = X(t) - X(t - 5) + X(t - 17.4) the dde argument might be X - Xtau + Xtau2, in which case this argument would be {Xtau: 5, Xtau2: 17.4}. Note that this means before calling this function, all three of X, Xtau, and Xtau2 should have been declared as symbolic variables: var("X, Xtau, Xtau2")

history: A Sage expression for the history, as a function of t (time)

tmax: Integrate from t = 0 to this time

timestep: The size of each time step (the smaller, the more accurate)

Output: A list of the values of the state variable (X in the example above) at each time step from t = 0 to tmax

Examples:

A version of Mackey-Glass, using a constant 0 history:

var("X, X_t")
sol = dde_solve(1 - X * X_t^5/(1 + X_t^5), X, {X_t: 8}, history=0, tmax=400, timestep=0.1)
list_plot(zip(srange(0, 400, 0.1), sol), plotjoined=True)

The same, but with an oscillating history function, giving weird transients that give way to a periodic steady state:

var("t, X, X_tau")
eq = 1 - X * X_tau^5/(1 + X_tau^5)
sol = dde_solve(eq, X, {X_tau: 8}, 5*cos(t), 400, 0.1)
list_plot(zip(srange(0, 400, 0.1), sol), plotjoined=True)

︡57fb1982-118d-4fdd-93b2-bedb6ff15384︡ ︠ac5084bc-0417-4eba-9efd-5bde129b5e61as︠ %auto def dde_solve(dde, statevar, delayedvars, history, tmax, timestep): # Check validity of delays. if min(delayedvars.values()) < 0: raise ValueError("This function will not work with negative delays. " "Consider consulting a fortune teller instead.") if any(val