在 C++/Fortran 中集成贝塞尔函数

Question

如何在 Fortran or/and C 中对包含从“0”到 "infinity" 的贝塞尔函数的方程进行积分？我在 matlab 中做过，但对于较大的输入和特定值后，贝塞尔函数给出了完全错误的结果（Matlab 中有限制）

Answer 1

您几乎可以 google 找到许多已经用 C 语言实现的 Bessel 函数。

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

最后，这些使用内置类型并将限制在它们可以表示的范围内（就像 MATLAB 一样）。充其量，使用双精度浮点表示法期望 15 位精度。因此，对于大数字，它们看起来是四舍五入的。例如。 1237846464123450000000000.00000

当然，Stack Overflow 上的其他人也对此进行了调查。

C++ Bessel function for complex numbers

Answer 2

贝塞尔函数的各种积分有大量解析结果（参见DLMF, Sect. 10.22），包括恰好在此范围内的定积分。你会过得更好，而且几乎肯定会更快更准确，努力将你的表达式重铸成可积分的东西并使用精确的结果。

Answer 3

上次我不得不处理此类事情时，最先进的技术是对零交叉点定义的区间进行简单积分。在大多数情况下，这是相对 stable 并且如果被积函数接近零，那么很容易做到。

作为尝试的起点，我添加了一些代码。当然，您需要进行收敛检测和错误检查。这不是生产代码，但我认为它可能为您提供了一个起点。它使用 gsl。

在我的 iMac 上，这段代码每次迭代大约需要 2 微秒。通过为间隔包含硬编码 table，它不会变得更快。

希望对您有所帮助。

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

在 C++/Fortran 中集成贝塞尔函数

Integration of Bessel functions in C++/Fortran

c

fortran

fortran90