conspire/math/integrate/ode/explicit/fixed_step/
mod.rs

1#[cfg(test)]
2mod test;
3
4use crate::math::{
5    Scalar, Tensor, TensorVec, Vector,
6    integrate::{Explicit, FixedStep, IntegrationError},
7};
8
9pub mod bogacki_shampine;
10pub mod dormand_prince;
11pub mod euler;
12pub mod heun;
13pub mod midpoint;
14pub mod ralston;
15pub mod verner_8;
16pub mod verner_9;
17
18/// Fixed-step explicit integrators for ordinary differential equations.
19pub trait FixedStepExplicit<Y, U>
20where
21    Self: Explicit<Y, U> + FixedStep,
22    Y: Tensor,
23    U: TensorVec<Item = Y>,
24{
25    fn integrate_fixed_step(
26        &self,
27        mut function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
28        time: &[Scalar],
29        initial_condition: Y,
30    ) -> Result<(Vector, U, U), IntegrationError> {
31        let t_0 = time[0];
32        let t_f = time[time.len() - 1];
33        let mut t_sol: Vector;
34        if time.len() < 2 {
35            return Err(IntegrationError::LengthTimeLessThanTwo);
36        } else if t_0 >= t_f {
37            return Err(IntegrationError::InitialTimeNotLessThanFinalTime);
38        } else if time.len() == 2 {
39            if self.dt() <= 0.0 || self.dt().is_nan() {
40                return Err(IntegrationError::TimeStepNotSet(
41                    time[0],
42                    time[1],
43                    format!("{self:?}"),
44                ));
45            } else {
46                let num_steps = ((t_f - t_0) / self.dt()).ceil() as usize;
47                t_sol = (0..num_steps)
48                    .map(|step| t_0 + (step as Scalar) * self.dt())
49                    .collect();
50                t_sol.push(t_f);
51            }
52        } else {
53            t_sol = time.iter().copied().collect();
54        }
55        let mut index = 0;
56        let mut t = t_0;
57        let mut dt;
58        let mut t_trial;
59        let mut k = vec![Y::default(); Self::SLOPES];
60        k[0] = function(t, &initial_condition)?;
61        let mut y = initial_condition.clone();
62        let mut y_sol = U::new();
63        y_sol.push(initial_condition.clone());
64        let mut dydt_sol = U::new();
65        dydt_sol.push(function(t, &y.clone())?);
66        let mut y_trial = Y::default();
67        while t < t_f {
68            t_trial = t_sol[index + 1];
69            dt = t_trial - t;
70            if let Err(error) = self.step(&mut function, &y, t, dt, &mut k, &mut y_trial) {
71                return Err(IntegrationError::Upstream(error, format!("{self:?}")));
72            } else {
73                t += dt;
74                y = y_trial.clone();
75                y_sol.push(y.clone());
76                dydt_sol.push(k[0].clone());
77                index += 1;
78            }
79        }
80        Ok((t_sol, y_sol, dydt_sol))
81    }
82    #[allow(clippy::too_many_arguments)]
83    fn step(
84        &self,
85        function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
86        y: &Y,
87        t: Scalar,
88        dt: Scalar,
89        k: &mut [Y],
90        y_trial: &mut Y,
91    ) -> Result<(), String>;
92}