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}
53
54impl Default for BogackiShampine {
55    fn default() -> Self {
56        Self {
57            abs_tol: ABS_TOL,
58            rel_tol: REL_TOL,
59            dt_beta: 0.9,
60            dt_expn: 3.0,
61        }
62    }
63}
64
65impl<Y, U> Explicit<Y, U> for BogackiShampine
66where
67    Self: InterpolateSolution<Y, U>,
68    Y: Tensor,
69    for<'a> &'a Y: Mul<TensorRank0, Output = Y> + Sub<&'a Y, Output = Y>,
70    U: TensorVec<Item = Y>,
71{
72    fn integrate(
73        &self,
74        mut function: impl FnMut(TensorRank0, &Y) -> Result<Y, IntegrationError>,
75        time: &[TensorRank0],
76        initial_condition: Y,
77    ) -> Result<(Vector, U, U), IntegrationError> {
78        let t_0 = time[0];
79        let t_f = time[time.len() - 1];
80        if time.len() < 2 {
81            return Err(IntegrationError::LengthTimeLessThanTwo);
82        } else if t_0 >= t_f {
83            return Err(IntegrationError::InitialTimeNotLessThanFinalTime);
84        }
85        let mut t = t_0;
86        let mut dt = t_f;
87        let mut e;
88        let mut k_1 = function(t, &initial_condition)?;
89        let mut k_2;
90        let mut k_3;
91        let mut k_4;
92        let mut t_sol = Vector::zero(0);
93        t_sol.push(t_0);
94        let mut y = initial_condition.clone();
95        let mut y_sol = U::zero(0);
96        y_sol.push(initial_condition.clone());
97        let mut dydt_sol = U::zero(0);
98        dydt_sol.push(k_1.clone());
99        let mut y_trial;
100        while t < t_f {
101            k_2 = function(t + 0.5 * dt, &(&k_1 * (0.5 * dt) + &y))?;
102            k_3 = function(t + 0.75 * dt, &(&k_2 * (0.75 * dt) + &y))?;
103            y_trial = (&k_1 * 2.0 + &k_2 * 3.0 + &k_3 * 4.0) * (dt / 9.0) + &y;
104            k_4 = function(t + dt, &y_trial)?;
105            e = ((&k_1 * -5.0 + k_2 * 6.0 + k_3 * 8.0 + &k_4 * -9.0) * (dt / 72.0)).norm_inf();
106            if e < self.abs_tol || e / y_trial.norm_inf() < self.rel_tol {
107                k_1 = k_4;
108                t += dt;
109                y = y_trial;
110                t_sol.push(t);
111                y_sol.push(y.clone());
112                dydt_sol.push(k_1.clone());
113            }
114            if e > 0.0 {
115                dt *= self.dt_beta * (self.abs_tol / e).powf(1.0 / self.dt_expn)
116            }
117            dt = dt.min(t_f - t)
118        }
119        if time.len() > 2 {
120            let t_int = Vector::new(time);
121            let (y_int, dydt_int) = self.interpolate(&t_int, &t_sol, &y_sol, function)?;
122            Ok((t_int, y_int, dydt_int))
123        } else {
124            Ok((t_sol, y_sol, dydt_sol))
125        }
126    }
127}
128
129impl<Y, U> InterpolateSolution<Y, U> for BogackiShampine
130where
131    Y: Tensor,
132    for<'a> &'a Y: Mul<TensorRank0, Output = Y> + Sub<&'a Y, Output = Y>,
133    U: TensorVec<Item = Y>,
134{
135    fn interpolate(
136        &self,
137        time: &Vector,
138        tp: &Vector,
139        yp: &U,
140        mut function: impl FnMut(TensorRank0, &Y) -> Result<Y, IntegrationError>,
141    ) -> Result<(U, U), IntegrationError> {
142        let mut dt;
143        let mut i;
144        let mut k_1;
145        let mut k_2;
146        let mut k_3;
147        let mut t;
148        let mut y;
149        let mut y_int = U::zero(0);
150        let mut dydt_int = U::zero(0);
151        let mut y_trial;
152        for time_k in time.iter() {
153            i = tp.iter().position(|tp_i| tp_i >= time_k).unwrap();
154            if time_k == &tp[i] {
155                t = tp[i];
156                y_trial = yp[i].clone();
157                dt = 0.0;
158            } else {
159                t = tp[i - 1];
160                y = yp[i - 1].clone();
161                dt = time_k - t;
162                k_1 = function(t, &y)?;
163                k_2 = function(t + 0.5 * dt, &(&k_1 * (0.5 * dt) + &y))?;
164                k_3 = function(t + 0.75 * dt, &(&k_2 * (0.75 * dt) + &y))?;
165                y_trial = (&k_1 * 2.0 + &k_2 * 3.0 + &k_3 * 4.0) * (dt / 9.0) + &y;
166            }
167            dydt_int.push(function(t + dt, &y_trial)?);
168            y_int.push(y_trial);
169        }
170        Ok((y_int, dydt_int))
171    }
172}
conspire/math/integrate/bogacki_shampine/mod.rs

conspire/math/integrate/bogacki_shampine/
mod.rs
