Интеграция осциллирующей (интерполированной) функции с GSL и python

Цифры в вашем примере не совсем понятны, и они сломают любую адаптивную схему. Позвольте мне объяснить, почему.

Вы пытаетесь интегрировать колебательную функцию с периодом 2*Pi от 0 до 10^10! Никакая адаптивная схема не сможет «увидеть» колебательное поведение на этом интервале, и они сойдутся с неправильным результатом (ложная сходимость)! Помните, что адаптивные схемы используют нисходящий подход. Они применяют какое-то правило ко всему интервалу, а затем делят этот интервал пополам и применяют одно и то же правило к каждому подразделу. После нескольких циклов (обычно 4 или 5) схема начинает проверять сходимость, сравнивая частичные результаты на последовательных шагах. В вашем примере схеме потребуется множество подразделений, чтобы наконец увидеть колебательное поведение, и это типичный случай, когда может произойти ложная сходимость!

Как интегрировать колебательную функцию на открытом интервале (a,\infinity)? Объяснение интегральная схема qawf вполне завершена. Интегрируйте функцию на подинтервалах, содержащих лишь несколько колебаний, и проверьте, как сходятся результаты, а затем экстраполируйте это!

Есть и другие цифры, которые не совсем понятны. Почему вам нужно сэмплировать sin(x)/(1+x^2) при каждом dx=1e-15? Любая разумная адаптивная схема может интегрировать sin(x) от 0 до 2Pi с ~10-20 точками выборки.

Правило Симпсона не сработало, потому что это не адаптивная схема. Код Python определит «dx» на основе предоставленного вами x-массива и будет использовать этот dx вплоть до 1e10! Однако я почти уверен, что ошибки округления в вашем коде довольно плохи, потому что вы выбрали dx ~ 1e-15.

РЕДАКТИРОВАТЬ 1 часть I: На самом деле проблема вызвана не только колебательным поведением подынтегральной функции. Учитывая, что конверт 1/x^2 сходится довольно быстро, ваша функция практически равна нулю, если x>>1. Итак, поскольку вы интегрируете эту огибающую в гигантский интервал [0,1e10], адаптивное интегрирование считает, что результат довольно мал, потому что он не может видеть небольшой (под)интервал, где функция не является незначительной. (вы можете подумать, что процедура интегрирования будет равномерно распределять точки оценки в близком интервале [0,1e10] - это не совсем верно для гауссовых интегралов, но это близко - так что вероятность того, что одна из этих точек попадет в интервал ~ [0,1e3], где подынтегральная функция не является незначительной, очень мала. После нескольких итераций процедура интегрирования получит, что ваш интеграл близок к нулю).

Изменить 1 часть II: я все еще думаю (после прочтения вашего комментария), что проблема в числах, которые вы подключили (или в оболочке python есть ошибка). Пожалуйста, попробуйте свой пример с разумными числами, как я сделал в следующем коде C++.

int main()  
{           
  const double omega = 1;
  auto g = [](double x)->double {return 1.0/(1.+x*x);};
  auto f = [&](double x)->double {return std::sin( omega * x) * g(x);};

  const int points_per_cycle  = 20;
  const int n_cycles = 10;
  const int size = n_cycles * points_per_cycle + 1;
  const double xmin = 0.0;
  const double xmax = n_cycles * (2 * M_PI);
  const double L = xmax-xmin;

  std::vector<double> x(size);
  std::vector<double> y(size);

  for (int i = 0; i <size; ++i) {
    x[i] = i * L/(size-1);
    y[i] = f(x[i]);
  }

  std::cout.precision(8); 
  // interpolation
  InterpolationGSL<std::vector<double>> GSLinterpol(x, y, GSLIT::cspline, false);
  // Integral of the interpolation
  std::cout << GSLinterpol.If((1+1e-12)*xmin, (1-1e-12)*xmax) << std::endl;

  // SECOND GSL INTEGRATION
  gsl_integration_workspace* w = gsl_integration_workspace_alloc (1000);

  gsl_integration_qawo_table* wf = gsl_integration_qawo_table_alloc 
    (omega, L, GSL_INTEG_SINE, 1000);

  int status = gsl_integration_qawo_table_set (wf, omega, L, GSL_INTEG_SINE);
  if(status) std::cerr<< "error: " << std::string(gsl_strerror (status)) << std::endl;

  double result;
  double abserr;

  std::function<double(double)> gg( std::cref(g) );
  GslFunction Fp(gg); 
  gsl_function *Fgsl = static_cast<gsl_function*>(&Fp);

  status = gsl_integration_qawo (Fgsl, xmin, 0.0, 1e-5, 1000, w, wf, &result, &abserr);

  if(status) std::cerr<< "error: " << std::string(gsl_strerror (status)) << std::endl;
  std::cout << result << std::endl;
  return 0;
}

Этот код использует мою gsl_function и interpolation обертки - так что вы можете найти код немного странным, но важно то, что он оценивает тот же интеграл, который вы упомянули, на разумном интервале, и результаты

Interpolation integral: 0.64631754
GSL integral: 0.64650827