precision loss with besselh, bessely for order n >16

[David Bateman]

David Bateman 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.
>
>   
You could even add a function like

function y = besselh (n, k, x)
  if (nargin < 2)
    error ("incorrect call to besselh");
  elseif (nargin == 2)
    x = k;
    k = 1;
  elseif (k != 1 && k != 2)
    error ("expecting K = 1 or 2")
  endif
 
  if (k == 1)
     y = bessel_Jn(n,x) + 1i * bessel_Yn(n,x);
  else
     y = bessel_Jn(n,x) - 1i * bessel_Yn(n,x);
  endif
endfunction

though I suspect it will have to be defined at the command line rather 
than a script to overload the oct-file containing besselh. Then you 
wouldn't need to modify your code.

Regards
David

-- 
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)