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

1#[cfg(test)]
2mod test;
3
4use crate::math::{
5    Scalar, Tensor, TensorVec, Vector,
6    integrate::{Explicit, IntegrationError, OdeIntegrator, VariableStep, VariableStepExplicit},
7    interpolate::InterpolateSolution,
8};
9use crate::{ABS_TOL, REL_TOL};
10use std::ops::{Mul, Sub};
11
12pub const C_2: Scalar = 0.03462;
13pub const C_3: Scalar = 0.097_024_350_638_780_44;
14pub const C_4: Scalar = 0.145_536_525_958_170_67;
15pub const C_5: Scalar = 0.561;
16pub const C_6: Scalar = 0.229_007_911_590_485;
17pub const C_7: Scalar = 0.544_992_088_409_515;
18pub const C_8: Scalar = 0.645;
19pub const C_9: Scalar = 0.48375;
20pub const C_10: Scalar = 0.06757;
21pub const C_11: Scalar = 0.2500;
22pub const C_12: Scalar = 0.659_065_061_873_099_9;
23pub const C_13: Scalar = 0.8206;
24pub const C_14: Scalar = 0.9012;
25
26pub const A_2_1: Scalar = 0.03462;
27pub const A_3_1: Scalar = -0.03893354388572875;
28pub const A_3_2: Scalar = 0.13595789452450918;
29pub const A_4_1: Scalar = 0.03638413148954267;
30pub const A_4_3: Scalar = 0.10915239446862801;
31pub const A_5_1: Scalar = 2.0257639143939694;
32pub const A_5_3: Scalar = -7.638023836496291;
33pub const A_5_4: Scalar = 6.173259922102322;
34pub const A_6_1: Scalar = 0.05112275589406061;
35pub const A_6_4: Scalar = 0.17708237945550218;
36pub const A_6_5: Scalar = 0.0008027762409222536;
37pub const A_7_1: Scalar = 0.13160063579752163;
38pub const A_7_4: Scalar = -0.2957276252669636;
39pub const A_7_5: Scalar = 0.08781378035642955;
40pub const A_7_6: Scalar = 0.6213052975225274;
41pub const A_8_1: Scalar = 0.07166666666666667;
42pub const A_8_6: Scalar = 0.33055335789153195;
43pub const A_8_7: Scalar = 0.2427799754418014;
44pub const A_9_1: Scalar = 0.071806640625;
45pub const A_9_6: Scalar = 0.3294380283228177;
46pub const A_9_7: Scalar = 0.1165190029271823;
47pub const A_9_8: Scalar = -0.034013671875;
48pub const A_10_1: Scalar = 0.04836757646340646;
49pub const A_10_6: Scalar = 0.03928989925676164;
50pub const A_10_7: Scalar = 0.10547409458903446;
51pub const A_10_8: Scalar = -0.021438652846483126;
52pub const A_10_9: Scalar = -0.10412291746271944;
53pub const A_11_1: Scalar = -0.026645614872014785;
54pub const A_11_6: Scalar = 0.03333333333333333;
55pub const A_11_7: Scalar = -0.1631072244872467;
56pub const A_11_8: Scalar = 0.03396081684127761;
57pub const A_11_9: Scalar = 0.1572319413814626;
58pub const A_11_10: Scalar = 0.21522674780318796;
59pub const A_12_1: Scalar = 0.03689009248708622;
60pub const A_12_6: Scalar = -0.1465181576725543;
61pub const A_12_7: Scalar = 0.2242577768172024;
62pub const A_12_8: Scalar = 0.02294405717066073;
63pub const A_12_9: Scalar = -0.0035850052905728597;
64pub const A_12_10: Scalar = 0.08669223316444385;
65pub const A_12_11: Scalar = 0.43838406519683376;
66pub const A_13_1: Scalar = -0.4866012215113341;
67pub const A_13_6: Scalar = -6.304602650282853;
68pub const A_13_7: Scalar = -0.2812456182894729;
69pub const A_13_8: Scalar = -2.679019236219849;
70pub const A_13_9: Scalar = 0.5188156639241577;
71pub const A_13_10: Scalar = 1.3653531876033418;
72pub const A_13_11: Scalar = 5.8850910885039465;
73pub const A_13_12: Scalar = 2.8028087862720628;
74pub const A_14_1: Scalar = 0.4185367457753472;
75pub const A_14_6: Scalar = 6.724547581906459;
76pub const A_14_7: Scalar = -0.42544428016461133;
77pub const A_14_8: Scalar = 3.3432791530012653;
78pub const A_14_9: Scalar = 0.6170816631175374;
79pub const A_14_10: Scalar = -0.9299661239399329;
80pub const A_14_11: Scalar = -6.099948804751011;
81pub const A_14_12: Scalar = -3.002206187889399;
82pub const A_14_13: Scalar = 0.2553202529443446;
83pub const A_15_1: Scalar = -0.7793740861228848;
84pub const A_15_6: Scalar = -13.937342538107776;
85pub const A_15_7: Scalar = 1.2520488533793563;
86pub const A_15_8: Scalar = -14.691500408016868;
87pub const A_15_9: Scalar = -0.494705058533141;
88pub const A_15_10: Scalar = 2.2429749091462368;
89pub const A_15_11: Scalar = 13.367893803828643;
90pub const A_15_12: Scalar = 14.396650486650687;
91pub const A_15_13: Scalar = -0.79758133317768;
92pub const A_15_14: Scalar = 0.4409353709534278;
93pub const A_16_1: Scalar = 2.0580513374668867;
94pub const A_16_6: Scalar = 22.357937727968032;
95pub const A_16_7: Scalar = 0.9094981099755646;
96pub const A_16_8: Scalar = 35.89110098240264;
97pub const A_16_9: Scalar = -3.442515027624454;
98pub const A_16_10: Scalar = -4.865481358036369;
99pub const A_16_11: Scalar = -18.909803813543427;
100pub const A_16_12: Scalar = -34.26354448030452;
101pub const A_16_13: Scalar = 1.2647565216956427;
102
103pub const B_1: Scalar = 0.014611976858423152;
104pub const B_8: Scalar = -0.3915211862331339;
105pub const B_9: Scalar = 0.23109325002895065;
106pub const B_10: Scalar = 0.12747667699928525;
107pub const B_11: Scalar = 0.2246434176204158;
108pub const B_12: Scalar = 0.5684352689748513;
109pub const B_13: Scalar = 0.058258715572158275;
110pub const B_14: Scalar = 0.13643174034822156;
111pub const B_15: Scalar = 0.030570139830827976;
112
113pub const D_1: Scalar = -0.005357988290444578;
114pub const D_8: Scalar = -2.583020491182464;
115pub const D_9: Scalar = 0.14252253154686625;
116pub const D_10: Scalar = 0.013420653512688676;
117pub const D_11: Scalar = -0.02867296291409493;
118pub const D_12: Scalar = 2.624999655215792;
119pub const D_13: Scalar = -0.2825509643291537;
120pub const D_14: Scalar = 0.13643174034822156;
121pub const D_15: Scalar = 0.030570139830827976;
122pub const D_16: Scalar = -0.04834231373823958;
123
124#[doc = include_str!("doc.md")]
125#[derive(Debug)]
126pub struct Verner9 {
127    /// Absolute error tolerance.
128    pub abs_tol: Scalar,
129    /// Relative error tolerance.
130    pub rel_tol: Scalar,
131    /// Multiplier for adaptive time steps.
132    pub dt_beta: Scalar,
133    /// Exponent for adaptive time steps.
134    pub dt_expn: Scalar,
135    /// Cut back factor for the time step.
136    pub dt_cut: Scalar,
137    /// Minimum value for the time step.
138    pub dt_min: Scalar,
139}
140
141impl Default for Verner9 {
142    fn default() -> Self {
143        Self {
144            abs_tol: ABS_TOL,
145            rel_tol: REL_TOL,
146            dt_beta: 0.9,
147            dt_expn: 9.0,
148            dt_cut: 0.5,
149            dt_min: ABS_TOL,
150        }
151    }
152}
153
154impl<Y, U> OdeIntegrator<Y, U> for Verner9
155where
156    Y: Tensor,
157    U: TensorVec<Item = Y>,
158{
159}
160
161impl VariableStep for Verner9 {
162    fn abs_tol(&self) -> Scalar {
163        self.abs_tol
164    }
165    fn rel_tol(&self) -> Scalar {
166        self.rel_tol
167    }
168    fn dt_beta(&self) -> Scalar {
169        self.dt_beta
170    }
171    fn dt_expn(&self) -> Scalar {
172        self.dt_expn
173    }
174    fn dt_cut(&self) -> Scalar {
175        self.dt_cut
176    }
177    fn dt_min(&self) -> Scalar {
178        self.dt_min
179    }
180}
181
182impl<Y, U> Explicit<Y, U> for Verner9
183where
184    Y: Tensor,
185    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
186    U: TensorVec<Item = Y>,
187{
188    const SLOPES: usize = 16;
189    fn integrate(
190        &self,
191        function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
192        time: &[Scalar],
193        initial_condition: Y,
194    ) -> Result<(Vector, U, U), IntegrationError> {
195        self.integrate_variable_step(function, time, initial_condition)
196    }
197}
198
199impl<Y, U> VariableStepExplicit<Y, U> for Verner9
200where
201    Self: Explicit<Y, U>,
202    Y: Tensor,
203    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
204    U: TensorVec<Item = Y>,
205{
206    fn error(dt: Scalar, k: &[Y]) -> Result<Scalar, String> {
207        Ok(((&k[0] * D_1
208            + &k[7] * D_8
209            + &k[8] * D_9
210            + &k[9] * D_10
211            + &k[10] * D_11
212            + &k[11] * D_12
213            + &k[12] * D_13
214            + &k[13] * D_14
215            + &k[14] * D_15
216            + &k[15] * D_16)
217            * dt)
218            .norm_inf())
219    }
220    fn slopes(
221        mut function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
222        y: &Y,
223        t: Scalar,
224        dt: Scalar,
225        k: &mut [Y],
226        y_trial: &mut Y,
227    ) -> Result<(), String> {
228        k[0] = function(t, y)?;
229        *y_trial = &k[0] * (A_2_1 * dt) + y;
230        k[1] = function(t + C_2 * dt, y_trial)?;
231        *y_trial = &k[0] * (A_3_1 * dt) + &k[1] * (A_3_2 * dt) + y;
232        k[2] = function(t + C_3 * dt, y_trial)?;
233        *y_trial = &k[0] * (A_4_1 * dt) + &k[2] * (A_4_3 * dt) + y;
234        k[3] = function(t + C_4 * dt, y_trial)?;
235        *y_trial = &k[0] * (A_5_1 * dt) + &k[2] * (A_5_3 * dt) + &k[3] * (A_5_4 * dt) + y;
236        k[4] = function(t + C_5 * dt, y_trial)?;
237        *y_trial = &k[0] * (A_6_1 * dt) + &k[3] * (A_6_4 * dt) + &k[4] * (A_6_5 * dt) + y;
238        k[5] = function(t + C_6 * dt, y_trial)?;
239        *y_trial = &k[0] * (A_7_1 * dt)
240            + &k[3] * (A_7_4 * dt)
241            + &k[4] * (A_7_5 * dt)
242            + &k[5] * (A_7_6 * dt)
243            + y;
244        k[6] = function(t + C_7 * dt, y_trial)?;
245        *y_trial = &k[0] * (A_8_1 * dt) + &k[5] * (A_8_6 * dt) + &k[6] * (A_8_7 * dt) + y;
246        k[7] = function(t + C_8 * dt, y_trial)?;
247        *y_trial = &k[0] * (A_9_1 * dt)
248            + &k[5] * (A_9_6 * dt)
249            + &k[6] * (A_9_7 * dt)
250            + &k[7] * (A_9_8 * dt)
251            + y;
252        k[8] = function(t + C_9 * dt, y_trial)?;
253        *y_trial = &k[0] * (A_10_1 * dt)
254            + &k[5] * (A_10_6 * dt)
255            + &k[6] * (A_10_7 * dt)
256            + &k[7] * (A_10_8 * dt)
257            + &k[8] * (A_10_9 * dt)
258            + y;
259        k[9] = function(t + C_10 * dt, y_trial)?;
260        *y_trial = &k[0] * (A_11_1 * dt)
261            + &k[5] * (A_11_6 * dt)
262            + &k[6] * (A_11_7 * dt)
263            + &k[7] * (A_11_8 * dt)
264            + &k[8] * (A_11_9 * dt)
265            + &k[9] * (A_11_10 * dt)
266            + y;
267        k[10] = function(t + C_11 * dt, y_trial)?;
268        *y_trial = &k[0] * (A_12_1 * dt)
269            + &k[5] * (A_12_6 * dt)
270            + &k[6] * (A_12_7 * dt)
271            + &k[7] * (A_12_8 * dt)
272            + &k[8] * (A_12_9 * dt)
273            + &k[9] * (A_12_10 * dt)
274            + &k[10] * (A_12_11 * dt)
275            + y;
276        k[11] = function(t + C_12 * dt, y_trial)?;
277        *y_trial = &k[0] * (A_13_1 * dt)
278            + &k[5] * (A_13_6 * dt)
279            + &k[6] * (A_13_7 * dt)
280            + &k[7] * (A_13_8 * dt)
281            + &k[8] * (A_13_9 * dt)
282            + &k[9] * (A_13_10 * dt)
283            + &k[10] * (A_13_11 * dt)
284            + &k[11] * (A_13_12 * dt)
285            + y;
286        k[12] = function(t + C_13 * dt, y_trial)?;
287        *y_trial = &k[0] * (A_14_1 * dt)
288            + &k[5] * (A_14_6 * dt)
289            + &k[6] * (A_14_7 * dt)
290            + &k[7] * (A_14_8 * dt)
291            + &k[8] * (A_14_9 * dt)
292            + &k[9] * (A_14_10 * dt)
293            + &k[10] * (A_14_11 * dt)
294            + &k[11] * (A_14_12 * dt)
295            + &k[12] * (A_14_13 * dt)
296            + y;
297        k[13] = function(t + C_14 * dt, y_trial)?;
298        *y_trial = &k[0] * (A_15_1 * dt)
299            + &k[5] * (A_15_6 * dt)
300            + &k[6] * (A_15_7 * dt)
301            + &k[7] * (A_15_8 * dt)
302            + &k[8] * (A_15_9 * dt)
303            + &k[9] * (A_15_10 * dt)
304            + &k[10] * (A_15_11 * dt)
305            + &k[11] * (A_15_12 * dt)
306            + &k[12] * (A_15_13 * dt)
307            + &k[13] * (A_15_14 * dt)
308            + y;
309        k[14] = function(t + dt, y_trial)?;
310        *y_trial = &k[0] * (A_16_1 * dt)
311            + &k[5] * (A_16_6 * dt)
312            + &k[6] * (A_16_7 * dt)
313            + &k[7] * (A_16_8 * dt)
314            + &k[8] * (A_16_9 * dt)
315            + &k[9] * (A_16_10 * dt)
316            + &k[10] * (A_16_11 * dt)
317            + &k[11] * (A_16_12 * dt)
318            + &k[12] * (A_16_13 * dt)
319            + y;
320        if k.len() == Self::SLOPES {
321            k[15] = function(t + dt, y_trial)?;
322        }
323        *y_trial = (&k[0] * B_1
324            + &k[7] * B_8
325            + &k[8] * B_9
326            + &k[9] * B_10
327            + &k[10] * B_11
328            + &k[11] * B_12
329            + &k[12] * B_13
330            + &k[13] * B_14
331            + &k[14] * B_15)
332            * dt
333            + y;
334        Ok(())
335    }
336}
337
338impl<Y, U> InterpolateSolution<Y, U> for Verner9
339where
340    Y: Tensor,
341    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
342    U: TensorVec<Item = Y>,
343{
344    fn interpolate(
345        &self,
346        time: &Vector,
347        tp: &Vector,
348        yp: &U,
349        function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
350    ) -> Result<(U, U), IntegrationError> {
351        Self::interpolate_variable_step(time, tp, yp, function)
352    }
353}
conspire/math/integrate/ode/explicit/variable_step/verner_9/mod.rs

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