Dear GSL developers,
I found a bug in the function gsl_sf_coulomb_wave_F_array which uses the
recurrence relation to obtain F for different l's.
To be short, when eta is large, and l is at least moderately large,
inaccuracy occurs in the implementation of the recurrence such that F
has unphysical spikes at certain x for smaller l's.
My observation is that F is determined directly for the largest l, and F
for smaller l's are determined from the recurrence relation.
When l is large, the number of recurrence is large to obtain F for small
l's.
Therefore, errors in F's get accumulated as we proceed to get F for
small l's.
Attached please find a test program.
The program uses two methods to obtain F for different l's and x's.
The first method obtains F directly, and the second method obtains F
from recurrence relation.
Please compile it using "g++ testCoulomb.cpp -lgsl -lgslcblas".
Now run the program using "./a.out > resultCoulomb".
We can now take a look at the output of the program using gnuplot.
After launching gnuplot, just type plot "resultCoulomb" using 1:2
every ::2::2000 to see the output for l=0 at different x's using
the first method.
And type plot "resultCoulomb" using 1:2 every ::2002::4000 to see
the output for l=0 at different x's using the second method.
Now you can find, F is fine using the first method, but F has unphysical
spikes.
You can see similar behaviors for other l's by "using 1:3", "using 1:4"
etc in gnuplot.
I can of course always stick to the first method, but the second method
is much faster.
So it would be great if you could fix this problem.
Thanks a lot!
--
Best regards,
Hongcheng Ni
Max-Planck-Institut für Physik komplexer Systeme
Nöthnitzer Str. 38, 01187 Dresden, Germany
#include <stdio.h>
#include <gsl/gsl_sf_coulomb.h>
int main( int argc, char *argv[] )
{
int l = 20;
double eta = -100.;
/* Method 1: obtain F directly */
{
printf("Method 1: obtain F directly\n");
gsl_sf_result F, Fp, G, Gp;
double exp_F, exp_G;
for(double rho=0.01; rho<20.; rho+=0.01)
{
printf("%f\t",rho);
for(int i=0; i<=l; i++)
{
gsl_sf_coulomb_wave_FG_e(eta, rho, i, 0, &F, &Fp, &G, &Gp, &exp_F, &exp_G);
printf("%e\t",F.val);
}
printf("\n");
}
}
/* Method 2: obtain F from recurrence */
{
printf("Method 2: obtain F from recurrence\n");
double F[1+l], exp_F[1+l];
for(double rho=0.01; rho<20.; rho+=0.01)
{
gsl_sf_coulomb_wave_F_array(0., l, eta, rho, F, exp_F);
printf("%f\t",rho);
for(int i=0; i<=l; i++)
{
printf("%e\t",F[i]);
}
printf("\n");
}
}
}
