mathhook_core/algebra/root_finding/
newton_raphson.rs

1//! Newton-Raphson method for root finding
2//!
3//! Implements Newton-Raphson method using numerical differentiation.
4//! Provides quadratic convergence near simple roots.
5//!
6//! # Algorithm
7//!
8//! Given initial guess x0:
9//! x(n+1) = x(n) - f(x(n)) / f'(x(n))
10//!
11//! # Convergence
12//!
13//! - Quadratic convergence near simple roots (error squares each iteration)
14//! - Requires good initial guess
15//! - May fail for multiple roots or if derivative is zero
16//! - Uses finite differences for derivative: f'(x) ≈ (f(x+h) - f(x))/h
17
18use super::{RootFinder, RootFindingConfig, RootResult};
19use crate::error::MathError;
20
21/// Newton-Raphson method with numerical differentiation
22pub struct NewtonRaphson {
23    /// Initial guess for the root
24    pub initial_guess: f64,
25}
26
27impl NewtonRaphson {
28    /// Create a new Newton-Raphson method with initial guess
29    ///
30    /// # Arguments
31    ///
32    /// * `initial_guess` - Starting point for iteration
33    ///
34    /// # Examples
35    ///
36    /// ```rust
37    /// use mathhook_core::algebra::root_finding::NewtonRaphson;
38    ///
39    /// let method = NewtonRaphson::new(1.5);
40    /// ```
41    #[inline]
42    pub fn new(initial_guess: f64) -> Self {
43        Self { initial_guess }
44    }
45
46    /// Compute numerical derivative using central difference
47    #[inline]
48    fn numerical_derivative<F>(&self, f: &F, x: f64, h: f64) -> f64
49    where
50        F: Fn(f64) -> f64,
51    {
52        (f(x + h) - f(x - h)) / (2.0 * h)
53    }
54
55    /// Check for convergence stagnation
56    #[inline]
57    fn is_stagnant(&self, x_new: f64, x_old: f64, tolerance: f64) -> bool {
58        (x_new - x_old).abs() < tolerance * x_old.abs().max(1.0)
59    }
60}
61
62impl RootFinder for NewtonRaphson {
63    fn find_root<F>(&self, f: F, config: &RootFindingConfig) -> Result<RootResult, MathError>
64    where
65        F: Fn(f64) -> f64,
66    {
67        let mut x = self.initial_guess;
68
69        for iteration in 0..config.max_iterations {
70            let fx = f(x);
71
72            if fx.abs() < config.tolerance {
73                return Ok(RootResult {
74                    root: x,
75                    iterations: iteration + 1,
76                    function_value: fx,
77                    converged: true,
78                });
79            }
80
81            let fpx = self.numerical_derivative(&f, x, config.derivative_h);
82
83            if fpx.abs() < 1e-14 {
84                return Err(MathError::ConvergenceFailed {
85                    reason: format!("Derivative too small at x = {}: f'(x) = {}", x, fpx),
86                });
87            }
88
89            let x_new = x - fx / fpx;
90
91            if x_new.is_nan() || x_new.is_infinite() {
92                return Err(MathError::ConvergenceFailed {
93                    reason: format!("Iteration produced invalid value: {}", x_new),
94                });
95            }
96
97            if self.is_stagnant(x_new, x, config.tolerance) {
98                let fx_new = f(x_new);
99                if fx_new.abs() < config.tolerance {
100                    return Ok(RootResult {
101                        root: x_new,
102                        iterations: iteration + 1,
103                        function_value: fx_new,
104                        converged: true,
105                    });
106                } else {
107                    return Err(MathError::ConvergenceFailed {
108                        reason: "Iteration stagnated without converging".to_owned(),
109                    });
110                }
111            }
112
113            x = x_new;
114        }
115
116        Err(MathError::MaxIterationsReached {
117            max_iterations: config.max_iterations,
118        })
119    }
120}
121
122#[cfg(test)]
123mod tests {
124    use super::*;
125
126    #[test]
127    fn test_newton_simple_linear() {
128        let method = NewtonRaphson::new(0.5);
129        let config = RootFindingConfig::default();
130
131        let result = method.find_root(|x| x - 1.0, &config).unwrap();
132
133        assert!((result.root - 1.0).abs() < 1e-9);
134        assert!(result.converged);
135    }
136
137    #[test]
138    fn test_newton_quadratic() {
139        let method = NewtonRaphson::new(1.5);
140        let config = RootFindingConfig::default();
141
142        let result = method.find_root(|x| x * x - 2.0, &config).unwrap();
143
144        assert!((result.root - 2.0_f64.sqrt()).abs() < 1e-10);
145        assert!(result.converged);
146        assert!(result.iterations < 10);
147    }
148
149    #[test]
150    fn test_newton_cubic() {
151        let method = NewtonRaphson::new(0.5);
152        let config = RootFindingConfig::default();
153
154        let result = method
155            .find_root(|x| x * x * x + x * x - 1.0, &config)
156            .unwrap();
157
158        let f_at_root = result.root.powi(3) + result.root.powi(2) - 1.0;
159        assert!(f_at_root.abs() < 1e-9);
160        assert!(result.converged);
161    }
162
163    #[test]
164    fn test_newton_transcendental() {
165        let method = NewtonRaphson::new(0.5);
166        let config = RootFindingConfig::default();
167
168        let result = method.find_root(|x| x.cos() - x, &config).unwrap();
169
170        let expected = 0.7390851332151607;
171        assert!((result.root - expected).abs() < 1e-10);
172        assert!(result.converged);
173    }
174
175    #[test]
176    fn test_newton_exponential() {
177        let method = NewtonRaphson::new(0.5);
178        let config = RootFindingConfig::default();
179
180        let result = method.find_root(|x| x.exp() - 2.0, &config).unwrap();
181
182        assert!((result.root - 2.0_f64.ln()).abs() < 1e-12);
183        assert!(result.converged);
184    }
185
186    #[test]
187    fn test_newton_sine() {
188        let method = NewtonRaphson::new(3.0);
189        let config = RootFindingConfig::default();
190
191        let result = method.find_root(|x| x.sin(), &config).unwrap();
192
193        assert!((result.root - std::f64::consts::PI).abs() < 1e-12);
194        assert!(result.converged);
195    }
196
197    #[test]
198    fn test_newton_polynomial() {
199        let method = NewtonRaphson::new(1.5);
200        let config = RootFindingConfig::default();
201
202        let result = method
203            .find_root(|x| x * x * x - 3.0 * x + 1.0, &config)
204            .unwrap();
205
206        let f_at_root = result.root.powi(3) - 3.0 * result.root + 1.0;
207        assert!(f_at_root.abs() < 1e-9);
208        assert!(result.converged);
209    }
210
211    #[test]
212    fn test_newton_high_degree_polynomial() {
213        let method = NewtonRaphson::new(1.2);
214        let config = RootFindingConfig::default();
215
216        let result = method.find_root(|x| x.powi(5) - x - 1.0, &config).unwrap();
217
218        let f_at_root = result.root.powi(5) - result.root - 1.0;
219        assert!(f_at_root.abs() < 1e-9);
220        assert!(result.converged);
221    }
222
223    #[test]
224    fn test_newton_zero_derivative_fails() {
225        let method = NewtonRaphson::new(0.0);
226        let config = RootFindingConfig::default();
227
228        let result = method.find_root(|x| x * x + 1.0, &config);
229
230        assert!(result.is_err());
231    }
232
233    #[test]
234    fn test_newton_fast_convergence() {
235        let method = NewtonRaphson::new(1.5);
236        let config = RootFindingConfig {
237            tolerance: 1e-14,
238            ..Default::default()
239        };
240
241        let result = method.find_root(|x| x * x - 2.0, &config).unwrap();
242
243        assert!(result.iterations < 6);
244        assert!((result.root - 2.0_f64.sqrt()).abs() < 1e-14);
245    }
246
247    #[test]
248    fn test_newton_different_initial_guesses() {
249        let config = RootFindingConfig::default();
250
251        let guesses = vec![0.5, 1.0, 1.5, 2.0];
252        for guess in guesses {
253            let method = NewtonRaphson::new(guess);
254            let result = method.find_root(|x| x * x - 2.0, &config).unwrap();
255
256            assert!((result.root - 2.0_f64.sqrt()).abs() < 1e-10);
257        }
258    }
259
260    #[test]
261    fn test_newton_negative_root() {
262        let method = NewtonRaphson::new(-1.5);
263        let config = RootFindingConfig::default();
264
265        let result = method.find_root(|x| x * x - 2.0, &config).unwrap();
266
267        assert!((result.root + 2.0_f64.sqrt()).abs() < 1e-10);
268        assert!(result.converged);
269    }
270
271    #[test]
272    fn test_newton_logarithmic() {
273        let method = NewtonRaphson::new(2.0);
274        let config = RootFindingConfig::default();
275
276        let result = method.find_root(|x| x.ln() - 1.0, &config).unwrap();
277
278        assert!((result.root - std::f64::consts::E).abs() < 1e-10);
279        assert!(result.converged);
280    }
281
282    #[test]
283    fn test_newton_trigonometric_combination() {
284        let method = NewtonRaphson::new(1.0);
285        let config = RootFindingConfig::default();
286
287        let result = method
288            .find_root(|x| x.sin() + x.cos() - 1.0, &config)
289            .unwrap();
290
291        let f_at_root = result.root.sin() + result.root.cos() - 1.0;
292        assert!(f_at_root.abs() < 1e-9);
293        assert!(result.converged);
294    }
295
296    #[test]
297    fn test_newton_rational_function() {
298        let method = NewtonRaphson::new(2.0);
299        let config = RootFindingConfig::default();
300
301        let result = method
302            .find_root(|x| (x * x - 4.0) / (x + 1.0), &config)
303            .unwrap();
304
305        assert!((result.root - 2.0).abs() < 1e-10 || (result.root + 2.0).abs() < 1e-10);
306        assert!(result.converged);
307    }
308
309    #[test]
310    fn test_newton_composite_function() {
311        let method = NewtonRaphson::new(1.5);
312        let config = RootFindingConfig::default();
313
314        let result = method.find_root(|x| (x * x).sin() - 0.5, &config).unwrap();
315
316        let f_at_root = (result.root * result.root).sin() - 0.5;
317        assert!(f_at_root.abs() < 1e-9);
318        assert!(result.converged);
319    }
320
321    #[test]
322    fn test_newton_near_zero() {
323        let method = NewtonRaphson::new(0.1);
324        let config = RootFindingConfig::default();
325
326        let result = method.find_root(|x| x * x * x, &config).unwrap();
327
328        assert!(result.root.abs() < 1e-3);
329        assert!(result.converged);
330    }
331
332    #[test]
333    fn test_newton_hyperbolic() {
334        let method = NewtonRaphson::new(1.0);
335        let config = RootFindingConfig::default();
336
337        let result = method.find_root(|x| x.sinh() - 1.0, &config).unwrap();
338
339        assert!((result.root - 1.0_f64.asinh()).abs() < 1e-10);
340        assert!(result.converged);
341    }
342
343    #[test]
344    fn test_newton_inverse_function() {
345        let method = NewtonRaphson::new(0.5);
346        let config = RootFindingConfig::default();
347
348        let result = method.find_root(|x| 1.0 / x - 2.0, &config).unwrap();
349
350        assert!((result.root - 0.5).abs() < 1e-10);
351        assert!(result.converged);
352    }
353
354    #[test]
355    fn test_newton_tolerance_control() {
356        let method = NewtonRaphson::new(1.5);
357
358        let config_loose = RootFindingConfig {
359            tolerance: 1e-6,
360            ..Default::default()
361        };
362        let result_loose = method.find_root(|x| x * x - 2.0, &config_loose).unwrap();
363
364        let config_tight = RootFindingConfig {
365            tolerance: 1e-14,
366            ..Default::default()
367        };
368        let result_tight = method.find_root(|x| x * x - 2.0, &config_tight).unwrap();
369
370        assert!(result_tight.function_value.abs() < result_loose.function_value.abs());
371    }
372}
mathhook_core/algebra/root_finding/newton_raphson.rs

mathhook_core/algebra/root_finding/
newton_raphson.rs
