This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
Laguerre Polynomial overflow...
- From: "Robert G. Brown" <rgb at phy dot duke dot edu>
- To: gsl-discuss at sources dot redhat dot com
- Date: Tue, 1 Oct 2002 14:40:48 -0400 (EDT)
- Subject: Laguerre Polynomial overflow...
Dear GSL listpersons:
I'm using gsl-1.0-1 prebuilt in Red Hat 7.2 (and have been using various
gsl routines from roughly 0.7-something on). Currently I'm testing an
ODE-based solution to the time-independent Schrodinger equation by
solving the O(3) simple harmonic oscillator problem, the radial
eigensolution to which is basically a laguerre polynomial times a few
other simple functions.
To compare the ODE solution to the exact solution, I've been evaluating
...
status = gsl_sf_laguerre_n_e(nsho, lsho+0.5 , rsho*rsho, &lgp);
if(status != GSL_SUCCESS) return(0.0);
ret = sqrt(2.0*fac(nsho)/fac(lsho + nsho + 1.0))*pow(rsho,(double)lsho)*
lgp.val*exp(-rsho*rsho*0.5)/norm;
fprintf(stderr,"sho3d( %2d, %2d, %2d, %.5e) = %.5e\n",nsho,lsho,msho,rsho,ret!
fprintf(stderr,"gsl_sf_laguerre_n( %2d, %2d, %.5e) = %.5e\n",nsho,lsho,rsho,lgp.val);
return(ret);
which works fine for nsho = 0 or 1, but when I hit nsho > 2, the routine
runs for some values of lsho and rsho but then dies on an overflow
beyond a certain point:
sho3d( 2, 0, 0, 2.21000e+00) = 7.38523e-02
gsl_sf_laguerre_n( 2, 0, 4.88410e+00) = 1.59197e+00
sho3d( 2, 0, 0, 2.22000e+00) = 7.70712e-02
gsl_sf_laguerre_n( 2, 0, 4.92840e+00) = 1.69856e+00
sho3d( 2, 0, 0, 2.23000e+00) = 8.02146e-02
gsl_sf_laguerre_n( 2, 0, 4.97290e+00) = 1.80762e+00
gsl: laguerre.c:96: ERROR: overflow
Abort
The interesting thing is that (as you can see) the solution is extremely
pedestrian here -- n is two, \ell is 0, and the argument is around 5.
The returned value is order unity. I can see no good reason for an
overflow to be generated externally, and presume that it is occuring
during the operation of e.g. the recursion relation that generates the
result internally. I also seem to generate no GSL error condition --
the code crashes even though I try to trap the overflow and continue.
For what it is worth, the solution is dead on up to that point,
validated somewhat backwards by the nearly exact coincidence of this
solution with the one I'm obtaining by explicit ODE solution with e.g.
rk45 (also from gsl elsewhere in the code).
If nobody on this list has encountered this or can suggest any fix
beyond going into the source and figuring out where and why the overflow
occurs and fixing it, I'm happy enough to do the latter, but given the
hassle of building from source and the possibility that it has already
been fixed in a more recent version, I thought I'd ask first.
rgb
Robert G. Brown http://www.phy.duke.edu/~rgb/
Duke University Dept. of Physics, Box 90305
Durham, N.C. 27708-0305
Phone: 1-919-660-2567 Fax: 919-660-2525 email:rgb@phy.duke.edu