conspire/math/integrate/bogacki_shampine/
mod.rs

1#[cfg(test)]
2mod test;
3
4use super::{
5    super::{Tensor, TensorRank0, TensorVec, Vector, interpolate::InterpolateSolution},
6    Explicit, IntegrationError,
7};
8use crate::{ABS_TOL, REL_TOL};
9use std::ops::{Mul, Sub};
10
11/// Explicit, three-stage, third-order, variable-step, Runge-Kutta method.[^cite]
12///
13/// [^cite]: P. Bogacki and L.F. Shampine, [Appl. Math. Lett. **2**, 321 (1989)](https://doi.org/10.1016/0893-9659(89)90079-7).
14///
15/// ```math
16/// \frac{dy}{dt} = f(t, y)
17/// ```
18/// ```math
19/// t_{n+1} = t_n + h
20/// ```
21/// ```math
22/// k_1 = f(t_n, y_n)
23/// ```
24/// ```math
25/// k_2 = f(t_n + \tfrac{1}{2} h, y_n + \tfrac{1}{2} h k_1)
26/// ```
27/// ```math
28/// k_3 = f(t_n + \tfrac{3}{4} h, y_n + \tfrac{3}{4} h k_2)
29/// ```
30/// ```math
31/// y_{n+1} = y_n + \frac{h}{9}\left(2k_1 + 3k_2 + 4k_3\right)
32/// ```
33/// ```math
34/// k_4 = f(t_{n+1}, y_{n+1})
35/// ```
36/// ```math
37/// e_{n+1} = \frac{h}{72}\left(-5k_1 + 6k_2 + 8k_3 - 9k_4\right)
38/// ```
39/// ```math
40/// h_{n+1} = \beta h \left(\frac{e_\mathrm{tol}}{e_{n+1}}\right)^{1/p}
41/// ```
42#[derive(Debug)]
43pub struct BogackiShampine {
44    /// Absolute error tolerance.
45    pub abs_tol: TensorRank0,
46    /// Relative error tolerance.
47    pub rel_tol: TensorRank0,
48    /// Multiplier for adaptive time steps.
49    pub dt_beta: TensorRank0,
50    /// Exponent for adaptive time steps.
51    pub dt_expn: TensorRank0,
52    /// Initial relative time step.
53    pub dt_init: TensorRank0,
54}
55
56impl Default for BogackiShampine {
57    fn default() -> Self {
58        Self {
59            abs_tol: ABS_TOL,
60            rel_tol: REL_TOL,
61            dt_beta: 0.9,
62            dt_expn: 3.0,
63            dt_init: 0.1,
64        }
65    }
66}
67
68impl<Y, U> Explicit<Y, U> for BogackiShampine
69where
70    Self: InterpolateSolution<Y, U>,
71    Y: Tensor,
72    for<'a> &'a Y: Mul<TensorRank0, Output = Y> + Sub<&'a Y, Output = Y>,
73    U: TensorVec<Item = Y>,
74{
75    fn integrate(
76        &self,
77        function: impl Fn(TensorRank0, &Y) -> Result<Y, IntegrationError>,
78        time: &[TensorRank0],
79        initial_condition: Y,
80    ) -> Result<(Vector, U, U), IntegrationError> {
81        if time.len() < 2 {
82            return Err(IntegrationError::LengthTimeLessThanTwo);
83        } else if time[0] >= time[time.len() - 1] {
84            return Err(IntegrationError::InitialTimeNotLessThanFinalTime);
85        }
86        let mut t = time[0];
87        let mut dt = self.dt_init * time[time.len() - 1];
88        let mut e;
89        let mut k_1 = function(t, &initial_condition)?;
90        let mut k_2;
91        let mut k_3;
92        let mut k_4;
93        let mut t_sol = Vector::zero(0);
94        t_sol.push(time[0]);
95        let mut y = initial_condition.clone();
96        let mut y_sol = U::zero(0);
97        y_sol.push(initial_condition.clone());
98        let mut dydt_sol = U::zero(0);
99        dydt_sol.push(k_1.clone());
100        let mut y_trial;
101        while t < time[time.len() - 1] {
102            k_2 = function(t + 0.5 * dt, &(&k_1 * (0.5 * dt) + &y))?;
103            k_3 = function(t + 0.75 * dt, &(&k_2 * (0.75 * dt) + &y))?;
104            y_trial = (&k_1 * 2.0 + &k_2 * 3.0 + &k_3 * 4.0) * (dt / 9.0) + &y;
105            k_4 = function(t + dt, &y_trial)?;
106            e = ((&k_1 * -5.0 + k_2 * 6.0 + k_3 * 8.0 + &k_4 * -9.0) * (dt / 72.0)).norm_inf();
107            if e < self.abs_tol || e / y_trial.norm_inf() < self.rel_tol {
108                k_1 = k_4;
109                t += dt;
110                y = y_trial;
111                t_sol.push(t);
112                y_sol.push(y.clone());
113                dydt_sol.push(k_1.clone());
114            }
115            dt *= self.dt_beta * (self.abs_tol / e).powf(1.0 / self.dt_expn);
116        }
117        if time.len() > 2 {
118            let t_int = Vector::new(time);
119            let (y_int, dydt_int) = self.interpolate(&t_int, &t_sol, &y_sol, function)?;
120            Ok((t_int, y_int, dydt_int))
121        } else {
122            Ok((t_sol, y_sol, dydt_sol))
123        }
124    }
125}
126
127impl<Y, U> InterpolateSolution<Y, U> for BogackiShampine
128where
129    Y: Tensor,
130    for<'a> &'a Y: Mul<TensorRank0, Output = Y> + Sub<&'a Y, Output = Y>,
131    U: TensorVec<Item = Y>,
132{
133    fn interpolate(
134        &self,
135        time: &Vector,
136        tp: &Vector,
137        yp: &U,
138        function: impl Fn(TensorRank0, &Y) -> Result<Y, IntegrationError>,
139    ) -> Result<(U, U), IntegrationError> {
140        let mut dt;
141        let mut i;
142        let mut k_1;
143        let mut k_2;
144        let mut k_3;
145        let mut t;
146        let mut y;
147        let mut y_int = U::zero(0);
148        let mut dydt_int = U::zero(0);
149        let mut y_trial;
150        for time_k in time.iter() {
151            i = tp.iter().position(|tp_i| tp_i > time_k).unwrap();
152            t = tp[i - 1];
153            y = yp[i - 1].clone();
154            dt = time_k - t;
155            k_1 = function(t, &y)?;
156            k_2 = function(t + 0.5 * dt, &(&k_1 * (0.5 * dt) + &y))?;
157            k_3 = function(t + 0.75 * dt, &(&k_2 * (0.75 * dt) + &y))?;
158            y_trial = (&k_1 * 2.0 + &k_2 * 3.0 + &k_3 * 4.0) * (dt / 9.0) + &y;
159            dydt_int.push(function(t + dt, &y_trial)?);
160            y_int.push(y_trial);
161        }
162        Ok((y_int, dydt_int))
163    }
164}
conspire/math/integrate/bogacki_shampine/mod.rs

conspire/math/integrate/bogacki_shampine/
mod.rs
