The dde module provides tools for solving delay differential equations (DDEs), which are differential equations where the derivative of the unknown function at a certain time depends on the solution at previous times. This module focuses on initial value problems for DDEs (DDE IVP builder).
- Defining a DDE
- The History Function
- Solving an Initial Value Problem (DDE IVP builder)
- Examples
- Notation
The DDE trait defines the delay differential equation dydt = f(t, y(t), y(t - lag_1), ..., y(t - lag_L)) and the lag functions.
diff(&self, t: T, y: &Y, yd: &[Y; L], dydt: &mut Y): Defines the differential equation.yis the current statey(t), andydis an array slice containing the delayed states[y(t - lag_1), ..., y(t - lag_L)].Lis a const generic parameter indicating the number of constant/state-dependent lags.lags(&self, t: T, y: &Y, lags: &mut [T; L]): Defines the actual time lag values[lag_1, ..., lag_L]at timetand statey. These lags determine the past timest - lag_iat which the solution is evaluated.
For event detection with precise zero-crossing detection, implement the separate Event trait:
config- Configure event detection parameters (direction filtering, termination count)event- Event functiong(t,y)whose zero crossings are detected using Brent-Dekker root finding
The Solout trait works similarly to how it does for ODEs, allowing customization of which points are saved in the solution.
use differential_equations::prelude::*;
use nalgebra::Vector3; // Using nalgebra for state vector
// Example: Breast Cancer Model with Delay
// du₁/dt = (v₀ / (1 + β₀ * u₃(t-τ)²)) * (p₀ - q₀) * u₁ - d₀ * u₁
// du₂/dt = (v₀ / (1 + β₀ * u₃(t-τ)²)) * (1 - p₀ + q₀) * u₁ + (v₁ / (1 + β₁ * u₃(t-τ)²)) * (p₁ - q₁) * u₂ - d₁ * u₂
// du₃/dt = (v₁ / (1 + β₁ * u₃(t-τ)²)) * (1 - p₁ + q₁) * u₂ - d₂ * u₃
struct BreastCancerModel {
p0: f64, q0: f64, v0: f64, d0: f64,
p1: f64, q1: f64, v1: f64, d1: f64, d2: f64,
beta0: f64, beta1: f64,
tau: f64, // The time delay
}
// L=1 because there is one delay term, Y is Vector3<f64>
impl DDE<1, f64, Vector3<f64>> for BreastCancerModel {
fn diff(&self, _t: f64, u: &Vector3<f64>, ud: &[Vector3<f64>; 1], dudt: &mut Vector3<f64>) {
// ud[0] corresponds to u(t - lags[0])
// ud[0][2] is u₃(t-τ)
let hist3 = ud[0][2];
let term0_common = self.v0 / (1.0 + self.beta0 * hist3.powi(2));
let term1_common = self.v1 / (1.0 + self.beta1 * hist3.powi(2));
dudt[0] = term0_common * (self.p0 - self.q0) * u[0] - self.d0 * u[0];
dudt[1] = term0_common * (1.0 - self.p0 + self.q0) * u[0]
+ term1_common * (self.p1 - self.q1) * u[1]
- self.d1 * u[1];
dudt[2] = term1_common * (1.0 - self.p1 + self.q1) * u[1] - self.d2 * u[2];
}
fn lags(&self, _t: f64, _y: &Vector3<f64>, lags: &mut [f64; 1]) {
// Define the constant lag
lags[0] = self.tau;
}
}
// Separate event detection implementation
impl Event<f64, Vector3<f64>> for BreastCancerModel {
fn config(&self) -> EventConfig {
EventConfig::new(CrossingDirection::Negative, Some(1)) // Terminate when u1 becomes negative
}
fn event(&self, _t: f64, y: &Vector3<f64>) -> f64 {
// Event function g(t,y) = y[0] - 0.0
// Zero crossing occurs when y[0] = 0.0 (first component becomes negative)
y[0]
}
}Generics <const L: usize, T, Y> are used: L is the number of discrete lags, T is the float type (e.g., f64), Y is the state vector type, and D is the type for event data (defaulting to String).
A crucial component for solving DDEs is the history function. This function, phi(t), provides the solution y(t) for all times t <= t0, where t0 is the initial time of the simulation. The solver uses this function to look up values of y at past times when evaluating delayed terms, especially at the beginning of the integration interval.
The history function is provided when creating a DDE IVP builder and should have the signature Fn(T) -> V.
// Example history function: y(t) = initial_state for t <= 0
// For the BreastCancerModel with state Vector3<f64>
let initial_state = Vector3::new(1.0, 1.0, 1.0);
let history_fn = |_t: f64| -> Vector3<f64> {
initial_state // Return the constant initial state for t <= t0
};The history function takes the time t (where t <= t0) and returns the historical value of y(t) as type Y.
The DDE IVP builder struct is used to set up and solve the DDE. It requires the DDE system, the time interval [t0, tf], the initial state y0 (which is y(t0)), and the history function phi.
fn main() {
let method = ExplicitRungeKutta::dopri5().rtol(1e-6).atol(1e-8); // Using the DDE45 (DOPRI5) solver
let t0 = 0.0;
let tf = 10.0;
let y0 = Vector3::new(1.0, 1.0, 1.0); // Initial state u(t0)
let system = BreastCancerModel {
p0: 0.2, q0: 0.3, v0: 1.0, d0: 5.0,
p1: 0.2, q1: 0.3, v1: 1.0, d1: 1.0, d2: 1.0,
beta0: 1.0, beta1: 1.0,
tau: 1.0, // Delay value
};
// History function: u(t) = y0 for t <= t0
let history_fn = |_t: f64| -> Vector3<f64> {
y0 // Return initial state for all t <= t0
};
let problem = IVP::dde(&system, t0, tf, y0, history_fn);
match problem
.even(0.5) // Example: Save solution every 0.5 time units
.event(&system) // Add event detection
.method(method)
.solve()
{
Ok(solution) => {
// Check if solver terminated due to event
if let Status::Interrupted = solution.status {
println!("Solver terminated due to event detection");
}
println!("Solution (Breast Cancer Model):");
for (t, u_vec) in solution.iter() {
println!("(t: {:.4}, u1: {:.4}, u2: {:.4}, u3: {:.4})", t, u_vec[0], u_vec[1], u_vec[2]);
}
println!("Function evaluations: {}", solution.evals.function);
println!("Steps: {}", solution.steps.total());
}
Err(e) => panic!("Error: {:?}", e),
}
}The solve method returns a Result<Solution, Status>, where Solution contains the time points, corresponding solution values, and solver statistics.
Solving Breast Cancer Model (tau=1) from t=0 to t=10...
Solver finished with status: Complete
Function evaluations: 198
Solver steps: 28
Accepted steps: 27
Rejected steps: 1
Number of output points: 21
t: 0.0000, u1: 1.0000, u2: 1.0000, u3: 1.0000
t: 0.5000, u1: 0.0801, u2: 0.6619, u3: 0.7841
t: 1.0000, u1: 0.0064, u2: 0.3972, u3: 0.5856
t: 1.5000, u1: 0.0005, u2: 0.2348, u3: 0.4282
t: 2.0000, u1: 0.0000, u2: 0.1377, u3: 0.3123
t: 2.5000, u1: 0.0000, u2: 0.0802, u3: 0.2253
t: 3.0000, u1: 0.0000, u2: 0.0466, u3: 0.1597
t: 3.5000, u1: -0.0000, u2: 0.0270, u3: 0.1111
t: 4.0000, u1: -0.0000, u2: 0.0156, u3: 0.0758
t: 4.5000, u1: 0.0000, u2: 0.0090, u3: 0.0510
t: 5.0000, u1: -0.0000, u2: 0.0052, u3: 0.0338
t: 5.5000, u1: 0.0000, u2: 0.0030, u3: 0.0222
t: 6.0000, u1: 0.0000, u2: 0.0017, u3: 0.0144
t: 6.5000, u1: 0.0000, u2: 0.0010, u3: 0.0093
t: 7.0000, u1: -0.0000, u2: 0.0006, u3: 0.0060
t: 7.5000, u1: -0.0000, u2: 0.0003, u3: 0.0038
t: 8.0000, u1: 0.0000, u2: 0.0002, u3: 0.0024
t: 8.5000, u1: 0.0000, u2: 0.0001, u3: 0.0015
t: 9.0000, u1: -0.0000, u2: 0.0001, u3: 0.0010
t: 9.5000, u1: -0.0000, u2: 0.0000, u3: 0.0006
t: 10.0000, u1: 0.0000, u2: 0.0000, u3: 0.0004
For more detailed examples of DDEs, refer to the examples/dde directory in the repository. These examples showcase different DDE systems, solver configurations, and usage patterns.
- Mackey-Glass Equation: A common benchmark DDE exhibiting chaotic behavior.
- Breast Cancer Model: A model of tumor growth with delays in the response of the immune system.
This library uses the following notation for DDEs:
t: The independent variable, typically time.y: The dependent state vectory(t).dydt: The derivative ofywith respect tot, i.e.,y'(t).lags: An array[lag_1, ..., lag_L]containing the time delay values.yd: An array[y(t - lag_1), ..., y(t - lag_L)]containing the state vectors at delayed times.k: Coefficients or intermediate stages within solver algorithms, often representing approximations to derivatives.
Consistency with this notation is encouraged for any future additions to the DDE module.