bacon_sci/roots/
polynomial.rs

1use crate::polynomial::Polynomial;
2use nalgebra::ComplexField;
3use num_complex::Complex;
4use num_traits::{FromPrimitive, Zero};
5
6/// Use Newton's method on a polynomial.
7///
8/// Given an initial guess of a polynomial root, use Newton's method to
9/// approximate within tol.
10///
11/// # Returns
12/// `Ok(root)` when a root has been found, `Err` on failure
13///
14/// # Params
15/// `initial` initial estimate of the root
16///
17/// `poly` the polynomial to solve the root for
18///
19/// `tol` The tolerance of relative error between iterations
20///
21/// `n_max` the maximum number of iterations
22///
23/// # Examples
24/// ```
25/// use bacon_sci::polynomial::Polynomial;
26/// use bacon_sci::roots::newton_polynomial;
27/// fn example() {
28///   let mut polynomial = Polynomial::new();
29///   polynomial.set_coefficient(2, 1.0);
30///   polynomial.set_coefficient(0, -1.0);
31///   let solution = newton_polynomial(0.5, &polynomial, 0.0001, 1000).unwrap();
32/// }
33/// ```
34pub fn newton_polynomial<N: ComplexField + FromPrimitive + Copy>(
35    initial: N,
36    poly: &Polynomial<N>,
37    tol: <N as ComplexField>::RealField,
38    n_max: usize,
39) -> Result<N, String>
40where
41    <N as ComplexField>::RealField: FromPrimitive + Copy,
42{
43    let mut n = 0;
44
45    let mut guess = initial;
46
47    let mut norm = guess.abs();
48    if norm <= tol {
49        return Ok(guess);
50    }
51
52    while n < n_max {
53        let (f_val, f_deriv_val) = poly.evaluate_derivative(guess);
54        let new_guess = guess - (f_val / f_deriv_val);
55        let new_norm = new_guess.abs();
56        if ((norm - new_norm) / norm).abs() <= tol || new_norm <= tol {
57            return Ok(new_guess);
58        }
59
60        norm = new_norm;
61        guess = new_guess;
62        n += 1;
63    }
64
65    Err("Newton_polynomial: maximum iterations exceeded".to_owned())
66}
67
68/// Use Mueller's method on a polynomial. Note this usually requires complex numbers.
69///
70/// Givin three initial guesses of a polynomial root, use Muller's method to
71/// approximate within tol.
72///
73/// # Returns
74/// `Ok(root)` when a root has been found, `Err` on failure
75///
76/// # Params
77/// `initial` Triplet of initial guesses
78///
79/// `poly` the polynomial to solve the root for
80///
81/// `tol` the tolerance of relative error between iterations
82///
83/// `n_max` Maximum number of iterations
84/// # Examples
85/// ```
86/// use bacon_sci::polynomial::Polynomial;
87/// use bacon_sci::roots::muller_polynomial;
88/// fn example() {
89///   let mut polynomial = Polynomial::new();
90///   polynomial.set_coefficient(2, 1.0);
91///   polynomial.set_coefficient(0, -1.0);
92///   let solution = muller_polynomial((0.0, 1.5, 2.0), &polynomial, 0.0001, 1000).unwrap();
93/// }
94/// ```
95pub fn muller_polynomial<N: ComplexField + FromPrimitive + Copy>(
96    initial: (N, N, N),
97    poly: &Polynomial<N>,
98    tol: <N as ComplexField>::RealField,
99    n_max: usize,
100) -> Result<Complex<<N as ComplexField>::RealField>, String>
101where
102    <N as ComplexField>::RealField: FromPrimitive + Copy,
103{
104    let poly = poly.make_complex();
105    let mut n = 0;
106    let mut poly_0 = Complex::<N::RealField>::new(initial.0.real(), initial.0.imaginary());
107    let mut poly_1 = Complex::<N::RealField>::new(initial.1.real(), initial.1.imaginary());
108    let mut poly_2 = Complex::<N::RealField>::new(initial.2.real(), initial.1.imaginary());
109    let mut h_1 = poly_1 - poly_0;
110    let mut h_2 = poly_2 - poly_1;
111    let poly_1_evaluated = poly.evaluate(poly_1);
112    let mut poly_2_evaluated = poly.evaluate(poly_2);
113    let mut delta_1 = (poly_1_evaluated - poly.evaluate(poly_0)) / h_1;
114    let mut delta_2 = (poly_2_evaluated - poly_1_evaluated) / h_2;
115    let mut delta = (delta_2 - delta_1) / (h_2 + h_1);
116
117    let negtwo = N::RealField::from_i32(-2).unwrap();
118    let four = N::RealField::from_i32(4).unwrap();
119
120    while n < n_max {
121        let b_coefficient = delta_2 + h_2 * delta;
122        let determinate = (b_coefficient.powi(2)
123            - Complex::<N::RealField>::new(four, N::RealField::zero()) * poly_2_evaluated * delta)
124            .sqrt();
125        let error = if (b_coefficient - determinate).abs() < (b_coefficient + determinate).abs() {
126            b_coefficient + determinate
127        } else {
128            b_coefficient - determinate
129        };
130        let step =
131            Complex::<N::RealField>::new(negtwo, N::RealField::zero()) * poly_2_evaluated / error;
132        let p = poly_2 + step;
133
134        if step.abs() <= tol {
135            return Ok(p);
136        }
137
138        poly_0 = poly_1;
139        poly_1 = poly_2;
140        poly_2 = p;
141        poly_2_evaluated = poly.evaluate(p);
142        h_1 = poly_1 - poly_0;
143        h_2 = poly_2 - poly_1;
144        let poly_1_evaluated = poly.evaluate(poly_1);
145        delta_1 = (poly_1_evaluated - poly.evaluate(poly_0)) / h_1;
146        delta_2 = (poly_2_evaluated - poly_1_evaluated) / h_2;
147        delta = (delta_2 - delta_1) / (h_1 + h_2);
148
149        n += 1;
150    }
151
152    Err("Muller: maximum iterations exceeded".to_owned())
153}
bacon_sci/roots/polynomial.rs

bacon_sci/roots/
polynomial.rs
