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

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