commit 59c27b110d1f07a0b453ab5295368a847e3f714a Author: Jérôme Date: Sat Mar 25 22:43:27 2023 +0100 does not work yet. problems with closure as argument diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..47b5aa9 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +/target +*.lock \ No newline at end of file diff --git a/Cargo.toml b/Cargo.toml new file mode 100644 index 0000000..1787ee1 --- /dev/null +++ b/Cargo.toml @@ -0,0 +1,8 @@ +[package] +name = "rust_numerical_solvers" +version = "0.1.0" +edition = "2021" + +# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html + +[dependencies] diff --git a/src/main.rs b/src/main.rs new file mode 100644 index 0000000..c774bb1 --- /dev/null +++ b/src/main.rs @@ -0,0 +1,31 @@ +mod univariate_solvers; + +/// x sin(x) +/// e + ────── +/// x +fn fct(x : f64) -> f64 { + if x != 0.0 { + return f64::sin(x)/x + f64::exp(x); + } else { + return 2.0; + } +} + +fn dfct(x : f64) -> f64 { + if x != 0.0 { + return f64::exp(x) + f64::cos(x)/x - f64::sin(x)/x.powi(2); + } else { + return 1.0; + } +} + +fn main() { + println!("Testing Rust numerical solvers."); + let x0 : f64 = 1.0; + let tol : f64 = 1e-10; + let max_iter : u32 = 100; + let x_mathematica = -3.26650043678562449167148755288; + let x_newton = univariate_solvers::newton_solve(&(fct as fn(f64) -> f64), &(dfct as fn(f64) -> f64), x0, tol, max_iter); + println!("Mathematica : x = {}", x_mathematica); + println!("Newton's method : x = {}", x_newton); +} diff --git a/src/univariate_solvers.rs b/src/univariate_solvers.rs new file mode 100644 index 0000000..2db983d --- /dev/null +++ b/src/univariate_solvers.rs @@ -0,0 +1,39 @@ +/// @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 max_iter maximum number of iterations +/// @return solution +pub fn newton_solve f64, F2: Fn(f64) -> f64>(f : &F, df : &F2, x0 : f64, tol : f64, max_iter : u32) -> f64 { + let mut x: f64 = x0; + let mut dx: f64; + let mut fx: f64; + let mut dfx: f64; + for _iter in 0..max_iter { + fx = f(x); + dfx = df(x); + if dfx == 0.0 { + dx = fx; + } else { + dx = fx/dfx; + } + x -= dx; + if f64::abs(dx) < tol { + break; + } + } + return x; +} + +/// @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 max_iter maximum number of iterations +/// @return solution +pub fn newton_solve_num f64>(f : &F, x0 : f64, tol : f64, dx_num : f64, max_iter : u32) -> f64 { + return newton_solve(f, &((|x: f64| (f(x + dx_num) - f(x - dx_num))/(2.0*dx_num)) as fn(f64) -> f64), x0, tol, max_iter); +} +