jerome/rust-numerical-solvers
1
0
Fork 0
Code Issues 2 Pull requests Projects Releases Packages Wiki Activity

Added optimized numerical Halley's method.

Browse source
This commit is contained in:
Jérôme 2023-03-26 00:13:37 +01:00
parent 176168924a
commit 35ea6fe20b
2 changed files with 39 additions and 1 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
src/main.rs
View file

@ -39,6 +39,7 @@ fn main() {
let x_newton = univariate_solvers::newton_solve(&(fct as fn(f64) -> f64), &(dfct as fn(f64) -> f64), x0, tol, max_iter);
let x_newton_num: f64 = univariate_solvers::newton_solve_num(&(fct as fn(f64) -> f64), x0, tol, dx_num, max_iter);
let x_halley: f64 = univariate_solvers::halley_solve(&(fct as fn(f64) -> f64), &(dfct as fn(f64) -> f64), &(ddfct as fn(f64) -> f64), x0, tol, max_iter, false).unwrap();
let x_halley_num: f64 = univariate_solvers::halley_solve_num(&(fct as fn(f64) -> f64), x0, tol, dx_num, max_iter, false).unwrap();
let x_bisection : f64 = univariate_solvers::bisection_solve(&(fct as fn(f64) -> f64), -5.0, 1.0, tol).unwrap();
let x_secant : f64 = univariate_solvers::secant_solve(&(fct as fn(f64) -> f64), -1.0, 1.0, tol, max_iter);
let x_ridder : f64 = univariate_solvers::ridder_solve(&(fct as fn(f64) -> f64), -5.0, 1.0, tol, max_iter).unwrap();
@ -46,6 +47,7 @@ fn main() {
println!("Newton's method : x = {}\tf(x) = {}", x_newton, fct(x_newton));
println!("Newton's method (num) : x = {}\tf(x) = {}", x_newton_num, fct(x_newton_num));
println!("Halley's method : x = {}\tf(x) = {}", x_halley, fct(x_halley));
println!("Halley's method (num) : x = {}\tf(x) = {}", x_halley_num, fct(x_halley_num));
println!("Bisection : x = {}\tf(x) = {}", x_bisection, fct(x_bisection));
println!("Secant method : x = {}\tf(x) = {}", x_secant, fct(x_secant));
println!("Ridder's method : x = {}\tf(x) = {}", x_ridder, fct(x_ridder));

38
src/univariate_solvers.rs
View file

@ -30,11 +30,12 @@ pub fn newton_solve<F, F2>(f : F, df : F2, x0 : f64, tol : f64, max_iter : u32)
/// @brief Newton's method for solving a function f(x) = 0
/// @param f function to solve
/// @param df derivative of function f
/// @param x0 initial guess
/// @param tol tolerance
/// @param dx_num numerical differentiation step size
/// @param max_iter maximum number of iterations
/// @return solution
/// @note This method uses numerical differentiation to compute the first and second derivatives.
pub fn newton_solve_num<F>(f : F, x0 : f64, tol : f64, dx_num : f64, max_iter : u32) -> f64
where F : Fn(f64) -> f64
{
@ -74,6 +75,41 @@ where F : Fn(f64) -> f64, F2 : Fn(f64) -> f64, F3 : Fn(f64) -> f64
return Err("Halley method did not converge after reaching the maximum number of iterations allowed.")
}
/// Halley's method for solving a function f(x) = 0
/// @param f function to solve
/// @param x0 initial guess
/// @param tol tolerance
/// @param dx_num numerical differentiation step size
/// @param max_iter maximum number of iterations
/// @return solution
/// @note This method is more efficient than Newton's method, but requires the second derivative of f.
/// @note This method uses numerical differentiation to compute the first and second derivatives.
pub fn halley_solve_num<F>(f: F, x0: f64, tol: f64, dx_num : f64, max_iter: u32, verbose: bool) -> Result<f64, &'static str>
where F : Fn(f64) -> f64
{
let mut x: f64 = x0;
let mut f_x: f64;
let mut df_x: f64;
let mut ddf_x: f64;
let mut fx_m_dx: f64;
let mut fx_p_dx: f64;
for _i in 0..max_iter {
f_x = f(x);
fx_m_dx = f(x - dx_num);
fx_p_dx = f(x + dx_num);
df_x = (fx_p_dx - fx_m_dx)/(2.0*dx_num);
ddf_x = (fx_p_dx - 2.0*f_x + fx_m_dx) / (dx_num.powi(2));
if verbose {
println!("x = {}, f(x) = {}, df(x) = {}, ddf(x) = {}", x, f_x, df_x, ddf_x);
}
if f64::abs(f_x) < tol {
return Ok(x);
}
x = x - 2.0*f_x*df_x / (2.0*df_x.powi(2) - f_x*ddf_x);
}
return Err("Halley method did not converge after reaching the maximum number of iterations allowed.")
}
// --------------------------------------------------------------------
// ------------------------ Bracketing methods ------------------------
// --------------------------------------------------------------------