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

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