Integration of a diverge function in c++ GSL Library

Question

I am trying to perform an Integral of Hypergeometric function 2F1(a,b,c,x) from -1 to 1 for some good values of $a,b,c$ (lets say $a=1,b=2,c=3$) .

I did it in Python and Mathematica and the answer is finite although it is divergent at the point $x=1$. Anyway, I need to do this integration as a part of my very long code in c++.

I am using GSL (thank you! creators of this amazing library) to do so. Since this function is divergent at the upper bound of my integration ($x=1$), I used the CQUAD to calculate the integral; however because the output of Hypergeometric 2F1 for $x=1$ is neither NaN nor Infinity, I get the following error:

gsl: hyperg_2F1.c:685: ERROR: domain error 
Default GSL error handler invoked.

Any of you guys can help me to resolve this problem?

Answer 1

I am not a frequent user of GSL, but I think I made it work by wrapping the hypergeometric function 2F1 provided by GSL in such a way, that it would return NaN for $x=1$. Then, the numerical integration will use this information for its purposes and the domain error will not be thrown.

So, (following the example provided in GSL docs) in function double f I compared the passed argument $x$ with a domain bound sing_value (don't judge the naming) in the floating-point sense. If they are sufficiently close, I return NaN, otherwise, I call the provided by GSL 2F1 function.

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_sf_hyperg.h>
#include <limits>


double f (double x, void * params) {
  double alpha = *(double *) params;
  double sing_value = 1.0;
  double diff = std::abs(x-sing_value);
  bool singular_point = (diff <= std::numeric_limits<double>::epsilon()*std::abs(x+sing_value)*2);
  // as per Kirill's comment, the comparison 
  // (diff < std::numeric_limits<double>::min()) 
  // is unnecessary for any practical sense.
  if (singular_point)
    return std::numeric_limits<double>::quiet_NaN();
  else
    return gsl_sf_hyperg_2F1(1,2,3,alpha*x);
}

int
main (void)
{
  gsl_integration_cquad_workspace * wcquad = gsl_integration_cquad_workspace_alloc (10000);

  double result, error;
  size_t nevals;
  double expected = log(16);
  double alpha = 1.0;

  gsl_function F;
  F.function = &f;
  F.params = &alpha;

  gsl_integration_cquad (&F, -1.0, 1.0, 0, 1e-14, wcquad, &result, &error,&nevals);

  printf ("result          = % .18f\n", result);
  printf ("exact result    = % .18f\n", expected);
  printf ("estimated error = % .18f\n", error);
  printf ("actual error    = % .18f\n", result - expected);
  printf ("# evaluations   = %zu\n", nevals);

  gsl_integration_cquad_workspace_free (wcquad);

  return 0;
}

According to Wolfram Alpha, $$ \int_{-1}^1 {}_2F_1(1,2,3,x)dx = \ln 16 $$

which seems to be confirmed by the code

result          =  2.772588722239701653
exact result    =  2.772588722239781145
estimated error =  0.000000000000023683
actual error    = -0.000000000000079492
# evaluations   = 2553

The code is not optimized in terms of tolerances, subdivisions and, especially, it being an awful mixture of C/C++ for no reason.

Tested on gcc 6.4 and GSL 2.4.2

Answer 2

Here is a slighly simpler solution to Anton's method. The GSL hyperg_2F1 function will raise a domain error flag as well as set the return value to NaN for input values it doesn't like. The integration routines know how to handle NaN return values, so all you need to do is turn off the GSL error handler. This can be accomplished as follows:

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_sf_hyperg.h>
#include <gsl/gsl_errno.h>

double f (double x, void * params) {
  double alpha = *(double *) params;
  return gsl_sf_hyperg_2F1(1,2,3,alpha*x);
}

int
main (void)
{
  gsl_integration_cquad_workspace * wcquad = gsl_integration_cquad_workspace_alloc (10000);

  double result, error;
  size_t nevals;
  double expected = log(16);
  double alpha = 1.0;

  gsl_set_error_handler_off();

  gsl_function F;
  F.function = &f;
  F.params = &alpha;

  gsl_integration_cquad (&F, -1.0, 1.0, 0, 1e-14, wcquad, &result, &error,&nevals);

  printf ("result          = % .18f\n", result);
  printf ("exact result    = % .18f\n", expected);
  printf ("estimated error = % .18f\n", error);
  printf ("actual error    = % .18f\n", result - expected);
  printf ("# evaluations   = %zu\n", nevals);

  gsl_integration_cquad_workspace_free (wcquad);

  return 0;
}

Here is the output I get:

result          =  2.772588722239768266
exact result    =  2.772588722239781145
estimated error =  0.000000000000305493
actual error    = -0.000000000000012879
# evaluations   = 1565