conspire/math/integrate/verner_8/

mod.rs

#[cfg(test)]
mod test;

use super::{
    super::{Scalar, Tensor, TensorVec, Vector, interpolate::InterpolateSolution},
    Explicit, IntegrationError,
};
use crate::{ABS_TOL, REL_TOL};
use std::ops::{Mul, Sub};

const C_2: Scalar = 0.05;
const C_3: Scalar = 0.1065625;
const C_4: Scalar = 0.15984375;
const C_5: Scalar = 0.39;
const C_6: Scalar = 0.465;
const C_7: Scalar = 0.155;
const C_8: Scalar = 0.943;
const C_9: Scalar = 0.901802041735857;
const C_10: Scalar = 0.909;
const C_11: Scalar = 0.94;

const A_2_1: Scalar = 0.05;
const A_3_1: Scalar = -0.0069931640625;
const A_3_2: Scalar = 0.1135556640625;
const A_4_1: Scalar = 0.0399609375;
const A_4_3: Scalar = 0.1198828125;
const A_5_1: Scalar = 0.36139756280045754;
const A_5_3: Scalar = -1.3415240667004928;
const A_5_4: Scalar = 1.3701265039000352;
const A_6_1: Scalar = 0.049047202797202795;
const A_6_4: Scalar = 0.23509720422144048;
const A_6_5: Scalar = 0.18085559298135673;
const A_7_1: Scalar = 0.06169289044289044;
const A_7_4: Scalar = 0.11236568314640277;
const A_7_5: Scalar = -0.03885046071451367;
const A_7_6: Scalar = 0.01979188712522046;
const A_8_1: Scalar = -1.767630240222327;
const A_8_4: Scalar = -62.5;
const A_8_5: Scalar = -6.061889377376669;
const A_8_6: Scalar = 5.6508231982227635;
const A_8_7: Scalar = 65.62169641937624;
const A_9_1: Scalar = -1.1809450665549708;
const A_9_4: Scalar = -41.50473441114321;
const A_9_5: Scalar = -4.434438319103725;
const A_9_6: Scalar = 4.260408188586133;
const A_9_7: Scalar = 43.75364022446172;
const A_9_8: Scalar = 0.00787142548991231;
const A_10_1: Scalar = -1.2814059994414884;
const A_10_4: Scalar = -45.047139960139866;
const A_10_5: Scalar = -4.731362069449576;
const A_10_6: Scalar = 4.514967016593808;
const A_10_7: Scalar = 47.44909557172985;
const A_10_8: Scalar = 0.01059228297111661;
const A_10_9: Scalar = -0.0057468422638446166;
const A_11_1: Scalar = -1.7244701342624853;
const A_11_4: Scalar = -60.92349008483054;
const A_11_5: Scalar = -5.951518376222392;
const A_11_6: Scalar = 5.556523730698456;
const A_11_7: Scalar = 63.98301198033305;
const A_11_8: Scalar = 0.014642028250414961;
const A_11_9: Scalar = 0.06460408772358203;
const A_11_10: Scalar = -0.0793032316900888;
const A_12_1: Scalar = -3.301622667747079;
const A_12_4: Scalar = -118.01127235975251;
const A_12_5: Scalar = -10.141422388456112;
const A_12_6: Scalar = 9.139311332232058;
const A_12_7: Scalar = 123.37594282840426;
const A_12_8: Scalar = 4.62324437887458;
const A_12_9: Scalar = -3.3832777380682018;
const A_12_10: Scalar = 4.527592100324618;
const A_12_11: Scalar = -5.828495485811623;
const A_13_1: Scalar = -3.039515033766309;
const A_13_4: Scalar = -109.26086808941763;
const A_13_5: Scalar = -9.290642497400293;
const A_13_6: Scalar = 8.43050498176491;
const A_13_7: Scalar = 114.20100103783314;
const A_13_8: Scalar = -0.9637271342145479;
const A_13_9: Scalar = -5.0348840888021895;
const A_13_10: Scalar = 5.958130824002923;

const B_1: Scalar = 0.04427989419007951;
const B_6: Scalar = 0.3541049391724449;
const B_7: Scalar = 0.24796921549564377;
const B_8: Scalar = -15.694202038838085;
const B_9: Scalar = 25.084064965558564;
const B_10: Scalar = -31.738367786260277;
const B_11: Scalar = 22.938283273988784;
const B_12: Scalar = -0.2361324633071542;

