1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125
use crate::List; use itertools::all; use na::RealField; use num_traits::NumCast; /// After so many iterations, consider the numerical method has failed to converge. pub const NUMBER_ITERATION_FAIL: usize = 1e6 as usize; /// Threshold that defines the convergence condition of the numerical Newton method. pub const NEWTON_METHOD_THRESHOLD: f64 = 1e-5; /// Newton's method algorithm. /// /// ## Definition /// /// Newton's method is an algorithm to find the approximated value of the exact solution of a /// variable in an equation that cannot be isolated from other terms. The algorithm consist in /// this, /// /// ```ignore /// loop: /// next_value = old_value - function(old_value, *args) / derivative(old_value, *args) /// if abs(next_value - old_value) < threshold: /// break /// ``` /// /// Newton's method function is expressed by moving all terms of the equation to one side to find /// `f(x) = 0`. The newton method's derivative is the derivative of this function. The `threshold` /// is the accepted residual between the estimated value and the exact solution. `old_value` is a /// initial guess of the solution. /// Newton's method is only guaranteed to converge if certain conditions are met (see /// [this](https://en.wikipedia.org/wiki/Newton%27s_method#Failure_analysis)). /// The maximum number of iterations is bounded to assume the algorithm could not converge. /// /// ## Usage /// /// In order to use this Newton's method implementaton, your need to: /// /// + create a struct to hold the arguments `*args` to be sent to both the function and the /// derivative of the Newton's method /// + add the trait [`NewtonMethodArguments`] to your struct /// + create the function of the Newton's method /// + create the derivative of the Newton's method /// /// This implementation uses [`List`] for input and output to allow component-wise computations. /// /// ## Example /// /// ``` /// use tool::{List, newton_method, NewtonMethodArguments, pows}; /// /// /// The struct that holds your arguments to be sent to the Newton's method's function and /// /// derivative. There you can define as many arguments as you want. /// pub struct MyArguments { /// some_value: f64, /// } /// /// /// We add the trait to your struct. /// impl NewtonMethodArguments for MyArguments {} /// /// /// Function of the Newton's method equation: /// /// /// /// f(x) = 100 - x ^ 4 - x * MyArguments.some_value = 0 /// pub fn newton_method_function( /// values: &List<f64>, /// args: &MyArguments, /// ) -> List<f64> { /// (-pows(values, 4)).add_scalar(100.) - values * args.some_value /// } /// /// /// Derivative of the Newton's method equation: /// /// /// /// f'(x) = - 4 * x ^ 3 - MyArguments.some_value /// pub fn newton_method_derivative( /// values: &List<f64>, /// args: &MyArguments, /// ) -> List<f64> { /// (- 4. * pows(values, 3)).add_scalar(-args.some_value) /// } /// /// // Initialize parameters. /// let initial_values = List::<f64>::zeros(3); /// let args = MyArguments { some_value: 12.0 }; /// /// // Calling the Newton's method. /// let results = newton_method( /// initial_values, /// newton_method_function, /// newton_method_derivative, /// args, /// ); /// ``` pub fn newton_method<T, A>( start_value: List<T>, newton_method_function: impl Fn(&List<T>, &A) -> List<T>, newton_method_derivative: impl Fn(&List<T>, &A) -> List<T>, newton_method_arguments: A, ) -> List<T> where T: RealField + NumCast, A: NewtonMethodArguments, { let mut new_value: List<T>; let mut current_value = start_value; let mut iteration = 0; 'convergence: loop { let func_res = newton_method_function(¤t_value, &newton_method_arguments); let deri_res = newton_method_derivative(¤t_value, &newton_method_arguments); new_value = ¤t_value - func_res.component_div(&deri_res); let all_converged = all((&new_value - ¤t_value).iter(), |residual| { residual.abs() < NumCast::from(NEWTON_METHOD_THRESHOLD).unwrap() }); if all_converged { break 'convergence; } else if iteration > NUMBER_ITERATION_FAIL { panic!("Newton Method could not converge."); } else { current_value = new_value; } iteration += 1; } new_value } /// Trait to be added to your custom struc holding Newton's method arguments. pub trait NewtonMethodArguments {}