Conseils sur la méthode pour intégrer les fonctions de Bessel en C++

Question

je dois intégrer une fonction oscillatoire de ce genre:

où j'ai une fonction Bessel qui oscille, alors que F est pas très oscillatoire. Je suis à la recherche de la méthode la plus précise/précise pour le faire en C++. Espérons que cela devrait être une implémentation déjà existante que je peux utiliser par exemple par une bibliothèque etc ...

La bibliothèque GSL pourrait être une option mais s'il vous plaît même dans ce cas pourriez-vous me recommander laquelle des nombreuses routines avaialble est le le plus utile pour moi?

EDIT: Je suis au courant de la présence de cette possibly duplicate question

mais je ne vois pas de réponse claire là. En outre, une bibliothèque Fortran ne serait pas bonne pour moi, à moins qu'il y ait une sorte de wrapper pour cela.

Answer 1

Dernière fois que j'ai eu à faire avec de telles choses, c'était l'état de l'art de faire une intégration simple des intervalles définis par les passages par zéro. C'est dans la plupart des cas relativement stable et si l'intégrand se rapproche de zéro rapide, facile à faire.

Comme point de départ pour jouer, j'ai inclus un peu de code. Bien sûr, vous devez travailler sur la détection de la convergence et la vérification des erreurs. Ce n'est pas un code de production mais je pensais que cela pourrait constituer un point de départ pour vous. C'est en utilisant gsl.

Sur mon iMac, ce code prend environ 2 μs par itération. Il ne deviendra pas plus rapide en incluant une table codée en dur pour les intervalles. J'espère que cela vous servira à quelque chose.

#include <iostream> 
#include <vector> 
#include <gsl/gsl_sf_bessel.h> 
#include <gsl/gsl_integration.h> 
#include <gsl/gsl_sf.h> 


double f (double x, void * params) { 
    double y = 1.0/(1.0 + x) * gsl_sf_bessel_J0 (x); 
    return y; 
} 

int main(int argc, const char * argv[]) { 

    double sum = 0; 
    double delta = 0.00001; 
    int max_steps = 1000; 
    gsl_integration_workspace * w = gsl_integration_workspace_alloc (max_steps); 

    gsl_function F; 
    F.function = &f; 
    F.params = 0; 

    double result, error; 
    double a,b; 
    for(int n=0; n < max_steps; n++) 
    { 
     if(n==0) 
     { 
      a = 0.0; 
      b = gsl_sf_bessel_zero_J0(1); 
     } 
     else 
     { 
      a = n; 
      b = gsl_sf_bessel_zero_J0(n+1); 
     } 
     gsl_integration_qag (&F, // function 
           besselj0_intervals[n], // from 
           besselj0_intervals[n+1], // to 
           0, // eps absolute 
           1e-4,// eps relative 
           max_steps, 
           GSL_INTEG_GAUSS15, 
           w, 
           &result, 
           &error); 
     sum += result; 

     std::cout << n << " " << result << " " << sum << "\n"; 

     if(abs(result) < delta) 
      break; 
    } 
    return 0; 
}