tool/core/
numerical_algorithms.rs

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

tool/core/
numerical_algorithms.rs
