Skip to content

Commit

Permalink
Rename RootFinder
Browse files Browse the repository at this point in the history
  • Loading branch information
cpmech committed Jun 26, 2024
1 parent e2d9ce6 commit b8eb8a3
Show file tree
Hide file tree
Showing 4 changed files with 41 additions and 45 deletions.
2 changes: 1 addition & 1 deletion russell_lab/examples/algo_min_and_root_solver_brent.rs
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ fn main() -> Result<(), StrError> {
println!("\n{}", stats);

// root
let solver = RootFinding::new();
let solver = RootFinder::new();
let (xo, stats) = solver.brent(0.3, 0.4, args, |x, _| {
Ok(1.0 / (1.0 - f64::exp(-2.0 * x) * f64::powi(f64::sin(5.0 * PI * x), 2)) - 1.5)
})?;
Expand Down
2 changes: 1 addition & 1 deletion russell_lab/examples/algo_root_finding_chebyshev.rs
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ fn main() -> Result<(), StrError> {
println!("N = {}", nn);

// find all roots in the interval
let solver = RootFinding::new();
let solver = RootFinder::new();
let roots = Vector::from(&solver.chebyshev(&interp)?);
let f_at_roots = roots.get_mapped(|x| f(x, args).unwrap());
println!("roots =\n{}", roots);
Expand Down
68 changes: 32 additions & 36 deletions russell_lab/src/algo/root_finder.rs
Original file line number Diff line number Diff line change
Expand Up @@ -2,26 +2,8 @@ use crate::StrError;
use crate::{mat_eigenvalues, InterpChebyshev, TOL_RANGE};
use crate::{Matrix, Vector};

/// Implements a root finding solver using Chebyshev interpolation
///
/// This struct depends on [InterpChebyshev], an interpolant using
/// Chebyshev-Gauss-Lobatto points.
///
/// It is essential that the interpolant best approximates the
/// data/function; otherwise, not all roots can be found.
///
/// The roots are the eigenvalues of the companion matrix.
///
/// # References
///
/// 1. Boyd JP (2002) Computing zeros on a real interval through Chebyshev expansion
/// and polynomial rootfinding, SIAM Journal of Numerical Analysis, 40(5):1666-1682
/// 2. Boyd JP (2013) Finding the zeros of a univariate equation: proxy rootfinders,
/// Chebyshev interpolation, and the companion matrix, SIAM Journal of Numerical
/// Analysis, 55(2):375-396.
/// 3. Boyd JP (2014) Solving Transcendental Equations: The Chebyshev Polynomial Proxy
/// and Other Numerical Rootfinders, Perturbation Series, and Oracles, SIAM, pp460
pub struct RootFinding {
/// Implements root finding algorithms
pub struct RootFinder {
/// Holds the tolerance to avoid division by zero with the trailing Chebyshev coefficient
///
/// Default = 1e-13
Expand Down Expand Up @@ -75,10 +57,10 @@ pub struct RootFinding {
h_cen: f64,
}

impl RootFinding {
impl RootFinder {
/// Allocates a new instance
pub fn new() -> Self {
RootFinding {
RootFinder {
tol_zero_an: 1e-13,
tol_abs_imaginary: 1.0e-8,
tol_abs_boundary: TOL_RANGE / 10.0,
Expand Down Expand Up @@ -108,6 +90,20 @@ impl RootFinding {
/// 1. It is essential that the interpolant best approximates the data/function;
/// otherwise, not all roots can be found.
///
/// # Method
///
/// The roots are the eigenvalues of the companion matrix as explained in the references.
///
/// # References
///
/// 1. Boyd JP (2002) Computing zeros on a real interval through Chebyshev expansion
/// and polynomial rootfinding, SIAM Journal of Numerical Analysis, 40(5):1666-1682
/// 2. Boyd JP (2013) Finding the zeros of a univariate equation: proxy rootfinders,
/// Chebyshev interpolation, and the companion matrix, SIAM Journal of Numerical
/// Analysis, 55(2):375-396.
/// 3. Boyd JP (2014) Solving Transcendental Equations: The Chebyshev Polynomial Proxy
/// and Other Numerical Rootfinders, Perturbation Series, and Oracles, SIAM, pp460
///
/// # Example
///
/// ```
Expand All @@ -125,7 +121,7 @@ impl RootFinding {
/// interp.set_function(nn, args, f)?;
///
/// // find all roots in the interval
/// let mut solver = RootFinding::new();
/// let mut solver = RootFinder::new();
/// let roots = solver.chebyshev(&interp)?;
/// array_approx_eq(&roots, &[-1.0, 1.0], 1e-15);
/// Ok(())
Expand Down Expand Up @@ -216,7 +212,7 @@ impl RootFinding {
/// interp.set_function(nn, args, f)?;
///
/// // find all roots in the interval
/// let mut solver = RootFinding::new();
/// let mut solver = RootFinder::new();
/// let mut roots = solver.chebyshev(&interp)?;
/// array_approx_eq(&roots, &[-0.5, 0.5], 1e-15); // inaccurate
///
Expand Down Expand Up @@ -288,7 +284,7 @@ impl RootFinding {

#[cfg(test)]
mod tests {
use super::RootFinding;
use super::RootFinder;
use crate::algo::NoArgs;
use crate::InterpChebyshev;
use crate::{approx_eq, array_approx_eq, get_test_functions};
Expand Down Expand Up @@ -368,7 +364,7 @@ mod tests {
let (xa, xb) = (-4.0, 4.0);
let nn = 2;
let interp = InterpChebyshev::new(nn, xa, xb).unwrap();
let solver = RootFinding::new();
let solver = RootFinder::new();
assert_eq!(
solver.chebyshev(&interp).err(),
Some("the interpolant must initialized first")
Expand All @@ -383,7 +379,7 @@ mod tests {
let args = &mut 0;
let mut interp = InterpChebyshev::new(nn, xa, xb).unwrap();
interp.set_function(nn, args, f).unwrap();
let solver = RootFinding::new();
let solver = RootFinder::new();
assert_eq!(
solver.chebyshev(&interp).err(),
Some("the trailing Chebyshev coefficient vanishes; try a smaller degree N")
Expand All @@ -403,7 +399,7 @@ mod tests {
interp.set_function(nn, args, f).unwrap();

// find roots
let solver = RootFinding::new();
let solver = RootFinder::new();
let roots = solver.chebyshev(&interp).unwrap();
let mut roots_refined = roots.clone();
solver.refine(&mut roots_refined, xa, xb, args, f).unwrap();
Expand Down Expand Up @@ -444,7 +440,7 @@ mod tests {
interp.set_function(nn, args, f).unwrap();

// find roots
let solver = RootFinding::new();
let solver = RootFinder::new();
let roots = solver.chebyshev(&interp).unwrap();
let mut roots_refined = roots.clone();
solver.refine(&mut roots_refined, xa, xb, args, f).unwrap();
Expand Down Expand Up @@ -483,7 +479,7 @@ mod tests {
interp.set_function(nn, args, f).unwrap();

// find roots
let solver = RootFinding::new();
let solver = RootFinder::new();
let roots = solver.chebyshev(&interp).unwrap();
let mut roots_refined = roots.clone();
solver.refine(&mut roots_refined, xa, xb, args, f).unwrap();
Expand Down Expand Up @@ -515,7 +511,7 @@ mod tests {
let args = &mut 0;
let _ = f(0.0, args);
let (xa, xb) = (-1.0, 1.0);
let mut solver = RootFinding::new();
let mut solver = RootFinder::new();
let mut roots = Vec::new();
assert_eq!(
solver.refine(&mut roots, xa, xb, args, f).err(),
Expand All @@ -542,7 +538,7 @@ mod tests {
interp.set_function(nn, args, f).unwrap();

// find roots
let solver = RootFinding::new();
let solver = RootFinder::new();
let roots = solver.chebyshev(&interp).unwrap();
let mut roots_refined = roots.clone();
solver.refine(&mut roots_refined, xa, xb, args, f).unwrap();
Expand Down Expand Up @@ -582,7 +578,7 @@ mod tests {
let (xa, xb) = test.range;
let mut interp = InterpChebyshev::new(nn_max, xa, xb).unwrap();
interp.adapt_function(tol, args, test.f).unwrap();
let solver = RootFinding::new();
let solver = RootFinder::new();
let roots = solver.chebyshev(&interp).unwrap();
let mut roots_refined = roots.clone();
if roots.len() > 0 {
Expand Down Expand Up @@ -638,7 +634,7 @@ mod tests {
interp.set_data(uu).unwrap();

// find all roots in the interval
let solver = RootFinding::new();
let solver = RootFinder::new();
let roots = &solver.chebyshev(&interp).unwrap();
let nroot = roots.len();
assert_eq!(nroot, 0)
Expand All @@ -657,7 +653,7 @@ mod tests {
interp.adapt_data(tol, uu).unwrap();

// find all roots in the interval
let solver = RootFinding::new();
let solver = RootFinder::new();
let roots = &solver.chebyshev(&interp).unwrap();
let nroot = roots.len();
assert_eq!(nroot, 0)
Expand All @@ -684,7 +680,7 @@ mod tests {
interp.adapt_data(tol, uu).unwrap();

// find all roots in the interval
let solver = RootFinding::new();
let solver = RootFinder::new();
let roots = solver.chebyshev(&interp).unwrap();
let nroot = roots.len();
assert_eq!(nroot, 1);
Expand Down
14 changes: 7 additions & 7 deletions russell_lab/src/algo/root_finder_brent.rs
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
use super::{RootFinding, Stats};
use super::{RootFinder, Stats};
use crate::StrError;

impl RootFinding {
impl RootFinder {
/// Employs Brent's method to find a single root of an equation
///
/// See: <https://mathworld.wolfram.com/BrentsMethod.html>
Expand Down Expand Up @@ -31,7 +31,7 @@ impl RootFinding {
///
/// fn main() -> Result<(), StrError> {
/// let args = &mut 0;
/// let solver = RootFinding::new();
/// let solver = RootFinder::new();
/// let (xa, xb) = (-4.0, 0.0);
/// let (xo, stats) = solver.brent(xa, xb, args, |x, _| Ok(4.0 - x * x))?;
/// println!("\nroot = {:?}", xo);
Expand Down Expand Up @@ -196,13 +196,13 @@ impl RootFinding {
mod tests {
use crate::algo::testing::get_test_functions;
use crate::approx_eq;
use crate::{NoArgs, RootFinding};
use crate::{NoArgs, RootFinder};

#[test]
fn brent_find_captures_errors_1() {
let f = |x, _: &mut NoArgs| Ok(x * x - 1.0);
let args = &mut 0;
let mut solver = RootFinding::new();
let mut solver = RootFinder::new();
assert_eq!(f(1.0, args).unwrap(), 0.0);
assert_eq!(
solver.brent(-0.5, -0.5, args, f).err(),
Expand Down Expand Up @@ -235,7 +235,7 @@ mod tests {
res
};
let args = &mut Args { count: 0, target: 0 };
let solver = RootFinding::new();
let solver = RootFinder::new();
// first function call
assert_eq!(solver.brent(-0.5, 2.0, args, f).err(), Some("stop"));
// second function call
Expand All @@ -251,7 +251,7 @@ mod tests {
#[test]
fn brent_find_works() {
let args = &mut 0;
let solver = RootFinding::new();
let solver = RootFinder::new();
for test in &get_test_functions() {
println!("\n===================================================================");
println!("\n{}", test.name);
Expand Down

0 comments on commit b8eb8a3

Please sign in to comment.