bracketingRootFinder.hpp

Go to the documentation of this file.

#ifndef JULIAN_BRACKETINGROOTFIDNER_HPP
#define JULIAN_BRACKETINGROOTFIDNER_HPP

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

namespace julian {

  enum class BracketingRootFinder {BISECTION,  
                   FALSEPOS,  
                   BRENT_DEKKER  
  };

  template<typename T>
  double  bracketingRootFinder(T f, double x_lo, double x_hi, BracketingRootFinder type, double precision = 1e-12, int max_iter = 100) {
    GslFunctionAdapter<T> pF(f);
    gsl_function *gsl_f = static_cast<gsl_function*>(&pF); 
    double root {0.0};
    
    int status;
    int iter = 0;

    const gsl_root_fsolver_type *type_of_solver;
    gsl_root_fsolver *s;
    
    switch (type){
    case BracketingRootFinder::BISECTION: type_of_solver = gsl_root_fsolver_bisection; break;
    case BracketingRootFinder::FALSEPOS: type_of_solver = gsl_root_fsolver_falsepos; break;
    case BracketingRootFinder::BRENT_DEKKER: type_of_solver = gsl_root_fsolver_brent; break;
    }
    
    s = gsl_root_fsolver_alloc (type_of_solver);
    gsl_root_fsolver_set (s, gsl_f, x_lo, x_hi);
    
    do {
      iter++;
      status = gsl_root_fsolver_iterate (s);
      root = gsl_root_fsolver_root (s);
      x_lo = gsl_root_fsolver_x_lower (s);
      x_hi = gsl_root_fsolver_x_upper (s);
      status = gsl_root_test_interval (x_lo, x_hi, 0, precision);
    } while (status == GSL_CONTINUE && iter < max_iter);
    
    gsl_root_fsolver_free (s);
    return root;
  }
}  // namespace julian
#endif