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

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

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