conspire/math/integrate/dae/explicit/variable_step/dormand_prince/
mod.rs

1use crate::math::{
2    Banded, Scalar, Tensor, TensorVec, Vector,
3    integrate::{
4        ExplicitDaeVariableStep, ExplicitDaeVariableStepFirstSameAsLast, IntegrationError,
5        ode::explicit::variable_step::dormand_prince::*,
6    },
7    optimize::{
8        EqualityConstraint, FirstOrderOptimization, FirstOrderRootFinding, SecondOrderOptimization,
9        ZerothOrderRootFinding,
10    },
11};
12use std::ops::{Mul, Sub};
13
14impl<Y, Z, U, V> ExplicitDaeVariableStep<Y, Z, U, V> for DormandPrince
15where
16    Self: ExplicitDaeZerothOrderRoot<Y, Z, U, V>,
17    Y: Tensor,
18    Z: Tensor,
19    U: TensorVec<Item = Y>,
20    V: TensorVec<Item = Z>,
21    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
22{
23    fn slopes_solve(
24        mut evolution: impl FnMut(Scalar, &Y, &Z) -> Result<Y, String>,
25        mut solution: impl FnMut(Scalar, &Y, &Z) -> Result<Z, String>,
26        y: &Y,
27        z: &Z,
28        t: Scalar,
29        dt: Scalar,
30        k: &mut [Y],
31        y_trial: &mut Y,
32        z_trial: &mut Z,
33    ) -> Result<(), String> {
34        *y_trial = &k[0] * (0.2 * dt) + y;
35        *z_trial = solution(t + 0.2 * dt, y_trial, z)?;
36        k[1] = evolution(t + 0.2 * dt, y_trial, z_trial)?;
37        *y_trial = &k[0] * (0.075 * dt) + &k[1] * (0.225 * dt) + y;
38        *z_trial = solution(t + 0.3 * dt, y_trial, z_trial)?;
39        k[2] = evolution(t + 0.3 * dt, y_trial, z_trial)?;
40        *y_trial = &k[0] * (C_44_45 * dt) - &k[1] * (C_56_15 * dt) + &k[2] * (C_32_9 * dt) + y;
41        *z_trial = solution(t + 0.8 * dt, y_trial, z_trial)?;
42        k[3] = evolution(t + 0.8 * dt, y_trial, z_trial)?;
43        *y_trial = &k[0] * (C_19372_6561 * dt) - &k[1] * (C_25360_2187 * dt)
44            + &k[2] * (C_64448_6561 * dt)
45            - &k[3] * (C_212_729 * dt)
46            + y;
47        *z_trial = solution(t + C_8_9 * dt, y_trial, z_trial)?;
48        k[4] = evolution(t + C_8_9 * dt, y_trial, z_trial)?;
49        *y_trial = &k[0] * (C_9017_3168 * dt) - &k[1] * (C_355_33 * dt)
50            + &k[2] * (C_46732_5247 * dt)
51            + &k[3] * (C_49_176 * dt)
52            - &k[4] * (C_5103_18656 * dt)
53            + y;
54        *z_trial = solution(t + dt, y_trial, z_trial)?;
55        k[5] = evolution(t + dt, y_trial, z_trial)?;
56        *y_trial = (&k[0] * C_35_384 + &k[2] * C_500_1113 + &k[3] * C_125_192
57            - &k[4] * C_2187_6784
58            + &k[5] * C_11_84)
59            * dt
60            + y;
61        *z_trial = solution(t + dt, y_trial, z_trial)?;
62        Ok(())
63    }
64    fn slopes_solve_and_error(
65        &self,
66        evolution: impl FnMut(Scalar, &Y, &Z) -> Result<Y, String>,
67        solution: impl FnMut(Scalar, &Y, &Z) -> Result<Z, String>,
68        y: &Y,
69        z: &Z,
70        t: Scalar,
71        dt: Scalar,
72        k: &mut [Y],
73        y_trial: &mut Y,
74        z_trial: &mut Z,
75    ) -> Result<Scalar, String> {
76        Self::slopes_solve_and_error_fsal(evolution, solution, y, z, t, dt, k, y_trial, z_trial)
77    }
78    fn step_solve(
79        &self,
80        _: impl FnMut(Scalar, &Y, &Z) -> Result<Y, String>,
81        y: &mut Y,
82        z: &mut Z,
83        t: &mut Scalar,
84        y_sol: &mut U,
85        z_sol: &mut V,
86        t_sol: &mut Vector,
87        dydt_sol: &mut U,
88        dt: &mut Scalar,
89        k: &mut [Y],
90        y_trial: &Y,
91        z_trial: &Z,
92        e: Scalar,
93    ) -> Result<(), String> {
94        self.step_solve_fsal(
95            y, z, t, y_sol, z_sol, t_sol, dydt_sol, dt, k, y_trial, z_trial, e,
96        )
97    }
98}
99
100impl<Y, Z, U, V> ExplicitDaeVariableStepFirstSameAsLast<Y, Z, U, V> for DormandPrince
101where
102    Y: Tensor,
103    Z: Tensor,
104    U: TensorVec<Item = Y>,
105    V: TensorVec<Item = Z>,
106    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
107{
108}
109
110super::implement_solvers!(DormandPrince);