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

1#[cfg(test)]
2mod test;
3
4use crate::math::{
5    Scalar, Tensor, TensorVec, Vector,
6    integrate::{
7        Explicit, IntegrationError, OdeIntegrator, VariableStep, VariableStepExplicit,
8        VariableStepExplicitFirstSameAsLast,
9    },
10    interpolate::InterpolateSolution,
11};
12use crate::{ABS_TOL, REL_TOL};
13use std::ops::{Mul, Sub};
14
15#[doc = include_str!("doc.md")]
16#[derive(Debug)]
17pub struct BogackiShampine {
18    /// Absolute error tolerance.
19    pub abs_tol: Scalar,
20    /// Relative error tolerance.
21    pub rel_tol: Scalar,
22    /// Multiplier for adaptive time steps.
23    pub dt_beta: Scalar,
24    /// Exponent for adaptive time steps.
25    pub dt_expn: Scalar,
26    /// Cut back factor for the time step.
27    pub dt_cut: Scalar,
28    /// Minimum value for the time step.
29    pub dt_min: Scalar,
30}
31
32impl Default for BogackiShampine {
33    fn default() -> Self {
34        Self {
35            abs_tol: ABS_TOL,
36            rel_tol: REL_TOL,
37            dt_beta: 0.9,
38            dt_expn: 3.0,
39            dt_cut: 0.5,
40            dt_min: ABS_TOL,
41        }
42    }
43}
44
45impl<Y, U> OdeIntegrator<Y, U> for BogackiShampine
46where
47    Y: Tensor,
48    U: TensorVec<Item = Y>,
49{
50}
51
52impl VariableStep for BogackiShampine {
53    fn abs_tol(&self) -> Scalar {
54        self.abs_tol
55    }
56    fn rel_tol(&self) -> Scalar {
57        self.rel_tol
58    }
59    fn dt_beta(&self) -> Scalar {
60        self.dt_beta
61    }
62    fn dt_expn(&self) -> Scalar {
63        self.dt_expn
64    }
65    fn dt_cut(&self) -> Scalar {
66        self.dt_cut
67    }
68    fn dt_min(&self) -> Scalar {
69        self.dt_min
70    }
71}
72
73impl<Y, U> Explicit<Y, U> for BogackiShampine
74where
75    Y: Tensor,
76    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
77    U: TensorVec<Item = Y>,
78{
79    const SLOPES: usize = 4;
80    fn integrate(
81        &self,
82        function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
83        time: &[Scalar],
84        initial_condition: Y,
85    ) -> Result<(Vector, U, U), IntegrationError> {
86        self.integrate_variable_step(function, time, initial_condition)
87    }
88}
89
90impl<Y, U> VariableStepExplicit<Y, U> for BogackiShampine
91where
92    Self: Explicit<Y, U>,
93    Y: Tensor,
94    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
95    U: TensorVec<Item = Y>,
96{
97    fn error(dt: Scalar, k: &[Y]) -> Result<Scalar, String> {
98        Ok(((&k[0] * -5.0 + &k[1] * 6.0 + &k[2] * 8.0 + &k[3] * -9.0) * (dt / 72.0)).norm_inf())
99    }
100    fn slopes(
101        mut function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
102        y: &Y,
103        t: Scalar,
104        dt: Scalar,
105        k: &mut [Y],
106        y_trial: &mut Y,
107    ) -> Result<(), String> {
108        *y_trial = &k[0] * (0.5 * dt) + y;
109        k[1] = function(t + 0.5 * dt, y_trial)?;
110        *y_trial = &k[1] * (0.75 * dt) + y;
111        k[2] = function(t + 0.75 * dt, y_trial)?;
112        *y_trial = (&k[0] * 2.0 + &k[1] * 3.0 + &k[2] * 4.0) * (dt / 9.0) + y;
113        Ok(())
114    }
115    fn slopes_and_error(
116        &self,
117        function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
118        y: &Y,
119        t: Scalar,
120        dt: Scalar,
121        k: &mut [Y],
122        y_trial: &mut Y,
123    ) -> Result<Scalar, String> {
124        Self::slopes_and_error_fsal(function, y, t, dt, k, y_trial)
125    }
126    fn step(
127        &self,
128        _function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
129        y: &mut Y,
130        t: &mut Scalar,
131        y_sol: &mut U,
132        t_sol: &mut Vector,
133        dydt_sol: &mut U,
134        dt: &mut Scalar,
135        k: &mut [Y],
136        y_trial: &Y,
137        e: Scalar,
138    ) -> Result<(), String> {
139        let dt_0 = *dt;
140        self.step_fsal(y, t, y_sol, t_sol, dydt_sol, dt, k, y_trial, e)?;
141        if e > 0.0 {
142            *dt = dt_0;
143            *dt *= self.dt_beta() * (self.abs_tol() / e).powf(1.0 / self.dt_expn())
144        }
145        Ok(()) // some temporary fixes to pass tests in fem that are barely failing
146    }
147}
148
149impl<Y, U> VariableStepExplicitFirstSameAsLast<Y, U> for BogackiShampine
150where
151    Y: Tensor,
152    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
153    U: TensorVec<Item = Y>,
154{
155}
156
157impl<Y, U> InterpolateSolution<Y, U> for BogackiShampine
158where
159    Y: Tensor,
160    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
161    U: TensorVec<Item = Y>,
162{
163    fn interpolate(
164        &self,
165        time: &Vector,
166        tp: &Vector,
167        yp: &U,
168        function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
169    ) -> Result<(U, U), IntegrationError> {
170        Self::interpolate_variable_step(time, tp, yp, function)
171    }
172}
conspire/math/integrate/ode/explicit/variable_step/bogacki_shampine/mod.rs

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