Delay and stochastic equations

import numpy as np
import sympy as sp
from odecast import var, Eq, solve, solve_sde

Delay differential equations

In a delay differential equation the rate of change depends on earlier values of the solution. y.delay(tau) stands for \(y(t - \tau)\), and y.d().delay(tau) for \(y'(t - \tau)\):

# Delayed logistic growth (Hutchinson's equation)
y = var("y")
sol = solve(Eq(y.d(), 1.8*y*(1 - y.delay(1))), ivp={y: 0.5}, tspan=(0, 50))

The delay may be a parameter, e.g. y.delay(tau) with params={tau: 1}, but must be constant in time and positive.

The equation needs the solution on the interval \([t_0 - \tau, t_0]\), the history. By default the solution is taken to be constant at its initial value before \(t_0\); history gives other values, as numbers or functions of \(t\):

sol = solve(Eq(y.d(), -y.delay(1)), ivp={y: 1}, tspan=(0, 10),
            history={y: lambda s: np.cos(s)})

The equations are solved with the method of steps. The interval is divided at \(t_0 + k\tau\) for every delay \(\tau\), where discontinuities of the derivatives of the solution propagate, and each piece is integrated with scipy.integrate.solve_ivp(). Delayed values come from the history or from the dense output of the pieces computed before.

Limitations: variables and derivatives below the highest order can be delayed, but not the highest derivative itself (neutral equations). Delay equations cannot be combined with events, sensitivities, boundary value problems or algebraic equations, and the analysis functions and the SymPy backend do not support them.

Stochastic differential equations

solve_sde() simulates equations with random forcing, in the Itô interpretation. The equations give the deterministic part (the drift); noise adds a term \(g\,dW\) with an independent Wiener process \(W\) for each entry:

# Geometric Brownian motion dY = mu Y dt + sigma Y dW
y = var("y")
sol = solve_sde(Eq(y.d(), 0.1*y), noise={y: 0.3*y}, ivp={y: 1.0},
                tspan=(0, 1), dt=1e-3, paths=1000, seed=1)

sol[y]          # all paths, shape (1000, len(sol.t))
sol.mean(y)     # ensemble mean as a function of t
sol.std(y)      # ensemble standard deviation
sol.quantile(y, [0.05, 0.95])
sol.plot(y)     # sample paths with mean and standard deviation

For a variable of higher order the noise enters the equation of its highest derivative. It is given for the state whose derivative that is; for the Langevin equation \(x'' + \gamma x' + \omega^2 x = \sigma\xi(t)\):

x = var("x")
sol = solve_sde(Eq(x.d(2) + 0.5*x.d() + 4*x, 0), noise={x.d(): 0.4},
                ivp={x: 0, x.d(): 0}, tspan=(0, 30), dt=1e-3, paths=500)

All paths are computed together with a fixed time step dt. Two methods are available:

"euler" (default)

Euler-Maruyama, strong order 1/2.

"milstein"

Milstein, strong order 1. It uses the derivative of each noise term with respect to its own state, computed symbolically.

sol.W contains the Wiener processes that drove the paths, so that results can be compared with exact solutions. seed makes the simulation reproducible.