[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)'.

My comments:

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


Feedbacks & comments are welcome.

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.


More information about the SciPy-user mailing list