[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [Bug-gsl] bug in gsl_sf_coulomb_wave_FG_e (gsl 1.14 compiled with gc

From: Alexey A. Illarionov
Subject: Re: [Bug-gsl] bug in gsl_sf_coulomb_wave_FG_e (gsl 1.14 compiled with gcc 4.5.0)
Date: Fri, 27 Aug 2010 13:36:08 -0400
User-agent: Mozilla/5.0 (X11; U; Linux i686; en-US; rv: Gecko/20100624 Lanikai/3.1

On 08/27/2010 06:47 AM, Brian Gough wrote:
At Wed, 25 Aug 2010 23:41:14 -0400,
Alexey A. Illarionov wrote:
It seems that I find a bug in gsl_sf_coulomb_wave_FG_e. For large values
of lambda and small values of x, with eta = 0 it produces nan values
without raising of the flag.

Here the sample cases (see attachment with example)
   gsl_sf_coulomb_wave_FG_e(0.,1.2693881947287221e-07, 37, 1, ...)
   gsl_sf_coulomb_wave_FG_e(0.,5.911135077240691e-07, 39, 1, ...)
   gsl_sf_coulomb_wave_FG_e(0.,6.118185507743884e-07, 40, 1, ...)
and lots of others.

Thanks for the bug report. I'm adding a test case for this.  Can you
confirm if the following values are correct (I'm not working with
these functions often):

   lam_F = 37.0;
   eta = 0.0;
   x = 1.2693881947287221e-07;
   k_G = 1;
   gsl_sf_coulomb_wave_FG_e(eta, x, lam_F, k_G,&F,&Fp,&G,&Gp,&Fe,&Ge);
   s = 0;
   message_buff[0] = 0;
   s += test_sf_check_result(message_buff,  F,  
-2.020327525317343313380459069e299, TEST_TOL3);
   s += test_sf_check_result(message_buff, Fp, 
5.729672862364037942289798904e307, TEST_TOL3);
   s += test_sf_check_result(message_buff,  G, 
4.397355593129492554472984282e299, TEST_TOL3);
   s += test_sf_check_result(message_buff, Gp,  - 
1.247095270068213211088454516e308, TEST_TOL3);

I've computed them using Pari as

   c(L,n) = { 2^L * exp(-Pi*n/2) * abs(gamma(L+1+I*n)) / gamma(2*L+2)}
   cofac(L,n,r) = { I * exp(-I*r)*r^(-L)/(gamma(2*L+2)*c(L,n)) }
   FplusIG(L,n,r) = { cofac(L,n,r) * 
intnum(t=0,[[1],1],exp(-t)*t^(L-I*n)*(t+2*I*r)^(L+I*n)) }

   FplusIG(37,0,1.2693881947287221e-07+x) # for F and F' (real part)
   FplusIG(36,0,1.2693881947287221e-07+x) # for G and G' (imaginary part)

Hi Brian,

It looks like that these values are incorrect since both F and Fp should be very small (correct me if I'm wrong since I'm not familiar with Pari and gsl test system). If you need verify the values than its better to use bessel function series.

For eta == 0. and integer lambda coulomb function are just bessel-ricatti functions, i.e. spherical_bessel functions times argument
F_l(x) = j_l (x) * x
G_l(x) = n_l (x) * x
F_l(x) = 1/(2l+1)!! x^(l+1) [1 - x^2/(2l+3)/2 + ...]
G_l(x) = (2l+1)!!/(2l+1) x^(-l) [1 + x^2/(2l-1)/2 + ...]
I think Abramowitz, Stegun should have the complete expressions.

Best regards, Alexey Alexandrovich Illarionov

reply via email to

[Prev in Thread] Current Thread [Next in Thread]