peroxide/numerical/
newton.rs

1use crate::numerical::utils::jacobian;
2use crate::structure::ad::*;
3use crate::traits::{
4    math::{Norm, Normed, Vector},
5    matrix::LinearAlgebra,
6    mutable::MutFP,
7};
8
9/// Newton-Raphson Method
10pub fn newton<F: Fn(&Vec<AD>) -> Vec<AD> + Copy>(init_cond: Vec<f64>, f: F, rtol: f64) -> Vec<f64> {
11    let mut x_next = init_cond;
12    let mut x = x_next.clone();
13    update(&mut x_next, f);
14    let mut err = x_next.sub_vec(&x).norm(Norm::L2);
15
16    while err >= rtol {
17        x = x_next.clone();
18        update(&mut x_next, f);
19        err = x_next.sub_vec(&x).norm(Norm::L2);
20    }
21
22    x_next
23}
24
25fn update<F: Fn(&Vec<AD>) -> Vec<AD> + Copy>(xs: &mut Vec<f64>, f: F) {
26    let j = jacobian(f, &xs);
27    let pinv_j = j.pseudo_inv();
28    //let fx = f(NumberVector::from_f64_vec(xs.clone())).to_f64_vec();
29    let xs_ad: Vec<AD> = xs.iter().map(|&t| AD::from(t)).collect();
30    let fx: Vec<f64> = f(&xs_ad).iter().map(|&t| t.into()).collect();
31
32    xs.mut_zip_with(|x, y| x - y, &(pinv_j * fx))
33}