Skip to main content

gridoxide_core/
roots.rs

1//! Bracketed root finding, shared by the event locators in
2//! [`continuation`](crate::continuation) and gridoxide's `dynamics`.
3
4/// Regula falsi with Illinois damping over `[lo, hi]`, given the bracketing
5/// values `f_lo < 0 <= f_hi`.
6///
7/// Each evaluation costs a full corrector in continuation and a full Newton
8/// solve in dynamics, so plain bisection would be wasteful; Illinois damping
9/// keeps the superlinear rate without losing the bracket the way unmodified
10/// regula falsi does when one endpoint stalls.
11pub struct Illinois {
12    pub lo: f64,
13    pub hi: f64,
14    pub f_lo: f64,
15    pub f_hi: f64,
16}
17
18impl Illinois {
19    /// The next trial point.
20    pub fn next(&self) -> f64 {
21        let denom = self.f_hi - self.f_lo;
22        if !denom.is_finite() || denom == 0.0 {
23            return 0.5 * (self.lo + self.hi);
24        }
25        let x = self.hi - self.f_hi * (self.hi - self.lo) / denom;
26        // Never leave the bracket, and never stall on an endpoint.
27        let (a, b) = (self.lo.min(self.hi), self.lo.max(self.hi));
28        let margin = 1e-3 * (b - a);
29        x.clamp(a + margin, b - margin)
30    }
31
32    /// Narrows the bracket with the value measured at `x`.
33    pub fn update(&mut self, x: f64, f_x: f64) {
34        if f_x >= 0.0 {
35            self.hi = x;
36            self.f_hi = f_x;
37            self.f_lo *= 0.5; // Illinois: halve the stale endpoint's weight
38        } else {
39            self.lo = x;
40            self.f_lo = f_x;
41            self.f_hi *= 0.5;
42        }
43    }
44
45    pub fn width(&self) -> f64 {
46        (self.hi - self.lo).abs()
47    }
48}