[SciPy-user] find roots for spherical bessel functions...
fred
fredantispam at free.fr
Wed Aug 16 17:04:03 CDT 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