1use crate::math::{
2 Derivative, Differentiable, Quantity, Scalar, Tensor, TensorVec,
3 integrate::{
4 BogackiShampine, ExplicitDaeVariableStepExplicit, ExplicitDaeVariableStepFirstSameAsLast,
5 FreeInterpolant, IntegrationError, Times,
6 },
7};
8use std::ops::{Div, Mul, Sub};
9
10impl<Y, Z, U, V, W, T> ExplicitDaeVariableStepExplicit<Y, Z, U, V, W, T> for BogackiShampine
11where
12 Self: ExplicitDaeVariableStepFirstSameAsLast<Y, Z, U, V, W, T>,
13 Y: Differentiable<T> + Div<Quantity<T>, Output = Derivative<Y, 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_and_error(
24 &self,
25 evolution: impl FnMut(Quantity<T>, &Y, &Z) -> Result<Derivative<Y, T>, String>,
26 solution: impl FnMut(Quantity<T>, &Y, &Z) -> Result<Z, String>,
27 y: &Y,
28 z: &Z,
29 t: Quantity<T>,
30 dt: Quantity<T>,
31 k: &mut [Derivative<Y, T>],
32 y_trial: &mut Y,
33 z_trial: &mut Z,
34 ) -> Result<Scalar, String> {
35 self.slopes_solve_and_error_fsal(evolution, solution, y, z, t, dt, k, y_trial, z_trial)
36 }
37 fn step_solve(
38 &self,
39 _: impl FnMut(Quantity<T>, &Y, &Z) -> Result<Derivative<Y, T>, String>,
40 y: &mut Y,
41 z: &mut Z,
42 t: &mut Quantity<T>,
43 y_sol: &mut U,
44 z_sol: &mut V,
45 t_sol: &mut Times<T>,
46 dydt_sol: &mut W,
47 k_sol: &mut Vec<W>,
48 dt: &mut Quantity<T>,
49 k: &mut [Derivative<Y, T>],
50 y_trial: &Y,
51 z_trial: &Z,
52 e: Scalar,
53 ) -> Result<(), String> {
54 self.step_solve_fsal(
55 y, z, t, y_sol, z_sol, t_sol, dydt_sol, k_sol, dt, k, y_trial, z_trial, e,
56 )
57 }
58 fn interpolate_explicit_dae_variable_step(
59 &self,
60 _evolution: impl FnMut(Quantity<T>, &Y, &Z) -> Result<Derivative<Y, T>, String>,
61 mut solution: impl FnMut(Quantity<T>, &Y, &Z) -> Result<Z, String>,
62 time: &Times<T>,
63 tp: &Times<T>,
64 yp: &U,
65 dydtp: &W,
66 _k_sol: &[W],
67 zp: &V,
68 ) -> Result<(U, W, V), IntegrationError> {
69 let (y_int, dydt_int) = Self::interpolate_free(time, tp, yp, dydtp);
70 let mut z_int = V::new();
71 for (idx, time_k) in time.iter().enumerate() {
72 let i = tp.iter().position(|tp_i| tp_i >= time_k).unwrap();
73 if time_k == &tp[i] {
74 z_int.push(zp[i].clone());
75 } else {
76 z_int.push(solution(*time_k, &y_int[idx], &zp[i - 1])?);
77 }
78 }
79 Ok((y_int, dydt_int, z_int))
80 }
81}
82
83impl<Y, Z, U, V, W, T> ExplicitDaeVariableStepFirstSameAsLast<Y, Z, U, V, W, T> for BogackiShampine
84where
85 Y: Differentiable<T> + Div<Quantity<T>, Output = Derivative<Y, T>> + Tensor,
86 Z: PartialEq + Tensor,
87 Derivative<Y, T>: Mul<Quantity<T>, Output = Y>,
88 U: TensorVec<Item = Y>,
89 V: TensorVec<Item = Z>,
90 W: TensorVec<Item = Derivative<Y, T>>,
91 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
92 for<'a> &'a Derivative<Y, T>:
93 Mul<Scalar, Output = Derivative<Y, T>> + Mul<Quantity<T>, Output = Y>,
94{
95}