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

1#[cfg(test)]
2mod test;
3
4use crate::math::Norm;
5use crate::math::{
6    Derivative, Differentiable, Quantity, Scalar, Tensor, TensorVec,
7    integrate::{
8        ButcherTableau, EmbeddedTableau, Explicit, IntegrationError, OdeIntegrator, Times,
9        VariableStep, VariableStepExplicit, VariableStepExplicitFirstSameAsLast,
10    },
11    interpolate::InterpolateSolution,
12};
13use crate::{ABS_TOL, REL_TOL};
14use std::ops::{Mul, Sub};
15
16pub(crate) const P_1_0: Scalar = 1.0;
17pub(crate) const P_1_1: Scalar = -8048581381.0 / 2820520608.0;
18pub(crate) const P_1_2: Scalar = 8663915743.0 / 2820520608.0;
19pub(crate) const P_1_3: Scalar = -12715105075.0 / 11282082432.0;
20pub(crate) const P_3_1: Scalar = 131558114200.0 / 32700410799.0;
21pub(crate) const P_3_2: Scalar = -68118460800.0 / 10900136933.0;
22pub(crate) const P_3_3: Scalar = 87487479700.0 / 32700410799.0;
23pub(crate) const P_4_1: Scalar = -1754552775.0 / 470086768.0;
24pub(crate) const P_4_2: Scalar = 14199869525.0 / 1410260304.0;
25pub(crate) const P_4_3: Scalar = -10690763975.0 / 1880347072.0;
26pub(crate) const P_5_1: Scalar = 127303824393.0 / 49829197408.0;
27pub(crate) const P_5_2: Scalar = -318862633887.0 / 49829197408.0;
28pub(crate) const P_5_3: Scalar = 701980252875.0 / 199316789632.0;
29pub(crate) const P_6_1: Scalar = -282668133.0 / 205662961.0;
30pub(crate) const P_6_2: Scalar = 2019193451.0 / 616988883.0;
31pub(crate) const P_6_3: Scalar = -1453857185.0 / 822651844.0;
32pub(crate) const P_7_1: Scalar = 40617522.0 / 29380423.0;
33pub(crate) const P_7_2: Scalar = -110615467.0 / 29380423.0;
34pub(crate) const P_7_3: Scalar = 69997945.0 / 29380423.0;
35
36/// The Dormand–Prince 5(4) tableau.
37#[derive(Debug)]
38pub struct Tableau;
39
40impl ButcherTableau for Tableau {
41    const STAGES: usize = 7;
42    const ORDER: Scalar = 5.0;
43    const A: &'static [&'static [Scalar]] = &[
44        &[],
45        &[0.2],
46        &[0.075, 0.225],
47        &[44.0 / 45.0, -56.0 / 15.0, 32.0 / 9.0],
48        &[
49            19372.0 / 6561.0,
50            -25360.0 / 2187.0,
51            64448.0 / 6561.0,
52            -212.0 / 729.0,
53        ],
54        &[
55            9017.0 / 3168.0,
56            -355.0 / 33.0,
57            46732.0 / 5247.0,
58            49.0 / 176.0,
59            -5103.0 / 18656.0,
60        ],
61        &[
62            35.0 / 384.0,
63            0.0,
64            500.0 / 1113.0,
65            125.0 / 192.0,
66            -2187.0 / 6784.0,
67            11.0 / 84.0,
68        ],
69    ];
70    const C: &'static [Scalar] = &[0.0, 0.2, 0.3, 0.8, 8.0 / 9.0, 1.0, 1.0];
71    const B: &'static [Scalar] = &[
72        35.0 / 384.0,
73        0.0,
74        500.0 / 1113.0,
75        125.0 / 192.0,
76        -2187.0 / 6784.0,
77        11.0 / 84.0,
78        0.0,
79    ];
80    const FSAL: bool = true;
81}
82
83impl EmbeddedTableau for Tableau {
84    const D: &'static [Scalar] = &[
85        71.0 / 57600.0,
86        0.0,
87        -71.0 / 16695.0,
88        71.0 / 1920.0,
89        -17253.0 / 339200.0,
90        22.0 / 525.0,
91        -0.025,
92    ];
93}
94
95#[doc = include_str!("doc.md")]
96#[derive(Debug)]
97pub struct DormandPrince {
98    /// Absolute error tolerance.
99    pub abs_tol: Scalar,
100    /// Relative error tolerance.
101    pub rel_tol: Scalar,
102    /// Multiplier for adaptive time steps.
103    pub dt_beta: Scalar,
104    /// Exponent for adaptive time steps.
105    pub dt_expn: Scalar,
106    /// Cut back factor for the time step.
107    pub dt_cut: Scalar,
108    /// Growth factor ceiling for the time step.
109    pub dt_grow: Scalar,
110    /// Minimum value for the time step.
111    pub dt_min: Scalar,
112    /// Norm type for error evaluation.
113    pub error_norm: Norm,
114}
115
116impl Default for DormandPrince {
117    fn default() -> Self {
118        Self {
119            abs_tol: ABS_TOL,
120            rel_tol: REL_TOL,
121            dt_beta: 0.9,
122            dt_expn: 5.0,
123            dt_cut: 0.5,
124            dt_grow: 5.0,
125            dt_min: ABS_TOL,
126            error_norm: Norm::Chebyshev,
127        }
128    }
129}
130
131impl<Y, U> OdeIntegrator<Y, U> for DormandPrince
132where
133    Y: Tensor,
134    U: TensorVec<Item = Y>,
135{
136}
137
138impl<T> VariableStep<T> for DormandPrince {
139    fn abs_tol(&self) -> Scalar {
140        self.abs_tol
141    }
142    fn rel_tol(&self) -> Scalar {
143        self.rel_tol
144    }
145    fn dt_beta(&self) -> Scalar {
146        self.dt_beta
147    }
148    fn dt_expn(&self) -> Scalar {
149        self.dt_expn
150    }
151    fn dt_cut(&self) -> Scalar {
152        self.dt_cut
153    }
154    fn dt_grow(&self) -> Scalar {
155        self.dt_grow
156    }
157    fn dt_min(&self) -> Quantity<T> {
158        Quantity::new(self.dt_min)
159    }
160    fn error_norm(&self) -> &Norm {
161        &self.error_norm
162    }
163}
164
165impl<Y, U, V, T> Explicit<Y, U, V, T> for DormandPrince
166where
167    Y: Differentiable<T> + Tensor,
168    Derivative<Y, T>: Mul<Quantity<T>, Output = Y>,
169    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
170    for<'a> &'a Derivative<Y, T>:
171        Mul<Scalar, Output = Derivative<Y, T>> + Mul<Quantity<T>, Output = Y>,
172    U: TensorVec<Item = Y>,
173    V: TensorVec<Item = Derivative<Y, T>>,
174{
175    const SLOPES: usize = 7;
176    fn integrate(
177        &self,
178        function: impl FnMut(Quantity<T>, &Y) -> Result<Derivative<Y, T>, String>,
179        time: &[Quantity<T>],
180        initial_condition: Y,
181    ) -> Result<(Times<T>, U, V), IntegrationError> {
182        self.integrate_variable_step(function, time, initial_condition)
183    }
184}
185
186impl<Y, U, V, T> VariableStepExplicit<Y, U, V, T> for DormandPrince
187where
188    Self: Explicit<Y, U, V, T>,
189    Y: Differentiable<T> + Tensor,
190    Derivative<Y, T>: Mul<Quantity<T>, Output = Y>,
191    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
192    for<'a> &'a Derivative<Y, T>:
193        Mul<Scalar, Output = Derivative<Y, T>> + Mul<Quantity<T>, Output = Y>,
194    U: TensorVec<Item = Y>,
195    V: TensorVec<Item = Derivative<Y, T>>,
196{
197    type Tableau = Tableau;
198    fn slopes_and_error(
199        &self,
200        function: impl FnMut(Quantity<T>, &Y) -> Result<Derivative<Y, T>, String>,
201        y: &Y,
202        t: Quantity<T>,
203        dt: Quantity<T>,
204        k: &mut [Derivative<Y, T>],
205        y_trial: &mut Y,
206    ) -> Result<Scalar, String> {
207        self.slopes_and_error_fsal(function, y, t, dt, k, y_trial)
208    }
209    fn step(
210        &self,
211        _function: impl FnMut(Quantity<T>, &Y) -> Result<Derivative<Y, T>, String>,
212        y: &mut Y,
213        t: &mut Quantity<T>,
214        y_sol: &mut U,
215        t_sol: &mut Times<T>,
216        dydt_sol: &mut V,
217        k_sol: &mut Vec<V>,
218        dt: &mut Quantity<T>,
219        k: &mut [Derivative<Y, T>],
220        y_trial: &Y,
221        e: Scalar,
222    ) -> Result<(), String> {
223        self.step_fsal(y, t, y_sol, t_sol, dydt_sol, k_sol, dt, k, y_trial, e)
224    }
225}
226
227impl<Y, U, V, T> VariableStepExplicitFirstSameAsLast<Y, U, V, T> for DormandPrince
228where
229    Y: Differentiable<T> + Tensor,
230    Derivative<Y, T>: Mul<Quantity<T>, Output = Y>,
231    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
232    for<'a> &'a Derivative<Y, T>:
233        Mul<Scalar, Output = Derivative<Y, T>> + Mul<Quantity<T>, Output = Y>,
234    U: TensorVec<Item = Y>,
235    V: TensorVec<Item = Derivative<Y, T>>,
236{
237}
238
239impl DormandPrince {
240    pub(crate) fn interpolate_free_dense<Y, U, V, T>(
241        time: &Times<T>,
242        tp: &Times<T>,
243        yp: &U,
244        dydtp: &V,
245        k_sol: &[V],
246    ) -> (U, V)
247    where
248        Y: Differentiable<T> + Tensor,
249        Derivative<Y, T>: Mul<Quantity<T>, Output = Y>,
250        for<'a> &'a Derivative<Y, T>:
251            Mul<Scalar, Output = Derivative<Y, T>> + Mul<Quantity<T>, Output = Y>,
252        U: TensorVec<Item = Y>,
253        V: TensorVec<Item = Derivative<Y, T>>,
254    {
255        let mut y_int = U::new();
256        let mut dydt_int = V::new();
257        for time_k in time.iter() {
258            let i = tp.iter().position(|tp_i| tp_i >= time_k).unwrap();
259            if time_k == &tp[i] {
260                y_int.push(yp[i].clone());
261                dydt_int.push(dydtp[i].clone());
262            } else {
263                let t_0 = tp[i - 1];
264                let h = tp[i] - t_0;
265                let theta = (*time_k - t_0).value() / h.value();
266                let theta2 = theta * theta;
267                let theta3 = theta2 * theta;
268                let theta4 = theta3 * theta;
269                let k = &k_sol[i - 1];
270                let c_1 = theta * P_1_0 + theta2 * P_1_1 + theta3 * P_1_2 + theta4 * P_1_3;
271                let c_3 = theta2 * P_3_1 + theta3 * P_3_2 + theta4 * P_3_3;
272                let c_4 = theta2 * P_4_1 + theta3 * P_4_2 + theta4 * P_4_3;
273                let c_5 = theta2 * P_5_1 + theta3 * P_5_2 + theta4 * P_5_3;
274                let c_6 = theta2 * P_6_1 + theta3 * P_6_2 + theta4 * P_6_3;
275                let c_7 = theta2 * P_7_1 + theta3 * P_7_2 + theta4 * P_7_3;
276                let dc_1 =
277                    P_1_0 + 2.0 * theta * P_1_1 + 3.0 * theta2 * P_1_2 + 4.0 * theta3 * P_1_3;
278                let dc_3 = 2.0 * theta * P_3_1 + 3.0 * theta2 * P_3_2 + 4.0 * theta3 * P_3_3;
279                let dc_4 = 2.0 * theta * P_4_1 + 3.0 * theta2 * P_4_2 + 4.0 * theta3 * P_4_3;
280                let dc_5 = 2.0 * theta * P_5_1 + 3.0 * theta2 * P_5_2 + 4.0 * theta3 * P_5_3;
281                let dc_6 = 2.0 * theta * P_6_1 + 3.0 * theta2 * P_6_2 + 4.0 * theta3 * P_6_3;
282                let dc_7 = 2.0 * theta * P_7_1 + 3.0 * theta2 * P_7_2 + 4.0 * theta3 * P_7_3;
283                let sum = &k[0] * c_1
284                    + &k[2] * c_3
285                    + &k[3] * c_4
286                    + &k[4] * c_5
287                    + &k[5] * c_6
288                    + &k[6] * c_7;
289                y_int.push(sum * h + &yp[i - 1]);
290                dydt_int.push(
291                    &k[0] * dc_1
292                        + &k[2] * dc_3
293                        + &k[3] * dc_4
294                        + &k[4] * dc_5
295                        + &k[5] * dc_6
296                        + &k[6] * dc_7,
297                );
298            }
299        }
300        (y_int, dydt_int)
301    }
302}
303
304impl<Y, U, V, T> InterpolateSolution<Y, U, V, T> for DormandPrince
305where
306    Y: Differentiable<T> + Tensor,
307    Derivative<Y, T>: Mul<Quantity<T>, Output = Y>,
308    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
309    for<'a> &'a Derivative<Y, T>:
310        Mul<Scalar, Output = Derivative<Y, T>> + Mul<Quantity<T>, Output = Y>,
311    U: TensorVec<Item = Y>,
312    V: TensorVec<Item = Derivative<Y, T>>,
313{
314    fn interpolate(
315        &self,
316        time: &Times<T>,
317        tp: &Times<T>,
318        yp: &U,
319        dydtp: &V,
320        k_sol: &[V],
321        _function: impl FnMut(Quantity<T>, &Y) -> Result<Derivative<Y, T>, String>,
322    ) -> Result<(U, V), IntegrationError> {
323        Ok(Self::interpolate_free_dense(time, tp, yp, dydtp, k_sol))
324    }
325}
conspire/math/integrate/ode/explicit/variable_step/dormand_prince/mod.rs

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