c - gsl 无限积分区间错误。发现不良的被积函数行为。如何解决？

Question

在尝试使用 C 中的 GSL 对无限区间 [0,inf) 进行数值积分后，我收到以下错误消息。

gsl: qags.c:553: ERROR: bad integrand behavior found in the integration interval
Default GSL error handler invoked.
Command terminated by signal 6

这是我正在集成的功能 $

double dI2dmu(double x, void * parametros){
  double *p,Ep,mu,M,T;  
  p=(double *) parametros;

  M=p[0];
  T=p[1];
  mu=p[2];

  Ep=sqrt(x*x+M*M);

  double fplus= -((exp((Ep - mu)/T)/(pow(1 + exp((Ep - mu)/T),2)*T) - exp((Ep + \
mu)/T)/(pow(1 + exp((Ep + mu)/T),2)*T))*pow(x,2))/(2.*Ep*pow(PI,2));
  return fplus;
}

以及集成过程的代码

 params[0]=0.007683; //M
 params[1]=0.284000;// T
 params[2]=0.1;   //mu

    gsl_function dI2mu_u; 
    dI2mu_u.function = &dI2dmu;
    dI2mu_u.params = &params;
    gsl_integration_qagiu (&dI2mu_u, 0, 0, 1e-7, 100000,
             w, &resultTest2, &error1Test2);

函数有以下几个方面:

在我看来，它的行为非常好。因此，我没有执行无限集成，而是执行了我认为可重新分区的上限的集成，例如:

  gsl_function G;
 G.function = &dI2dmu;
 G.params = &params;

 gsl_integration_qags (&G, 0, 1e2*A, 0, 1e-7, 100000,
                    w, &result1, &error1); 

得到与Mathematica结果一致的结果进行无限积分

result definite up to 10*A        =  0.005065263943958745
result up to infinity             =  nan
Mathematica result up to infinity =  0.005065260000000000

但是 GSL 无限积分一直是“nan”。有任何想法吗？在此先感谢您的帮助。

Answer 1

正如@yonatan zuleta ochoa 正确指出的那样，问题出在 exp(t)/pow(exp(t)+1,2) 中。 . exp(t)可以溢出 ieee754 DBL_MAX对于 t 的值低至 nextafter(log(DBL_MAX), INFINITY) ，即 ~7.09783e2 .

当 exp(t) == INFINITY ,

exp(t)/pow(exp(t)+1,2) == ∞/pow(∞+1,2) == ∞/∞ == NAN

Yonatan 提出的解决方案是使用对数，可以按如下方式完成:

exp(t)/pow(exp(t)+1,2) == exp(log(exp(t)) - log(pow(exp(t)+1,2)))
                       == exp(t - 2*log(exp(t)+1))
                       == exp(t - 2*log1p(exp(t))) //<math.h> function avoiding loss of precision for log(exp(t)+1)) if exp(t) << 1.0

这是一个完全合理的方法，避免了 NAN高达非常高的值 t .但是，在您的代码中，t == (Ep ± mu)/T可以是INFINITY如果abs(T) < 1.0对于 x 的值靠近DBL_MAX , 即使x 不是无穷大。在这种情况下，减法 t - 2*log1p(exp(t))变成 ∞ - ∞ ，即 NAN再次。

另一种方法是替换 exp(x)/pow(exp(x)+1,2)与 1.0/(pow(exp(x)+1,2)*pow(exp(x), -1))将分母和分子除以 exp(x) (对于任何有限的 x 都不为零)。这简化为 1.0/(exp(x)+exp(-x)+2.0) .

这里是避免NAN的函数的实现对于 x 的值直至并包括 DBL_MAX :

static double auxfun4(double a, double b, double c, double d)
{
  return 1.0/(a*b+2.0+c*d);
}
double dI2dmu(double x, void * parametros)
{
  double *p = (double *) parametros;
  double invT = 1.0/p[1];
  double Ep = hypot(x, p[0]);
  double muexp = exp(p[2]*invT);
  double Epexp = exp(Ep*invT);
  double muinv = 1.0/muexp;
  double Epinv = 1.0/Epexp;
  double subterm = auxfun4(Epexp, muinv, Epinv, muexp);
  subterm -= auxfun4(Epexp, muexp, Epinv, muinv);
  double fminus = subterm*(x/Ep)*invT*(0.5/(M_PI*M_PI))*x;;
  return -fminus;
}

此实现还使用了 hypot(x,M) , 而不是 sqrt(x*x, M*M) , 并避免计算 x*x通过重新排列乘法/除法的顺序到组 x/Ep一起。自 hypot(x,M)将是 abs(x)对于 abs(x) >> abs(M) , 术语 x/Ep方法 1.0对于大 x .