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