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 bogacki_shampine;
12pub mod dormand_prince;
13pub mod verner_8;
14pub mod verner_9;
15
16/// Variable-step explicit integrators for ordinary differential equations.
17pub trait VariableStepExplicit<Y, U>
18where
19    Self: InterpolateSolution<Y, U> + Explicit<Y, U> + VariableStep,
20    Y: Tensor,
21    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
22    U: TensorVec<Item = Y>,
23{
24    fn integrate_variable_step(
25        &self,
26        mut function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
27        time: &[Scalar],
28        initial_condition: Y,
29    ) -> Result<(Vector, U, U), IntegrationError> {
30        let t_0 = time[0];
31        let t_f = time[time.len() - 1];
32        if time.len() < 2 {
33            return Err(IntegrationError::LengthTimeLessThanTwo);
34        } else if t_0 >= t_f {
35            return Err(IntegrationError::InitialTimeNotLessThanFinalTime);
36        }
37        let mut t = t_0;
38        let mut dt = t_f - t_0;
39        let mut k = vec![Y::default(); Self::SLOPES];
40        k[0] = function(t, &initial_condition)?;
41        let mut t_sol = Vector::new();
42        t_sol.push(t_0);
43        let mut y = initial_condition.clone();
44        let mut y_sol = U::new();
45        y_sol.push(initial_condition.clone());
46        let mut dydt_sol = U::new();
47        dydt_sol.push(k[0].clone());
48        let mut y_trial = Y::default();
49        while t < t_f {
50            match self.slopes_and_error(&mut function, &y, t, dt, &mut k, &mut y_trial) {
51                Ok(e) => {
52                    if let Err(error) = self.step(
53                        &mut function,
54                        &mut y,
55                        &mut t,
56                        &mut y_sol,
57                        &mut t_sol,
58                        &mut dydt_sol,
59                        &mut dt,
60                        &mut k,
61                        &y_trial,
62                        e,
63                    ) {
64                        dt *= self.dt_cut();
65                        if dt < self.dt_min() {
66                            return Err(IntegrationError::MinimumStepSizeUpstream(
67                                self.dt_min(),
68                                error,
69                                format!("{self:?}"),
70                            ));
71                        }
72                    } else {
73                        dt = dt.min(t_f - t);
74                        if dt < self.dt_min() && t < t_f {
75                            return Err(IntegrationError::MinimumStepSizeReached(
76                                self.dt_min(),
77                                format!("{self:?}"),
78                            ));
79                        }
80                    }
81                }
82                Err(error) => {
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                }
92            }
93        }
94        if time.len() > 2 {
95            let t_int = Vector::from(time);
96            let (y_int, dydt_int) = self.interpolate(&t_int, &t_sol, &y_sol, function)?;
97            Ok((t_int, y_int, dydt_int))
98        } else {
99            Ok((t_sol, y_sol, dydt_sol))
100        }
101    }
102    fn interpolate_variable_step(
103        time: &Vector,
104        tp: &Vector,
105        yp: &U,
106        mut function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
107    ) -> Result<(U, U), IntegrationError> {
108        let mut dt;
109        let mut i;
110        let mut k = vec![Y::default(); Self::SLOPES];
111        let mut t;
112        let mut y;
113        let mut y_int = U::new();
114        let mut dydt_int = U::new();
115        let mut y_trial = Y::default();
116        for time_k in time.iter() {
117            i = tp.iter().position(|tp_i| tp_i >= time_k).unwrap();
118            if time_k == &tp[i] {
119                t = tp[i];
120                y_trial = yp[i].clone();
121                dt = 0.0;
122            } else {
123                t = tp[i - 1];
124                y = &yp[i - 1];
125                dt = time_k - t;
126                k[0] = function(t, y)?;
127                Self::slopes(&mut function, y, t, dt, &mut k, &mut y_trial)?;
128            }
129            dydt_int.push(function(t + dt, &y_trial)?);
130            y_int.push(y_trial.clone());
131        }
132        Ok((y_int, dydt_int))
133    }
134    fn error(dt: Scalar, k: &[Y]) -> Result<Scalar, String>;
135    fn slopes(
136        function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
137        y: &Y,
138        t: Scalar,
139        dt: Scalar,
140        k: &mut [Y],
141        y_trial: &mut Y,
142    ) -> Result<(), String>;
143    fn slopes_and_error(
144        &self,
145        mut function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
146        y: &Y,
147        t: Scalar,
148        dt: Scalar,
149        k: &mut [Y],
150        y_trial: &mut Y,
151    ) -> Result<Scalar, String> {
152        Self::slopes(&mut function, y, t, dt, k, y_trial)?;
153        Self::error(dt, k)
154    }
155    #[allow(clippy::too_many_arguments)]
156    fn step(
157        &self,
158        mut function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
159        y: &mut Y,
160        t: &mut Scalar,
161        y_sol: &mut U,
162        t_sol: &mut Vector,
163        dydt_sol: &mut U,
164        dt: &mut Scalar,
165        _k: &mut [Y],
166        y_trial: &Y,
167        e: Scalar,
168    ) -> Result<(), String> {
169        if e < self.abs_tol() || e / y_trial.norm_inf() < self.rel_tol() {
170            *t += *dt;
171            *y = y_trial.clone();
172            t_sol.push(*t);
173            y_sol.push(y.clone());
174            dydt_sol.push(function(*t, y)?);
175        }
176        self.time_step(e, dt);
177        Ok(())
178    }
179    /// Provides the adaptive time step as a function of the error.
180    ///
181    /// ```math
182    /// h_{n+1} = \beta h \left(\frac{e_\mathrm{tol}}{e_{n+1}}\right)^{1/p}
183    /// ```
184    fn time_step(&self, error: Scalar, dt: &mut Scalar) {
185        if error > 0.0 {
186            *dt *= (self.dt_beta() * (self.abs_tol() / error).powf(1.0 / self.dt_expn()))
187                .max(self.dt_cut())
188        }
189    }
190}
191
192/// First-same-as-last property for explicit ordinary differential equation integrators.
193pub trait VariableStepExplicitFirstSameAsLast<Y, U>
194where
195    Self: VariableStepExplicit<Y, U>,
196    Y: Tensor,
197    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
198    U: TensorVec<Item = Y>,
199{
200    fn slopes_and_error_fsal(
201        mut function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
202        y: &Y,
203        t: Scalar,
204        dt: Scalar,
205        k: &mut [Y],
206        y_trial: &mut Y,
207    ) -> Result<Scalar, String> {
208        Self::slopes(&mut function, y, t, dt, k, y_trial)?;
209        k[Self::SLOPES - 1] = function(t + dt, y_trial)?;
210        Self::error(dt, k)
211    }
212    #[allow(clippy::too_many_arguments)]
213    fn step_fsal(
214        &self,
215        y: &mut Y,
216        t: &mut Scalar,
217        y_sol: &mut U,
218        t_sol: &mut Vector,
219        dydt_sol: &mut U,
220        dt: &mut Scalar,
221        k: &mut [Y],
222        y_trial: &Y,
223        e: Scalar,
224    ) -> Result<(), String> {
225        if e < self.abs_tol() || e / y_trial.norm_inf() < self.rel_tol() {
226            k[0] = k[Self::SLOPES - 1].clone();
227            *t += *dt;
228            *y = y_trial.clone();
229            t_sol.push(*t);
230            y_sol.push(y.clone());
231            dydt_sol.push(k[0].clone());
232        }
233        self.time_step(e, dt);
234        Ok(())
235    }
236}
conspire/math/integrate/ode/explicit/variable_step/mod.rs

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