无限积分间隔中的 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

我认为这里的问题是，与 Mathematica 不同，C 在计算中不使用任意精度。然后，在计算 Exp [Ep] 时，数值计算会溢出。

现在，GSL 使用变换 x = (1-t)/t 映射到区间 (0,1]。 因此，对于 t

Exp[ ( Ep(x) - \Mu)/T ] / { 1 + Exp[( Ep(x) - \Mu )/T] }^2

使用 A/B = Exp[ Ln A - Ln B]，您可以获得更好的数值行为。

我会尝试，如果我有好的结果，然后我会告诉你。

解决方案

正如我之前所说，您必须注意不确定形式所带来的问题。所以，让我们用对数版本写出有问题的术语：

  double dIdmu(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  -2.0*log(1.0 + exp((Ep - mu)/T) ) ) -  exp( (Ep + mu)/T  -2.0*log(1.0 + exp((Ep + mu)/T) ) ) )  * pow(x,2)   /  (2.*  T * Ep*pow(M_PI,2));

return fplus;
        }

还有这个主要功能

    int main()
{
  double params[3];

  double resultTest2, error1Test2;

  gsl_integration_workspace * w 
    = gsl_integration_workspace_alloc (10000);

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

    gsl_function dI2mu_u; 
    dI2mu_u.function = &dIdmu;
    dI2mu_u.params = &params;
    gsl_integration_qagiu (&dI2mu_u, 0.0, 1e-7, 1e-7, 10000, w, &resultTest2, &error1Test2);


    printf("%e\n", resultTest2);
    gsl_integration_workspace_free ( w);

    return 0;
}

你会得到答案： -5.065288e-03。 我很好奇……这就是我在 Mathematica 中定义函数的方式

所以比较答案：

• GSL -5.065288e-03
• Mathematica -0.005065287633739702

Answer 2

正如@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 值。但是，在您的代码中，如果abs(T) < 1.0 的值接近DBL_MAX，t == (Ep ± mu)/T 可以是INFINITY，即使x 是非无穷大。在这种情况下，减法t - 2*log1p(exp(t))变成∞ - ∞，又是NAN。

另一种方法是用1.0/(pow(exp(x)+1,2)*pow(exp(x), -1)) 替换exp(x)/pow(exp(x)+1,2)，方法是将分母和分子都除以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/Ep 组合在一起来避免计算x*x。因为hypot(x,M) 将是abs(x) 对应abs(x) >> abs(M)，所以术语x/Ep 接近1.0 对应大x。