#[cfg(test)]
mod test;
use crate::math::Norm;
use crate::math::{
Derivative, Differentiate, Quantity, Scalar, Tensor, TensorVec,
integrate::{
Explicit, FreeInterpolant, IntegrationError, OdeIntegrator, Times, VariableStep,
VariableStepExplicit, VariableStepExplicitFirstSameAsLast,
},
interpolate::InterpolateSolution,
};
use crate::{ABS_TOL, REL_TOL};
use std::ops::{Div, Mul, Sub};
#[doc = include_str!("doc.md")]
#[derive(Debug)]
pub struct BogackiShampine {
pub abs_tol: Scalar,
pub rel_tol: Scalar,
pub dt_beta: Scalar,
pub dt_expn: Scalar,
pub dt_cut: Scalar,
pub dt_min: Scalar,
pub error_norm: Norm,
}
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_cut: 0.5,
dt_min: ABS_TOL,
error_norm: Norm::Chebyshev,
}
}
}
impl<Y, U> OdeIntegrator<Y, U> for BogackiShampine
where
Y: Tensor,
U: TensorVec<Item = Y>,
{
}
impl<T> VariableStep<T> for BogackiShampine {
fn abs_tol(&self) -> Scalar {
self.abs_tol
}
fn rel_tol(&self) -> Scalar {
self.rel_tol
}
fn dt_beta(&self) -> Scalar {
self.dt_beta
}
fn dt_expn(&self) -> Scalar {
self.dt_expn
}
fn dt_cut(&self) -> Scalar {
self.dt_cut
}
fn dt_min(&self) -> Quantity<T> {
Quantity::new(self.dt_min)
}
fn error_norm(&self) -> &Norm {
&self.error_norm
}
}
impl<Y, U, V, T> Explicit<Y, U, V, T> for BogackiShampine
where
Y: Differentiate<T> + Div<Quantity<T>, Output = Derivative<Y, T>> + Tensor,
Derivative<Y, T>: Mul<Quantity<T>, Output = Y>,
for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
for<'a> &'a Derivative<Y, T>:
Mul<Scalar, Output = Derivative<Y, T>> + Mul<Quantity<T>, Output = Y>,
U: TensorVec<Item = Y>,
V: TensorVec<Item = Derivative<Y, T>>,
{
const SLOPES: usize = 4;
fn integrate(
&self,
function: impl FnMut(Quantity<T>, &Y) -> Result<Derivative<Y, T>, String>,
time: &[Quantity<T>],
initial_condition: Y,
) -> Result<(Times<T>, U, V), IntegrationError> {
self.integrate_variable_step(function, time, initial_condition)
}
}
impl<Y, U, V, T> VariableStepExplicit<Y, U, V, T> for BogackiShampine
where
Self: Explicit<Y, U, V, T>,
Y: Differentiate<T> + Div<Quantity<T>, Output = Derivative<Y, T>> + Tensor,
Derivative<Y, T>: Mul<Quantity<T>, Output = Y>,
for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
for<'a> &'a Derivative<Y, T>:
Mul<Scalar, Output = Derivative<Y, T>> + Mul<Quantity<T>, Output = Y>,
U: TensorVec<Item = Y>,
V: TensorVec<Item = Derivative<Y, T>>,
{
fn error(&self, dt: Quantity<T>, k: &[Derivative<Y, T>]) -> Result<Scalar, String> {
Ok(self
.error_norm
.measure(&((&k[0] * -5.0 + &k[1] * 6.0 + &k[2] * 8.0 + &k[3] * -9.0) * (dt / 72.0))))
}
fn slopes(
mut function: impl FnMut(Quantity<T>, &Y) -> Result<Derivative<Y, T>, String>,
y: &Y,
t: Quantity<T>,
dt: Quantity<T>,
k: &mut [Derivative<Y, T>],
y_trial: &mut Y,
) -> Result<(), String> {
*y_trial = &k[0] * (0.5 * dt) + y;
k[1] = function(t + 0.5 * dt, y_trial)?;
*y_trial = &k[1] * (0.75 * dt) + y;
k[2] = function(t + 0.75 * dt, y_trial)?;
*y_trial = (&k[0] * 2.0 + &k[1] * 3.0 + &k[2] * 4.0) * (dt / 9.0) + y;
Ok(())
}
fn slopes_and_error(
&self,
function: impl FnMut(Quantity<T>, &Y) -> Result<Derivative<Y, T>, String>,
y: &Y,
t: Quantity<T>,
dt: Quantity<T>,
k: &mut [Derivative<Y, T>],
y_trial: &mut Y,
) -> Result<Scalar, String> {
self.slopes_and_error_fsal(function, y, t, dt, k, y_trial)
}
fn step(
&self,
_function: impl FnMut(Quantity<T>, &Y) -> Result<Derivative<Y, T>, String>,
y: &mut Y,
t: &mut Quantity<T>,
y_sol: &mut U,
t_sol: &mut Times<T>,
dydt_sol: &mut V,
k_sol: &mut Vec<V>,
dt: &mut Quantity<T>,
k: &mut [Derivative<Y, T>],
y_trial: &Y,
e: Scalar,
) -> Result<(), String> {
let dt_0 = *dt;
self.step_fsal(y, t, y_sol, t_sol, dydt_sol, k_sol, dt, k, y_trial, e)?;
if e > 0.0 {
let (beta, tol, expn) = (
VariableStep::<T>::dt_beta(self),
VariableStep::<T>::abs_tol(self),
VariableStep::<T>::dt_expn(self),
);
*dt = dt_0;
*dt *= beta * (tol / e).powf(1.0 / expn)
}
Ok(()) }
}
impl<Y, U, V, T> VariableStepExplicitFirstSameAsLast<Y, U, V, T> for BogackiShampine
where
Y: Differentiate<T> + Div<Quantity<T>, Output = Derivative<Y, T>> + Tensor,
Derivative<Y, T>: Mul<Quantity<T>, Output = Y>,
for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
for<'a> &'a Derivative<Y, T>:
Mul<Scalar, Output = Derivative<Y, T>> + Mul<Quantity<T>, Output = Y>,
U: TensorVec<Item = Y>,
V: TensorVec<Item = Derivative<Y, T>>,
{
}
impl<Y, U, V, T> FreeInterpolant<Y, U, V, T> for BogackiShampine
where
Y: Differentiate<T> + Div<Quantity<T>, Output = Derivative<Y, T>> + Tensor,
Derivative<Y, T>: Mul<Quantity<T>, Output = Y>,
for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
for<'a> &'a Derivative<Y, T>:
Mul<Scalar, Output = Derivative<Y, T>> + Mul<Quantity<T>, Output = Y>,
U: TensorVec<Item = Y>,
V: TensorVec<Item = Derivative<Y, T>>,
{
}
impl<Y, U, V, T> InterpolateSolution<Y, U, V, T> for BogackiShampine
where
Y: Differentiate<T> + Div<Quantity<T>, Output = Derivative<Y, T>> + Tensor,
Derivative<Y, T>: Mul<Quantity<T>, Output = Y>,
for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
for<'a> &'a Derivative<Y, T>:
Mul<Scalar, Output = Derivative<Y, T>> + Mul<Quantity<T>, Output = Y>,
U: TensorVec<Item = Y>,
V: TensorVec<Item = Derivative<Y, T>>,
{
fn interpolate(
&self,
time: &Times<T>,
tp: &Times<T>,
yp: &U,
dydtp: &V,
_k_sol: &[V],
_function: impl FnMut(Quantity<T>, &Y) -> Result<Derivative<Y, T>, String>,
) -> Result<(U, V), IntegrationError> {
Ok(Self::interpolate_free(time, tp, yp, dydtp))
}
}