Skip to content
Snippets Groups Projects
Commit 2864548d authored by Thomas Kennedy's avatar Thomas Kennedy
Browse files

Add iterator version of Newton's Solver

parent 6bae3bdd
No related branches found
No related tags found
No related merge requests found
use super::errors::InvariantError;
trait IterableStateAdapter {
fn next_state(&mut self) -> Result<f64, InvariantError>;
fn get_current_guess(&self) -> Option<f64>;
}
pub struct NewtonSolver<'a> {
pub x_n: Option<f64>,
f: Box<dyn 'a + Fn(f64) -> f64>,
df: Box<dyn 'a + Fn(f64) -> f64>,
}
impl<'a> NewtonSolver<'a> {
pub fn new(
x_0: f64,
math_f: &'a impl Fn(f64) -> f64,
math_df: &'a impl Fn(f64) -> f64,
) -> Self {
NewtonSolver {
x_n: Some(x_0),
f: Box::new(math_f),
df: Box::new(math_df),
}
}
pub fn iter(&mut self) -> &mut Self {
self
}
}
impl<'a> IterableStateAdapter for NewtonSolver<'a> {
fn next_state(&mut self) -> Result<f64, InvariantError> {
let x_n = self.x_n.unwrap();
let f = &self.f;
let df = &self.df;
let next_x_n = x_n - (f(x_n) / df(x_n));
match next_x_n.is_finite() {
false => Err(InvariantError::from("$df(x_n) == 0$")),
true => {
self.x_n = Some(next_x_n);
Ok(next_x_n)
}
}
}
fn get_current_guess(&self) -> Option<f64> {
self.x_n
}
}
impl<'a> std::iter::Iterator for NewtonSolver<'a> {
type Item = f64;
fn next(&mut self) -> Option<Self::Item> {
let current_guess = self.get_current_guess();
if let Some(x_n) = current_guess {
match self.next_state() {
Ok(x_n) => self.x_n = Some(x_n),
Err(_) => self.x_n = None,
}
return Some(x_n);
}
None
}
}
#[cfg(test)]
mod tests {
use super::*;
use hamcrest2::prelude::*;
#[test]
fn test_newton_solver_some() {
let f = |x: f64| -> f64 { x.powf(2.0) };
let df = |x: f64| -> f64 { 2.0 * x };
let mut solver = NewtonSolver::new(-1.0, &f, &df);
assert_eq!(solver.next(), Some(-1.0));
assert_eq!(solver.next(), Some(-0.5));
assert_eq!(solver.next(), Some(-0.25));
assert_eq!(solver.next(), Some(-0.125));
assert_eq!(solver.next(), Some(-0.0625));
}
#[test]
fn test_newton_solver_none() {
let f = |x: f64| -> f64 { x.powf(2.0) };
let df = |x: f64| -> f64 { 2.0 * x };
let mut solver = NewtonSolver::new(0.0, &f, &df);
assert_eq!(solver.next(), Some(0.0));
assert_eq!(solver.next(), None);
}
}
......@@ -3,4 +3,5 @@
extern crate hamcrest2;
pub mod errors;
pub mod iterators;
pub mod loop_monoliths;
......@@ -5,7 +5,7 @@ const MAX_ITERATIONS: u64 = 100;
const ONE_HALF: f64 = 0.5;
/// Compute a solution to f using the bisection method
///
/// ```ignore
/// a0 = a
/// b0 = b
///
......@@ -64,7 +64,7 @@ pub fn bisection(
}
/// Compute a solution to f using the false position method
///
/// ```ignore
/// a0 = a
/// b0 = b;
///
......@@ -111,7 +111,7 @@ pub fn regula_falsi(
}
/// Compute a solution to f using the secant method
///
/// ```ignore
/// x_{n-1} = a0 = a
/// x_n = b0 = b;
///
......@@ -152,7 +152,7 @@ pub fn secant(
}
/// Compute a solution to f using Newton's method.
///
/// ```ignore
/// Compute x_{n+1} = x_n - (f(x_n) / df(x_n)) until
/// |x_{n+1} - x_n| <= eps
pub fn newton(
......
use std::env;
use nonlinear_equation_solvers::loop_monoliths;
use nonlinear_equation_solvers::iterators;
/// 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)
}
/// 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)
}
fn main() {
let (a, b) = __handle_cli_args();
......@@ -64,73 +132,41 @@ fn main() {
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);
//--------------------------------------------------------------------------
iter_example_1(a, &math_f, &math_df);
println!();
iter_example_2(a, &math_f, &math_df);
}
/// 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
};
fn iter_example_1(
a: f64,
math_f: &impl Fn(f64) -> f64,
math_df: &impl Fn(f64) -> f64,
) {
let mut newton_solver = iterators::NewtonSolver::new(a, &math_f, &math_df);
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
});
for (idx, x_n) in newton_solver.enumerate().take_while(|(idx, _)| idx < &10) {
println!("{:>3}: {:>16.10}", idx, x_n);
}
}
(f, df)
fn iter_example_2(
a: f64,
math_f: &impl Fn(f64) -> f64,
math_df: &impl Fn(f64) -> f64,
) {
let mut newton_solver = iterators::NewtonSolver::new(a, &math_f, &math_df);
let first_few_guesses = newton_solver
.enumerate()
.take_while(|(idx, _)| idx < &10)
.map(|(_, guess)| guess)
.collect::<Vec<_>>();
for x_n in first_few_guesses {
println!("{:>16.10}", x_n);
}
}
#[cfg(test)]
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment