precision loss with besselh, bessely for order n >16

[Jaroslav Hajek]

On Sat, Jun 6, 2009 at 10:00 PM, David Bateman<dbateman at dbateman.org> wrote:
> 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.
>

boost::math is another option. The slight problem with GSL is that it
is double precision only; for bessel functions that is probably no
problem but for a number of possible other applications it matters. It
is unfortunate, IMHO. Boost is much more generic and also more
actively developed.


-- 
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz