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

Question

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

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

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

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

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

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

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

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

F = f (двойной x, void *params)

и я, вероятно, могу использовать адаптивную интеграцию 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?

Answer 1

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 через.

Извиняюсь за любые опечатки.

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