Integration of Bessel functions in C++/Fortran

Question

How can I integrate an equation including bessel functions numerically from "0" to "infinity" in Fortran or/and C? I did in matlab, but it's not true for larger inputs and after a specific values , the bessel functions give completely wrong results(there is a restriction in Matlab)

Answer 1

There's a large number of analytic results for various integrals of the Bessel functions (see DLMF, Sect. 10.22), including definite integrals over precisely this range. You'd be much better off, and almost certainly faster and more accurate, trying hard to recast your expression into something that's integrable and using an exact result.

Answer 2

Last time I had to do with such things, it was state of the art to do simple integration of the intervals defined by the zero crossings. That is in most cases relatively stable and if the integrand is approaching zero reasonable fast easy to do.

As a starting point for playing around I´ve included a bit of code. Of course you need to work on the convergence detection and error checking. This is no production code but I thought maybe it provides a starting point for you. Its using gsl.

On my iMac this code takes about 2 µs per iteration. It will not become faster by including a hardcoded table for the intervals.

I hope this is of some use for you.

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

Answer 3

You can pretty much google and find lots of Bessel functions implemented in C already.

http://www.atnf.csiro.au/computing/software/gipsy/sub/bessel.c

http://jean-pierre.moreau.pagesperso-orange.fr/c_bessel.html

https://msdn.microsoft.com/en-us/library/h7zkk1bz.aspx

In the end, these use the built in types and will be limited to the ranges they can represent (just as MATLAB is). At best, expect 15 digits of precision using double precision floating point representation. So, for large numbers, they will appear to be rounded. eg. 1237846464123450000000000.00000

And, of course, others on Stack Overflow have looked into it.

C++ Bessel function for complex numbers