const D_1: Scalar = -0.00003272103901028138;
const D_6: Scalar = -0.0005046250618777704;
const D_7: Scalar = 0.0001211723589784759;
const D_8: Scalar = -20.142336771313868;
const D_9: Scalar = 5.2371785994398286;
const D_10: Scalar = -8.156744408794658;
const D_11: Scalar = 22.938283273988784;
const D_12: Scalar = -0.2361324633071542;
const D_13: Scalar = 0.36016794372897754;

/// Explicit, thirteen-stage, eighth-order, variable-step, Runge-Kutta method.[^cite]
///
/// [^cite]: J.H. Verner, [Numer. Algor. **53**, 383 (2010)](https://doi.org/10.1007/s11075-009-9290-3).
///
/// ```math
/// \frac{dy}{dt} = f(t, y)
/// ```
/// ```math
/// t_{n+1} = t_n + h
/// ```
/// ```math
/// k_1 = f(t_n, y_n)
/// ```
/// ```math
/// \cdots
/// ```
/// ```math
/// h_{n+1} = \beta h \left(\frac{e_\mathrm{tol}}{e_{n+1}}\right)^{1/p}
/// ```
#[derive(Debug)]
pub struct Verner8 {
    /// Absolute error tolerance.
    pub abs_tol: Scalar,
    /// Relative error tolerance.
    pub rel_tol: Scalar,
    /// Multiplier for adaptive time steps.
    pub dt_beta: Scalar,
    /// Exponent for adaptive time steps.
    pub dt_expn: Scalar,
}

impl Default for Verner8 {
    fn default() -> Self {
        Self {
            abs_tol: ABS_TOL,
            rel_tol: REL_TOL,
            dt_beta: 0.9,
            dt_expn: 8.0,
        }
    }
}

