From b4140bf2ec18f6743cb90cba7c67d26b03018757 Mon Sep 17 00:00:00 2001
From: "Thomas J. Kennedy" <tkennedy@cs.odu.edu>
Date: Sat, 20 Jun 2020 13:40:41 -0400
Subject: [PATCH] Port Python Non-Linear Solvers to Rust
---
.../{ => Python}/run_solvers.py | 0
.../{ => Python}/solvers.py | 0
NonLinearEquationSolvers/Rust/Cargo.lock | 149 +++++++++++++
NonLinearEquationSolvers/Rust/Cargo.toml | 10 +
NonLinearEquationSolvers/Rust/rustfmt.toml | 20 ++
NonLinearEquationSolvers/Rust/src/lib.rs | 5 +
.../Rust/src/loop_monoliths.rs | 196 ++++++++++++++++++
NonLinearEquationSolvers/Rust/src/main.rs | 159 ++++++++++++++
8 files changed, 539 insertions(+)
rename NonLinearEquationSolvers/{ => Python}/run_solvers.py (100%)
rename NonLinearEquationSolvers/{ => Python}/solvers.py (100%)
create mode 100644 NonLinearEquationSolvers/Rust/Cargo.lock
create mode 100644 NonLinearEquationSolvers/Rust/Cargo.toml
create mode 100644 NonLinearEquationSolvers/Rust/rustfmt.toml
create mode 100644 NonLinearEquationSolvers/Rust/src/lib.rs
create mode 100644 NonLinearEquationSolvers/Rust/src/loop_monoliths.rs
create mode 100644 NonLinearEquationSolvers/Rust/src/main.rs
diff --git a/NonLinearEquationSolvers/run_solvers.py b/NonLinearEquationSolvers/Python/run_solvers.py
similarity index 100%
rename from NonLinearEquationSolvers/run_solvers.py
rename to NonLinearEquationSolvers/Python/run_solvers.py
diff --git a/NonLinearEquationSolvers/solvers.py b/NonLinearEquationSolvers/Python/solvers.py
similarity index 100%
rename from NonLinearEquationSolvers/solvers.py
rename to NonLinearEquationSolvers/Python/solvers.py
diff --git a/NonLinearEquationSolvers/Rust/Cargo.lock b/NonLinearEquationSolvers/Rust/Cargo.lock
new file mode 100644
index 0000000..0305a80
--- /dev/null
+++ b/NonLinearEquationSolvers/Rust/Cargo.lock
@@ -0,0 +1,149 @@
+# This file is automatically @generated by Cargo.
+# It is not intended for manual editing.
+[[package]]
+name = "aho-corasick"
+version = "0.7.10"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "8716408b8bc624ed7f65d223ddb9ac2d044c0547b6fa4b0d554f3a9540496ada"
+dependencies = [
+ "memchr",
+]
+
+[[package]]
+name = "autocfg"
+version = "1.0.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "f8aac770f1885fd7e387acedd76065302551364496e46b3dd00860b2f8359b9d"
+
+[[package]]
+name = "hamcrest2"
+version = "0.3.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "49f837c62de05dc9cc71ff6486cd85de8856a330395ae338a04bfcefe5e91075"
+dependencies = [
+ "num",
+ "regex",
+]
+
+[[package]]
+name = "lazy_static"
+version = "1.4.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "e2abad23fbc42b3700f2f279844dc832adb2b2eb069b2df918f455c4e18cc646"
+
+[[package]]
+name = "memchr"
+version = "2.3.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "3728d817d99e5ac407411fa471ff9800a778d88a24685968b36824eaf4bee400"
+
+[[package]]
+name = "nonlinear_equation_solvers"
+version = "0.1.0"
+dependencies = [
+ "hamcrest2",
+]
+
+[[package]]
+name = "num"
+version = "0.2.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "b8536030f9fea7127f841b45bb6243b27255787fb4eb83958aa1ef9d2fdc0c36"
+dependencies = [
+ "num-bigint",
+ "num-complex",
+ "num-integer",
+ "num-iter",
+ "num-rational",
+ "num-traits",
+]
+
+[[package]]
+name = "num-bigint"
+version = "0.2.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "090c7f9998ee0ff65aa5b723e4009f7b217707f1fb5ea551329cc4d6231fb304"
+dependencies = [
+ "autocfg",
+ "num-integer",
+ "num-traits",
+]
+
+[[package]]
+name = "num-complex"
+version = "0.2.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "b6b19411a9719e753aff12e5187b74d60d3dc449ec3f4dc21e3989c3f554bc95"
+dependencies = [
+ "autocfg",
+ "num-traits",
+]
+
+[[package]]
+name = "num-integer"
+version = "0.1.43"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "8d59457e662d541ba17869cf51cf177c0b5f0cbf476c66bdc90bf1edac4f875b"
+dependencies = [
+ "autocfg",
+ "num-traits",
+]
+
+[[package]]
+name = "num-iter"
+version = "0.1.41"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "7a6e6b7c748f995c4c29c5f5ae0248536e04a5739927c74ec0fa564805094b9f"
+dependencies = [
+ "autocfg",
+ "num-integer",
+ "num-traits",
+]
+
+[[package]]
+name = "num-rational"
+version = "0.2.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "5c000134b5dbf44adc5cb772486d335293351644b801551abe8f75c84cfa4aef"
+dependencies = [
+ "autocfg",
+ "num-bigint",
+ "num-integer",
+ "num-traits",
+]
+
+[[package]]
+name = "num-traits"
+version = "0.2.12"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "ac267bcc07f48ee5f8935ab0d24f316fb722d7a1292e2913f0cc196b29ffd611"
+dependencies = [
+ "autocfg",
+]
+
+[[package]]
+name = "regex"
+version = "1.3.9"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "9c3780fcf44b193bc4d09f36d2a3c87b251da4a046c87795a0d35f4f927ad8e6"
+dependencies = [
+ "aho-corasick",
+ "memchr",
+ "regex-syntax",
+ "thread_local",
+]
+
+[[package]]
+name = "regex-syntax"
+version = "0.6.18"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "26412eb97c6b088a6997e05f69403a802a92d520de2f8e63c2b65f9e0f47c4e8"
+
+[[package]]
+name = "thread_local"
+version = "1.0.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "d40c6d1b69745a6ec6fb1ca717914848da4b44ae29d9b3080cbee91d72a69b14"
+dependencies = [
+ "lazy_static",
+]
diff --git a/NonLinearEquationSolvers/Rust/Cargo.toml b/NonLinearEquationSolvers/Rust/Cargo.toml
new file mode 100644
index 0000000..0452cdd
--- /dev/null
+++ b/NonLinearEquationSolvers/Rust/Cargo.toml
@@ -0,0 +1,10 @@
+[package]
+name = "nonlinear_equation_solvers"
+version = "0.1.0"
+authors = ["Thomas J. Kennedy <tkennedy@cs.odu.edu>"]
+edition = "2018"
+
+# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
+
+[dev-dependencies]
+hamcrest2 = "*"
diff --git a/NonLinearEquationSolvers/Rust/rustfmt.toml b/NonLinearEquationSolvers/Rust/rustfmt.toml
new file mode 100644
index 0000000..cab6f59
--- /dev/null
+++ b/NonLinearEquationSolvers/Rust/rustfmt.toml
@@ -0,0 +1,20 @@
+# Commented out options are only available in the 'nightly' channel
+
+# control_brace_style = "ClosingNextLine"
+# merge_imports = true
+newline_style = "Unix"
+# normalize_comments = true
+remove_nested_parens = true
+# reorder_impl_items = true
+reorder_imports = false
+# trailing_semicolon = true
+
+tab_spaces = 4
+hard_tabs = false
+max_width = 80
+# wrap_comments = true
+# report_todo = "Always"
+# report_fixme = "Always"
+# space_after_colon = true
+
+# blank_lines_lower_bound = 1
diff --git a/NonLinearEquationSolvers/Rust/src/lib.rs b/NonLinearEquationSolvers/Rust/src/lib.rs
new file mode 100644
index 0000000..027ffb5
--- /dev/null
+++ b/NonLinearEquationSolvers/Rust/src/lib.rs
@@ -0,0 +1,5 @@
+#[cfg(test)]
+#[macro_use]
+extern crate hamcrest2;
+
+pub mod loop_monoliths;
diff --git a/NonLinearEquationSolvers/Rust/src/loop_monoliths.rs b/NonLinearEquationSolvers/Rust/src/loop_monoliths.rs
new file mode 100644
index 0000000..9a39f01
--- /dev/null
+++ b/NonLinearEquationSolvers/Rust/src/loop_monoliths.rs
@@ -0,0 +1,196 @@
+const EPSILON: f64 = 10e-8;
+const MAX_ITERATIONS: u64 = 100;
+const ONE_HALF: f64 = 0.5;
+
+#[derive(Debug, Clone)]
+pub struct InvariantError {
+ pub msg: String,
+}
+
+impl InvariantError {
+ fn new(msg: String) -> Self {
+ InvariantError { msg }
+ }
+}
+
+/// Compute a solution to f using the bisection method
+///
+/// a0 = a
+/// b0 = b
+///
+/// for n = 1, 2, ... do
+/// xn = 0.5 (a_{n-1} + b_{n-1});
+/// if f(xn) < 0 then
+/// an = xn, bn = bn-1
+/// end if;
+///
+/// if f(xn) > 0 then
+/// an = an-1, bn = xn
+/// end if;
+///
+/// if bn - an < eps then
+/// break;
+/// end if;
+///
+/// end for
+pub fn bisection(
+ f: &impl Fn(f64) -> f64,
+ a: f64,
+ b: f64,
+) -> Result<(u64, f64), InvariantError> {
+ let mut a_n = a;
+ let mut b_n = b;
+
+ let mut x_n = 0_f64;
+
+ for n in 1..MAX_ITERATIONS {
+ if f(b_n) < 0.0 {
+ return Err(InvariantError::new(format!(
+ "$f(b_{} = {}) < 0$",
+ (n - 1),
+ b_n
+ )));
+ }
+
+ x_n = ONE_HALF * (a_n + b_n);
+
+ let f_of_x_n = f(x_n);
+
+ if f_of_x_n < 0.0 {
+ a_n = x_n;
+ // b_n = b_n-1 # unchanged
+ }
+
+ if f_of_x_n >= 0.0 {
+ // a_n = a_n-1 # unchanged
+ b_n = x_n;
+ }
+
+ // Stop Condition
+ if (b_n - a_n).abs() < EPSILON {
+ return Ok((n, x_n));
+ }
+ }
+
+ Ok((MAX_ITERATIONS, x_n))
+}
+
+/// Compute a solution to f using the false position method
+///
+/// a0 = a
+/// b0 = b;
+///
+/// for n = 1, 2, ... do
+/// xn = an - ((an - bn)/(f(an) - f(bn)))f(an);
+///
+/// if f(xn) · f(an) > 0 then
+/// an+1 = xn, bn+1 = bn
+/// else
+/// an+1 = an, bn+1 = xn
+/// end if
+/// end for
+pub fn regula_falsi(
+ f: &impl Fn(f64) -> f64,
+ a: f64,
+ b: f64,
+) -> Result<(u64, f64), InvariantError> {
+ let mut a_n = a;
+ let mut b_n = b;
+
+ let mut x_n = 0_f64;
+
+ for n in 1..MAX_ITERATIONS {
+ x_n = a_n - ((a_n - b_n) / (f(a_n) - f(b_n))) * f(a_n);
+
+ if !x_n.is_finite() {
+ return Err(InvariantError::new(
+ format!("$f(a_n) - f(b_n) == 0$",),
+ ));
+ }
+
+ if f(x_n) * f(a_n) > 0.0 {
+ a_n = x_n;
+ // b_n - No change
+ } else {
+ // a_n - No Change
+ b_n = x_n;
+ }
+
+ if (b_n - a_n).abs() < EPSILON {
+ return Ok((n, x_n));
+ }
+ }
+
+ Ok((MAX_ITERATIONS, x_n))
+}
+
+/// Compute a solution to f using the secant method
+///
+/// x_{n-1} = a0 = a
+/// x_n = b0 = b;
+///
+/// for n = 1, 2, ... do
+/// xn = xn - ((xn - x_n-1)/(f(xn) - f(x_n-1)))f(xn);
+///
+/// if f(xn) · f(an) > 0 then
+/// an+1 = xn, bn+1 = bn
+/// else
+/// an+1 = an, bn+1 = xn
+/// end if
+/// end for
+pub fn secant(
+ f: &impl Fn(f64) -> f64,
+ x_n_minus_1_in: f64,
+ x_n_in: f64,
+) -> Result<(u64, f64), InvariantError> {
+ let mut x_n_minus_1 = x_n_minus_1_in;
+ let mut x_n = x_n_in;
+
+ for n in 2..MAX_ITERATIONS {
+ let next_x_n =
+ x_n - ((x_n - x_n_minus_1) / (f(x_n) - f(x_n_minus_1))) * f(x_n);
+
+ if !next_x_n.is_finite() {
+ return Err(InvariantError::new(format!(
+ "$f(x_n) - f(x_nm1) == 0$",
+ )));
+ }
+
+ x_n_minus_1 = x_n;
+ x_n = next_x_n;
+
+ if (x_n - x_n_minus_1).abs() < EPSILON {
+ return Ok((n, x_n));
+ }
+ }
+
+ Ok((MAX_ITERATIONS, x_n))
+}
+
+/// Compute a solution to f using Newton's method.
+///
+/// Compute x_{n+1} = x_n - (f(x_n) / df(x_n)) until
+/// |x_{n+1} - x_n| <= eps
+pub fn newton(
+ f: &impl Fn(f64) -> f64,
+ df: &impl Fn(f64) -> f64,
+ x_0: f64,
+) -> Result<(u64, f64), InvariantError> {
+ let mut x_n = x_0;
+
+ for n in 1..MAX_ITERATIONS {
+ let next_x_n = x_n - (f(x_n) / df(x_n));
+
+ if !next_x_n.is_finite() {
+ return Err(InvariantError::new(format!("$df(x_n) == 0$")));
+ }
+
+ if (x_n - next_x_n).abs() < EPSILON {
+ return Ok((n, x_n));
+ }
+
+ x_n = next_x_n;
+ }
+
+ Ok((MAX_ITERATIONS, x_n))
+}
diff --git a/NonLinearEquationSolvers/Rust/src/main.rs b/NonLinearEquationSolvers/Rust/src/main.rs
new file mode 100644
index 0000000..f2c91c3
--- /dev/null
+++ b/NonLinearEquationSolvers/Rust/src/main.rs
@@ -0,0 +1,159 @@
+use std::env;
+use nonlinear_equation_solvers::loop_monoliths;
+
+fn main() {
+ let (a, b) = __handle_cli_args();
+
+ let (math_f, math_df) = __build_f_df();
+
+ println!("# Bisection");
+ match loop_monoliths::bisection(&math_f, a, b) {
+ Ok((_num_iter, solution)) => {
+ let fx = math_f(solution);
+
+ print_solution(solution, fx);
+ }
+ Err(err) => {
+ println!();
+ println!("## Method Failed");
+ println!("{:?}", err);
+ }
+ }
+
+ println!();
+ println!("# Regula Falsi (False Position)");
+ match loop_monoliths::regula_falsi(&math_f, a, b) {
+ Ok((_num_iter, solution)) => {
+ let fx = math_f(solution);
+
+ print_solution(solution, fx);
+ }
+ Err(err) => {
+ println!();
+ println!("## Method Failed");
+ println!("{:?}", err);
+ }
+ }
+
+ println!();
+ println!("# Secant");
+ match loop_monoliths::secant(&math_f, a, b) {
+ Ok((_num_iter, solution)) => {
+ let fx = math_f(solution);
+
+ print_solution(solution, fx);
+ }
+ Err(err) => {
+ println!();
+ println!("## Method Failed");
+ println!("{:?}", err);
+ }
+ }
+
+ println!();
+ println!("# Newton's Method");
+ match loop_monoliths::newton(&math_f, &math_df, a) {
+ Ok((_num_iter, solution)) => {
+ let fx = math_f(solution);
+
+ print_solution(solution, fx);
+ }
+ Err(err) => {
+ println!();
+ println!("## Method Failed");
+ println!("{:?}", err);
+ }
+ }
+}
+
+/// Handle positional command line argument validation.
+///
+/// # Return
+/// Tuple in the form (lower limit, upper limit)
+fn __handle_cli_args() -> (f64, f64) {
+ let args: Vec<String> = env::args().collect();
+
+ if args.len() < 3 {
+ eprintln!("Usage:\n {} a b", args[0]);
+ std::process::exit(1);
+ }
+
+ let limit_a = args[1]
+ .parse::<f64>()
+ .unwrap_or_else(|_err| std::process::exit(2));
+
+ let limit_b = args[2]
+ .parse::<f64>()
+ .unwrap_or_else(|_err| std::process::exit(3));
+
+ (limit_a, limit_b)
+}
+
+/// Print the solution (x) and f(x) under a subsection heading.
+///
+/// :param solution: approximate solution
+/// :param fx_solution: f(solution)
+fn print_solution(solution: f64, fx_solution: f64) {
+ println!();
+ println!("## Solution");
+ println!("$x={:20.8}", solution);
+ println!();
+ println!("$f(x)={:20.8}", fx_solution);
+ println!();
+}
+
+/// Wrapper function that returns two math style functions:
+///
+/// - f - a function in the form Real -> Real
+/// - df - the derivative of f in the form Real -> Real
+///
+/// A where clause does not work...
+///
+/// <https://stackoverflow.com/questions/47514930/what-are-the-differences-between-an-impl-trait-argument-and-generic-function-par>
+///
+/// Rust lambdas/closures syntax is tied with C++'s shenanigans for frustration
+///
+fn __build_f_df() -> (impl Fn(f64) -> f64, impl Fn(f64) -> f64) {
+ let f = |x: f64| -> f64 {
+ // (x ** 2) - 1
+ // Fraction(math.cos(x))
+ // Fraction(math.log(x)) // can not have negative operand (Real)
+ // x**5 - (7 * x)**2
+ x.powf(2.0) - (3.0 * x) - 4.0
+ };
+
+ let df = Box::new(|x: f64| -> f64 {
+ // 2 * x
+ // Fraction(-1 * math.sin(x))
+ // Fraction(numerator=1, denominator=x)
+ // 5 * (x ** 4) - (14 * x)
+ 2.0 * x - 3.0
+ });
+
+ (f, df)
+}
+
+#[cfg(test)]
+mod tests {
+ use super::*;
+
+ use hamcrest2::prelude::*;
+
+ #[test]
+ fn test_math_f() {
+ let (math_f, _) = __build_f_df();
+
+ assert_that!(math_f(0.0), close_to(-4.0, 1e-8));
+ assert_that!(math_f(1.0), close_to(-6.0, 1e-8));
+ assert_that!(math_f(2.0), close_to(-6.0, 1e-8));
+ }
+
+ #[test]
+ fn test_math_df() {
+ let (_, math_df) = __build_f_df();
+
+ assert_that!(math_df(0.0), close_to(-3.0, 1e-8));
+ assert_that!(math_df(1.0), close_to(-1.0, 1e-8));
+ assert_that!(math_df(2.0), close_to(1.0, 1e-8));
+ }
+}
--
GitLab