2
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ Do you have tried QAGS instead of CQUAD? $\endgroup$ – gammatester Apr 30 '18 at 18:37
  • $\begingroup$ Yes I have tried QAG, QAGS , QAGP and CQUAD. $\endgroup$ – Kamran Salehi Vaziri Apr 30 '18 at 19:54
  • $\begingroup$ Can you use the tanh-sinh quadrature in Boost.Math? $\endgroup$ – user14717 May 1 '18 at 2:39
3
$\begingroup$

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

$\endgroup$
  • 2
    $\begingroup$ Nitpick (I only mention this because I think it's a common source of confusion, I hope it's okay): so long as abs(sing_value)>1e-290, the comparison diff < std::numeric_limits<double>::min() is always false and can be omitted. The smallest positive real double is almost counter-intuitively small. The comparison is in most cases equivalent to an exact comparison with zero, diff == 0. $\endgroup$ – Kirill May 2 '18 at 11:15
  • $\begingroup$ @Kirill that's actually a very nice detail. I never (years!!!!!) bothered to get a sense of what I am comparing to, just did it based on what I think is right as opposed to reality. $\endgroup$ – Anton Menshov May 2 '18 at 14:57
2
$\begingroup$

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
$\endgroup$
  • $\begingroup$ nice trick. I wonder, why is the number of evaluations in your case so much lower than in mine with the same accuracy parameters. $\endgroup$ – Anton Menshov Sep 27 '18 at 17:03
  • 1
    $\begingroup$ I noticed that too and don't have an answer. If I had to guess, the hyperg_2F1 function may return valid results even in cases where your singular_point test returns NaN. $\endgroup$ – vibe Sep 27 '18 at 20:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.