conspire/math/integrate/ode/implicit/midpoint/
mod.rs

1#[cfg(test)]
2mod test;
3
4use crate::math::{
5    Derivative, Differentiable, Quantity, Scalar, Tensor, TensorArray, TensorVec,
6    integrate::{
7        FixedStep, ImplicitFirstOrder, ImplicitZerothOrder, IntegrationError, OdeIntegrator,
8    },
9};
10use std::{
11    fmt::Debug,
12    ops::{Add, Mul, Sub},
13};
14
15#[doc = include_str!("doc.md")]
16#[derive(Debug, Default)]
17pub struct Midpoint {
18    /// Fixed value for the time step.
19    dt: Scalar,
20}
21
22impl<Y, U> OdeIntegrator<Y, U> for Midpoint
23where
24    Y: Tensor,
25    U: TensorVec<Item = Y>,
26{
27}
28
29impl<T> FixedStep<T> for Midpoint {
30    fn dt(&self) -> Quantity<T> {
31        Quantity::new(self.dt)
32    }
33}
34
35impl<Y, U, V, T> ImplicitZerothOrder<Y, U, V, T> for Midpoint
36where
37    Y: Differentiable<T> + Tensor,
38    Derivative<Y, T>: Mul<Quantity<T>, Output = Y>,
39    for<'a> &'a Y: Add<&'a Y, Output = Y> + Sub<&'a Y, Output = Y>,
40    U: TensorVec<Item = Y>,
41    V: TensorVec<Item = Derivative<Y, T>>,
42{
43    fn residual(
44        &self,
45        mut function: impl FnMut(Quantity<T>, &Y) -> Result<Derivative<Y, T>, IntegrationError>,
46        t: Quantity<T>,
47        y: &Y,
48        _t_trial: Quantity<T>,
49        y_trial: &Y,
50        dt: Quantity<T>,
51    ) -> Result<Y, String> {
52        Ok(y_trial - y - function(t + 0.5 * dt, &((y + y_trial) * 0.5))? * dt)
53    }
54}
55
56impl<Y, J, U, V, T> ImplicitFirstOrder<Y, J, U, V, T> for Midpoint
57where
58    Y: Differentiable<T> + Tensor,
59    Derivative<Y, T>: Mul<Quantity<T>, Output = Y>,
60    J: Differentiable<T> + Tensor + TensorArray,
61    Derivative<J, T>: Mul<Quantity<T>, Output = J>,
62    for<'a> &'a Y: Add<&'a Y, Output = Y> + Sub<&'a Y, Output = Y>,
63    U: TensorVec<Item = Y>,
64    V: TensorVec<Item = Derivative<Y, T>>,
65{
66    fn hessian(
67        &self,
68        mut jacobian: impl FnMut(Quantity<T>, &Y) -> Result<Derivative<J, T>, IntegrationError>,
69        t: Quantity<T>,
70        y: &Y,
71        _t_trial: Quantity<T>,
72        y_trial: &Y,
73        dt: Quantity<T>,
74    ) -> Result<J, String> {
75        Ok(J::identity() - jacobian(t + 0.5 * dt, &((y + y_trial) * 0.5))? * (dt * 0.5))
76    }
77}
conspire/math/integrate/ode/implicit/midpoint/mod.rs

conspire/math/integrate/ode/implicit/midpoint/
mod.rs
