src/newton.rs

/*!
Newton method in 2D.
*/

use alg_tools::types::*;

/// Approximately solves $f(x)=b$ using Newton's method in 1D.
///
/// The function `g` should return $(f'(x), f(x))$.
/// The initial iterate will be `x`, and exactly `iters` iterations are taken.
#[inline]
pub fn newton_sym1x1<F : Float>(
    g : impl Fn(F) -> (F, F),
    mut x : F,
    iters : usize
) -> F {
    for _i in 0..iters {
        let (a, b) = g(x);
        x -= b / a
    }
    x
}

/// Approximately solves $f(x)=b$ using Newton's method in "D.
///
/// The function `g` should return $(∇f(x), f(x))$.
/// The Hessian $A=∇f(x)$ should be symmetric and given in the form $[A_{11}, A_{12}, A_{22}]$.
/// The initial iterate will be `[x1, x2]`, and exactly `iters` iterations are taken.
#[inline]
pub fn newton_sym2x2<F : Float>(
    g : impl Fn(F, F) -> ([F; 3], [F; 2]),
    [mut x1, mut x2] : [F; 2],
    iters : usize
) -> [F; 2] {
    for _i in 0..iters {
        let ([a11, a12, a22], [b1, b2]) = g(x1, x2);
        let q = a11 * a22 - a12 * a12;
        x1 -= (a22 * b1 - a12 * b2) / q;
        x2 -= (a11 * b2 - a12 * b1) / q;
    }
    [x1, x2]
}