ordinary_diffeq/
controller.rs

1#[derive(Debug, Clone, Copy, PartialEq)]
2pub enum TryStep {
3    Accepted(f64, f64),
4    NotYetAccepted(f64),
5}
6
7impl TryStep {
8    pub fn extract(&self) -> f64 {
9        match self {
10            TryStep::Accepted(h, _) => *h,
11            TryStep::NotYetAccepted(h) => *h,
12        }
13    }
14
15    pub fn is_accepted(&self) -> bool {
16        matches!(self, TryStep::Accepted(_, _))
17    }
18
19    pub fn reset(&mut self) -> Result<TryStep, &str> {
20        match self {
21            TryStep::Accepted(_, h) => Ok(TryStep::NotYetAccepted(*h)),
22            TryStep::NotYetAccepted(_) => Err("Cannot reset a NotYetAccepted TryStep"),
23        }
24    }
25}
26
27pub trait Controller<const D: usize> {
28    fn determine_step(&mut self, h: f64, err: f64) -> TryStep;
29}
30
31#[derive(Debug, Clone, Copy)]
32pub struct PIController {
33    pub alpha: f64,
34    pub beta: f64,
35    pub factor_c1: f64,
36    pub factor_c2: f64,
37    pub factor_old: f64,
38    pub h_max: f64,
39    pub safety_factor: f64,
40    pub next_step_guess: TryStep,
41}
42
43impl<const D: usize> Controller<D> for PIController {
44    /// Determines if the previously run step size and error were valid or not. Either way, it also
45    /// returns what the next step size should be
46    fn determine_step(&mut self, prev_step: f64, err: f64) -> TryStep {
47        let factor_11 = err.powf(self.alpha);
48        let factor = self.factor_c2.max(
49            self.factor_c1
50                .min(factor_11 * self.factor_old.powf(-self.beta) / self.safety_factor),
51        );
52        if err <= 1.0 {
53            let mut h = prev_step / factor;
54            // Accept the stepsize and provide what the next step size should be
55            self.factor_old = err.max(1.0e-4);
56            if h.abs() > self.h_max {
57                // If the step goes past the maximum allowed, though, we shrink it
58                h = self.h_max.copysign(h);
59            }
60            TryStep::Accepted(prev_step, h)
61        } else {
62            // Reject the stepsize and propose a smaller one for the current step
63            TryStep::NotYetAccepted(prev_step / (self.factor_c1.min(factor_11 / self.safety_factor)))
64        }
65    }
66}
67
68impl PIController {
69    pub fn new(
70        alpha: f64,
71        beta: f64,
72        max_factor: f64,
73        min_factor: f64,
74        h_max: f64,
75        safety_factor: f64,
76        initial_h: f64,
77    ) -> Self {
78        Self {
79            alpha,
80            beta,
81            factor_c1: 1.0 / min_factor,
82            factor_c2: 1.0 / max_factor,
83            factor_old: 1.0e-4,
84            h_max: h_max.abs(),
85            safety_factor,
86            next_step_guess: TryStep::NotYetAccepted(initial_h),
87        }
88    }
89}
90
91impl Default for PIController {
92    fn default() -> Self {
93        Self::new(0.17, 0.04, 10.0, 0.2, 100000.0, 0.9, 1e-4)
94    }
95}
96
97#[cfg(test)]
98mod tests {
99    use super::*;
100
101    #[test]
102    fn test_controller_creation() {
103        let controller = PIController::new(0.17, 0.04, 10.0, 0.2, 10.0, 0.9, 1e-4);
104
105        assert!(controller.alpha == 0.17);
106        assert!(controller.beta == 0.04);
107        assert!(controller.factor_c1 == 1.0 / 0.2);
108        assert!(controller.factor_c2 == 1.0 / 10.0);
109        assert!(controller.factor_old == 1.0e-4);
110        assert!(controller.h_max == 10.0);
111        assert!(controller.safety_factor == 0.9);
112        assert!(controller.next_step_guess == TryStep::NotYetAccepted(1e-4));
113    }
114}
ordinary_diffeq/controller.rs

ordinary_diffeq/
controller.rs
