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

1#[cfg(test)]
2mod test;
3
4use crate::math::{
5    Scalar, Tensor, TensorVec, Vector,
6    integrate::{Explicit, IntegrationError, VariableStep},
7    interpolate::InterpolateSolution,
8};
9use std::ops::{Mul, Sub};
10
11pub mod internal_variables;
12
13pub mod bogacki_shampine;
14pub mod dormand_prince;
15pub mod verner_8;
16pub mod verner_9;
17
18/// Variable-step explicit ordinary differential equation solvers.
19pub trait VariableStepExplicit<Y, U>
20where
21    Self: InterpolateSolution<Y, U> + Explicit<Y, U> + VariableStep,
22    Y: Tensor,
23    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
24    U: TensorVec<Item = Y>,
25{
26    fn integrate_variable_step(
27        &self,
28        mut function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
29        time: &[Scalar],
30        initial_condition: Y,
31    ) -> Result<(Vector, U, U), IntegrationError> {
32        let t_0 = time[0];
33        let t_f = time[time.len() - 1];
34        if time.len() < 2 {
35            return Err(IntegrationError::LengthTimeLessThanTwo);
36        } else if t_0 >= t_f {
37            return Err(IntegrationError::InitialTimeNotLessThanFinalTime);
38        }
39        let mut t = t_0;
40        let mut dt = t_f - t_0;
41        let mut k = vec![Y::default(); Self::SLOPES];
42        k[0] = function(t, &initial_condition)?;
43        let mut t_sol = Vector::new();
44        t_sol.push(t_0);
45        let mut y = initial_condition.clone();
46        let mut y_sol = U::new();
47        y_sol.push(initial_condition.clone());
48        let mut dydt_sol = U::new();
49        dydt_sol.push(k[0].clone());
50        let mut y_trial = Y::default();
51        while t < t_f {
52            match self.slopes(&mut function, &y, t, dt, &mut k, &mut y_trial) {
53                Ok(e) => {
54                    if let Err(error) = self.step(
55                        &mut function,
56                        &mut y,
57                        &mut t,
58                        &mut y_sol,
59                        &mut t_sol,
60                        &mut dydt_sol,
61                        &mut dt,
62                        &mut k,
63                        &y_trial,
64                        e,
65                    ) {
66                        dt *= self.dt_cut();
67                        if dt < self.dt_min() {
68                            return Err(IntegrationError::MinimumStepSizeUpstream(
69                                self.dt_min(),
70                                error,
71                                format!("{self:?}"),
72                            ));
73                        }
74                    } else {
75                        dt = dt.min(t_f - t);
76                        if dt < self.dt_min() && t < t_f {
77                            return Err(IntegrationError::MinimumStepSizeReached(
78                                self.dt_min(),
79                                format!("{self:?}"),
80                            ));
81                        }
82                    }
83                }
84                Err(error) => {
85                    dt *= self.dt_cut();
86                    if dt < self.dt_min() {
87                        return Err(IntegrationError::MinimumStepSizeUpstream(
88                            self.dt_min(),
89                            error,
90                            format!("{self:?}"),
91                        ));
92                    }
93                }
94            }
95        }
96        if time.len() > 2 {
97            let t_int = Vector::from(time);
98            let (y_int, dydt_int) = self.interpolate(&t_int, &t_sol, &y_sol, function)?;
99            Ok((t_int, y_int, dydt_int))
100        } else {
101            Ok((t_sol, y_sol, dydt_sol))
102        }
103    }
104    fn slopes(
105        &self,
106        function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
107        y: &Y,
108        t: Scalar,
109        dt: Scalar,
110        k: &mut [Y],
111        y_trial: &mut Y,
112    ) -> Result<Scalar, String>;
113    #[allow(clippy::too_many_arguments)]
114    fn step(
115        &self,
116        function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
117        y: &mut Y,
118        t: &mut Scalar,
119        y_sol: &mut U,
120        t_sol: &mut Vector,
121        dydt_sol: &mut U,
122        dt: &mut Scalar,
123        k: &mut [Y],
124        y_trial: &Y,
125        e: Scalar,
126    ) -> Result<(), String>;
127    /// Provides the adaptive time step as a function of the error.
128    ///
129    /// ```math
130    /// h_{n+1} = \beta h \left(\frac{e_\mathrm{tol}}{e_{n+1}}\right)^{1/p}
131    /// ```
132    fn time_step(&self, error: Scalar, dt: &mut Scalar) {
133        if error > 0.0 {
134            *dt *= (self.dt_beta() * (self.abs_tol() / error).powf(1.0 / self.dt_expn()))
135                .max(self.dt_cut())
136        }
137    }
138}
conspire/math/integrate/ode/explicit/variable_step/mod.rs

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