precision loss with besselh, bessely for order n >16

[David Bateman]

jessica.alfonsi at unipd.it wrote:
> Hello,
> I'm running some MatLab codes for Mie theory found in a book, which I
> ported very easily to Octave 3.0.1 on Ubuntu Jaunty 9.04. These codes rely
> on the computation of special Bessel functions, in particular I
> experienced problems (NaN, Inf + Infi) with besselh and bessely due to
> precision loss for order higher than 16. With lower orders the codes
> produced the correct results as with MatLab.
>
> See for example the following command in Octave:
>
> besselh(50,[0.1:0.1:2.0]')
>
> which gives an array with all elements
> Inf + Infi
>
> Note also that the same code on MatLab 2009a (now expired) trial run
> without issues for order up to 50. It would be very important for my
> research to be able to run these codes in Octave for orders of besselh as
> close as those obtained with MatLab. Is there any workaround to improve
> this limit ?
>
> Thank you for your attention
>
> Best
>
> Jessica
>   

Yes, the issue is that Octave uses the rather old AMOS library for its 
Bessel functions and this library overflows for large arguments, thus 
the infinite return value.. There has been discussion in the past of 
replacing this code with the code from GSL, but that iss a larger task 
than it sounds, as if  we use GSL for  the bessel functions it makes 
sense using it for a large number of other elementary functions.

Now is there a workaround? Well in fact yes there is by using the 
octave-forge package GSL. This package is a simplistic binding to most 
of the GSL functions. If you have this octave-forge package loaded, then 
a Hankel function of the first kind as you defined is given by

 bessel_Jn(50,[0.1:0.1:2.0]') + 1i * bessel_Yn(50,[0.1:0.1:2.0]')

for your example, and for a hankel function of the second kind it would be

 bessel_Jn(50,[0.1:0.1:2.0]') - 1i * bessel_Yn(50,[0.1:0.1:2.0]')

and this works fine.. If some one was willing to do the working of 
importing GSL into Octave I think that they would have support from the 
rest of the Octave community..

D.

-- 
David Bateman                                dbateman at dbateman.org
35 rue Gambetta                              +33 1 46 04 02 18 (Home)
92100 Boulogne-Billancourt FRANCE            +33 6 72 01 06 33 (Mob)