rustyqlib/core/solvers/
solver_1d.rs

1//! The shared solver types and the pluggable method dispatch.
2//!
3//! [`Solver1d`] is the configuration (tolerance + iteration cap) every
4//! algorithm takes; its inherent methods delegate to the per-algorithm
5//! files for direct, statically-chosen calls. [`Problem`] + [`Method`] are
6//! the pluggable layer: describe the objective once, then swap algorithms
7//! with an enum value.
8
9use super::{bisection, halley, newton_raphson, newton_safeguarded, secant};
10use crate::core::errors::RustyQLibError;
11
12/// Result of a 1-D root search.
13#[derive(Debug, Clone, Copy)]
14pub struct Root {
15    /// The best estimate of the root.
16    pub x: f64,
17    /// Iterations actually taken.
18    pub iterations: usize,
19    /// True when `|f(x)| <= tol` was reached.
20    pub converged: bool,
21}
22
23/// A 1-D root finder: residual tolerance plus an iteration cap.
24#[derive(Debug, Clone, Copy)]
25pub struct Solver1d {
26    /// Convergence tolerance on `|f(x)|`.
27    pub tol: f64,
28    /// Maximum number of iterations.
29    pub max_iter: usize,
30}
31
32impl Default for Solver1d {
33    fn default() -> Self {
34        Self { tol: 1e-12, max_iter: 100 }
35    }
36}
37
38impl Solver1d {
39    pub fn new(tol: f64, max_iter: usize) -> Self {
40        Self { tol, max_iter }
41    }
42
43    // ── direct calls (method fixed at the call site) ────────────────────
44
45    /// [Bisection](bisection::bisection) on a sign-changing bracket.
46    pub fn bisection(&self, f: impl Fn(f64) -> f64, lo: f64, hi: f64) -> Result<Root, RustyQLibError> {
47        bisection::bisection(self, f, lo, hi)
48    }
49
50    /// [Newton-Raphson](newton_raphson::newton_raphson) with an analytic
51    /// derivative.
52    pub fn newton_raphson(
53        &self,
54        f: impl Fn(f64) -> f64,
55        df: impl Fn(f64) -> f64,
56        x0: f64,
57    ) -> Root {
58        newton_raphson::newton_raphson(self, f, df, x0)
59    }
60
61    /// [Secant](secant::secant) from two starting points.
62    pub fn secant(&self, f: impl Fn(f64) -> f64, x0: f64, x1: f64) -> Root {
63        secant::secant(self, f, x0, x1)
64    }
65
66    /// [Halley](halley::halley) with analytic first and second derivatives.
67    pub fn halley(
68        &self,
69        f: impl Fn(f64) -> f64,
70        df: impl Fn(f64) -> f64,
71        d2f: impl Fn(f64) -> f64,
72        x0: f64,
73    ) -> Root {
74        halley::halley(self, f, df, d2f, x0)
75    }
76
77    /// [Bracket-safeguarded Newton](newton_safeguarded::newton_safeguarded).
78    pub fn newton_safeguarded(
79        &self,
80        f: impl Fn(f64) -> f64,
81        df: impl Fn(f64) -> f64,
82        lo: f64,
83        hi: f64,
84        x0: f64,
85    ) -> Root {
86        newton_safeguarded::newton_safeguarded(self, f, df, lo, hi, x0)
87    }
88
89    // ── pluggable dispatch (method chosen at run time) ──────────────────
90
91    /// Solve `problem` with the chosen [`Method`].
92    ///
93    /// Derivatives the problem does not carry are approximated by central
94    /// finite differences, so Newton/Halley run on derivative-free
95    /// problems too. `Bisection` and `NewtonSafeguarded` require a
96    /// bracket and error without one.
97    pub fn solve(&self, method: Method, problem: &Problem) -> Result<Root, RustyQLibError> {
98        let f = |x: f64| (problem.f)(x);
99        let df = |x: f64| match problem.df {
100            Some(df) => df(x),
101            None => numeric_derivative(problem.f, x),
102        };
103        let d2f = |x: f64| match problem.d2f {
104            Some(d2f) => d2f(x),
105            None => numeric_second_derivative(problem.f, x),
106        };
107        let bracket = |name: &str| {
108            problem
109                .bracket
110                .ok_or_else(|| RustyQLibError::NumericalError(format!("{name} needs a bracket: use Problem::with_bracket")))
111        };
112        match method {
113            Method::Bisection => {
114                let (lo, hi) = bracket("bisection")?;
115                self.bisection(f, lo, hi)
116            }
117            Method::NewtonRaphson => Ok(self.newton_raphson(f, df, problem.x0)),
118            Method::Secant => {
119                // second start: the far bracket end when one is given,
120                // otherwise a small relative step from x0
121                let x1 = match problem.bracket {
122                    Some((lo, hi)) => {
123                        if (problem.x0 - lo).abs() > (problem.x0 - hi).abs() { lo } else { hi }
124                    }
125                    None => problem.x0 + 1e-4 * (1.0 + problem.x0.abs()),
126                };
127                Ok(self.secant(f, problem.x0, x1))
128            }
129            Method::Halley => Ok(self.halley(f, df, d2f, problem.x0)),
130            Method::NewtonSafeguarded => {
131                let (lo, hi) = bracket("newton_safeguarded")?;
132                Ok(self.newton_safeguarded(f, df, lo, hi, problem.x0))
133            }
134        }
135    }
136}
137
138/// The pluggable algorithm choice for [`Solver1d::solve`].
139#[derive(Debug, Clone, Copy, PartialEq, Eq)]
140pub enum Method {
141    Bisection,
142    NewtonRaphson,
143    Secant,
144    Halley,
145    NewtonSafeguarded,
146}
147
148/// A 1-D root-finding problem: the objective, whatever derivatives are
149/// available, a starting point, and an optional bracket.
150pub struct Problem<'a> {
151    pub f: &'a dyn Fn(f64) -> f64,
152    pub df: Option<&'a dyn Fn(f64) -> f64>,
153    pub d2f: Option<&'a dyn Fn(f64) -> f64>,
154    pub x0: f64,
155    pub bracket: Option<(f64, f64)>,
156}
157
158impl<'a> Problem<'a> {
159    pub fn new(f: &'a dyn Fn(f64) -> f64, x0: f64) -> Self {
160        Self { f, df: None, d2f: None, x0, bracket: None }
161    }
162
163    pub fn with_derivative(mut self, df: &'a dyn Fn(f64) -> f64) -> Self {
164        self.df = Some(df);
165        self
166    }
167
168    pub fn with_second_derivative(mut self, d2f: &'a dyn Fn(f64) -> f64) -> Self {
169        self.d2f = Some(d2f);
170        self
171    }
172
173    pub fn with_bracket(mut self, lo: f64, hi: f64) -> Self {
174        self.bracket = Some((lo.min(hi), lo.max(hi)));
175        self
176    }
177}
178
179fn numeric_derivative(f: &dyn Fn(f64) -> f64, x: f64) -> f64 {
180    let h = 1e-6 * (1.0 + x.abs());
181    (f(x + h) - f(x - h)) / (2.0 * h)
182}
183
184fn numeric_second_derivative(f: &dyn Fn(f64) -> f64, x: f64) -> f64 {
185    let h = 1e-4 * (1.0 + x.abs());
186    (f(x + h) - 2.0 * f(x) + f(x - h)) / (h * h)
187}
188
189#[cfg(test)]
190mod tests {
191    use super::*;
192
193    const SQRT2: f64 = std::f64::consts::SQRT_2;
194
195    fn f(x: f64) -> f64 {
196        x * x - 2.0
197    }
198    fn df(x: f64) -> f64 {
199        2.0 * x
200    }
201
202    #[test]
203    fn every_method_is_pluggable_on_one_problem() {
204        let obj = |x: f64| f(x);
205        let problem = Problem::new(&obj, 1.0).with_bracket(0.0, 2.0);
206        let solver = Solver1d::default();
207        for method in [
208            Method::Bisection,
209            Method::NewtonRaphson,
210            Method::Secant,
211            Method::Halley,
212            Method::NewtonSafeguarded,
213        ] {
214            let root = solver.solve(method, &problem).unwrap();
215            assert!(
216                root.converged && (root.x - SQRT2).abs() < 1e-9,
217                "{method:?}: {root:?}"
218            );
219        }
220    }
221
222    #[test]
223    fn numeric_derivative_fallback_matches_analytic() {
224        // same problem with and without the analytic derivative: Newton
225        // must land on the same root either way
226        let obj = |x: f64| f(x);
227        let d = |x: f64| df(x);
228        let with = Problem::new(&obj, 1.0).with_derivative(&d);
229        let without = Problem::new(&obj, 1.0);
230        let solver = Solver1d::default();
231        let ra = solver.solve(Method::NewtonRaphson, &with).unwrap();
232        let rn = solver.solve(Method::NewtonRaphson, &without).unwrap();
233        assert!(ra.converged && rn.converged);
234        assert!((ra.x - rn.x).abs() < 1e-9, "{} vs {}", ra.x, rn.x);
235    }
236
237    #[test]
238    fn bracketed_methods_error_without_a_bracket() {
239        let obj = |x: f64| f(x);
240        let problem = Problem::new(&obj, 1.0);
241        let solver = Solver1d::default();
242        assert!(solver.solve(Method::Bisection, &problem).is_err());
243        assert!(solver.solve(Method::NewtonSafeguarded, &problem).is_err());
244        // non-bracketed methods still work
245        assert!(solver.solve(Method::Secant, &problem).unwrap().converged);
246    }
247
248    #[test]
249    fn iteration_counts_are_reported() {
250        let r = Solver1d::default().newton_raphson(f, df, 1.0);
251        assert!(r.iterations > 0 && r.iterations < 10);
252    }
253}
rustyqlib/core/solvers/solver_1d.rs

rustyqlib/core/solvers/
solver_1d.rs
