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

fred fredantispam at free.fr
Mon Aug 21 18:08:06 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, so I have found a first solution, I call "recursive", because the
intervals [a,b] of Jn depend of zeros
of Jn-1 (thanks to you and A&S ;-)

This is useful when you want to find out _all_ zeros for 1<m,n<100 for
example.
But when you only want the nth zero of the mth Jn, I think I can find an
easier algorithm (I'm currently working
on it), although the recursive method is quite fast (a couple of seconds
for m=n=100).

I post here an example of using my code, with pylab, in order to display
zeros of Jn and (rJn)'.

1) For m >= 30 & n >= 15, one can see that the curves begin to be quite
"noisy".
I don't understand/know why. Any idea ?

2) I have defined Jn as spherical Bessel, because:
a) You can't display it in an easy way, knowing that sph_jn() does not
accept array as arg ;
b) I can't get solving sph_jn = 0, since brentq/fsolve can only accept
n as a second argument.
But in sph_jn(), n is the first argument (this is what I have understood).

from scipy import arange, pi, sqrt, zeros
from scipy.special import jv, jvp
from scipy.optimize import brentq
from sys import argv
from pylab import *

def Jn(r,n):
return (sqrt(pi/(2*r))*jv(n+0.5,r))
def Jn_zeros(n,nt):
zerosj = zeros((n+1, nt))
zerosj[0] = arange(1,nt+1)*pi
points = arange(1,nt+n+1)*pi
racines = zeros(nt+n)
for i in range(1,n+1):
for j in range(nt+n-i):
foo = brentq(Jn, points[j], points[j+1], (i,))
racines[j] = foo
points = racines
zerosj[i][:nt] = racines[:nt]
return (zerosj)

def rJnp(r,n):
return (0.5*sqrt(pi/(2*r))*jv(n+0.5,r) + sqrt(pi*r/2)*jvp(n+0.5,r))
def rJnp_zeros(n,nt):
zerosj = zeros((n+1, nt))
zerosj[0] = (2.*arange(1,nt+1)-1)*pi/2
points = (2.*arange(1,nt+n+1)-1)*pi/2
racines = zeros(nt+n)
for i in range(1,n+1):
for j in range(nt+n-i):
foo = brentq(rJnp, points[j], points[j+1], (i,))
racines[j] = foo
points = racines
zerosj[i][:nt] = racines[:nt]
return (zerosj)

n = int(argv[1])
nt = int(argv[2])

dr = 0.01

jnz = Jn_zeros(n,nt)[n]
r1 = arange(0,jnz[len(jnz)-1],dr)

jnzp = rJnp_zeros(n,nt)[n]
r2 = arange(0,jnzp[len(jnzp)-1],dr)

grid(True)
plot(r1,Jn(r1,n),'b', jnz,zeros(len(jnz)),'bo', r2,rJnp(r2,n),'r',
jnzp,zeros(len(jnzp)),'rd')
legend((r'\$j_{'+str(n)+'}(r)\$','', r'\$(rj_{'+str(n)+'}(r))\'\$',''))
show()

Cheers,

--
Fred.