大变量的高阶贝塞尔函数计算答案

【问题标题】：high order bessel function computation with large variables大变量的高阶贝塞尔函数计算
【发布时间】：2015-12-05 05:04:06
【问题描述】：

我的工作涉及计算大变量值的高阶贝塞尔函数。在 MATLAB 中，这已毫无问题地完成。但是，为了扩大问题的规模，我已经调整为使用 MPI 编写 C++ 代码。当然，生成贝塞尔函数的步骤是通过调用一些库来完成的。为了具体说明问题，让我考虑一下这个非常具体的错误。

在matlab中，假设我想计算$J_46341(86840.0)$，并且

matlab 给我：besselj(46341,86840)=0.001309896212292

不过，调用一个简单的测试示例

gsl_sf_bessel_Jn_e 返回“错误：NaN”

我已经检查了订单 46340，matlab 和 gsl 在可接受的精度内返回相同的答案 0.00292895。 GSL 中的多一步会导致 NaN 错误，而 matlab 仍然保留了良好的准确数字答案。

我确实尝试使用递归关系来生成更高阶的值，从一个不那么小的阶，比如从 20000 及以上的阶开始，但是，这只会延迟 NaN 错误而不能完全解决问题。

将注意力转移到其他可用的软件库上，我尝试了 NAG，但令我完全失望的是，

nag_bessel_j_alpha (s18ekc) 的约束为 abs(nl)

，换句话说，它只能计算到 101 阶，这显然不符合我的研究兴趣。

所以，我的问题很简单：

是否有更可靠的库方法来获得高阶贝塞尔大 x 的函数值？

渐近地，贝塞尔函数接近 0，如果尾部接近下溢限制，我肯定可以将这些值设置为零。然而，NaN 问题似乎发生在强烈振荡曲线和渐近衰减尾部之间。

【问题讨论】：

为什么不自己计算函数呢？
使用贝塞尔函数的递归关系不会帮助我更进一步。事实上，与使用 GSL 库例程相比，随着订单的增加，它的分解速度要快得多。
@macduf 你肯定不是手动的意思。

标签： c++ matlab gsl bessel-functions

【解决方案1】：

问题解决了。感谢您的社区工作，您的知识和贡献让我感到惊讶！！！

请看这里， how to call fortran routines from C++?

https://mathoverflow.net/questions/225121/computation-of-high-order-bessel-function-at-large-variable-value

MATLAB、R、Python 和 JuliaLang/openspecfun 均以 Donald E. Amos 博士（桑迪亚国家实验室）的原始 fortran 源代码为基础，引用论文：

D. E. Amos, "A subroutine package for Bessel functions of a complex
argument and nonnegative order", Sandia National Laboratory Report,
SAND85-1018, May, 1985.
D. E. Amos, "A portable package for Bessel functions of a complex
argument and nonnegative order", Trans. Math. Software, 1986.

现在称为 ACM 收集的 Amos 算法 644。

http://dl.acm.org/citation.cfm?id=212078
http://dl.acm.org/citation.cfm?id=1268783
http://dl.acm.org/citation.cfm?id=98299

然而，托管在 netlib 上的源代码并非没有错误，而且可能不是最新的，

http://netlib.sandia.gov/master/index.html
http://netlib.sandia.gov/amos/

虽然 openspecfun 采用的版本运行良好，

https://github.com/JuliaLang/openspecfun

【讨论】：