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    Tensor, TensorArray, TensorRank0, TensorVec, Vector, interpolate::InterpolateSolution,
26    optimize::OptimizeError,
27};
28use crate::defeat_message;
29use std::{
30    fmt::{self, Debug, Display, Formatter},
31    ops::{Div, Mul, Sub},
32};
33
34/// Base trait for ordinary differential equation solvers.
35pub trait OdeSolver<Y, U>
36where
37    Self: Debug,
38    Y: Tensor,
39    U: TensorVec<Item = Y>,
40{
41}
42
43impl<A, Y, U> OdeSolver<Y, U> for A
44where
45    A: Debug,
46    Y: Tensor,
47    U: TensorVec<Item = Y>,
48{
49}
50
51/// Base trait for explicit ordinary differential equation solvers.
52pub trait Explicit<Y, U>: OdeSolver<Y, U>
53where
54    Self: InterpolateSolution<Y, U>,
55    Y: Tensor,
56    for<'a> &'a Y: Mul<TensorRank0, Output = Y> + Sub<&'a Y, Output = Y>,
57    U: TensorVec<Item = Y>,
58{
59    /// Solves an initial value problem by explicitly integrating a system of ordinary differential equations.
60    ///
61    /// ```math
62    /// \frac{dy}{dt} = f(t, y),\quad y(t_0) = y_0
63    /// ```
64    fn integrate(
65        &self,
66        function: impl Fn(TensorRank0, &Y) -> Result<Y, IntegrationError>,
67        time: &[TensorRank0],
68        initial_condition: Y,
69    ) -> Result<(Vector, U, U), IntegrationError>;
70}
71
72/// Base trait for implicit ordinary differential equation solvers.
73pub trait Implicit<Y, J, U>: OdeSolver<Y, U>
74where
75    Self: InterpolateSolution<Y, U>,
76    Y: Tensor + Div<J, Output = Y>,
77    for<'a> &'a Y: Mul<TensorRank0, Output = Y> + Sub<&'a Y, Output = Y>,
78    J: Tensor + TensorArray,
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        jacobian: impl Fn(TensorRank0, &Y) -> Result<J, IntegrationError>,
90        time: &[TensorRank0],
91        initial_condition: Y,
92    ) -> Result<(Vector, U, U), IntegrationError>;
93}
94
95/// Possible errors encountered when integrating.
96pub enum IntegrationError {
97    InitialTimeNotLessThanFinalTime,
98    LengthTimeLessThanTwo,
99}
100
101impl Debug for IntegrationError {
102    fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
103        let error = match self {
104            Self::InitialTimeNotLessThanFinalTime => {
105                "\x1b[1;91mThe initial time must precede the final time.".to_string()
106            }
107            Self::LengthTimeLessThanTwo => {
108                "\x1b[1;91mThe time must contain at least two entries.".to_string()
109            }
110        };
111        write!(f, "\n{}\n\x1b[0;2;31m{}\x1b[0m\n", error, defeat_message())
112    }
113}
114
115impl Display for IntegrationError {
116    fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
117        let error = match self {
118            Self::InitialTimeNotLessThanFinalTime => {
119                "\x1b[1;91mThe initial time must precede the final time.".to_string()
120            }
121            Self::LengthTimeLessThanTwo => {
122                "\x1b[1;91mThe time must contain at least two entries.".to_string()
123            }
124        };
125        write!(f, "{}\x1b[0m", error)
126    }
127}
128
129impl From<OptimizeError> for IntegrationError {
130    fn from(_error: OptimizeError) -> Self {
131        todo!()
132    }
133}