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

1use crate::math::{
2    Derivative, Differentiate, Quantity, Scalar, Tensor, TensorVec,
3    integrate::{
4        ExplicitDaeVariableStepExplicit, ExplicitDaeVariableStepFirstSameAsLast, IntegrationError,
5        Times, ode::explicit::variable_step::dormand_prince::*,
6    },
7};
8use std::ops::{Mul, Sub};
9
10impl<Y, Z, U, V, W, T> ExplicitDaeVariableStepExplicit<Y, Z, U, V, W, T> for DormandPrince
11where
12    Self: ExplicitDaeVariableStepFirstSameAsLast<Y, Z, U, V, W, T>,
13    Y: Differentiate<T> + Tensor,
14    Z: PartialEq + Tensor,
15    Derivative<Y, T>: Mul<Quantity<T>, Output = Y>,
16    U: TensorVec<Item = Y>,
17    V: TensorVec<Item = Z>,
18    W: TensorVec<Item = Derivative<Y, T>>,
19    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
20    for<'a> &'a Derivative<Y, T>:
21        Mul<Scalar, Output = Derivative<Y, T>> + Mul<Quantity<T>, Output = Y>,
22{
23    fn slopes_solve(
24        mut evolution: impl FnMut(Quantity<T>, &Y, &Z) -> Result<Derivative<Y, T>, String>,
25        mut solution: impl FnMut(Quantity<T>, &Y, &Z) -> Result<Z, String>,
26        y: &Y,
27        z: &Z,
28        t: Quantity<T>,
29        dt: Quantity<T>,
30        k: &mut [Derivative<Y, T>],
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(Quantity<T>, &Y, &Z) -> Result<Derivative<Y, T>, String>,
67        solution: impl FnMut(Quantity<T>, &Y, &Z) -> Result<Z, String>,
68        y: &Y,
69        z: &Z,
70        t: Quantity<T>,
71        dt: Quantity<T>,
72        k: &mut [Derivative<Y, T>],
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(Quantity<T>, &Y, &Z) -> Result<Derivative<Y, T>, String>,
81        y: &mut Y,
82        z: &mut Z,
83        t: &mut Quantity<T>,
84        y_sol: &mut U,
85        z_sol: &mut V,
86        t_sol: &mut Times<T>,
87        dydt_sol: &mut W,
88        k_sol: &mut Vec<W>,
89        dt: &mut Quantity<T>,
90        k: &mut [Derivative<Y, T>],
91        y_trial: &Y,
92        z_trial: &Z,
93        e: Scalar,
94    ) -> Result<(), String> {
95        self.step_solve_fsal(
96            y, z, t, y_sol, z_sol, t_sol, dydt_sol, k_sol, dt, k, y_trial, z_trial, e,
97        )
98    }
99    #[allow(clippy::too_many_arguments)]
100    fn interpolate_explicit_dae_variable_step(
101        &self,
102        _evolution: impl FnMut(Quantity<T>, &Y, &Z) -> Result<Derivative<Y, T>, String>,
103        mut solution: impl FnMut(Quantity<T>, &Y, &Z) -> Result<Z, String>,
104        time: &Times<T>,
105        tp: &Times<T>,
106        yp: &U,
107        dydtp: &W,
108        k_sol: &[W],
109        zp: &V,
110    ) -> Result<(U, W, V), IntegrationError> {
111        let (y_int, dydt_int) = Self::interpolate_free_dense(time, tp, yp, dydtp, k_sol);
112        let mut z_int = V::new();
113        for (idx, time_k) in time.iter().enumerate() {
114            let i = tp.iter().position(|tp_i| tp_i >= time_k).unwrap();
115            if time_k == &tp[i] {
116                z_int.push(zp[i].clone());
117            } else {
118                z_int.push(solution(*time_k, &y_int[idx], &zp[i - 1])?);
119            }
120        }
121        Ok((y_int, dydt_int, z_int))
122    }
123}
124
125impl<Y, Z, U, V, W, T> ExplicitDaeVariableStepFirstSameAsLast<Y, Z, U, V, W, T> for DormandPrince
126where
127    Y: Differentiate<T> + Tensor,
128    Z: PartialEq + Tensor,
129    Derivative<Y, T>: Mul<Quantity<T>, Output = Y>,
130    U: TensorVec<Item = Y>,
131    V: TensorVec<Item = Z>,
132    W: TensorVec<Item = Derivative<Y, T>>,
133    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
134    for<'a> &'a Derivative<Y, T>:
135        Mul<Scalar, Output = Derivative<Y, T>> + Mul<Quantity<T>, Output = Y>,
136{
137}
conspire/math/integrate/dae/explicit/variable_step/dormand_prince/mod.rs

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