rustyqlib/core/solvers/
bisection.rs

1//! Bisection: derivative-free and unconditionally convergent on a
2//! sign-changing bracket, at a linear rate.
3
4use super::solver_1d::{Root, Solver1d};
5use crate::core::errors::RustyQLibError;
6
7/// Bisection on a sign-changing bracket `[lo, hi]`. Errors when the
8/// bracket does not straddle a root.
9pub fn bisection(
10    cfg: &Solver1d,
11    f: impl Fn(f64) -> f64,
12    lo: f64,
13    hi: f64,
14) -> Result<Root, RustyQLibError> {
15    let (mut lo, mut hi) = (lo.min(hi), lo.max(hi));
16    let f_lo = f(lo);
17    if f_lo.abs() <= cfg.tol {
18        return Ok(Root { x: lo, iterations: 0, converged: true });
19    }
20    let f_hi = f(hi);
21    if f_hi.abs() <= cfg.tol {
22        return Ok(Root { x: hi, iterations: 0, converged: true });
23    }
24    if f_lo * f_hi > 0.0 {
25        return Err(RustyQLibError::NumericalError(format!(
26            "bisection needs a sign change: f({lo}) = {f_lo}, f({hi}) = {f_hi}"
27        )));
28    }
29    let mut x = 0.5 * (lo + hi);
30    for i in 1..=cfg.max_iter {
31        let fx = f(x);
32        if fx.abs() <= cfg.tol {
33            return Ok(Root { x, iterations: i, converged: true });
34        }
35        if fx * f_lo < 0.0 {
36            hi = x;
37        } else {
38            lo = x;
39        }
40        x = 0.5 * (lo + hi);
41        if hi - lo <= f64::EPSILON * (1.0 + x.abs()) {
42            return Ok(Root { x, iterations: i, converged: f(x).abs() <= cfg.tol });
43        }
44    }
45    Ok(Root { x, iterations: cfg.max_iter, converged: false })
46}
47
48#[cfg(test)]
49mod tests {
50    use super::*;
51
52    #[test]
53    fn finds_sqrt_two() {
54        let r = bisection(&Solver1d::default(), |x| x * x - 2.0, 0.0, 2.0).unwrap();
55        assert!(r.converged && (r.x - std::f64::consts::SQRT_2).abs() < 1e-9, "{r:?}");
56    }
57
58    #[test]
59    fn rejects_a_bracket_without_sign_change() {
60        assert!(bisection(&Solver1d::default(), |x| x * x - 2.0, 2.0, 3.0).is_err());
61    }
62}
rustyqlib/core/solvers/bisection.rs

rustyqlib/core/solvers/
bisection.rs
