find roots for spherical bessel functions...
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.
On 16/08/06, fred <fredantispam@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? 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. 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. 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.
A. M. Archibald a écrit :
On 16/08/06, fred <fredantispam@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. Ok. (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. idem.
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, etc). So I was thinking it will be so easy with scipy. Thanks. -- Fred.
On 17/08/06, fred <fredantispam@free.fr> wrote:
PS : with fsolve's octave routine, I have no such problem (find out x0, etc). So I was thinking it will be so easy with scipy.
I am mystified by this. What initial values are you providing to fsolve, and how are you sure that they're giving you the right roots? (Si la langue est une problème, vous pouvez me contacter en Français hors la liste) A. M. Archibald
A. M. Archibald a écrit :
On 16/08/06, fred <fredantispam@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.
fred a écrit :
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) In fact, my second idea was to interpolate intevals [a,b] (which should contain zero) for a given Jn, and then for Jn for all n. Both should be linear interpolation.
But I can't get it working fine because linear interpolation does not match very well [a,b] values :-( So... I don't know. Nevertheless, I have at least one solution to find out zeros of spherical Bessel, which works fine. Thanks to all. Cheers, -- Fred.
participants (2)
-
A. M. Archibald -
fred