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

fred fredantispam at free.fr
Wed Aug 16 18:04:03 EDT 2006


Hi all,

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.

So, I decided to use optimize func, such as fsolve() or brentq().

My problem is to find out x0, the "seed" to start to find the root, for
fsolve(), or
the bounds [a,b] for brentq().

I tried x0 = alpha + beta*n, where alpha \approx 3.5 & beta \approx 1,
but there is always an order for which it fails : for brentq(), you must
have f(a)*f(b) < 0.

Is there anyone who has any idea to estimate x0 or [a,b] for the nth
root of the mth spherical Bessel function
(same question for (rJn(r))')

Here's my code.

def Jn(r,n):
  return (sqrt(pi/(2*r))*jv(n+0.5,r))
def Jn_zeros(n,nt):
  r0 = 3.5 + n
  val = zeros(nt)
  i = 0

  while (i < nt):
    foo = brentq(Jn, r0, r0+pi, (n,))
    val[i] = foo
    i = i + 1
    r0 = foo + pi
  return (val)

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):
  r0 = 3.5 + n - pi/2
  val = zeros(nt)
  i = 0

  while (i < nt):
    foo = brentq(rJnp, r0, r0+pi, (n,))
    val[i] = foo
    i = i + 1
    r0 = foo + pi
  return (val)

Thanks in advance.

Cheers,

-- 
Fred.



More information about the SciPy-User mailing list