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

Question

Я пытаюсь интегрировать сильно колебательные данные с функцией qawo из научной библиотеки GSL и pygsl в python. Поскольку я имею дело с данными, я думал, что интерполированная функция могла бы работать, но GSL дает мне неверный результат !!, Позвольте мне объяснить, взяв в качестве примера функцию Sin(x)/(1+x²).

Следующий код работает нормально:

import pygsl
from  pygsl import integrate

def f1(x,y):
    return 1./(1 + x**2)
sys = integrate.gsl_function(f1, None)

w = integrate.workspace(100)
cyclew = integrate.workspace(1000000)

table = integrate.qawo_table(1, 10000, integrate.SINE, 100)
flag, result, error = integrate.qawo(sys, 0, 1e-8, 1e-8, 100, w, table)

дает 0,626761 как следует. Но если мы смоделировали точки данных, используя вышеуказанную функцию...

xarr = np.linspace(0,1e15,1e30)  
yarr = np.sin(xarr)/(1.+xarr**2)

interp = interpol.interp1d(xarr,yarr)

def fgsl(x,y):
    return interp(x)

syst = integrate.gsl_function(fgsl, None)

w = integrate.workspace(1000)
cyclew = integrate.workspace(100000000)

table = integrate.qawo_table(1, 1e10, integrate.SINE, 100)

flag, result, error = integrate.qawo(syst, 0, 1e-15, 1e-15, 100, w, table)

что дает совершенно неверный результат: 4.2426e-21

Более того, если мы интегрировали yarr с функцией simps:

import scipy.integrate as ints
res = ints.simps(yarr,xarr)

дает довольно хорошее приближение: 0,64676099.

Просто предположим, что я не могу использовать правило Симпсона. Кто-нибудь знает, как я могу использовать функцию интерполяции с gsl? или как я могу изменить код, чтобы сделать интеграцию?

Заранее спасибо.

Answer 1

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

Вы пытаетесь интегрировать колебательную функцию с периодом 2* Пи от 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 integral: 0.64631754
GSL integral: 0.64650827