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

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

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