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