conspire/math/integrate/ode/explicit/
mod.rs

1#[cfg(test)]
2mod test;
3
4use crate::{
5    math::{
6        Derivative, Differentiable, Quantity, Scalar, Tensor, TensorVec,
7        integrate::{IntegrationError, OdeIntegrator, StabilityRegion, Times},
8    },
9    units::Time,
10};
11use std::f64::consts::FRAC_PI_2;
12
13pub(crate) mod fixed_step;
14pub(crate) mod variable_step;
15
16/// The fastest time scale of a right-hand side function, which bounds a stable time step.
17///
18/// The time scale is the reciprocal of the spectral radius of the Jacobian.
19#[derive(Clone, Copy, Debug)]
20pub enum Spectrum<T = Time> {
21    /// The eigenvalues are real and negative, as in diffusion.
22    Real(Quantity<T>),
23    /// The eigenvalues are imaginary, as in waves.
24    Imaginary(Quantity<T>),
25    /// The eigenvalues lie along a ray at the given angle from the negative real axis,
26    /// as in damped waves, between $`0`$ (diffusion) and $`\pi/2`$ (waves).
27    Complex(Quantity<T>, Scalar),
28    /// The eigenvalues are real and negative at the first time scale and imaginary at the
29    /// second, as in diffusion coupled to waves.
30    Mixed(Quantity<T>, Quantity<T>),
31}
32
33impl<T> Spectrum<T> {
34    pub(crate) fn limit(self, stability: &mut Stability) -> Quantity<T> {
35        match self {
36            Self::Real(scale) => scale * stability.extent(0.0),
37            Self::Imaginary(scale) => scale * stability.extent(FRAC_PI_2),
38            Self::Complex(scale, angle) => scale * stability.extent(angle),
39            Self::Mixed(real, imaginary) => {
40                let real = real * stability.extent(0.0);
41                let imaginary = imaginary * stability.extent(FRAC_PI_2);
42                if real < imaginary { real } else { imaginary }
43            }
44        }
45    }
46}
47
48pub(crate) struct Stability {
49    region: StabilityRegion,
50    angle: Scalar,
51    extent: Scalar,
52}
53
54impl Stability {
55    pub(crate) fn new(region: StabilityRegion) -> Self {
56        Self {
57            region,
58            angle: Scalar::NAN,
59            extent: 0.0,
60        }
61    }
62    fn extent(&mut self, angle: Scalar) -> Scalar {
63        if angle != self.angle {
64            self.angle = angle;
65            self.extent = self.region.extent(angle);
66        }
67        self.extent
68    }
69}
70
71pub(crate) fn check_safety(safety: Scalar) -> Result<(), IntegrationError> {
72    if safety > 0.0 && safety <= 1.0 {
73        Ok(())
74    } else {
75        Err(IntegrationError::InvalidSafetyFactor(safety))
76    }
77}
78
79/// Explicit integrators for ordinary differential equations.
80pub trait Explicit<Y, U, V, T = Time>
81where
82    Self: OdeIntegrator<Y, U>,
83    Y: Differentiable<T> + Tensor,
84    U: TensorVec<Item = Y>,
85    V: TensorVec<Item = Derivative<Y, T>>,
86{
87    const SLOPES: usize;
88    #[doc = include_str!("doc.md")]
89    fn integrate(
90        &self,
91        function: impl FnMut(Quantity<T>, &Y) -> Result<Derivative<Y, T>, String>,
92        time: &[Quantity<T>],
93        initial_condition: Y,
94    ) -> Result<(Times<T>, U, V), IntegrationError>;
95    #[doc = include_str!("bounded.md")]
96    fn integrate_bounded(
97        &self,
98        function: impl FnMut(Quantity<T>, &Y) -> Result<Derivative<Y, T>, String>,
99        bound: impl FnMut(Quantity<T>, &Y) -> Result<Spectrum<T>, String>,
100        safety: Scalar,
101        time: &[Quantity<T>],
102        initial_condition: Y,
103    ) -> Result<(Times<T>, U, V), IntegrationError>;
104}
conspire/math/integrate/ode/explicit/mod.rs

conspire/math/integrate/ode/explicit/
mod.rs
