conspire/math/integrate/dormand_prince/
mod.rs

1#[cfg(test)]
2mod test;
3
4use super::{
5    super::{Scalar, Tensor, TensorVec, Vector, interpolate::InterpolateSolution},
6    Explicit, IntegrationError, OdeSolver,
7};
8use crate::{ABS_TOL, REL_TOL};
9use std::ops::{Mul, Sub};
10
11const C_44_45: Scalar = 44.0 / 45.0;
12const C_56_15: Scalar = 56.0 / 15.0;
13const C_32_9: Scalar = 32.0 / 9.0;
14const C_8_9: Scalar = 8.0 / 9.0;
15const C_19372_6561: Scalar = 19372.0 / 6561.0;
16const C_25360_2187: Scalar = 25360.0 / 2187.0;
17const C_64448_6561: Scalar = 64448.0 / 6561.0;
18const C_212_729: Scalar = 212.0 / 729.0;
19const C_9017_3168: Scalar = 9017.0 / 3168.0;
20const C_355_33: Scalar = 355.0 / 33.0;
21const C_46732_5247: Scalar = 46732.0 / 5247.0;
22const C_49_176: Scalar = 49.0 / 176.0;
23const C_5103_18656: Scalar = 5103.0 / 18656.0;
24const C_35_384: Scalar = 35.0 / 384.0;
25const C_500_1113: Scalar = 500.0 / 1113.0;
26const C_125_192: Scalar = 125.0 / 192.0;
27const C_2187_6784: Scalar = 2187.0 / 6784.0;
28const C_11_84: Scalar = 11.0 / 84.0;
29const C_71_57600: Scalar = 71.0 / 57600.0;
30const C_71_16695: Scalar = 71.0 / 16695.0;
31const C_71_1920: Scalar = 71.0 / 1920.0;
32const C_17253_339200: Scalar = 17253.0 / 339200.0;
33const C_22_525: Scalar = 22.0 / 525.0;
34
35#[doc = include_str!("doc.md")]
36#[derive(Debug)]
37pub struct DormandPrince {
38    /// Absolute error tolerance.
39    pub abs_tol: Scalar,
40    /// Relative error tolerance.
41    pub rel_tol: Scalar,
42    /// Multiplier for adaptive time steps.
43    pub dt_beta: Scalar,
44    /// Exponent for adaptive time steps.
45    pub dt_expn: Scalar,
46    /// Cut back factor for the time step.
47    pub dt_cut: Scalar,
48    /// Minimum value for the time step.
49    pub dt_min: Scalar,
50}
51
52impl Default for DormandPrince {
53    fn default() -> Self {
54        Self {
55            abs_tol: ABS_TOL,
56            rel_tol: REL_TOL,
57            dt_beta: 0.9,
58            dt_expn: 5.0,
59            dt_cut: 0.5,
60            dt_min: ABS_TOL,
61        }
62    }
63}
64
65impl<Y, U> OdeSolver<Y, U> for DormandPrince
66where
67    Y: Tensor,
68    U: TensorVec<Item = Y>,
69{
70    fn abs_tol(&self) -> Scalar {
71        self.abs_tol
72    }
73    fn dt_cut(&self) -> Scalar {
74        self.dt_cut
75    }
76    fn dt_min(&self) -> Scalar {
77        self.dt_min
78    }
79}
80
81impl<Y, U> Explicit<Y, U> for DormandPrince
82where
83    Self: OdeSolver<Y, U>,
84    Y: Tensor,
85    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
86    U: TensorVec<Item = Y>,
87{
88    const SLOPES: usize = 7;
89    fn dt_beta(&self) -> Scalar {
90        self.dt_beta
91    }
92    fn dt_expn(&self) -> Scalar {
93        self.dt_expn
94    }
95    fn slopes(
96        &self,
97        mut function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
98        y: &Y,
99        t: Scalar,
100        dt: Scalar,
101        k: &mut [Y],
102        y_trial: &mut Y,
103    ) -> Result<Scalar, String> {
104        *y_trial = &k[0] * (0.2 * dt) + y;
105        k[1] = function(t + 0.2 * dt, y_trial)?;
106        *y_trial = &k[0] * (0.075 * dt) + &k[1] * (0.225 * dt) + y;
107        k[2] = function(t + 0.3 * dt, y_trial)?;
108        *y_trial = &k[0] * (C_44_45 * dt) - &k[1] * (C_56_15 * dt) + &k[2] * (C_32_9 * dt) + y;
109        k[3] = function(t + 0.8 * dt, y_trial)?;
110        *y_trial = &k[0] * (C_19372_6561 * dt) - &k[1] * (C_25360_2187 * dt)
111            + &k[2] * (C_64448_6561 * dt)
112            - &k[3] * (C_212_729 * dt)
113            + y;
114        k[4] = function(t + C_8_9 * dt, y_trial)?;
115        *y_trial = &k[0] * (C_9017_3168 * dt) - &k[1] * (C_355_33 * dt)
116            + &k[2] * (C_46732_5247 * dt)
117            + &k[3] * (C_49_176 * dt)
118            - &k[4] * (C_5103_18656 * dt)
119            + y;
120        k[5] = function(t + dt, y_trial)?;
121        *y_trial = (&k[0] * C_35_384 + &k[2] * C_500_1113 + &k[3] * C_125_192
122            - &k[4] * C_2187_6784
123            + &k[5] * C_11_84)
124            * dt
125            + y;
126        k[6] = function(t + dt, y_trial)?;
127        Ok(
128            ((&k[0] * C_71_57600 - &k[2] * C_71_16695 + &k[3] * C_71_1920
129                - &k[4] * C_17253_339200
130                + &k[5] * C_22_525
131                - &k[6] * 0.025)
132                * dt)
133                .norm_inf(),
134        )
135    }
136    fn step(
137        &self,
138        _function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
139        y: &mut Y,
140        t: &mut Scalar,
141        y_sol: &mut U,
142        t_sol: &mut Vector,
143        dydt_sol: &mut U,
144        dt: &mut Scalar,
145        k: &mut [Y],
146        y_trial: &Y,
147        e: Scalar,
148    ) -> Result<(), String> {
149        if e < self.abs_tol || e / y_trial.norm_inf() < self.rel_tol {
150            k[0] = k[6].clone();
151            *t += *dt;
152            *y = y_trial.clone();
153            t_sol.push(*t);
154            y_sol.push(y.clone());
155            dydt_sol.push(k[0].clone());
156        }
157        self.time_step(e, dt);
158        Ok(())
159    }
160}
161
162impl<Y, U> InterpolateSolution<Y, U> for DormandPrince
163where
164    Y: Tensor,
165    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
166    U: TensorVec<Item = Y>,
167{
168    fn interpolate(
169        &self,
170        time: &Vector,
171        tp: &Vector,
172        yp: &U,
173        mut function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
174    ) -> Result<(U, U), IntegrationError> {
175        let mut dt;
176        let mut i;
177        let mut k_1;
178        let mut k_2;
179        let mut k_3;
180        let mut k_4;
181        let mut k_5;
182        let mut k_6;
183        let mut t;
184        let mut y;
185        let mut y_int = U::new();
186        let mut dydt_int = U::new();
187        let mut y_trial;
188        for time_k in time.iter() {
189            i = tp.iter().position(|tp_i| tp_i >= time_k).unwrap();
190            if time_k == &tp[i] {
191                t = tp[i];
192                y_trial = yp[i].clone();
193                dt = 0.0;
194            } else {
195                t = tp[i - 1];
196                y = yp[i - 1].clone();
197                dt = time_k - t;
198                k_1 = function(t, &y)?;
199                y_trial = &k_1 * (0.2 * dt) + &y;
200                k_2 = function(t + 0.2 * dt, &y_trial)?;
201                y_trial = &k_1 * (0.075 * dt) + &k_2 * (0.225 * dt) + &y;
202                k_3 = function(t + 0.3 * dt, &y_trial)?;
203                y_trial = &k_1 * (C_44_45 * dt) - &k_2 * (C_56_15 * dt) + &k_3 * (C_32_9 * dt) + &y;
204                k_4 = function(t + 0.8 * dt, &y_trial)?;
205                y_trial = &k_1 * (C_19372_6561 * dt) - &k_2 * (C_25360_2187 * dt)
206                    + &k_3 * (C_64448_6561 * dt)
207                    - &k_4 * (C_212_729 * dt)
208                    + &y;
209                k_5 = function(t + C_8_9 * dt, &y_trial)?;
210                y_trial = &k_1 * (C_9017_3168 * dt) - &k_2 * (C_355_33 * dt)
211                    + &k_3 * (C_46732_5247 * dt)
212                    + &k_4 * (C_49_176 * dt)
213                    - &k_5 * (C_5103_18656 * dt)
214                    + &y;
215                k_6 = function(t + dt, &y_trial)?;
216                y_trial = (&k_1 * C_35_384 + &k_3 * C_500_1113 + &k_4 * C_125_192
217                    - &k_5 * C_2187_6784
218                    + &k_6 * C_11_84)
219                    * dt
220                    + &y;
221            }
222            dydt_int.push(function(t + dt, &y_trial)?);
223            y_int.push(y_trial);
224        }
225        Ok((y_int, dydt_int))
226    }
227}
conspire/math/integrate/dormand_prince/mod.rs

conspire/math/integrate/dormand_prince/
mod.rs
