diff --git a/src/main.rs b/src/main.rs index 55ccc3f..f10bd10 100644 --- a/src/main.rs +++ b/src/main.rs @@ -29,8 +29,10 @@ 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_bisection : f64 = univariate_solvers::bisection_solve(&(fct as fn(f64) -> f64), -5.0, 1.0, tol).unwrap(); - println!("Mathematica : x = {}", x_mathematica); - println!("Newton's method : x = {}", x_newton); - println!("Newton's method (num) : x = {}", x_newton_num); - println!("Bisection : x = {}", x_bisection); + let x_secant : f64 = univariate_solvers::secant_solve(&(fct as fn(f64) -> f64), -1.0, 1.0, tol, max_iter); + println!("Mathematica : x = {}\tf(x) = {}", x_mathematica, fct(x_mathematica)); + 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!("Bisection : x = {}\tf(x) = {}", x_bisection, fct(x_bisection)); + println!("Secant method : x = {}\tf(x) = {}", x_secant, fct(x_secant)); } diff --git a/src/univariate_solvers.rs b/src/univariate_solvers.rs index 189eca1..91a97c1 100644 --- a/src/univariate_solvers.rs +++ b/src/univariate_solvers.rs @@ -47,6 +47,13 @@ where F : Fn(f64) -> f64 // ------------------------ Bracketing methods ------------------------ // -------------------------------------------------------------------- +/// @brief Bisection method for solving a function f(x) = 0 +/// @param f function to solve +/// @param a left bracket +/// @param b right bracket +/// @param tol tolerance +/// @return solution +/// @note The interval [a, b] must bracket the root, meaning f(a) and f(b) must be of a different sign. pub fn bisection_solve(f : F, mut a : f64, mut b : f64, tol : f64) -> Result where F : Fn(f64) -> f64 { @@ -69,4 +76,35 @@ where F : Fn(f64) -> f64 } } return Result::Ok((a + b)/2.0); +} + +/// @brief Secant method for solving a function f(x) = 0 +/// @param f function to solve +/// @param a left bracket +/// @param b right bracket +/// @param tol tolerance +/// @param max_iter maximum number of iterations +/// @return solution +/// @note The interval [a, b] does not have to bracket the root. +/// @note The secant method is not guaranteed to converge. +pub fn secant_solve(f : F, mut a : f64, mut b : f64, tol : f64, max_iter : u32) -> f64 +where F : Fn(f64) -> f64 +{ + let mut c: f64 = (a + b)/2.0; + let mut fa: f64 = f(a); + let mut fb: f64 = f(b); + let mut fc: f64; + for _iter in 0..max_iter { + // c is x[n], a is x[n-1], b is x[n-2] + c = a - fa*(a - b)/(fa - fb); + fc = f(c); + b = a; + fb = fa; + a = c; + fa = fc; + if (b - a).abs() < tol { + break; + } + } + return c; } \ No newline at end of file