вычисление функции Бесселя высокого порядка с большими переменными

Моя работа связана с вычислением функции Бесселя высокого порядка при большом значении переменной. В MATLAB это было сделано без проблем. Однако, чтобы масштабировать проблему, я настроился на написание кода на C++ с помощью MPI. Конечно, этап создания функции Бесселя выполняется путем вызова некоторых библиотек. Чтобы конкретизировать проблему, позвольте мне рассмотреть эту очень конкретную ошибку.

Предположим, в 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

, другими словами, он может вычислять только до порядка 101, и это явно не в моих интересах.

Итак, мой вопрос довольно прост:

Есть ли более надежный библиотечный подход для получения значения функции Бесселя высокого порядка для больших x?

Асимптотически функция Бесселя приближается к 0, я могу с уверенностью установить эти значения равными нулю, если хвост приближается к пределу потери значимости. Однако проблема NaN, по-видимому, возникает где-то между сильно осциллирующей кривой и асимптотически затухающим хвостом.


person Franklin Dong    schedule 02.12.2015    source источник
comment
почему бы не вычислить функцию самостоятельно?   -  person marsei    schedule 03.12.2015
comment
использование рекурсивного отношения функции Бесселя не помогло бы мне двигаться дальше. На самом деле, с увеличением количества заказов он выходит из строя гораздо быстрее, чем при использовании библиотечных подпрограмм GSL.   -  person Franklin Dong    schedule 03.12.2015
comment
@macduf, конечно, ты не имеешь в виду вручную.   -  person erip    schedule 03.12.2015


Ответы (1)


Задача решена. Спасибо за работу сообщества, и это действительно поразило меня вашими знаниями и вкладом!!!

См. здесь, , как вызывать подпрограммы fortran из C++?

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

MATLAB, R, Python и JuliaLang/openspecfun основаны на оригинальном исходном коде фортрана доктора Дональда Э. Амоса (национальная лаборатория Сандии), цитируемая статья:

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.

Теперь известный как Алгоритм Амоса 644, собранный ACM.

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
person Franklin Dong    schedule 04.12.2015