conspire/math/integrate/dae/explicit/variable_step/dormand_prince/
mod.rs1use 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}