c++ - GSL 集成行为异常

Question

我正在使用 gsl_integration_qagi 例程对 (-infty, +infty) 进行集成。我的期望是，在沿 x 轴平移时，积分结果(曲线下面积)不应改变。然而，这不是，我的观察者。我在某处犯错了吗？代码附在下面:

变量 offset 基本上创建了翻译。对于 0、10.0、20.0 的偏移值，该区域保持相同(如预期)，但在偏移 ~ 40.0 后突然下降到零

double offset=200.0;

double f (double x, void * params) {
   double alpha = *(double *) params;
   x += offset;
   double f = exp(-x*x);
   return f;
}

int main(int argc, const char * argv[])
{

    gsl_integration_workspace * w
     = gsl_integration_workspace_alloc (1000);

    double result, error;
    double expected = -4.0;
    double alpha = 1.0;

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

    gsl_integration_qagi (&F, 0, 0.001, 1000,
                          w, &result, &error);

    printf ("result = % .18f\n", result);

    return 0;
}

提前致谢，尼克尔

Answer 1

您正在尝试使用 qagi 集成一个支持非常有限的功能，这很糟糕。积分完全错过被积函数的可能性很大。为什么？

Qagi 使用 15 点高斯规则。这大约意味着它将在以下固定点(第一次迭代)评估函数

  const double center = 0.5 * (a + b);
  const double half_length = 0.5 * (b - a);
  const double abscissa = half_length * xgk[jtw];
  const double fval1 = GSL_FN_EVAL (f, center - abscissa);
  const double fval2 = GSL_FN_EVAL (f, center + abscissa);  

在哪里

  static const double xgk[8] =    /* abscissae of the 15-point kronrod rule */
  {
     0.991455371120812639206854697526329,
     0.949107912342758524526189684047851,
     0.864864423359769072789712788640926,
     0.741531185599394439863864773280788,
     0.586087235467691130294144838258730,
     0.405845151377397166906606412076961,
     0.207784955007898467600689403773245,
     0.000000000000000000000000000000000
 };

(这直接取自 GSL 代码)。然后，根据 GSL 从这些点获得的值，它可以进一步划分特定区域并再次应用此规则。

从非线性变换 x = (1-t)/t(这是 gsl 应用于将 [-infinity, infinity] 映射到 (0-1] 区间的变换)，我们可以假设 x = 0 映射到 t = 1。此外，评估点之一是 t = half_length (1.0 + 0.991455371120812639206854697526329) ~ 1。那么，当offset is zero 非常小(这就是为什么使用 qagi 在 x=0 处集成合理的函数中心没有问题)。但是，当您将 x 平移一个偏移量时，您适合整个函数(支持非常有限)在 30 个评估点中的两个之间。在这种情况下，GSL 完全错过您的功能并返回零。

简单总结:GSL 试图在第一次迭代中仅使用 30 个点来分析整个 [-infinity, infinity] 区间。它错过以任意 x 为中心的支持非常有限的功能的可能性非常高!只有当你的功能支持非常大时才使用 qagi!