conspire/math/integrate/bogacki_shampine/
mod.rs#[cfg(test)]
mod test;
use super::{
super::{
interpolate::InterpolateSolution, Tensor, TensorArray, TensorRank0, TensorVec, Vector,
},
Explicit, IntegrationError,
};
use crate::{ABS_TOL, REL_TOL};
use std::ops::{Mul, Sub};
#[derive(Debug)]
pub struct BogackiShampine {
pub abs_tol: TensorRank0,
pub rel_tol: TensorRank0,
pub dt_beta: TensorRank0,
pub dt_expn: TensorRank0,
pub dt_init: TensorRank0,
}
impl Default for BogackiShampine {
fn default() -> Self {
Self {
abs_tol: ABS_TOL,
rel_tol: REL_TOL,
dt_beta: 0.9,
dt_expn: 3.0,
dt_init: 0.1,
}
}
}
impl<Y, U> Explicit<Y, U> for BogackiShampine
where
Self: InterpolateSolution<Y, U>,
Y: Tensor + TensorArray,
for<'a> &'a Y: Mul<TensorRank0, Output = Y> + Sub<&'a Y, Output = Y>,
U: TensorVec<Item = Y>,
{
fn integrate(
&self,
function: impl Fn(&TensorRank0, &Y) -> Y,
time: &[TensorRank0],
initial_condition: Y,
) -> Result<(Vector, U), IntegrationError> {
if time.len() < 2 {
return Err(IntegrationError::LengthTimeLessThanTwo);
} else if time[0] >= time[time.len() - 1] {
return Err(IntegrationError::InitialTimeNotLessThanFinalTime);
}
let mut t = time[0];
let mut dt = self.dt_init * time[time.len() - 1];
let mut e;
let mut k_1 = function(&t, &initial_condition);
let mut k_2;
let mut k_3;
let mut k_4;
let mut t_sol = Vector::zero(0);
t_sol.push(time[0]);
let mut y = initial_condition.clone();
let mut y_sol = U::zero(0);
y_sol.push(initial_condition.clone());
let mut y_trial;
while t < time[time.len() - 1] {
k_2 = function(&(t + 0.5 * dt), &(&k_1 * (0.5 * dt) + &y));
k_3 = function(&(t + 0.75 * dt), &(&k_2 * (0.75 * dt) + &y));
y_trial = (&k_1 * 2.0 + &k_2 * 3.0 + &k_3 * 4.0) * (dt / 9.0) + &y;
k_4 = function(&(t + dt), &y_trial);
e = ((&k_1 * -5.0 + k_2 * 6.0 + k_3 * 8.0 + &k_4 * -9.0) * (dt / 72.0)).norm();
if e < self.abs_tol || e / y_trial.norm() < self.rel_tol {
k_1 = k_4;
t += dt;
y = y_trial;
t_sol.push(t);
y_sol.push(y.clone());
}
dt *= self.dt_beta * (self.abs_tol / e).powf(1.0 / self.dt_expn);
}
if time.len() > 2 {
let t_int = Vector::new(time);
let y_int = self.interpolate(&t_int, &t_sol, &y_sol, function);
Ok((t_int, y_int))
} else {
Ok((t_sol, y_sol))
}
}
}
impl<Y, U> InterpolateSolution<Y, U> for BogackiShampine
where
Y: Tensor + TensorArray,
for<'a> &'a Y: Mul<TensorRank0, Output = Y> + Sub<&'a Y, Output = Y>,
U: TensorVec<Item = Y>,
{
fn interpolate(
&self,
time: &Vector,
tp: &Vector,
yp: &U,
function: impl Fn(&TensorRank0, &Y) -> Y,
) -> U {
let mut dt = 0.0;
let mut i = 0;
let mut k_1 = Y::zero();
let mut k_2 = Y::zero();
let mut k_3 = Y::zero();
let mut t = 0.0;
let mut y = Y::zero();
time.iter()
.map(|time_k| {
i = tp.iter().position(|tp_i| tp_i > time_k).unwrap();
t = tp[i - 1];
y = yp[i - 1].clone();
dt = time_k - t;
k_1 = function(&t, &y);
k_2 = function(&(t + 0.5 * dt), &(&k_1 * (0.5 * dt) + &y));
k_3 = function(&(t + 0.75 * dt), &(&k_2 * (0.75 * dt) + &y));
(&k_1 * 2.0 + &k_2 * 3.0 + &k_3 * 4.0) * (dt / 9.0) + &y
})
.collect()
}
}