2

I have tried boost::math::cyl_bessel_k(x,y) * exp(y). In most cases, this is equal to Matlab's scaled besselk(x,y,1). But in some cases (e.g., x=1, y=2000) when both besselk(x,y)=0 and boost::math::cyl_bessel_k(x,y)=0, Matlab's scaled version besselk(x,y,1) gives me different values varies around 10^-3. But boost::math::cyl_bessel_k(x,y) * exp(y) returns -nan. I'd like to find an equivalent statement for Matlab's besselk(x,y,1). How can I handle this?

horchler
  • 18,384
  • 4
  • 37
  • 73
taha
  • 172
  • 8

1 Answers1

1

I'm not seeing anything in Boost that does what you need (though you might be able to implement it yourself by using lower-level functions). As you've found out, the scaled Bessel functions are not computed simply by multiplying exp(z). The GSL appears to have incorporated this functionality, e.g., gsl_sf_bessel_Knu_scaled. For an "exact equivalent," you might look at the paper and code by Amos, e.g., CBESK. Both Matlab and Octave appear to use this implementation. However, the code is written in Fortran, so you'd need to translate it or put a wrapper around it (this project appears to have done that so it might be useful – there are others out there too).

You may also be able to use Matlab's Coder and codegen to output something as well.

horchler
  • 18,384
  • 4
  • 37
  • 73
  • `gsl_sf_bessel_Knu_scaled` is working correctly except for first parameter is negative. When first parameter is negative, I get an error as follow: `gsl: bessel_Knu.c:42: ERROR: domain error Default GSL error handler invoked. Aborted (core dumped)`. Also, I put this line for error handling `gsl_set_error_handler_off();`, but this time the output of `gsl_sf_bessel_Knu_scaled` is again `nan`. How can I handle this? – taha Jul 08 '14 at 02:09
  • I found that if first parameter is negative, result will be the same with `gsl_sf_bessel_Knu_scaled(abs(firstparam),y)`. So it works. Thank you very much! – taha Jul 09 '14 at 01:12