Интеграция GSL ведет себя странно

Question

Я делаю интеграцию через (-infty, +infty), используя процедуру gsl_integration_qagi. Я ожидаю, что при перемещении по оси 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 = α

    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, Тогда шансы, что интегратор пропустит вашу функцию, когда смещение равно нулю, весьма малы (поэтому нет проблем интегрировать разумный центр функций в точке x = 0 с помощью qagi). Однако, когда вы переводите x по смещению, вы помещаете всю функцию (которая имеет очень ограниченную поддержку) между двумя из 30 точек оценки. В этом случае GSL полностью пропускает вашу функцию и возвращает ноль.

Резюме простыми словами: GSL пытается проанализировать весь интервал [-infinity, infinity], используя только 30 точек в первой итерации. Вероятность того, что он пропустит функцию с очень ограниченной поддержкой с центром в произвольном x, очень высока! Используйте qagi, только если поддержка вашей функции очень велика!