derivativeRootFinder.hpp

Go to the documentation of this file.

#ifndef JULIAN_DERIVATIVEROOTFINDER_HPP
#define JULIAN_DERIVATIVEROOTFINDER_HPP

#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>
#include <mathematics/numericalAlgorithms/GslFunctionAdapter.hpp>
#include <mathematics/numericalAlgorithms/GslFunctionFdfAdapter.hpp>
#include <gsl/gsl_deriv.h>

namespace julian {

  enum class DerivativeRootFinder {NEWTON_RAPHSON , 
                   SECANT,  
                   STEFFENSEN 
  };

  template<typename T, typename U>
  double derivativeRootFinder(T f, U df, double x_0, DerivativeRootFinder type,double precision = 1e-12, int max_iter = 100) {
    GslFunctionFdfAdapter<decltype(f),decltype(df)> Fdf(f,df);
    int status;
    int iter = 0;
    double x = x_0;
    const gsl_root_fdfsolver_type *Type = gsl_root_fdfsolver_newton;
    gsl_root_fdfsolver *s;
    
    switch (type) {
    case DerivativeRootFinder::NEWTON_RAPHSON: Type = gsl_root_fdfsolver_newton; break;
    case DerivativeRootFinder::SECANT: Type = gsl_root_fdfsolver_secant; break;
    case DerivativeRootFinder::STEFFENSEN: Type = gsl_root_fdfsolver_steffenson; break;
    }
    
    s = gsl_root_fdfsolver_alloc (Type);
    gsl_root_fdfsolver_set (s, &Fdf, x_0);

    do {
      iter++;
      status = gsl_root_fdfsolver_iterate (s);
      x_0 = x;
      x = gsl_root_fdfsolver_root (s);
      status = gsl_root_test_delta (x, x_0, 0, precision);
    } while (status == GSL_CONTINUE && iter < max_iter);
    
    gsl_root_fdfsolver_free (s);
    return x;
  }

  template<typename T>
  double derivativeRootFinder(T f, double x_0, DerivativeRootFinder type,double precision = 1e-12, int max_iter = 100, double h = 1e-6) {
    GslFunctionAdapter<decltype(f)> pF(f);
    gsl_function *gsl_f = static_cast<gsl_function*>(&pF); 
    
    auto df = [&](double x)->double{
      double result, abserr;
      gsl_deriv_forward (gsl_f, x, h, &result, &abserr);
      return result;
    };
    
    GslFunctionFdfAdapter<decltype(f),decltype(df)> Fdf(f,df);
    
    int status;
    int iter = 0;
    double x = x_0;
    const gsl_root_fdfsolver_type *Type;
    gsl_root_fdfsolver *s;
    
    switch (type) {
    case DerivativeRootFinder::NEWTON_RAPHSON: Type = gsl_root_fdfsolver_newton; break;
    case DerivativeRootFinder::SECANT : Type = gsl_root_fdfsolver_secant; break;
    case DerivativeRootFinder::STEFFENSEN : Type = gsl_root_fdfsolver_steffenson; break;
    }
    
    s = gsl_root_fdfsolver_alloc (Type);
    gsl_root_fdfsolver_set (s, &Fdf, x_0);
    
    do {
      iter++;
      status = gsl_root_fdfsolver_iterate (s);
      x_0 = x;
      x = gsl_root_fdfsolver_root (s);
      status = gsl_root_test_delta (x, x_0, 0, precision);
    } while (status == GSL_CONTINUE && iter < max_iter);

    gsl_root_fdfsolver_free (s);
    return x;
  }
}  // namespace julian
#endif