conspire/math/integrate/dormand_prince/
mod.rs

1#[cfg(test)]
2mod test;
3
4use super::{
5    super::{Tensor, TensorRank0, TensorVec, Vector, interpolate::InterpolateSolution},
6    Explicit, IntegrationError,
7};
8use crate::{ABS_TOL, REL_TOL};
9use std::ops::{Mul, Sub};
10
11const C_44_45: TensorRank0 = 44.0 / 45.0;
12const C_56_15: TensorRank0 = 56.0 / 15.0;
13const C_32_9: TensorRank0 = 32.0 / 9.0;
14const C_8_9: TensorRank0 = 8.0 / 9.0;
15const C_19372_6561: TensorRank0 = 19372.0 / 6561.0;
16const C_25360_2187: TensorRank0 = 25360.0 / 2187.0;
17const C_64448_6561: TensorRank0 = 64448.0 / 6561.0;
18const C_212_729: TensorRank0 = 212.0 / 729.0;
19const C_9017_3168: TensorRank0 = 9017.0 / 3168.0;
20const C_355_33: TensorRank0 = 355.0 / 33.0;
21const C_46732_5247: TensorRank0 = 46732.0 / 5247.0;
22const C_49_176: TensorRank0 = 49.0 / 176.0;
23const C_5103_18656: TensorRank0 = 5103.0 / 18656.0;
24const C_35_384: TensorRank0 = 35.0 / 384.0;
25const C_500_1113: TensorRank0 = 500.0 / 1113.0;
26const C_125_192: TensorRank0 = 125.0 / 192.0;
27const C_2187_6784: TensorRank0 = 2187.0 / 6784.0;
28const C_11_84: TensorRank0 = 11.0 / 84.0;
29const C_71_57600: TensorRank0 = 71.0 / 57600.0;
30const C_71_16695: TensorRank0 = 71.0 / 16695.0;
31const C_71_1920: TensorRank0 = 71.0 / 1920.0;
32const C_17253_339200: TensorRank0 = 17253.0 / 339200.0;
33const C_22_525: TensorRank0 = 22.0 / 525.0;
34
35/// Explicit, six-stage, fifth-order, variable-step, Runge-Kutta method.[^cite]
36///
37/// [^cite]: J.R. Dormand and P.J. Prince, [J. Comput. Appl. Math. **6**, 19 (1980)](https://doi.org/10.1016/0771-050X(80)90013-3).
38///
39/// ```math
40/// \frac{dy}{dt} = f(t, y)
41/// ```
42/// ```math
43/// t_{n+1} = t_n + h
44/// ```
45/// ```math
46/// k_1 = f(t_n, y_n)
47/// ```
48/// ```math
49/// k_2 = f(t_n + \tfrac{1}{5} h, y_n + \tfrac{1}{5} h k_1)
50/// ```
51/// ```math
52/// k_3 = f(t_n + \tfrac{3}{10} h, y_n + \tfrac{3}{40} h k_1 + \tfrac{9}{40} h k_2)
53/// ```
54/// ```math
55/// k_4 = f(t_n + \tfrac{4}{5} h, y_n + \tfrac{44}{45} h k_1 - \tfrac{56}{15} h k_2 + \tfrac{32}{9} h k_3)
56/// ```
57/// ```math
58/// k_5 = f(t_n + \tfrac{8}{9} h, y_n + \tfrac{19372}{6561} h k_1 - \tfrac{25360}{2187} h k_2 + \tfrac{64448}{6561} h k_3 - \tfrac{212}{729} h k_4)
59/// ```
60/// ```math
61/// k_6 = f(t_n + h, y_n + \tfrac{9017}{3168} h k_1 - \tfrac{355}{33} h k_2 - \tfrac{46732}{5247} h k_3 + \tfrac{49}{176} h k_4 - \tfrac{5103}{18656} h k_5)
62/// ```
63/// ```math
64/// y_{n+1} = y_n + h\left(\frac{35}{384}\,k_1 + \frac{500}{1113}\,k_3 + \frac{125}{192}\,k_4 - \frac{2187}{6784}\,k_5 + \frac{11}{84}\,k_6\right)
65/// ```
66/// ```math
67/// k_7 = f(t_{n+1}, y_{n+1})
68/// ```
69/// ```math
70/// e_{n+1} = \frac{h}{5}\left(\frac{71}{11520}\,k_1 - \frac{71}{3339}\,k_3 + \frac{71}{384}\,k_4 - \frac{17253}{67840}\,k_5 + \frac{22}{105}\,k_6 - \frac{1}{8}\,k_7\right)
71/// ```
72/// ```math
73/// h_{n+1} = \beta h \left(\frac{e_\mathrm{tol}}{e_{n+1}}\right)^{1/p}
74/// ```
75#[derive(Debug)]
76pub struct DormandPrince {
77    /// Absolute error tolerance.
78    pub abs_tol: TensorRank0,
79    /// Relative error tolerance.
80    pub rel_tol: TensorRank0,
81    /// Multiplier for adaptive time steps.
82    pub dt_beta: TensorRank0,
83    /// Exponent for adaptive time steps.
84    pub dt_expn: TensorRank0,
85    /// Initial relative time step.
86    pub dt_init: TensorRank0,
87}
88
89impl Default for DormandPrince {
90    fn default() -> Self {
91        Self {
92            abs_tol: ABS_TOL,
93            rel_tol: REL_TOL,
94            dt_beta: 0.9,
95            dt_expn: 5.0,
96            dt_init: 0.1,
97        }
98    }
99}
100
101impl<Y, U> Explicit<Y, U> for DormandPrince
102where
103    Self: InterpolateSolution<Y, U>,
104    Y: Tensor,
105    for<'a> &'a Y: Mul<TensorRank0, Output = Y> + Sub<&'a Y, Output = Y>,
106    U: TensorVec<Item = Y>,
107{
108    fn integrate(
109        &self,
110        function: impl Fn(TensorRank0, &Y) -> Result<Y, IntegrationError>,
111        time: &[TensorRank0],
112        initial_condition: Y,
113    ) -> Result<(Vector, U, U), IntegrationError> {
114        if time.len() < 2 {
115            return Err(IntegrationError::LengthTimeLessThanTwo);
116        } else if time[0] >= time[time.len() - 1] {
117            return Err(IntegrationError::InitialTimeNotLessThanFinalTime);
118        }
119        let mut t = time[0];
120        let mut dt = self.dt_init * time[time.len() - 1];
121        let mut e;
122        let mut k_1 = function(t, &initial_condition)?;
123        let mut k_2;
124        let mut k_3;
125        let mut k_4;
126        let mut k_5;
127        let mut k_6;
128        let mut k_7;
129        let mut t_sol = Vector::zero(0);
130        t_sol.push(time[0]);
131        let mut y = initial_condition.clone();
132        let mut y_sol = U::zero(0);
133        y_sol.push(initial_condition.clone());
134        let mut dydt_sol = U::zero(0);
135        dydt_sol.push(k_1.clone());
136        let mut y_trial;
137        while t < time[time.len() - 1] {
138            k_2 = function(t + 0.2 * dt, &(&k_1 * (0.2 * dt) + &y))?;
139            k_3 = function(
140                t + 0.3 * dt,
141                &(&k_1 * (0.075 * dt) + &k_2 * (0.225 * dt) + &y),
142            )?;
143            k_4 = function(
144                t + 0.8 * dt,
145                &(&k_1 * (C_44_45 * dt) - &k_2 * (C_56_15 * dt) + &k_3 * (C_32_9 * dt) + &y),
146            )?;
147            k_5 = function(
148                t + C_8_9 * dt,
149                &(&k_1 * (C_19372_6561 * dt) - &k_2 * (C_25360_2187 * dt)
150                    + &k_3 * (C_64448_6561 * dt)
151                    - &k_4 * (C_212_729 * dt)
152                    + &y),
153            )?;
154            k_6 = function(
155                t + dt,
156                &(&k_1 * (C_9017_3168 * dt) - &k_2 * (C_355_33 * dt)
157                    + &k_3 * (C_46732_5247 * dt)
158                    + &k_4 * (C_49_176 * dt)
159                    - &k_5 * (C_5103_18656 * dt)
160                    + &y),
161            )?;
162            y_trial = (&k_1 * C_35_384 + &k_3 * C_500_1113 + &k_4 * C_125_192 - &k_5 * C_2187_6784
163                + &k_6 * C_11_84)
164                * dt
165                + &y;
166            k_7 = function(t + dt, &y_trial)?;
167            e = ((&k_1 * C_71_57600 - k_3 * C_71_16695 + k_4 * C_71_1920 - k_5 * C_17253_339200
168                + k_6 * C_22_525
169                - &k_7 * 0.025)
170                * dt)
171                .norm_inf();
172            if e < self.abs_tol || e / y_trial.norm_inf() < self.rel_tol {
173                k_1 = k_7;
174                t += dt;
175                y = y_trial;
176                t_sol.push(t);
177                y_sol.push(y.clone());
178                dydt_sol.push(k_1.clone());
179            }
180            dt *= self.dt_beta * (self.abs_tol / e).powf(1.0 / self.dt_expn);
181        }
182        if time.len() > 2 {
183            let t_int = Vector::new(time);
184            let (y_int, dydt_int) = self.interpolate(&t_int, &t_sol, &y_sol, function)?;
185            Ok((t_int, y_int, dydt_int))
186        } else {
187            Ok((t_sol, y_sol, dydt_sol))
188        }
189    }
190}
191
192impl<Y, U> InterpolateSolution<Y, U> for DormandPrince
193where
194    Y: Tensor,
195    for<'a> &'a Y: Mul<TensorRank0, Output = Y> + Sub<&'a Y, Output = Y>,
196    U: TensorVec<Item = Y>,
197{
198    fn interpolate(
199        &self,
200        time: &Vector,
201        tp: &Vector,
202        yp: &U,
203        function: impl Fn(TensorRank0, &Y) -> Result<Y, IntegrationError>,
204    ) -> Result<(U, U), IntegrationError> {
205        let mut dt;
206        let mut i;
207        let mut k_1;
208        let mut k_2;
209        let mut k_3;
210        let mut k_4;
211        let mut k_5;
212        let mut k_6;
213        let mut t;
214        let mut y;
215        let mut y_int = U::zero(0);
216        let mut dydt_int = U::zero(0);
217        let mut y_trial;
218        for time_k in time.iter() {
219            i = tp.iter().position(|tp_i| tp_i > time_k).unwrap();
220            t = tp[i - 1];
221            y = yp[i - 1].clone();
222            dt = time_k - t;
223            k_1 = function(t, &y)?;
224            k_2 = function(t + 0.2 * dt, &(&k_1 * (0.2 * dt) + &y))?;
225            k_3 = function(
226                t + 0.3 * dt,
227                &(&k_1 * (0.075 * dt) + &k_2 * (0.225 * dt) + &y),
228            )?;
229            k_4 = function(
230                t + 0.8 * dt,
231                &(&k_1 * (C_44_45 * dt) - &k_2 * (C_56_15 * dt) + &k_3 * (C_32_9 * dt) + &y),
232            )?;
233            k_5 = function(
234                t + C_8_9 * dt,
235                &(&k_1 * (C_19372_6561 * dt) - &k_2 * (C_25360_2187 * dt)
236                    + &k_3 * (C_64448_6561 * dt)
237                    - &k_4 * (C_212_729 * dt)
238                    + &y),
239            )?;
240            k_6 = function(
241                t + dt,
242                &(&k_1 * (C_9017_3168 * dt) - &k_2 * (C_355_33 * dt)
243                    + &k_3 * (C_46732_5247 * dt)
244                    + &k_4 * (C_49_176 * dt)
245                    - &k_5 * (C_5103_18656 * dt)
246                    + &y),
247            )?;
248            y_trial = (&k_1 * C_35_384 + &k_3 * C_500_1113 + &k_4 * C_125_192 - &k_5 * C_2187_6784
249                + &k_6 * C_11_84)
250                * dt
251                + &y;
252            dydt_int.push(function(t + dt, &y_trial)?);
253            y_int.push(y_trial);
254        }
255        Ok((y_int, dydt_int))
256    }
257}
conspire/math/integrate/dormand_prince/mod.rs

conspire/math/integrate/dormand_prince/
mod.rs
