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

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