conspire/math/integrate/backward_euler/
mod.rs1#[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#[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}