[SciPy-user] find roots for spherical bessel functions...

fred fredantispam at free.fr
Thu Aug 17 02:26:15 CDT 2006

A. M. Archibald a écrit :
> On 16/08/06, fred <fredantispam at free.fr> wrote:
>>I'm trying to find out the first n roots of the first m spherical bessel
>>functions Jn(r)
>>(and for the derivative of (r*Jn(r))), for 1<m,n<100.
> This is only about 20,000 values, so I recommend precomputing a table
> (which will be slow, but only necessary once). Probably the tables
> exist somewhere, but whether they're freely available (can you
> copyright those numbers?) in digital form... perhaps you should
> publish them (and the code, so we can check!) when you succeed?
Ok, only I succeed ;-)

> Looking at Abramowitz and Stegun (
> http://www.math.sfu.ca/~cbm/aands/page_370.htm ), you can see that the
> zeros of the m+1st interlace with the zeros of the mth. So if you can
> find the first 200 zeros of the first spherical Bessel function, all
> the rest fall immediately using (say) brentq.
(PS: A & S online ? Great ! :-)

> As for the first, well, you could probably eyeball it, see what the
> shortest space between zeros is for the first two hundred, and just
> take values at less than half that spacing. That should get you
> brackets for all of them. (You could also use some clever estimates of
> the sizes of first and second derivatives to make sure you didn't miss
> any roots, if you wanted to do it more automatically.) Sturm theory
> would probably also let you bracket roots between roots of
> integer-weight Bessel functions, if you needed a random-access
> function.
Ok. I will look that.
> Your second function, r*Jn'(r)+Jn(r) will be more painful, but some
> clever reasoning about the signs and zeros of Jn and Jn' (and Jn'',
> which you can get from the differential equation) should let you
> bracket the roots without undue difficulty.

> Good luck, and please do make the table available online,
> A. M. Archibald
> P.S. scipy.special can compute spherical Bessel functions directly.
Yes, I know that.
The problem I have with this function (as associated Legendre polys, etc),
is that its argument can only be a scalar. And I need to pass it arrays.
So I have to define spherical Bessel func with jv(), which accepts
arrays as arg.

PS : with fsolve's octave routine, I have no such problem (find out x0,
So I was thinking it will be so easy with scipy.



More information about the SciPy-user mailing list