Comment intégrer std :: valarray <double> avec gsl?

Question

Je suis relativement nouveau au C++, mais j'ai quelques (rares) codages et expériences numériques.

Je sais que cette question est publiée de temps en temps, comment intégrer un tableau. Dans MATLAB vous pouvez forcer votre tableau à être une fonction (j'ai oublié comment, mais je sais que je l'ai déjà fait) et l'envoyer aux intégrateurs intégrés, donc ma question est de savoir comment le faire en C++.

J'ai cette intégrale:

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

q est juste à double const, z est la variable d'intégration, mais A (z) est un tableau (je vais l'appeler actualfunction à partir de maintenant) qui a même nombre de points en tant qu'axe z dans mon code. Les limites d'intégration sont z [0] et z [nz-1].

J'ai calculé cette intégrale en utilisant la règle de trapèze, et pour l'axe z de 5000 points cela prend 0,06 sec. Mon problème est que ce calcul se produit à peu près 300 * 30 * 20 fois (j'en ai 3 pour des boucles), et ce 0.06 sec se développe très rapidement à 3 heures de la simulation. Et tout le goulot d'étranglement de mon code est cette intégration (je peux évidemment accélérer en réduisant z, mais ce n'est pas le point.)

Je sais que les fonctions de bibliothèque sont généralement bien meilleures que celles écrites par l'utilisateur. Je sais aussi que je ne peux pas utiliser quelque chose de plus simple comme la règle de Simpson, parce que l'intégrand est fortement oscillatoire, et je veux éviter ma propre implémentation d'un algoritm numérique complexe.

GSL a besoin d'une fonction sous forme:

F = f (double x, void * params)

et je peux probablement utiliser l'intégration adaptative QAWO de GSL, mais comment puis-je faire ma fonction en forme qui transforme mon tableau en fonction?

Je pense que quelque chose:

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]; 
} 

Est-ce quelque chose comme cela possible? Je doute que l'algorithme numérique utilise la même différence spatiale que la fonction actuelle, devrais-je alors faire une interpolation de la fonction réelle ou quelque chose?

Y at-il quelque chose de mieux que 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)}; } 

utilisation:

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]; 
}); 

puis

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

et vous pouvez passer pfun et pvoid à travers.

Toutes nos excuses pour les fautes de frappe. L'idée est d'écrire un lambda qui fait ce que nous voulons. Puis lambda_to_pfunc l'enveloppe pour que nous puissions le passer comme un void* et un pointeur de fonction vers les API de style C.