differential_equations/dde/
problem.rs

1//! Initial Value Problem Struct and Constructors for Delay Differential Equations (DDEs)
2
3use crate::{
4    dde::{DDE, numerical_method::DelayNumericalMethod, solve_dde},
5    error::Error,
6    interpolate::Interpolation,
7    solout::*,
8    solution::Solution,
9    traits::{Real, State},
10};
11
12/// Initial Value Problem for Delay Differential Equations (DDEs)
13///
14/// The Initial Value Problem takes the form:
15/// y'(t) = f(t, y(t), y(t-tau1), ...),  t ∈ [t0, tf],  y(t) = phi(t) for t <= t0
16///
17/// # Overview
18///
19/// The DDEProblem struct provides a simple interface for solving DDEs:
20///
21/// # Example
22///
23/// ```
24/// use differential_equations::prelude::*;
25/// use nalgebra::{Vector2, vector};
26///
27/// let mut rkf45 = ExplicitRungeKutta::rkf45()
28///    .rtol(1e-6)
29///    .atol(1e-6);
30///
31/// let t0 = 0.0;
32/// let tf = 10.0;
33/// let y0 = vector![1.0, 0.0];
34/// let phi = |t| {
35///     y0
36/// };
37/// struct Example;
38/// impl DDE<1, f64, Vector2<f64>> for Example {
39///     fn diff(&self, t: f64, y: &Vector2<f64>, yd: &[Vector2<f64>; 1], dydt: &mut Vector2<f64>) {
40///        dydt[0] = y[1] + yd[0][0];
41///        dydt[1] = -y[0] + yd[0][1];
42///     }
43///
44///     fn lags(&self, _t: f64, _y: &Vector2<f64>, lags: &mut [f64; 1]) {
45///        lags[0] = 1.0; // Fixed delay of 1.0
46///     }
47/// }
48/// let solution = DDEProblem::new(Example, t0, tf, y0, phi).solve(&mut rkf45).unwrap();
49///
50/// let (t, y) = solution.last().unwrap();
51/// println!("Solution: ({}, {})", t, y);
52/// ```
53///
54/// # Fields
55///
56/// * `dde` - DDE implementing the differential equation
57/// * `t0` - Initial time
58/// * `tf` - Final time
59/// * `y0` - Initial state vector at `t0`
60/// * `phi` - Initial history function `phi(t)` for `t <= t0`
61///
62/// # Basic Usage
63///
64/// * `new(dde, t0, tf, y0, phi)` - Create a new DDEProblem
65/// * `solve(&mut solver)` - Solve using default output (solver step points)
66///
67/// # Output Control Methods
68///
69/// These methods configure how solution points are generated and returned:
70///
71/// * `even(dt)` - Generate evenly spaced output points with interval `dt`
72/// * `dense(n)` - Include `n` interpolated points between each solver step
73/// * `t_eval(points)` - Evaluate solution at specific time points
74/// * `solout(custom_solout)` - Use a custom output handler
75///
76/// Each returns a solver configuration that can be executed with `.solve(&mut solver)`.
77#[derive(Clone, Debug)] // Note: Clone requires F and H to be Clone
78pub struct DDEProblem<'a, const L: usize, T, Y, F, H>
79where
80    T: Real,
81    Y: State<T>,
82    F: DDE<L, T, Y>,
83    H: Fn(T) -> Y,
84{
85    // Initial Value Problem Fields
86    pub dde: &'a F, // DDE containing the Differential Equation and Optional Terminate Function.
87    pub t0: T,      // Initial Time.
88    pub tf: T,      // Final Time.
89    pub y0: Y,      // Initial State Vector.
90    pub phi: H,     // Initial History Function.
91}
92
93impl<'a, const L: usize, T, Y, F, H> DDEProblem<'a, L, T, Y, F, H>
94where
95    T: Real,
96    Y: State<T>,
97    F: DDE<L, T, Y>,
98    H: Fn(T) -> Y + Clone, // Require Clone for H if DDEProblem needs Clone
99{
100    /// Create a new Initial Value Problem for DDEs
101    ///
102    /// # Arguments
103    /// * `dde`  - DDE containing the Differential Equation and Optional Terminate Function.
104    /// * `t0`   - Initial Time.
105    /// * `tf`   - Final Time.
106    /// * `y0`   - Initial State Vector at `t0`.
107    /// * `phi`  - Initial History Function `phi(t)` for `t <= t0`.
108    ///
109    /// # Returns
110    /// * DDEProblem Problem ready to be solved.
111    ///
112    pub fn new(dde: &'a F, t0: T, tf: T, y0: Y, phi: H) -> Self {
113        DDEProblem {
114            dde,
115            t0,
116            tf,
117            y0,
118            phi,
119        }
120    }
121
122    /// Solve the DDEProblem using a default solout, e.g. outputting solutions at calculated steps.
123    ///
124    /// # Returns
125    /// * `Result<Solution<T, Y>, Error<T, Y>>` - `Ok(Solution)` if successful or interrupted by events, `Err(Error)` if errors or issues are encountered.
126    ///
127    pub fn solve<S>(&self, solver: &mut S) -> Result<Solution<T, Y>, Error<T, Y>>
128    where
129        S: DelayNumericalMethod<L, T, Y, H> + Interpolation<T, Y>,
130        H: Clone, // phi needs to be cloneable for solve_dde
131    {
132        let mut default_solout = DefaultSolout::new(); // Default solout implementation
133        solve_dde(
134            solver,
135            self.dde,
136            self.t0,
137            self.tf,
138            &self.y0,
139            self.phi.clone(),
140            &mut default_solout,
141        )
142    }
143
144    /// Returns a DDEProblem DelayNumericalMethod with the provided solout function for outputting points.
145    ///
146    /// # Returns
147    /// * DDEProblem DelayNumericalMethod with the provided solout function ready for .solve() method.
148    ///
149    pub fn solout<O: Solout<T, Y>>(
150        &'a self,
151        solout: &'a mut O,
152    ) -> DDEProblemMutRefSoloutPair<'a, L, T, Y, F, H, O> {
153        DDEProblemMutRefSoloutPair::new(self, solout)
154    }
155
156    /// Uses the an Even Solout implementation to output evenly spaced points between the initial and final time.
157    /// Note that this does not include the solution of the calculated steps.
158    ///
159    /// # Arguments
160    /// * `dt` - Interval between each output point.
161    ///
162    /// # Returns
163    /// * DDEProblem DelayNumericalMethod with Even Solout function ready for .solve() method.
164    ///
165    pub fn even(&self, dt: T) -> DDEProblemSoloutPair<'_, L, T, Y, F, H, EvenSolout<T>> {
166        let even_solout = EvenSolout::new(dt, self.t0, self.tf); // Even solout implementation
167        DDEProblemSoloutPair::new(self, even_solout)
168    }
169
170    /// Uses the Dense Output method to output n number of interpolation points between each step.
171    /// Note this includes the solution of the calculated steps.
172    ///
173    /// # Arguments
174    /// * `n` - Number of interpolation points between each step.
175    ///
176    /// # Returns
177    /// * DDEProblem DelayNumericalMethod with Dense Output function ready for .solve() method.
178    ///
179    pub fn dense(&self, n: usize) -> DDEProblemSoloutPair<'_, L, T, Y, F, H, DenseSolout> {
180        let dense_solout = DenseSolout::new(n); // Dense solout implementation
181        DDEProblemSoloutPair::new(self, dense_solout)
182    }
183
184    /// Uses the provided time points for evaluation instead of the default method.
185    /// Note this does not include the solution of the calculated steps.
186    ///
187    /// # Arguments
188    /// * `points` - Custom output points.
189    ///
190    /// # Returns
191    /// * DDEProblem DelayNumericalMethod with Custom Time Evaluation function ready for .solve() method.
192    ///
193    pub fn t_eval(
194        &self,
195        points: impl AsRef<[T]>,
196    ) -> DDEProblemSoloutPair<'_, L, T, Y, F, H, TEvalSolout<T>> {
197        let t_eval_solout = TEvalSolout::new(points, self.t0, self.tf); // Custom time evaluation solout implementation
198        DDEProblemSoloutPair::new(self, t_eval_solout)
199    }
200
201    /// Uses the CrossingSolout method to output points when a specific component crosses a threshold.
202    /// Note this does not include the solution of the calculated steps.
203    ///
204    /// # Arguments
205    /// * `component_idx` - Index of the component to monitor for crossing.
206    /// * `threshhold` - Value to cross.
207    /// * `direction` - Direction of crossing (positive or negative).
208    ///
209    /// # Returns
210    /// * DDEProblem DelayNumericalMethod with CrossingSolout function ready for .solve() method.
211    ///
212    pub fn crossing(
213        &self,
214        component_idx: usize,
215        threshhold: T,
216        direction: CrossingDirection,
217    ) -> DDEProblemSoloutPair<'_, L, T, Y, F, H, CrossingSolout<T>> {
218        let crossing_solout =
219            CrossingSolout::new(component_idx, threshhold).with_direction(direction); // Crossing solout implementation
220        DDEProblemSoloutPair::new(self, crossing_solout)
221    }
222
223    /// Uses the HyperplaneCrossingSolout method to output points when a specific hyperplane is crossed.
224    /// Note this does not include the solution of the calculated steps.
225    ///
226    /// # Arguments
227    /// * `point` - Point on the hyperplane.
228    /// * `normal` - Normal vector of the hyperplane.
229    /// * `extractor` - Function to extract the component from the state vector.
230    /// * `direction` - Direction of crossing (positive or negative).
231    ///
232    /// # Returns
233    /// * DDEProblem DelayNumericalMethod with HyperplaneCrossingSolout function ready for .solve() method.
234    ///
235    pub fn hyperplane_crossing<Y1>(
236        &self,
237        point: Y1,
238        normal: Y1,
239        extractor: fn(&Y) -> Y1,
240        direction: CrossingDirection,
241    ) -> DDEProblemSoloutPair<'_, L, T, Y, F, H, HyperplaneCrossingSolout<T, Y1, Y>>
242    where
243        Y1: State<T>,
244    {
245        let solout =
246            HyperplaneCrossingSolout::new(point, normal, extractor).with_direction(direction);
247
248        DDEProblemSoloutPair::new(self, solout)
249    }
250
251    /// Uses an `EventSolout` to detect zero crossings g(t,y)=0 of a user-defined event function.
252    pub fn event<E>(
253        &self,
254        event: &'a E,
255    ) -> DDEProblemSoloutPair<'_, L, T, Y, F, H, EventSolout<'_, T, Y, E>>
256    where
257        E: Event<T, Y>,
258    {
259        let solout = EventSolout::new(event, self.t0, self.tf);
260        DDEProblemSoloutPair::new(self, solout)
261    }
262}
263
264/// DDEProblemMutRefSoloutPair serves as an intermediate between the DDEProblem struct and a custom solout provided by the user.
265pub struct DDEProblemMutRefSoloutPair<'a, const L: usize, T, Y, F, H, O>
266where
267    T: Real,
268    Y: State<T>,
269    F: DDE<L, T, Y>,
270    H: Fn(T) -> Y,
271    O: Solout<T, Y>,
272{
273    pub problem: &'a DDEProblem<'a, L, T, Y, F, H>, // Reference to the DDEProblem struct
274    pub solout: &'a mut O,                          // Reference to the solout implementation
275}
276
277impl<'a, const L: usize, T, Y, F, H, O> DDEProblemMutRefSoloutPair<'a, L, T, Y, F, H, O>
278where
279    T: Real,
280    Y: State<T>,
281    F: DDE<L, T, Y>,
282    H: Fn(T) -> Y + Clone, // Require Clone for H
283    O: Solout<T, Y>,
284{
285    /// Create a new DDEProblemMutRefSoloutPair
286    ///
287    /// # Arguments
288    /// * `problem` - Reference to the DDEProblem struct
289    /// * `solout` - Mutable reference to the solout implementation
290    ///
291    pub fn new(problem: &'a DDEProblem<L, T, Y, F, H>, solout: &'a mut O) -> Self {
292        DDEProblemMutRefSoloutPair { problem, solout }
293    }
294
295    /// Solve the DDEProblem using the provided solout
296    ///
297    /// # Arguments
298    /// * `solver` - DelayNumericalMethod to use for solving the DDEProblem
299    ///
300    /// # Returns
301    /// * `Result<Solution<T, Y>, Error<T, Y>>` - `Ok(Solution)` if successful or interrupted by events, `Err(Error)` if errors or issues are encountered.
302    ///
303    pub fn solve<S>(&mut self, solver: &mut S) -> Result<Solution<T, Y>, Error<T, Y>>
304    where
305        S: DelayNumericalMethod<L, T, Y, H> + Interpolation<T, Y>,
306    {
307        solve_dde(
308            solver,
309            self.problem.dde,
310            self.problem.t0,
311            self.problem.tf,
312            &self.problem.y0,
313            self.problem.phi.clone(),
314            self.solout,
315        )
316    }
317}
318
319/// DDEProblemSoloutPair serves as an intermediate between the DDEProblem struct and solve_dde when a predefined solout is used.
320#[derive(Clone, Debug)]
321pub struct DDEProblemSoloutPair<'a, const L: usize, T, Y, F, H, O>
322where
323    T: Real,
324    Y: State<T>,
325    F: DDE<L, T, Y>,
326    H: Fn(T) -> Y,
327    O: Solout<T, Y>,
328{
329    pub problem: &'a DDEProblem<'a, L, T, Y, F, H>,
330    pub solout: O,
331}
332
333impl<'a, const L: usize, T, Y, F, H, O> DDEProblemSoloutPair<'a, L, T, Y, F, H, O>
334where
335    T: Real,
336    Y: State<T>,
337    F: DDE<L, T, Y>,
338    H: Fn(T) -> Y + Clone,
339    O: Solout<T, Y>,
340{
341    /// Create a new DDEProblemSoloutPair
342    ///
343    /// # Arguments
344    /// * `problem` - Reference to the DDEProblem struct
345    /// * `solout` - Solout implementation
346    ///
347    pub fn new(problem: &'a DDEProblem<L, T, Y, F, H>, solout: O) -> Self {
348        DDEProblemSoloutPair { problem, solout }
349    }
350
351    /// Solve the DDEProblem using the provided solout
352    ///
353    /// # Arguments
354    /// * `solver` - DelayNumericalMethod to use for solving the DDEProblem
355    ///
356    /// # Returns
357    /// * `Result<Solution<T, Y>, Error<T, Y>>` - `Ok(Solution)` if successful or interrupted by events, `Err(Error)` if errors or issues are encountered.
358    ///
359    pub fn solve<S>(mut self, solver: &mut S) -> Result<Solution<T, Y>, Error<T, Y>>
360    where
361        S: DelayNumericalMethod<L, T, Y, H> + Interpolation<T, Y>,
362    {
363        solve_dde(
364            solver,
365            self.problem.dde,
366            self.problem.t0,
367            self.problem.tf,
368            &self.problem.y0,
369            self.problem.phi.clone(),
370            &mut self.solout,
371        )
372    }
373
374    /// Wrap current solout with event detection while preserving original output strategy.
375    pub fn event<E>(
376        self,
377        event: &'a E,
378    ) -> DDEProblemSoloutPair<'a, L, T, Y, F, H, EventWrappedSolout<'a, T, Y, O, E>>
379    where
380        E: Event<T, Y>,
381    {
382        let wrapped = EventWrappedSolout::new(self.solout, event, self.problem.t0, self.problem.tf);
383        DDEProblemSoloutPair {
384            problem: self.problem,
385            solout: wrapped,
386        }
387    }
388}
differential_equations/dde/problem.rs

differential_equations/dde/
problem.rs