impl<Y, U> Explicit<Y, U> for Verner8
where
    Self: InterpolateSolution<Y, U>,
    Y: Tensor,
    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
    U: TensorVec<Item = Y>,
{
    const SLOPES: usize = 13;
    fn slopes(
        &self,
        mut function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
        y: &Y,
        t: &Scalar,
        dt: &Scalar,
        k: &mut [Y],
        y_trial: &mut Y,
    ) -> Result<Scalar, String> {
        k[0] = function(*t, y)?;
        k[1] = function(t + C_2 * dt, &(&k[0] * (A_2_1 * dt) + y))?;
        k[2] = function(
            t + C_3 * dt,
            &(&k[0] * (A_3_1 * dt) + &k[1] * (A_3_2 * dt) + y),
        )?;
        k[3] = function(
            t + C_4 * dt,
            &(&k[0] * (A_4_1 * dt) + &k[2] * (A_4_3 * dt) + y),
        )?;
        k[4] = function(
            t + C_5 * dt,
            &(&k[0] * (A_5_1 * dt) + &k[2] * (A_5_3 * dt) + &k[3] * (A_5_4 * dt) + y),
        )?;
        k[5] = function(
            t + C_6 * dt,
            &(&k[0] * (A_6_1 * dt) + &k[3] * (A_6_4 * dt) + &k[4] * (A_6_5 * dt) + y),
        )?;
        k[6] = function(
            t + C_7 * dt,
            &(&k[0] * (A_7_1 * dt)
                + &k[3] * (A_7_4 * dt)
                + &k[4] * (A_7_5 * dt)
                + &k[5] * (A_7_6 * dt)
                + y),
        )?;
        k[7] = function(
            t + C_8 * dt,
            &(&k[0] * (A_8_1 * dt)
                + &k[3] * (A_8_4 * dt)
                + &k[4] * (A_8_5 * dt)
                + &k[5] * (A_8_6 * dt)
                + &k[6] * (A_8_7 * dt)
                + y),
        )?;
        k[8] = function(
            t + C_9 * dt,
            &(&k[0] * (A_9_1 * dt)
                + &k[3] * (A_9_4 * dt)
                + &k[4] * (A_9_5 * dt)
                + &k[5] * (A_9_6 * dt)
                + &k[6] * (A_9_7 * dt)
                + &k[7] * (A_9_8 * dt)
                + y),
        )?;
        k[9] = function(
            t + C_10 * dt,
            &(&k[0] * (A_10_1 * dt)
                + &k[3] * (A_10_4 * dt)
                + &k[4] * (A_10_5 * dt)
                + &k[5] * (A_10_6 * dt)
                + &k[6] * (A_10_7 * dt)
                + &k[7] * (A_10_8 * dt)
                + &k[8] * (A_10_9 * dt)
                + y),
        )?;
        k[10] = function(
            t + C_11 * dt,
            &(&k[0] * (A_11_1 * dt)
                + &k[3] * (A_11_4 * dt)
                + &k[4] * (A_11_5 * dt)
                + &k[5] * (A_11_6 * dt)
                + &k[6] * (A_11_7 * dt)
                + &k[7] * (A_11_8 * dt)
                + &k[8] * (A_11_9 * dt)
                + &k[9] * (A_11_10 * dt)
                + y),
        )?;
        k[11] = function(
            t + dt,
            &(&k[0] * (A_12_1 * dt)
                + &k[3] * (A_12_4 * dt)
                + &k[4] * (A_12_5 * dt)
                + &k[5] * (A_12_6 * dt)
                + &k[6] * (A_12_7 * dt)
                + &k[7] * (A_12_8 * dt)
                + &k[8] * (A_12_9 * dt)
                + &k[9] * (A_12_10 * dt)
                + &k[10] * (A_12_11 * dt)
                + y),
        )?;
        *y_trial = (&k[0] * B_1
            + &k[5] * B_6
            + &k[6] * B_7
            + &k[7] * B_8
            + &k[8] * B_9
            + &k[9] * B_10
            + &k[10] * B_11
            + &k[11] * B_12)
            * *dt
            + y;
        k[12] = function(
            t + dt,
            &(&k[0] * (A_13_1 * dt)
                + &k[3] * (A_13_4 * dt)
                + &k[4] * (A_13_5 * dt)
                + &k[5] * (A_13_6 * dt)
                + &k[6] * (A_13_7 * dt)
                + &k[7] * (A_13_8 * dt)
                + &k[8] * (A_13_9 * dt)
                + &k[9] * (A_13_10 * dt)
                + y),
        )?;
        Ok(((&k[0] * D_1
            + &k[5] * D_6
            + &k[6] * D_7
            + &k[7] * D_8
            + &k[8] * D_9
            + &k[9] * D_10
            + &k[10] * D_11
            + &k[11] * D_12
            + &k[12] * D_13)
            * *dt)
            .norm_inf())
    }
    fn step(
        &self,
        mut function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
        y: &mut Y,
        t: &mut Scalar,
        y_sol: &mut U,
        t_sol: &mut Vector,
        dydt_sol: &mut U,
        dt: &mut Scalar,
        _k: &mut [Y],
        y_trial: &Y,
        e: &Scalar,
    ) -> Result<(), String> {
        if e < &self.abs_tol || e / y_trial.norm_inf() < self.rel_tol {
            *t += *dt;
            *y = y_trial.clone();
            t_sol.push(*t);
            y_sol.push(y.clone());
            dydt_sol.push(function(*t, y)?);
        }
        if e > &0.0 {
            *dt *= self.dt_beta * (self.abs_tol / e).powf(1.0 / self.dt_expn)
        }
        Ok(())
    }
}

