integrate.hpp

Go to the documentation of this file.

#ifndef JULIAN_INTEGRATION_HPP
#define JULIAN_INTEGRATION_HPP

#include <mathematics/numericalAlgorithms/GslFunctionAdapter.hpp>
#include <gsl/gsl_integration.h>

namespace julian {

  enum class Integration {NON_ADAPTIVE_GAUSS_KRONROD,  
              ADAPTIVE_GAUSS_KRONROD,  
              ADAPTIVE_GAUSS_KRONROD_WITH_SINGULARITIES,  
              ADAPTIVE_GAUSS_KRONROD_INFINITE_RANGE,  
              ADAPTIVE_GAUSS_KRONROD_UPPER_INFINITE,  
              ADAPTIVE_GAUSS_KRONROD_LOWER_INFINITE  
  }; 

  template<typename F>
  double integrate(F f, double x_lo, double x_hi, double precision_abs,double precision_rel, int n, Integration t) {
    GslFunctionAdapter<F> pF(f);
    gsl_function *gsl_f = static_cast<gsl_function*>(&pF); 
    
    double ret;
    double error;
    size_t number_of_eval;
    
    gsl_integration_workspace * workspace = gsl_integration_workspace_alloc (n);
    
    switch (t) {
    case Integration::NON_ADAPTIVE_GAUSS_KRONROD:
      gsl_integration_qng(gsl_f, x_lo, x_hi, precision_abs, precision_rel, &ret, &error, &number_of_eval);break;
    case Integration::ADAPTIVE_GAUSS_KRONROD:
      gsl_integration_qag(gsl_f, x_lo, x_hi, precision_abs, precision_rel, n, 4, workspace, &ret, &error);break; // 4 means it is a 41 point rule
    case Integration::ADAPTIVE_GAUSS_KRONROD_WITH_SINGULARITIES:
      gsl_integration_qags (gsl_f, x_lo, x_hi, precision_abs, precision_rel, n, workspace, &ret, &error);break;
    case Integration::ADAPTIVE_GAUSS_KRONROD_INFINITE_RANGE:
      gsl_integration_qagi (gsl_f, precision_abs, precision_rel, n, workspace, &ret, &error);break;
    case Integration::ADAPTIVE_GAUSS_KRONROD_LOWER_INFINITE:
      gsl_integration_qagil(gsl_f, x_lo, precision_abs, precision_rel, n, workspace, &ret, &error );break;
    case Integration::ADAPTIVE_GAUSS_KRONROD_UPPER_INFINITE:
      gsl_integration_qagiu(gsl_f, x_hi, precision_abs, precision_rel, n, workspace, &ret, &error );break;    
    }
    
    gsl_integration_workspace_free (workspace);
    return ret;
  }
}  // namespace julian
  
#endif /* INTEGRATION_HPP */