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(crate) mod bogacki_shampine;
12pub(crate) mod dormand_prince;
13pub(crate) mod verner_8;
14pub(crate) 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 k_sol: Vec<U> = Vec::new();
49        let mut y_trial = Y::default();
50        while t < t_f {
51            match self.slopes_and_error(&mut function, &y, t, dt, &mut k, &mut y_trial) {
52                Ok(e) => {
53                    if let Err(error) = self.step(
54                        &mut function,
55                        &mut y,
56                        &mut t,
57                        &mut y_sol,
58                        &mut t_sol,
59                        &mut dydt_sol,
60                        &mut k_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) =
99                self.interpolate(&t_int, &t_sol, &y_sol, &dydt_sol, &k_sol, function)?;
100            Ok((t_int, y_int, dydt_int))
101        } else {
102            Ok((t_sol, y_sol, dydt_sol))
103        }
104    }
105    fn interpolate_variable_step(
106        time: &Vector,
107        tp: &Vector,
108        yp: &U,
109        mut function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
110    ) -> Result<(U, U), IntegrationError> {
111        let mut dt;
112        let mut i;
113        let mut k = vec![Y::default(); Self::SLOPES];
114        let mut t;
115        let mut y;
116        let mut y_int = U::new();
117        let mut dydt_int = U::new();
118        let mut y_trial = Y::default();
119        for time_k in time.iter() {
120            i = tp.iter().position(|tp_i| tp_i >= time_k).unwrap();
121            if time_k == &tp[i] {
122                t = tp[i];
123                y_trial = yp[i].clone();
124                dt = 0.0;
125            } else {
126                t = tp[i - 1];
127                y = &yp[i - 1];
128                dt = time_k - t;
129                k[0] = function(t, y)?;
130                Self::slopes(&mut function, y, t, dt, &mut k, &mut y_trial)?;
131            }
132            dydt_int.push(function(t + dt, &y_trial)?);
133            y_int.push(y_trial.clone());
134        }
135        Ok((y_int, dydt_int))
136    }
137    fn error(&self, dt: Scalar, k: &[Y]) -> Result<Scalar, String>;
138    fn slopes(
139        function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
140        y: &Y,
141        t: Scalar,
142        dt: Scalar,
143        k: &mut [Y],
144        y_trial: &mut Y,
145    ) -> Result<(), String>;
146    fn slopes_and_error(
147        &self,
148        mut function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
149        y: &Y,
150        t: Scalar,
151        dt: Scalar,
152        k: &mut [Y],
153        y_trial: &mut Y,
154    ) -> Result<Scalar, String> {
155        Self::slopes(&mut function, y, t, dt, k, y_trial)?;
156        self.error(dt, k)
157    }
158    #[allow(clippy::too_many_arguments)]
159    fn step(
160        &self,
161        mut function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
162        y: &mut Y,
163        t: &mut Scalar,
164        y_sol: &mut U,
165        t_sol: &mut Vector,
166        dydt_sol: &mut U,
167        k_sol: &mut Vec<U>,
168        dt: &mut Scalar,
169        k: &mut [Y],
170        y_trial: &Y,
171        e: Scalar,
172    ) -> Result<(), String> {
173        if e < self.abs_tol() || e < self.rel_tol() * self.norm().apply(y_trial) {
174            k_sol.push(k.iter().cloned().collect());
175            *t += *dt;
176            *y = y_trial.clone();
177            t_sol.push(*t);
178            y_sol.push(y.clone());
179            dydt_sol.push(function(*t, y)?);
180        }
181        self.time_step(e, dt);
182        Ok(())
183    }
184    /// Provides the adaptive time step as a function of the error.
185    ///
186    /// ```math
187    /// h_{n+1} = \beta h \left(\frac{e_\mathrm{tol}}{e_{n+1}}\right)^{1/p}
188    /// ```
189    fn time_step(&self, error: Scalar, dt: &mut Scalar) {
190        if error > 0.0 {
191            *dt *= (self.dt_beta() * (self.abs_tol() / error).powf(1.0 / self.dt_expn()))
192                .max(self.dt_cut())
193        }
194    }
195}
196
197/// Free (dense-output) interpolant for explicit ordinary differential equation integrators.
198///
199/// Uses cubic Hermite interpolation over the accepted-step values and derivatives already
200/// computed during integration, so it requires no additional evaluations of the right-hand side
201/// function.
202pub trait FreeInterpolant<Y, U>
203where
204    Self: VariableStepExplicit<Y, U>,
205    Y: Tensor,
206    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
207    U: TensorVec<Item = Y>,
208{
209    fn interpolate_free(time: &Vector, tp: &Vector, yp: &U, dydtp: &U) -> (U, U) {
210        let mut y_int = U::new();
211        let mut dydt_int = U::new();
212        for time_k in time.iter() {
213            let i = tp.iter().position(|tp_i| tp_i >= time_k).unwrap();
214            if time_k == &tp[i] {
215                y_int.push(yp[i].clone());
216                dydt_int.push(dydtp[i].clone());
217            } else {
218                let t_0 = tp[i - 1];
219                let h = tp[i] - t_0;
220                let theta = (time_k - t_0) / h;
221                let theta2 = theta * theta;
222                let theta3 = theta2 * theta;
223                let h00 = 2.0 * theta3 - 3.0 * theta2 + 1.0;
224                let h10 = theta3 - 2.0 * theta2 + theta;
225                let h01 = -2.0 * theta3 + 3.0 * theta2;
226                let h11 = theta3 - theta2;
227                let dh00 = 6.0 * theta2 - 6.0 * theta;
228                let dh10 = 3.0 * theta2 - 4.0 * theta + 1.0;
229                let dh01 = -6.0 * theta2 + 6.0 * theta;
230                let dh11 = 3.0 * theta2 - 2.0 * theta;
231                y_int.push(
232                    &yp[i - 1] * h00
233                        + &dydtp[i - 1] * (h10 * h)
234                        + &yp[i] * h01
235                        + &dydtp[i] * (h11 * h),
236                );
237                dydt_int.push(
238                    &yp[i - 1] * (dh00 / h)
239                        + &dydtp[i - 1] * dh10
240                        + &yp[i] * (dh01 / h)
241                        + &dydtp[i] * dh11,
242                );
243            }
244        }
245        (y_int, dydt_int)
246    }
247}
248
249/// First-same-as-last property for explicit ordinary differential equation integrators.
250pub trait VariableStepExplicitFirstSameAsLast<Y, U>
251where
252    Self: VariableStepExplicit<Y, U>,
253    Y: Tensor,
254    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
255    U: TensorVec<Item = Y>,
256{
257    fn slopes_and_error_fsal(
258        &self,
259        mut function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
260        y: &Y,
261        t: Scalar,
262        dt: Scalar,
263        k: &mut [Y],
264        y_trial: &mut Y,
265    ) -> Result<Scalar, String> {
266        Self::slopes(&mut function, y, t, dt, k, y_trial)?;
267        k[Self::SLOPES - 1] = function(t + dt, y_trial)?;
268        self.error(dt, k)
269    }
270    #[allow(clippy::too_many_arguments)]
271    fn step_fsal(
272        &self,
273        y: &mut Y,
274        t: &mut Scalar,
275        y_sol: &mut U,
276        t_sol: &mut Vector,
277        dydt_sol: &mut U,
278        k_sol: &mut Vec<U>,
279        dt: &mut Scalar,
280        k: &mut [Y],
281        y_trial: &Y,
282        e: Scalar,
283    ) -> Result<(), String> {
284        if e < self.abs_tol() || e < self.rel_tol() * self.norm().apply(y_trial) {
285            k_sol.push(k.iter().cloned().collect());
286            k[0] = k[Self::SLOPES - 1].clone();
287            *t += *dt;
288            *y = y_trial.clone();
289            t_sol.push(*t);
290            y_sol.push(y.clone());
291            dydt_sol.push(k[0].clone());
292        }
293        self.time_step(e, dt);
294        Ok(())
295    }
296}
conspire/math/integrate/ode/explicit/variable_step/mod.rs

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