conspire/math/integrate/ode/explicit/variable_step/verner_8/
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.05;
17pub(crate) const C_3: Scalar = 0.1065625;
18pub(crate) const C_4: Scalar = 0.15984375;
19pub(crate) const C_5: Scalar = 0.39;
20pub(crate) const C_6: Scalar = 0.465;
21pub(crate) const C_7: Scalar = 0.155;
22pub(crate) const C_8: Scalar = 0.943;
23pub(crate) const C_9: Scalar = 0.901802041735857;
24pub(crate) const C_10: Scalar = 0.909;
25pub(crate) const C_11: Scalar = 0.94;
26
27pub(crate) const A_2_1: Scalar = 0.05;
28pub(crate) const A_3_1: Scalar = -0.0069931640625;
29pub(crate) const A_3_2: Scalar = 0.1135556640625;
30pub(crate) const A_4_1: Scalar = 0.0399609375;
31pub(crate) const A_4_3: Scalar = 0.1198828125;
32pub(crate) const A_5_1: Scalar = 0.36139756280045754;
33pub(crate) const A_5_3: Scalar = -1.3415240667004928;
34pub(crate) const A_5_4: Scalar = 1.3701265039000352;
35pub(crate) const A_6_1: Scalar = 0.049047202797202795;
36pub(crate) const A_6_4: Scalar = 0.23509720422144048;
37pub(crate) const A_6_5: Scalar = 0.18085559298135673;
38pub(crate) const A_7_1: Scalar = 0.06169289044289044;
39pub(crate) const A_7_4: Scalar = 0.11236568314640277;
40pub(crate) const A_7_5: Scalar = -0.03885046071451367;
41pub(crate) const A_7_6: Scalar = 0.01979188712522046;
42pub(crate) const A_8_1: Scalar = -1.767630240222327;
43pub(crate) const A_8_4: Scalar = -62.5;
44pub(crate) const A_8_5: Scalar = -6.061889377376669;
45pub(crate) const A_8_6: Scalar = 5.6508231982227635;
46pub(crate) const A_8_7: Scalar = 65.62169641937624;
47pub(crate) const A_9_1: Scalar = -1.1809450665549708;
48pub(crate) const A_9_4: Scalar = -41.50473441114321;
49pub(crate) const A_9_5: Scalar = -4.434438319103725;
50pub(crate) const A_9_6: Scalar = 4.260408188586133;
51pub(crate) const A_9_7: Scalar = 43.75364022446172;
52pub(crate) const A_9_8: Scalar = 0.00787142548991231;
53pub(crate) const A_10_1: Scalar = -1.2814059994414884;
54pub(crate) const A_10_4: Scalar = -45.047139960139866;
55pub(crate) const A_10_5: Scalar = -4.731362069449576;
56pub(crate) const A_10_6: Scalar = 4.514967016593808;
57pub(crate) const A_10_7: Scalar = 47.44909557172985;
58pub(crate) const A_10_8: Scalar = 0.01059228297111661;
59pub(crate) const A_10_9: Scalar = -0.0057468422638446166;
60pub(crate) const A_11_1: Scalar = -1.7244701342624853;
61pub(crate) const A_11_4: Scalar = -60.92349008483054;
62pub(crate) const A_11_5: Scalar = -5.951518376222392;
63pub(crate) const A_11_6: Scalar = 5.556523730698456;
64pub(crate) const A_11_7: Scalar = 63.98301198033305;
65pub(crate) const A_11_8: Scalar = 0.014642028250414961;
66pub(crate) const A_11_9: Scalar = 0.06460408772358203;
67pub(crate) const A_11_10: Scalar = -0.0793032316900888;
68pub(crate) const A_12_1: Scalar = -3.301622667747079;
69pub(crate) const A_12_4: Scalar = -118.01127235975251;
70pub(crate) const A_12_5: Scalar = -10.141422388456112;
71pub(crate) const A_12_6: Scalar = 9.139311332232058;
72pub(crate) const A_12_7: Scalar = 123.37594282840426;
73pub(crate) const A_12_8: Scalar = 4.62324437887458;
74pub(crate) const A_12_9: Scalar = -3.3832777380682018;
75pub(crate) const A_12_10: Scalar = 4.527592100324618;
76pub(crate) const A_12_11: Scalar = -5.828495485811623;
77pub(crate) const A_13_1: Scalar = -3.039515033766309;
78pub(crate) const A_13_4: Scalar = -109.26086808941763;
79pub(crate) const A_13_5: Scalar = -9.290642497400293;
80pub(crate) const A_13_6: Scalar = 8.43050498176491;
81pub(crate) const A_13_7: Scalar = 114.20100103783314;
82pub(crate) const A_13_8: Scalar = -0.9637271342145479;
83pub(crate) const A_13_9: Scalar = -5.0348840888021895;
84pub(crate) const A_13_10: Scalar = 5.958130824002923;
85
86pub(crate) const B_1: Scalar = 0.04427989419007951;
87pub(crate) const B_6: Scalar = 0.3541049391724449;
88pub(crate) const B_7: Scalar = 0.24796921549564377;
89pub(crate) const B_8: Scalar = -15.694202038838085;
90pub(crate) const B_9: Scalar = 25.084064965558564;
91pub(crate) const B_10: Scalar = -31.738367786260277;
92pub(crate) const B_11: Scalar = 22.938283273988784;
93pub(crate) const B_12: Scalar = -0.2361324633071542;
94
95pub(crate) const D_1: Scalar = -0.00003272103901028138;
96pub(crate) const D_6: Scalar = -0.0005046250618777704;
97pub(crate) const D_7: Scalar = 0.0001211723589784759;
98pub(crate) const D_8: Scalar = -20.142336771313868;
99pub(crate) const D_9: Scalar = 5.2371785994398286;
100pub(crate) const D_10: Scalar = -8.156744408794658;
101pub(crate) const D_11: Scalar = 22.938283273988784;
102pub(crate) const D_12: Scalar = -0.2361324633071542;
103pub(crate) const D_13: Scalar = 0.36016794372897754;
104
105/// The Verner 8(7) tableau.
106#[derive(Debug)]
107pub struct Tableau;
108
109impl ButcherTableau for Tableau {
110    const STAGES: usize = 13;
111    const ORDER: Scalar = 8.0;
112    #[rustfmt::skip]
113    const A: &'static [&'static [Scalar]] = &[
114        &[],
115        &[A_2_1],
116        &[A_3_1, A_3_2],
117        &[A_4_1, 0.0, A_4_3],
118        &[A_5_1, 0.0, A_5_3, A_5_4],
119        &[A_6_1, 0.0, 0.0, A_6_4, A_6_5],
120        &[A_7_1, 0.0, 0.0, A_7_4, A_7_5, A_7_6],
121        &[A_8_1, 0.0, 0.0, A_8_4, A_8_5, A_8_6, A_8_7],
122        &[A_9_1, 0.0, 0.0, A_9_4, A_9_5, A_9_6, A_9_7, A_9_8],
123        &[A_10_1, 0.0, 0.0, A_10_4, A_10_5, A_10_6, A_10_7, A_10_8, A_10_9],
124        &[A_11_1, 0.0, 0.0, A_11_4, A_11_5, A_11_6, A_11_7, A_11_8, A_11_9, A_11_10],
125        &[A_12_1, 0.0, 0.0, A_12_4, A_12_5, A_12_6, A_12_7, A_12_8, A_12_9, A_12_10, A_12_11],
126        &[A_13_1, 0.0, 0.0, A_13_4, A_13_5, A_13_6, A_13_7, A_13_8, A_13_9, A_13_10, 0.0, 0.0],
127    ];
128    const C: &'static [Scalar] = &[
129        0.0, C_2, C_3, C_4, C_5, C_6, C_7, C_8, C_9, C_10, C_11, 1.0, 1.0,
130    ];
131    #[rustfmt::skip]
132    const B: &'static [Scalar] =
133        &[B_1, 0.0, 0.0, 0.0, 0.0, B_6, B_7, B_8, B_9, B_10, B_11, B_12, 0.0];
134}
135
136impl EmbeddedTableau for Tableau {
137    #[rustfmt::skip]
138    const D: &'static [Scalar] =
139        &[D_1, 0.0, 0.0, 0.0, 0.0, D_6, D_7, D_8, D_9, D_10, D_11, D_12, D_13];
140}
141
142#[doc = include_str!("doc.md")]
143#[derive(Debug)]
144pub struct Verner8 {
145    /// Absolute error tolerance.
146    pub abs_tol: Scalar,
147    /// Relative error tolerance.
148    pub rel_tol: Scalar,
149    /// Multiplier for adaptive time steps.
150    pub dt_beta: Scalar,
151    /// Exponent for adaptive time steps.
152    pub dt_expn: Scalar,
153    /// Cut back factor for the time step.
154    pub dt_cut: Scalar,
155    /// Growth factor ceiling for the time step.
156    pub dt_grow: Scalar,
157    /// Minimum value for the time step.
158    pub dt_min: Scalar,
159    /// Norm type for error evaluation.
160    pub error_norm: Norm,
161}
162
163impl Default for Verner8 {
164    fn default() -> Self {
165        Self {
166            abs_tol: ABS_TOL,
167            rel_tol: REL_TOL,
168            dt_beta: 0.9,
169            dt_expn: 8.0,
170            dt_cut: 0.5,
171            dt_grow: 5.0,
172            dt_min: ABS_TOL,
173            error_norm: Norm::Chebyshev,
174        }
175    }
176}
177
178impl<Y, U> OdeIntegrator<Y, U> for Verner8
179where
180    Y: Tensor,
181    U: TensorVec<Item = Y>,
182{
183}
184
185impl<T> VariableStep<T> for Verner8 {
186    fn abs_tol(&self) -> Scalar {
187        self.abs_tol
188    }
189    fn rel_tol(&self) -> Scalar {
190        self.rel_tol
191    }
192    fn dt_beta(&self) -> Scalar {
193        self.dt_beta
194    }
195    fn dt_expn(&self) -> Scalar {
196        self.dt_expn
197    }
198    fn dt_cut(&self) -> Scalar {
199        self.dt_cut
200    }
201    fn dt_grow(&self) -> Scalar {
202        self.dt_grow
203    }
204    fn dt_min(&self) -> Quantity<T> {
205        Quantity::new(self.dt_min)
206    }
207    fn error_norm(&self) -> &Norm {
208        &self.error_norm
209    }
210}
211
212impl<Y, U, V, T> Explicit<Y, U, V, T> for Verner8
213where
214    Y: Differentiable<T> + Tensor,
215    Derivative<Y, T>: Mul<Quantity<T>, Output = Y>,
216    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
217    for<'a> &'a Derivative<Y, T>:
218        Mul<Scalar, Output = Derivative<Y, T>> + Mul<Quantity<T>, Output = Y>,
219    U: TensorVec<Item = Y>,
220    V: TensorVec<Item = Derivative<Y, T>>,
221{
222    const SLOPES: usize = 13;
223    fn integrate(
224        &self,
225        function: impl FnMut(Quantity<T>, &Y) -> Result<Derivative<Y, T>, String>,
226        time: &[Quantity<T>],
227        initial_condition: Y,
228    ) -> Result<(Times<T>, U, V), IntegrationError> {
229        self.integrate_variable_step(function, time, initial_condition)
230    }
231}
232
233impl<Y, U, V, T> VariableStepExplicit<Y, U, V, T> for Verner8
234where
235    Self: Explicit<Y, U, V, T>,
236    Y: Differentiable<T> + Tensor,
237    Derivative<Y, T>: Mul<Quantity<T>, Output = Y>,
238    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
239    for<'a> &'a Derivative<Y, T>:
240        Mul<Scalar, Output = Derivative<Y, T>> + Mul<Quantity<T>, Output = Y>,
241    U: TensorVec<Item = Y>,
242    V: TensorVec<Item = Derivative<Y, T>>,
243{
244    type Tableau = Tableau;
245}
246
247impl<Y, U, V, T> InterpolateSolution<Y, U, V, T> for Verner8
248where
249    Y: Differentiable<T> + Tensor,
250    Derivative<Y, T>: Mul<Quantity<T>, Output = Y>,
251    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
252    for<'a> &'a Derivative<Y, T>:
253        Mul<Scalar, Output = Derivative<Y, T>> + Mul<Quantity<T>, Output = Y>,
254    U: TensorVec<Item = Y>,
255    V: TensorVec<Item = Derivative<Y, T>>,
256{
257    fn interpolate(
258        &self,
259        time: &Times<T>,
260        tp: &Times<T>,
261        yp: &U,
262        _dydtp: &V,
263        _k_sol: &[V],
264        function: impl FnMut(Quantity<T>, &Y) -> Result<Derivative<Y, T>, String>,
265    ) -> Result<(U, V), IntegrationError> {
266        Self::interpolate_variable_step(time, tp, yp, function)
267    }
268}
conspire/math/integrate/ode/explicit/variable_step/verner_8/mod.rs

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