conspire/math/integrate/backward_euler/
mod.rs

1#[cfg(test)]
2mod test;
3
4use super::{
5    super::{
6        Hessian, Jacobian, Solution, TensorArray, TensorRank0, 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<TensorRank0, 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(TensorRank0, &Y) -> Result<Y, IntegrationError>,
34        time: &[TensorRank0],
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 = solver.root(
61                |y_trial: &Y| Ok(y_trial - &y - &(&function(t_trial, y_trial)? * dt)),
62                y.clone(),
63                EqualityConstraint::None,
64            )?;
65            t = t_trial;
66            y = y_trial;
67            t_sol.push(t);
68            y_sol.push(y.clone());
69            dydt_sol.push(function(t, &y)?);
70            index += 1;
71        }
72        Ok((t_sol, y_sol, dydt_sol))
73    }
74}
75
76impl<Y, J, U> ImplicitFirstOrder<Y, J, U> for BackwardEuler
77where
78    Self: InterpolateSolution<Y, U>,
79    Y: Jacobian + Solution + Div<J, Output = Y>,
80    for<'a> &'a Y: Mul<TensorRank0, Output = Y> + Sub<&'a Y, Output = Y>,
81    J: Hessian + TensorArray,
82    U: TensorVec<Item = Y>,
83    Vector: From<Y>,
84{
85    fn integrate(
86        &self,
87        function: impl Fn(TensorRank0, &Y) -> Result<Y, IntegrationError>,
88        jacobian: impl Fn(TensorRank0, &Y) -> Result<J, IntegrationError>,
89        time: &[TensorRank0],
90        initial_condition: Y,
91        solver: impl FirstOrderRootFinding<Y, J, Y>,
92    ) -> Result<(Vector, U, U), IntegrationError> {
93        let t_0 = time[0];
94        let t_f = time[time.len() - 1];
95        if time.len() < 2 {
96            return Err(IntegrationError::LengthTimeLessThanTwo);
97        } else if t_0 >= t_f {
98            return Err(IntegrationError::InitialTimeNotLessThanFinalTime);
99        }
100        let mut index = 0;
101        let mut t = t_0;
102        let mut dt;
103        let identity = J::identity();
104        let mut t_sol = Vector::zero(0);
105        t_sol.push(t_0);
106        let mut t_trial;
107        let mut y = initial_condition.clone();
108        let mut y_sol = U::zero(0);
109        y_sol.push(initial_condition.clone());
110        let mut dydt_sol = U::zero(0);
111        dydt_sol.push(function(t, &y.clone())?);
112        let mut y_trial;
113        while t < t_f {
114            t_trial = time[index + 1];
115            dt = t_trial - t;
116            y_trial = solver.root(
117                |y_trial: &Y| Ok(y_trial - &y - &(&function(t_trial, y_trial)? * dt)),
118                |y_trial: &Y| Ok(jacobian(t_trial, y_trial)? * -dt + &identity),
119                y.clone(),
120                EqualityConstraint::None,
121            )?;
122            t = t_trial;
123            y = y_trial;
124            t_sol.push(t);
125            y_sol.push(y.clone());
126            dydt_sol.push(function(t, &y)?);
127            index += 1;
128        }
129        Ok((t_sol, y_sol, dydt_sol))
130    }
131}
132
133impl<Y, U> InterpolateSolution<Y, U> for BackwardEuler
134where
135    Y: Jacobian + Solution + TensorArray,
136    for<'a> &'a Y: Mul<TensorRank0, Output = Y> + Sub<&'a Y, Output = Y>,
137    U: TensorVec<Item = Y>,
138    Vector: From<Y>,
139{
140    fn interpolate(
141        &self,
142        _time: &Vector,
143        _tp: &Vector,
144        _yp: &U,
145        mut _function: impl FnMut(TensorRank0, &Y) -> Result<Y, IntegrationError>,
146    ) -> Result<(U, U), IntegrationError> {
147        unimplemented!()
148    }
149}
conspire/math/integrate/backward_euler/mod.rs

conspire/math/integrate/backward_euler/
mod.rs
