ddeint is a small toolkit for numerically solving delay differential equations (DDEs), i.e. ODEs whose right-hand side depends on the state at past times t - delay. It wraps integrators from the ode_solvers crate, plus a piecewise-linear history interpolator (DdeVar/DdeVars) that is updated after every accepted integration step so that delayed lookups always reflect the actual solved trajectory rather than only the initial history.

Three integrators are available, all with the same shape (new/step/ integrate, a results: Vec<(f64, DVector<f64>)> field) so switching between them is close to a drop-in change:

For all three, step_size controls how often the delayed history is updated and a result is recorded. For DdeRK4 that's also the literal integration step; for the two adaptive solvers, the actual number of internal steps between two recorded points is chosen by their own error control, and the state at exactly t + step_size is obtained from their dense-output interpolant rather than by counting steps. See examples/integrator_comparison.rs for all three solving the same DDE side by side (and a note on why the higher-order integrator isn't automatically the more accurate one here).

use ddeint::Interpolator;

let x = vec![0.0, 1.0, 2.0];
let y = vec![0.0, 1.0, 4.0];
let interpolator = Interpolator::new(x, y, /* fill_value = */ 0.0);

assert_eq!(interpolator.interpolate(0.5), Some(0.5));

A System that needs delayed state has to implement DdeSystem in addition to ode_solvers's System trait; DdeSystem::record is called by DdeRK4 after every step and is where the system should push the newly solved state into its own history storage (typically an Rc<RefCell<DdeVars>> shared with the closure/struct that reads it back in system()):

use std::cell::RefCell;
use std::rc::Rc;
use std::sync::Arc;
use nalgebra::DVector;
use ddeint::{DdeRK4, DdeSystem, DdeVars};

#[derive(Clone)]
struct DelayedPredatorPrey {
    delay: f64,
    history: Rc<RefCell<DdeVars>>,
}

impl ode_solvers::System<f64, DVector<f64>> for DelayedPredatorPrey {
    fn system(&self, t: f64, y: &DVector<f64>, dy: &mut DVector<f64>) {
        let history = self.history.borrow();
        let prey_delayed = history[0].instance_value(t - self.delay);
        let predator_delayed = history[1].instance_value(t - self.delay);

        dy[0] = 0.5 * y[0] * (1.0 - predator_delayed);
        dy[1] = -0.5 * y[1] * (1.0 - prey_delayed);
    }
}

impl DdeSystem for DelayedPredatorPrey {
    fn record(&self, t: f64, y: &DVector<f64>) {
        let mut history = self.history.borrow_mut();
        history[0].update(t, y[0]);
        history[1].update(t, y[1]);
    }
}

See src/lotka_delay.rs for the full, runnable version of this example (ddeint::solve_lotka_volterra), also wired into main.rs.

Runnable examples live under examples/. Each one writes its result to a CSV file (in the directory cargo is run from) in addition to printing a short summary:

cargo run --release --example three_var_dde          # 3-variable, 2-delay system (Willé & Baker's classic example)
cargo run --release --example sine_delay              # single delay, checked against its exact closed-form solution
cargo run --release --example variable_delay          # a delay that is itself a function of time, not a constant
cargo run --release --example integrator_comparison   # DdeRK4 vs DdeDopri5 vs DdeDop853 on the same DDE

The exact-solution and convergence checks these examples use for verification are also captured as cargo test-run regression tests in tests/analytic_verification.rs.

This project was inspired by Zulko's Python ddeint, a lightweight DDE solver built on top of scipy.integrate.ode. Many thanks to Zulko for the original design and the history-interpolation (ddeVar) idea this crate is built around.

One thing this crate does differently: in the original ddeint, the whole state is a single vector-valued history function, so every state variable is looked up with the same delay. This port's DdeVars instead keeps a separate history per state variable, so different variables can use different delay values within the same system -- see examples/three_var_dde.rs, where y1 is delayed by 1.0 and y2 by 0.2 in the same system.