[SciPy-User] Strange behavior of sph_jn

[Aronne Merrelli]

On Tue, Apr 10, 2012 at 8:49 AM, ZIEGLER Yann (ETU EOT)
<yann.ziegler at etu.unistra.fr> wrote:
> Hi,
>
> It seems that the bug I have with sph_jn for some values is due to scipy
> itself. I finded someone having the same kind of trouble on
> projects.scipy.org/scipy (the website seems to have some technical problems
> at this time).


Hi,

I believe that the SciPy implementation is from this reference (it is
mentioned in the header file of the fortran file in SciPy):

"Computation of Special Functions", 1996, John Wiley & Sons, Inc.
Shanjie Zhang and Jianming Jin

The fortran subroutine called from the python level sph_jn is 'SPHJ';
I don't understand the implementation very well (it is very cryptic to
my untrained eye, and I don't have access to the book), but at the
core of it is an iterative algorithm that just breaks down for the
large arguments that you were testing. I think it is running into the
max floating point exponent.

So, I don't think there is any SciPy bug exactly, it is just this
particular numerical implementation of the spherical bessel function
doesn't work for those large arguments.  Since you don't need the
derivatives, then your own simple function relating it to the bessel
function should work fine.

Aronne