Как интегрировать std::valarray‹double› с gsl?

Я относительно новичок в C++, но у меня есть некоторый (недостаточный) опыт кодирования и числовых вычислений.

Я знаю, что этот вопрос время от времени публикуется, как вы интегрируете массив. В MATLAB вы можете заставить свой массив быть функцией (я забыл, как это сделать, но я знаю, что делал это раньше) и отправить его встроенным интеграторам, поэтому мой вопрос заключается в том, как вы делаете это в C++.

У меня есть этот интеграл:

I = integral(A(z)*sin(qz)*dz)

q — просто двойная константа, z — интегрирующая переменная, но A(z) — это массив (с этого момента я буду называть его фактической функцией), который имеет то же количество точек, что и ось z в моем коде. Интегрирующие границы: z[0] и z[nz-1].

Я вычислил этот интеграл, используя правило трапеций, и для оси Z из 5000 точек это занимает 0,06 секунды. Моя проблема в том, что этот расчет происходит примерно 300 * 30 * 20 раз (у меня 3 цикла for), и эти 0,06 секунды очень быстро вырастают до 3 часов симуляции. И все узкое место моего кода — вот эта интеграция (я, очевидно, могу ускорить за счет уменьшения z, но не в этом дело).

Я знаю, что библиотечные функции обычно намного лучше, чем написанные пользователем. Я также знаю, что не могу использовать что-то более простое, как правило Симпсона, потому что подынтегральная функция сильно колеблется, и я хочу избежать собственной реализации какого-то сложного числового алгоритма.

GSL нужна функция в форме:

F = f(двойной х, пустота *параметры)

и я, вероятно, могу использовать адаптивную интеграцию QAWO из gsl, но как мне сделать свою функцию в форме, которая превращает мой массив в функцию?

Я думаю что-то вроде:

F(double z, void *params)
{
  std::valarray<double> actualfunction = *(std::valarray<double> *) params;
  double dz = *(double *) params; // Pretty sure this is wrong
  unsigned int actual_index = z / dz; // crazy assumption (my z[0] was 0)
  return actualfunction[actual_index];
 }

Возможно ли что-то подобное? Я сомневаюсь, что численный алгоритм будет использовать ту же пространственную разницу, что и фактическая функция, должен ли я тогда как-то выполнять интерполяцию фактической функции или что-то в этом роде?

Есть ли что-то лучше, чем gsl?


c++ integration gsl valarray
person Aleksandar Demić    schedule 06.07.2017    source источник

Ответы (1)


arrow_upward
0
arrow_downward 

template<class F>
struct func_ptr_helper {
  F f;
  void* pvoid(){ return std::addressof(f); }
  template<class R, class...Args>
  using before_ptr=R(*)(void*,Args...);
  template<class R, class...Args>
  using after_ptr=R(*)(Args...,void*);
  template<class R, class...Args>
  static before_ptr<R,Args...> before_func() {
    return [](void* p, Args...args)->R{
      return (*static_cast<F*>(p))(std::forward<Args>(args)...);
    };
  }
  template<class R, class...Args>
  static after_ptr<R,Args...> after_func() {
    return [](Args...args, void* p)->R{
      return (*static_cast<F*>(p))(std::forward<Args>(args)...);
    };
  }
};
template<class F>
func_ptr_helper<F> lambda_to_pfunc( F f ){ return {std::move(f)}; }

использовать:

auto f = lambda_to_pfunc([&actualfunction, &dz](double z){
  unsigned int actual_index = z / dz; // crazy assumption (my z[0] was 0)
  return actualfunction[actual_index];
});

потом

void* pvoid - f.pvoid();
void(*pfun)(double, void*) = f.after_func();

и вы можете пройти pfun и pvoid.

Приносим извинения за любые опечатки.

Идея в том, что мы пишем лямбду, которая делает то, что мы хотим. Затем lambda_to_pfunc оборачивает его, чтобы мы могли передать его как void* и указатель функции в API-интерфейсы в стиле C.

person Yakk - Adam Nevraumont    schedule 06.07.2017