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

1use crate::math::{
2    Scalar, Tensor, TensorVec, Vector, assert_eq_within_tols,
3    integrate::{Explicit, IntegrationError, VariableStep},
4    interpolate::InterpolateSolutionInternalVariables,
5};
6use std::ops::{Mul, Sub};
7
8/// Variable-step explicit ordinary differential equation solvers with internal variables.
9pub trait VariableStepExplicitInternalVariables<Y, Z, U, V>
10where
11    Self: InterpolateSolutionInternalVariables<Y, Z, U, V> + Explicit<Y, U> + VariableStep,
12    Y: Tensor,
13    Z: Tensor,
14    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
15    U: TensorVec<Item = Y>,
16    V: TensorVec<Item = Z>,
17{
18    fn integrate_and_evaluate_variable_step(
19        &self,
20        mut function: impl FnMut(Scalar, &Y, &Z) -> Result<Y, String>,
21        mut evaluate: impl FnMut(Scalar, &Y, &Z) -> Result<Z, String>,
22        time: &[Scalar],
23        initial_condition: Y,
24        initial_evaluation: Z,
25    ) -> Result<(Vector, U, U, V), IntegrationError> {
26        let t_0 = time[0];
27        let t_f = time[time.len() - 1];
28        if time.len() < 2 {
29            return Err(IntegrationError::LengthTimeLessThanTwo);
30        } else if t_0 >= t_f {
31            return Err(IntegrationError::InitialTimeNotLessThanFinalTime);
32        }
33        let mut t = t_0;
34        let mut dt = t_f - t_0;
35        let mut t_sol = Vector::new();
36        t_sol.push(t_0);
37        let mut y = initial_condition;
38        let mut z = initial_evaluation;
39        if assert_eq_within_tols(&evaluate(t, &y, &z)?, &z).is_err() {
40            return Err(IntegrationError::InconsistentInitialConditions);
41        }
42        let mut k = vec![Y::default(); Self::SLOPES];
43        k[0] = function(t, &y, &z)?;
44        let mut y_sol = U::new();
45        y_sol.push(y.clone());
46        let mut z_sol = V::new();
47        z_sol.push(z.clone());
48        let mut dydt_sol = U::new();
49        dydt_sol.push(k[0].clone());
50        let mut y_trial = Y::default();
51        let mut z_trial = Z::default();
52        while t < t_f {
53            match self.slopes(
54                &mut function,
55                &mut evaluate,
56                &y,
57                &z,
58                t,
59                dt,
60                &mut k,
61                &mut y_trial,
62                &mut z_trial,
63            ) {
64                Ok(e) => {
65                    if let Some(error) = self
66                        .step(
67                            &mut function,
68                            &mut y,
69                            &mut z,
70                            &mut t,
71                            &mut y_sol,
72                            &mut z_sol,
73                            &mut t_sol,
74                            &mut dydt_sol,
75                            &mut dt,
76                            &mut k,
77                            &y_trial,
78                            &z_trial,
79                            e,
80                        )
81                        .err()
82                    {
83                        dt *= self.dt_cut();
84                        if dt < self.dt_min() {
85                            return Err(IntegrationError::MinimumStepSizeUpstream(
86                                self.dt_min(),
87                                error,
88                                format!("{self:?}"),
89                            ));
90                        }
91                    } else {
92                        dt = dt.min(t_f - t);
93                        if dt < self.dt_min() && t < t_f {
94                            return Err(IntegrationError::MinimumStepSizeReached(
95                                self.dt_min(),
96                                format!("{self:?}"),
97                            ));
98                        }
99                    }
100                }
101                Err(error) => {
102                    dt *= self.dt_cut();
103                    if dt < self.dt_min() {
104                        return Err(IntegrationError::MinimumStepSizeUpstream(
105                            self.dt_min(),
106                            error,
107                            format!("{self:?}"),
108                        ));
109                    }
110                }
111            }
112        }
113        if time.len() > 2 {
114            let t_int = Vector::from(time);
115            let (y_int, dydt_int, z_int) =
116                self.interpolate(&t_int, &t_sol, &y_sol, &z_sol, function, evaluate)?;
117            Ok((t_int, y_int, dydt_int, z_int))
118        } else {
119            Ok((t_sol, y_sol, dydt_sol, z_sol))
120        }
121    }
122    #[allow(clippy::too_many_arguments)]
123    fn slopes(
124        &self,
125        function: impl FnMut(Scalar, &Y, &Z) -> Result<Y, String>,
126        evaluate: impl FnMut(Scalar, &Y, &Z) -> Result<Z, String>,
127        y: &Y,
128        z: &Z,
129        t: Scalar,
130        dt: Scalar,
131        k: &mut [Y],
132        y_trial: &mut Y,
133        z_trial: &mut Z,
134    ) -> Result<Scalar, String>;
135    #[allow(clippy::too_many_arguments)]
136    fn step(
137        &self,
138        function: impl FnMut(Scalar, &Y, &Z) -> Result<Y, String>,
139        y: &mut Y,
140        z: &mut Z,
141        t: &mut Scalar,
142        y_sol: &mut U,
143        z_sol: &mut V,
144        t_sol: &mut Vector,
145        dydt_sol: &mut U,
146        dt: &mut Scalar,
147        k: &mut [Y],
148        y_trial: &Y,
149        z_trial: &Z,
150        e: Scalar,
151    ) -> Result<(), String>;
152}
conspire/math/integrate/ode/explicit/variable_step/internal_variables/mod.rs

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