[SciPy-User] Strange behavior of sph_jn
ZIEGLER Yann (ETU EOT)
Wed Apr 11 03:36:34 CDT 2012
> On Tue, Apr 10, 2012 at 8:49 AM, ZIEGLER Yann (ETU EOT)
> <email@example.com> 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).
> 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.
You're probably right, it is coherent with what I read in some docs :
for big values it is efficient enough to use the asymptotic approximation
of these functions... but I'm not convinced anyway returning inf or NaN
is the behavior one could expect by calling sph_jn which must be defined
for all real or complex numbers. So, you're right: I was mistaken in
speaking about a bug, but for me it's a pitfall that should be removed
from scipy... at least for people like me who don't know all the ins and
outs of special functions computation.
> 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.
In fact, I needed the derivatives but, fortunately, there are some useful
recurrence relations, such as :
Thanks for your comment!
PS : Sorry, I'm not familiar with mailing lists and I've created a new
message when I answered previously. I hope it will be all right now!
More information about the SciPy-User