conspire/math/integrate/backward_euler/
mod.rs

1#[cfg(test)]
2mod test;
3
4use super::{
5    super::{
6        Hessian, Jacobian, Scalar, Solution, TensorArray, TensorVec, Vector,
7        interpolate::InterpolateSolution,
8        optimize::{EqualityConstraint, FirstOrderRootFinding, ZerothOrderRootFinding},
9    },
10    ImplicitFirstOrder, ImplicitZerothOrder, IntegrationError,
11};
12use std::{
13    fmt::Debug,
14    ops::{Div, Mul, Sub},
15};
16
17/// Implicit, single-stage, first-order, fixed-step, Runge-Kutta method.[^cite]
18///
19/// [^cite]: Also known as the backward Euler method.
20#[derive(Debug, Default)]
21pub struct BackwardEuler {}
22
23impl<Y, U> ImplicitZerothOrder<Y, U> for BackwardEuler
24where
25    Self: InterpolateSolution<Y, U>,
26    Y: Solution,
27    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
28    U: TensorVec<Item = Y>,
29    Vector: From<Y>,
30{
31    fn integrate(
32        &self,
33        function: impl Fn(Scalar, &Y) -> Result<Y, IntegrationError>,
34        time: &[Scalar],
35        initial_condition: Y,
36        solver: impl ZerothOrderRootFinding<Y>,
37    ) -> Result<(Vector, U, U), IntegrationError> {
38        let t_0 = time[0];
39        let t_f = time[time.len() - 1];
40        if time.len() < 2 {
41            return Err(IntegrationError::LengthTimeLessThanTwo);
42        } else if t_0 >= t_f {
43            return Err(IntegrationError::InitialTimeNotLessThanFinalTime);
44        }
45        let mut index = 0;
46        let mut t = t_0;
47        let mut dt;
48        let mut t_sol = Vector::zero(0);
49        t_sol.push(t_0);
50        let mut t_trial;
51        let mut y = initial_condition.clone();
52        let mut y_sol = U::zero(0);
53        y_sol.push(initial_condition.clone());
54        let mut dydt_sol = U::zero(0);
55        dydt_sol.push(function(t, &y.clone())?);
56        let mut y_trial;
57        while t < t_f {
58            t_trial = time[index + 1];
59            dt = t_trial - t;
60            y_trial = match solver.root(
61                |y_trial: &Y| Ok(y_trial - &y - &(&function(t_trial, y_trial)? * dt)),
62                y.clone(),
63                EqualityConstraint::None,
64            ) {
65                Ok(solution) => solution,
66                Err(error) => {
67                    return Err(IntegrationError::Upstream(
68                        format!("{error:?}"),
69                        format!("{self:?}"),
70                    ));
71                }
72            };
73            t = t_trial;
74            y = y_trial;
75            t_sol.push(t);
76            y_sol.push(y.clone());
77            dydt_sol.push(function(t, &y)?);
78            index += 1;
79        }
80        Ok((t_sol, y_sol, dydt_sol))
81    }
82}
83
84impl<Y, J, U> ImplicitFirstOrder<Y, J, U> for BackwardEuler
85where
86    Self: InterpolateSolution<Y, U>,
87    Y: Jacobian + Solution + Div<J, Output = Y>,
88    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
89    J: Hessian + TensorArray,
90    U: TensorVec<Item = Y>,
91    Vector: From<Y>,
92{
93    fn integrate(
94        &self,
95        function: impl Fn(Scalar, &Y) -> Result<Y, IntegrationError>,
96        jacobian: impl Fn(Scalar, &Y) -> Result<J, IntegrationError>,
97        time: &[Scalar],
98        initial_condition: Y,
99        solver: impl FirstOrderRootFinding<Y, J, Y>,
100    ) -> Result<(Vector, U, U), IntegrationError> {
101        let t_0 = time[0];
102        let t_f = time[time.len() - 1];
103        if time.len() < 2 {
104            return Err(IntegrationError::LengthTimeLessThanTwo);
105        } else if t_0 >= t_f {
106            return Err(IntegrationError::InitialTimeNotLessThanFinalTime);
107        }
108        let mut index = 0;
109        let mut t = t_0;
110        let mut dt;
111        let identity = J::identity();
112        let mut t_sol = Vector::zero(0);
113        t_sol.push(t_0);
114        let mut t_trial;
115        let mut y = initial_condition.clone();
116        let mut y_sol = U::zero(0);
117        y_sol.push(initial_condition.clone());
118        let mut dydt_sol = U::zero(0);
119        dydt_sol.push(function(t, &y.clone())?);
120        let mut y_trial;
121        while t < t_f {
122            t_trial = time[index + 1];
123            dt = t_trial - t;
124            y_trial = match solver.root(
125                |y_trial: &Y| Ok(y_trial - &y - &(&function(t_trial, y_trial)? * dt)),
126                |y_trial: &Y| Ok(jacobian(t_trial, y_trial)? * -dt + &identity),
127                y.clone(),
128                EqualityConstraint::None,
129            ) {
130                Ok(solution) => solution,
131                Err(error) => {
132                    return Err(IntegrationError::Upstream(
133                        format!("{error:?}"),
134                        format!("{self:?}"),
135                    ));
136                }
137            };
138            t = t_trial;
139            y = y_trial;
140            t_sol.push(t);
141            y_sol.push(y.clone());
142            dydt_sol.push(function(t, &y)?);
143            index += 1;
144        }
145        Ok((t_sol, y_sol, dydt_sol))
146    }
147}
148
149impl<Y, U> InterpolateSolution<Y, U> for BackwardEuler
150where
151    Y: Jacobian + Solution + TensorArray,
152    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
153    U: TensorVec<Item = Y>,
154    Vector: From<Y>,
155{
156    fn interpolate(
157        &self,
158        _time: &Vector,
159        _tp: &Vector,
160        _yp: &U,
161        mut _function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
162    ) -> Result<(U, U), IntegrationError> {
163        unimplemented!()
164    }
165}
conspire/math/integrate/backward_euler/mod.rs

conspire/math/integrate/backward_euler/
mod.rs
