Laguerre Polynomial overflow...

Laguerre Polynomial overflow...

[Robert G. Brown]

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

Re: Laguerre Polynomial overflow...

[Gert Van den Eynde]

[-- Attachment #1: Type: text/plain, Size: 197 bytes --]


Dear Robert,

I cannot reproduce this using gsl-1.2 (my own compilation) on a SuSE 8.0 linux machine. Could you try to compile and run the attached program and see if that fails too?

Bye,

Gert


[-- Attachment #2: lag.c --]
[-- Type: text/plain, Size: 729 bytes --]

#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf_laguerre.h>

int main () {

    double x = 0.0, step = 0.2;
    int    k = 0, l = 0;

    gsl_sf_result res;

    for (k = 0; k < 4; k ++) {
        for (x = 4.0; x < 6.1; x += step)  {
            gsl_sf_laguerre_n_e(k,0.5,x,&res);
            printf("L_%d(%.1f) = %.8f        error = %.8f\n",k,x,res.val,res.err);
        }
        printf("\n");
    }


    gsl_sf_laguerre_n_e(2,0.5,4.97290e+00,&res);
    printf("L_%d(%.5f) = %.8f        error = %.8f\n",2,4.97290e+00,res.val,res.err);

    gsl_sf_laguerre_n_e(3,0.5,4.97290e+00,&res);
    printf("L_%d(%.5f) = %.8f        error = %.8f\n",3,4.97290e+00,res.val,res.err);

    return 0;
}

Re: Laguerre Polynomial overflow...

[Robert G. Brown]

On Wed, 2 Oct 2002, Gert Van den Eynde wrote:

> 
> Dear Robert,

> I cannot reproduce this using gsl-1.2 (my own compilation) on a SuSE
> 8.0 linux machine. Could you try to compile and run the attached program
> and see if that fails too?

Before starting to mess with the internal code I upgraded to the gsl 1.1
I had available in the RH 7.3 RPM set.  It already worked at that point.
Obviously little brownies were at work removing odd convergence/overflow
bugs in between...;-)

The moral, I suppose, is that I should make sure I'm using the latest
version before filing a bug report.  I'm a Bad Dog...

I appreciate it, though, and apologize for wasting your time.

I'll definitely compile and work through 1.2 before worrying about my
NEXT gsl problem -- the lack of (visible) adaptive stepsizes in the gsl
ode solvers.  Thus far I've had to solve on a very fine grid to get
decent accuracy (bound state eigensolutions integrated out from the
origin spend a bit of time crossing the axis where they oscillate and
eventually become stiff), but from the manual it looks like some
adaptive stepsize routines have already made their way into the later
revisions.  I can't wait to try them out.

BTW, the GSL is a fabulous toolset and I don't mean to sound critical at
all.  NR sucks (yes, I've read the rant, and then there is the
strangulation noose -- I mean "copyright").  I'd like to contribute, if
I ever have the chance.  For example, I have (somewhere, I'd have to dig
it out) a C translation of Forsythe, Malcom and Moler's "QUANC8" (8
panel Newton-Cotes quadrature) if that is of interest.  I'm not really a
computer scientist, but I've done a lot of numerical computations over
the last decade or two.  The other thing I might be able to contribute
is at least simple genetic optimization algorithms and (if I could
figure out how to frame them) neural network algorithms, as this is one
of the things I do.

   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

Re: Laguerre Polynomial overflow...

[Brian Gough]

Robert G. Brown writes:
 > The moral, I suppose, is that I should make sure I'm using the latest
 > version before filing a bug report. 

The GNU Project always provides a NEWS file for each package, which
lists the updates in each version -- check there first, before
installing anything.

Brian

Re: Laguerre Polynomial overflow...

[Brian Gough]

Robert G. Brown writes:
 > BTW, the GSL is a fabulous toolset and I don't mean to sound critical at
 > all.  NR sucks (yes, I've read the rant, and then there is the
 > strangulation noose -- I mean "copyright").  I'd like to contribute, if
 > I ever have the chance.  For example, I have (somewhere, I'd have to dig
 > it out) a C translation of Forsythe, Malcom and Moler's "QUANC8" (8
 > panel Newton-Cotes quadrature) if that is of interest.  I'm not really a
 > computer scientist, but I've done a lot of numerical computations over
 > the last decade or two.  The other thing I might be able to contribute
 > is at least simple genetic optimization algorithms and (if I could
 > figure out how to frame them) neural network algorithms, as this is one
 > of the things I do.

Thanks.  Information for developers is available in the TODO and
HACKING files in gsl cvs at http://sources.redhat.com/gsl

The design document doc/gsl-design.texi at the same location is also
useful reading.

Brian