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

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