GSL 集成行为奇怪答案

【问题标题】：GSL integration behaves strangeGSL 集成行为奇怪
【发布时间】：2013-08-28 14:56:22
【问题描述】：

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

变量偏移量基本上创建了平移。对于 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 = &alpha;

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

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

    return 0;
}

提前致谢，尼基尔

【问题讨论】：

你能解决这个问题吗？我复制了您的问题，但由于对函数的有限支持和 qagi 使用的低阶 (15) 高斯近似，GSL 错过了被积函数，我想不出其他解释。
嗨，我认为你是对的。问题是由于两个原因：（a）您正确指出的支持有限，但更重要的是，（b）由于 x= (1-t)/t 变换导致的非线性收缩因子，它收缩了 a因子 1/t^2... 在 t = 1000 时，由于 10^-6 的收缩，支撑基本上看起来像一条狭窄的垂直线:(...我还没有找到一个好的方法来获得摆脱这个问题。我的主要问题是两个（或更多）由 x = 1000+ 分隔的高斯需要积分......所以我需要一些绝对好的方法来处理积分，没有太多损失
这些功能的大致支持你都知道，所以不用qagi了！从 x_median - epsilon 到 x_median + epsilon 的简单 qag 适合，其中 epsilon 是一个合理的数字。对每个高斯分别执行此操作，然后对最终结果求和。没有适用于所有情况的通用自适应积分器。您需要向集成商提供有关该功能的任何“额外”信息（限制集成范围是提供这些额外信息的一种方式）
但这是一个很好的数字问题。我给了 +1 票
但如果你真的想要一个非常强大且速度慢的积分器，请尝试新 GSL 1.15/1.16 版本中的 CQUAD 双自适应

标签： c++ gsl

【解决方案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。那么，当偏移量为零时，积分器错过您的函数的机会非常小（这就是为什么使用 qagi 在 x=0 处积分合理的函数没有问题的原因）。但是，当您将 x 平移一个偏移量时，您会将整个函数（其支持非常有限）置于 30 个评估点中的两个之间。在这种情况下，GSL 完全错过了您的函数并返回零。

简而言之：GSL 试图在第一次迭代中仅使用 30 个点来分析整个 [-infinity, infinity] 区间。它错过以任意 x 为中心的支持非常有限的函数的可能性非常高！仅当您的功能支持非常大时才使用 qagi！

【讨论】：

我可以问一个愚蠢的问题吗？在 cmets 和您的帖子中，您使用了“有限支持”一词。我想知道它指的是什么？
支持是函数非零的区间/区间集合。