1

私の仕事には、大きな変数値での高次ベッセル関数の計算が含まれます。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 でもう 1 ステップ実行すると NaN エラーが発生しますが、matlab では正確な数値の回答が保持されます。

再帰関係を使用して、それほど小さくないオーダー、たとえば 20000 以上のオーダーからより高いオーダーの値を生成しようとしましたが、これは問題を完全に解決せずに NaN エラーを遅らせるだけです。

他の利用可能なソフトウェア ライブラリに注意を向けて、NAG を試してみましたが、まったくがっかりしました。

nag_bessel_j_alpha (s18ekc) には abs(nl)<=101 の制約があります

つまり、101 のオーダーまでしか計算できず、明らかに私の研究対象ではありません。

だから、私の質問はかなり簡単です：

大きなxの高次ベッセル関数値を取得するためのより信頼できるライブラリアプローチはありますか?

漸近的に、ベッセル関数は 0 に近づきます。テールがアンダーフローの限界に近づいている場合、これらの値を確実にゼロに設定できます。ただし、NaN の問題は、強く振動する曲線と漸近的に減衰するテールの間で発生するようです。

4

1 に答える 1

1

問題が解決しました。コミュニティ活動に感謝します。皆さんの知識と貢献に本当に驚かされました!!!

こちらを参照してください 。C++ から fortran ルーチンを呼び出す方法は?

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 Algorithm 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
于 2015-12-04T21:04:06.410 に答える