impl<Y, U> InterpolateSolution<Y, U> for Verner8
where
    Y: Tensor,
    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
    U: TensorVec<Item = Y>,
{
    fn interpolate(
        &self,
        time: &Vector,
        tp: &Vector,
        yp: &U,
        mut function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
    ) -> Result<(U, U), IntegrationError> {
        let mut dt;
        let mut i;
        let mut k_1;
        let mut k_2;
        let mut k_3;
        let mut k_4;
        let mut k_5;
        let mut k_6;
        let mut k_7;
        let mut k_8;
        let mut k_9;
        let mut k_10;
        let mut k_11;
        let mut k_12;
        let mut t;
        let mut y;
        let mut y_int = U::new();
        let mut dydt_int = U::new();
        let mut y_trial;
        for time_k in time.iter() {
            i = tp.iter().position(|tp_i| tp_i >= time_k).unwrap();
            if time_k == &tp[i] {
                t = tp[i];
                y_trial = yp[i].clone();
                dt = 0.0;
            } else {
                t = tp[i - 1];
                y = yp[i - 1].clone();
                dt = time_k - t;
                k_1 = function(t, &y)?;
                k_2 = function(t + C_2 * dt, &(&k_1 * (A_2_1 * dt) + &y))?;
                k_3 = function(
                    t + C_3 * dt,
                    &(&k_1 * (A_3_1 * dt) + &k_2 * (A_3_2 * dt) + &y),
                )?;
                k_4 = function(
                    t + C_4 * dt,
                    &(&k_1 * (A_4_1 * dt) + &k_3 * (A_4_3 * dt) + &y),
                )?;
                k_5 = function(
                    t + C_5 * dt,
                    &(&k_1 * (A_5_1 * dt) + &k_3 * (A_5_3 * dt) + &k_4 * (A_5_4 * dt) + &y),
                )?;
                k_6 = function(
                    t + C_6 * dt,
                    &(&k_1 * (A_6_1 * dt) + &k_4 * (A_6_4 * dt) + &k_5 * (A_6_5 * dt) + &y),
                )?;
                k_7 = function(
                    t + C_7 * dt,
                    &(&k_1 * (A_7_1 * dt)
                        + &k_4 * (A_7_4 * dt)
                        + &k_5 * (A_7_5 * dt)
                        + &k_6 * (A_7_6 * dt)
                        + &y),
                )?;
                k_8 = function(
                    t + C_8 * dt,
                    &(&k_1 * (A_8_1 * dt)
                        + &k_4 * (A_8_4 * dt)
                        + &k_5 * (A_8_5 * dt)
                        + &k_6 * (A_8_6 * dt)
                        + &k_7 * (A_8_7 * dt)
                        + &y),
                )?;
                k_9 = function(
                    t + C_9 * dt,
                    &(&k_1 * (A_9_1 * dt)
                        + &k_4 * (A_9_4 * dt)
                        + &k_5 * (A_9_5 * dt)
                        + &k_6 * (A_9_6 * dt)
                        + &k_7 * (A_9_7 * dt)
                        + &k_8 * (A_9_8 * dt)
                        + &y),
                )?;
                k_10 = function(
                    t + C_10 * dt,
                    &(&k_1 * (A_10_1 * dt)
                        + &k_4 * (A_10_4 * dt)
                        + &k_5 * (A_10_5 * dt)
                        + &k_6 * (A_10_6 * dt)
                        + &k_7 * (A_10_7 * dt)
                        + &k_8 * (A_10_8 * dt)
                        + &k_9 * (A_10_9 * dt)
                        + &y),
                )?;
                k_11 = function(
                    t + C_11 * dt,
                    &(&k_1 * (A_11_1 * dt)
                        + &k_4 * (A_11_4 * dt)
                        + &k_5 * (A_11_5 * dt)
                        + &k_6 * (A_11_6 * dt)
                        + &k_7 * (A_11_7 * dt)
                        + &k_8 * (A_11_8 * dt)
                        + &k_9 * (A_11_9 * dt)
                        + &k_10 * (A_11_10 * dt)
                        + &y),
                )?;
                k_12 = function(
                    t + dt,
                    &(&k_1 * (A_12_1 * dt)
                        + &k_4 * (A_12_4 * dt)
                        + &k_5 * (A_12_5 * dt)
                        + &k_6 * (A_12_6 * dt)
                        + &k_7 * (A_12_7 * dt)
                        + &k_8 * (A_12_8 * dt)
                        + &k_9 * (A_12_9 * dt)
                        + &k_10 * (A_12_10 * dt)
                        + &k_11 * (A_12_11 * dt)
                        + &y),
                )?;
                y_trial = (&k_1 * B_1
                    + &k_6 * B_6
                    + &k_7 * B_7
                    + &k_8 * B_8
                    + &k_9 * B_9
                    + &k_10 * B_10
                    + &k_11 * B_11
                    + &k_12 * B_12)
                    * dt
                    + &y;
            }
            dydt_int.push(function(t + dt, &y_trial)?);
            y_int.push(y_trial);
        }
        Ok((y_int, dydt_int))
    }
}