我正在尝试执行2F1(a,b,c,x)从-1到1的超几何函数积分以获得一些好的值（让我们说) .
我是在 Python 和 Mathematica 中完成的，答案是有限的，尽管它在这一点上是不同的. 无论如何，我需要将此集成作为我在 C++ 中很长的代码的一部分。
我正在使用 GSL（谢谢！这个神奇库的创建者）这样做。因为这个函数在我的积分上限是发散的（)，我用CQUAD来计算积分；然而，因为超几何的2F1输出既不是NaN也不是Infinity，我收到以下错误：
gsl: hyperg_2F1.c:685: ERROR: domain error
Default GSL error handler invoked.
你们中的任何人都可以帮我解决这个问题吗？
我不是 GSL 的常客，但我认为我通过包装2F1GSL 提供的超几何函数使其工作，它会NaN返回. 然后，数值积分将使用此信息用于其目的，并且domain error不会被抛出。
因此，（按照 GSL 文档中提供的示例）在函数中double f我比较了传递的参数sing_value在浮点意义上具有域绑定（不要判断命名）。如果它们足够接近，我返回NaN，否则，我调用 GSL 提供的2F1函数。
#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 = α
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;
}
根据 Wolfram Alpha 的说法，
这似乎已被代码确认
result = 2.772588722239701653
exact result = 2.772588722239781145
estimated error = 0.000000000000023683
actual error = -0.000000000000079492
# evaluations = 2553
该代码没有在容差、细分方面进行优化，尤其是，它无缘无故地混合了 C/C++。
在 gcc 6.4 和 GSL 2.4.2 上测试
这是安东方法的一个稍微简单的解决方案。GSLhyperg_2F1函数将引发域错误标志，并将返回值设置NaN为它不喜欢的输入值。集成例程知道如何处理NaN返回值，因此您需要做的就是关闭 GSL 错误处理程序。这可以通过以下方式完成：
#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 = α
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;
}
这是我得到的输出：
result = 2.772588722239768266
exact result = 2.772588722239781145
estimated error = 0.000000000000305493
actual error = -0.000000000000012879
# evaluations = 1565