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

bacon_sci_1/roots/
polynomial.rs
