conspire/math/integrate/
mod.rs

1#[cfg(test)]
2mod test;
3
4mod backward_euler;
5mod bogacki_shampine;
6mod dormand_prince;
7mod verner_8;
8mod verner_9;
9
10pub use backward_euler::BackwardEuler;
11pub use bogacki_shampine::BogackiShampine;
12pub use dormand_prince::DormandPrince;
13pub use verner_8::Verner8;
14pub use verner_9::Verner9;
15
16pub type Ode1be = BackwardEuler;
17pub type Ode23 = BogackiShampine;
18pub type Ode45 = DormandPrince;
19pub type Ode78 = Verner8;
20pub type Ode89 = Verner9;
21
22// consider symplectic integrators for dynamics eventually
23
24use super::{
25    Solution, Tensor, TensorArray, TensorRank0, TensorVec, TestError, Vector,
26    interpolate::InterpolateSolution,
27    optimize::{FirstOrderRootFinding, OptimizeError, ZerothOrderRootFinding},
28};
29use crate::defeat_message;
30use std::{
31    fmt::{self, Debug, Display, Formatter},
32    ops::{Div, Mul, Sub},
33};
34
35/// Base trait for ordinary differential equation solvers.
36pub trait OdeSolver<Y, U>
37where
38    Self: Debug,
39    Y: Tensor,
40    U: TensorVec<Item = Y>,
41{
42}
43
44impl<A, Y, U> OdeSolver<Y, U> for A
45where
46    A: Debug,
47    Y: Tensor,
48    U: TensorVec<Item = Y>,
49{
50}
51
52/// Base trait for explicit ordinary differential equation solvers.
53pub trait Explicit<Y, U>: OdeSolver<Y, U>
54where
55    Self: InterpolateSolution<Y, U>,
56    Y: Tensor,
57    for<'a> &'a Y: Mul<TensorRank0, Output = Y> + Sub<&'a Y, Output = Y>,
58    U: TensorVec<Item = Y>,
59{
60    /// Solves an initial value problem by explicitly integrating a system of ordinary differential equations.
61    ///
62    /// ```math
63    /// \frac{dy}{dt} = f(t, y),\quad y(t_0) = y_0
64    /// ```
65    fn integrate(
66        &self,
67        function: impl FnMut(TensorRank0, &Y) -> Result<Y, IntegrationError>,
68        time: &[TensorRank0],
69        initial_condition: Y,
70    ) -> Result<(Vector, U, U), IntegrationError>;
71}
72
73/// Base trait for zeroth-order implicit ordinary differential equation solvers.
74pub trait ImplicitZerothOrder<Y, U>: OdeSolver<Y, U>
75where
76    Self: InterpolateSolution<Y, U>,
77    Y: Solution,
78    for<'a> &'a Y: Mul<TensorRank0, Output = Y> + Sub<&'a Y, Output = Y>,
79    U: TensorVec<Item = Y>,
80{
81    /// Solves an initial value problem by implicitly integrating a system of ordinary differential equations.
82    ///
83    /// ```math
84    /// \frac{dy}{dt} = f(t, y),\quad y(t_0) = y_0,\quad \frac{\partial f}{\partial y} = J(t, y)
85    /// ```
86    fn integrate(
87        &self,
88        function: impl Fn(TensorRank0, &Y) -> Result<Y, IntegrationError>,
89        time: &[TensorRank0],
90        initial_condition: Y,
91        solver: impl ZerothOrderRootFinding<Y>,
92    ) -> Result<(Vector, U, U), IntegrationError>;
93}
94
95/// Base trait for first-order implicit ordinary differential equation solvers.
96pub trait ImplicitFirstOrder<Y, J, U>: OdeSolver<Y, U>
97where
98    Self: InterpolateSolution<Y, U>,
99    Y: Solution + Div<J, Output = Y>,
100    for<'a> &'a Y: Mul<TensorRank0, Output = Y> + Sub<&'a Y, Output = Y>,
101    J: Tensor + TensorArray,
102    U: TensorVec<Item = Y>,
103{
104    /// Solves an initial value problem by implicitly integrating a system of ordinary differential equations.
105    ///
106    /// ```math
107    /// \frac{dy}{dt} = f(t, y),\quad y(t_0) = y_0,\quad \frac{\partial f}{\partial y} = J(t, y)
108    /// ```
109    fn integrate(
110        &self,
111        function: impl Fn(TensorRank0, &Y) -> Result<Y, IntegrationError>,
112        jacobian: impl Fn(TensorRank0, &Y) -> Result<J, IntegrationError>,
113        time: &[TensorRank0],
114        initial_condition: Y,
115        solver: impl FirstOrderRootFinding<Y, J, Y>,
116    ) -> Result<(Vector, U, U), IntegrationError>;
117}
118
119/// Possible errors encountered when integrating.
120pub enum IntegrationError {
121    InitialTimeNotLessThanFinalTime,
122    LengthTimeLessThanTwo,
123}
124
125impl Debug for IntegrationError {
126    fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
127        let error = match self {
128            Self::InitialTimeNotLessThanFinalTime => {
129                "\x1b[1;91mThe initial time must precede the final time.".to_string()
130            }
131            Self::LengthTimeLessThanTwo => {
132                "\x1b[1;91mThe time must contain at least two entries.".to_string()
133            }
134        };
135        write!(f, "\n{}\n\x1b[0;2;31m{}\x1b[0m\n", error, defeat_message())
136    }
137}
138
139impl Display for IntegrationError {
140    fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
141        let error = match self {
142            Self::InitialTimeNotLessThanFinalTime => {
143                "\x1b[1;91mThe initial time must precede the final time.".to_string()
144            }
145            Self::LengthTimeLessThanTwo => {
146                "\x1b[1;91mThe time must contain at least two entries.".to_string()
147            }
148        };
149        write!(f, "{error}\x1b[0m")
150    }
151}
152
153impl From<OptimizeError> for IntegrationError {
154    fn from(_error: OptimizeError) -> Self {
155        todo!()
156    }
157}
158
159impl From<IntegrationError> for TestError {
160    fn from(error: IntegrationError) -> Self {
161        TestError {
162            message: error.to_string(),
163        }
164    }
165}