Каков точный эквивалент Matlab besselk (x, y, 1) в С++?

Я пробовал boost::math::cyl_bessel_k(x,y) * exp(y). В большинстве случаев это равно масштабируемому Matlab besselk(x,y,1). Но в некоторых случаях (например, x=1, y=2000), когда и besselk(x,y)=0, и boost::math::cyl_bessel_k(x,y)=0, масштабированная версия Matlab besselk(x,y,1) дает мне разные значения, варьирующиеся вокруг 10^-3. Но boost::math::cyl_bessel_k(x,y) * exp(y) возвращает -nan. Я хотел бы найти эквивалентное утверждение для Matlab besselk(x,y,1). Как я могу справиться с этим?


person taha    schedule 07.07.2014    source источник


Ответы (1)


Я не вижу в Boost ничего, что делало бы то, что вам нужно (хотя вы могли бы реализовать это самостоятельно, используя функции более низкого уровня). Как вы узнали, масштабированные функции Бесселя не вычисляются простым умножением exp(z). GSL, похоже, включает эту функциональность, например, gsl_sf_bessel_Knu_scaled. Для «точного эквивалента» вы можете посмотреть статью и код Амоса, например, CBESK. И Matlab, и Octave используют эту реализацию. Однако код написан на Фортране, поэтому вам нужно будет его перевести или поместить в оболочку (этот проект, кажется, сделал это, поэтому он может быть полезен — есть и другие).

Вы также можете использовать Matlab Coder и codegen для вывода чего-либо в виде хорошо.

person horchler    schedule 07.07.2014
comment
gsl_sf_bessel_Knu_scaled работает правильно, за исключением того, что первый параметр отрицательный. Когда первый параметр отрицательный, я получаю следующее сообщение об ошибке: gsl: bessel_Knu.c:42: ERROR: domain error Default GSL error handler invoked. Aborted (core dumped). Кроме того, я поставил эту строку для обработки ошибок gsl_set_error_handler_off();, но на этот раз вывод gsl_sf_bessel_Knu_scaled снова nan. Как я могу справиться с этим? - person taha; 08.07.2014
comment
Я обнаружил, что если первый параметр отрицательный, результат будет таким же, как и у gsl_sf_bessel_Knu_scaled(abs(firstparam),y). Так что это работает. Большое спасибо! - person taha; 09.07.2014