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}