Labels

Friday, August 2, 2024

GSL integration

//https://stackoverflow.com/questions/47038457/gsl-integration-within-a-function-c

#include 
#include 
#include 
#include 

namespace details
{
  extern "C" {
  double gsl_function_cpp_wrapper_helper(double z, void *cpp_f)
  {
    const auto p_cpp_f =
        reinterpret_cast *>(cpp_f);
    return (*p_cpp_f)(z);
  }
  }
}

class gsl_function_cpp_wrapper
{
  using Cpp_F = std::function;

 public:
  operator const gsl_function_struct *() const { return &gsl_f_; };

  gsl_function_cpp_wrapper(Cpp_F &&cpp_f) : cpp_f_(std::move(cpp_f))
  {
    gsl_f_.function = details::gsl_function_cpp_wrapper_helper;
    gsl_f_.params = &cpp_f_;
  }

  gsl_function_cpp_wrapper(double(f)(double))
      : gsl_function_cpp_wrapper(Cpp_F(f))
  {
  }

  template 
  gsl_function_cpp_wrapper(double (OBJ::*method)(double) const, const OBJ *obj)
      : gsl_function_cpp_wrapper(
            Cpp_F(std::bind(method, obj, std::placeholders::_1)))
  {
  }
  template 
  gsl_function_cpp_wrapper(double (OBJ::*method)(double), OBJ *obj)
      : gsl_function_cpp_wrapper(
            Cpp_F(std::bind(method, obj, std::placeholders::_1)))
  {
  }

 protected:
  Cpp_F cpp_f_;
  gsl_function_struct gsl_f_;
};

//----------------

class Class_Example
{
 public:
  double f(double x) const { return std::log(alpha_ * x) / std::sqrt(x); }

 protected:
  double alpha_ = 1;
};

double free_function_example(double x)
{
  const double alpha = 1;
  return std::log(alpha * x) / std::sqrt(x);
}

//----------------

int main()
{
  double result, error;
  gsl_integration_workspace *const w = gsl_integration_workspace_alloc(1000);
  assert(w != nullptr);

  //----------------

  Class_Example class_example;

  gsl_function_cpp_wrapper wrapper_1(&Class_Example::f, &class_example);

  gsl_integration_qags(wrapper_1, 0, 1, 0, 1e-7, 1000, w, &result, &error);

  std::printf("result          = % .18f\n", result);
  std::printf("estimated error = % .18f\n", error);

  //----------------

  gsl_function_cpp_wrapper wrapper_2(free_function_example);

  gsl_integration_qags(wrapper_2, 0, 1, 0, 1e-7, 1000, w, &result, &error);

  std::printf("result          = % .18f\n", result);
  std::printf("estimated error = % .18f\n", error);

  //----------------

  gsl_integration_workspace_free(w);

  return EXIT_SUCCESS;
}

No comments:

Post a Comment