multi_roots.h

Go to the documentation of this file.

 
 
#ifndef THEORETICA_MULTI_ROOTS_H
#define THEORETICA_MULTI_ROOTS_H
 
#include "../autodiff/autodiff.h"
 
 
namespace theoretica {
 
 
    template<unsigned int N>
    inline vec<real, N> multiroot_newton(
        autodiff::dvec_t<N>(*f)(autodiff::dvec_t<N>),
        vec<real, N> guess = vec<real, N>(0),
        real tolerance = OPTIMIZATION_MINGRAD_TOLERANCE,
        unsigned int max_iter = OPTIMIZATION_MINGRAD_ITER) {
 
        // Current position
        vec<real, N> x = guess;
 
        // Number of iterations
        unsigned int iter = 0;
 
        // The current value of f(x)
        vec<real, N> f_x = multidual<N>::extract_real(
            f(multidual<N>::make_argument(x))
        );
 
        // Jacobian matrix
        mat<real, N, N> J;
 
        while(f_x.sqr_norm() > square(tolerance) && iter <= max_iter) {
 
            // Compute the function value and Jacobian
            // using automatic differentiation
            multidual<N>::extract(
                f(multidual<N>::make_argument(x)), f_x, J
            );
 
            // Update the current best guess
            x -= algebra::solve(J, f_x);
            iter++;
        }
 
        if(iter > max_iter) {
            TH_MATH_ERROR("multi_root_newton", iter, NO_ALGO_CONVERGENCE);
            return vec<real, N>(nan());
        }
 
        return x;
    }
 
}
 
#endif