[Bug-gsl] gsl_integration_qagiu failure

[Manoj Rajagopalan]

GSL version: 1.4
Obtained from: http://sources.redhat.com/gsl/ (one of the US mirrors)
Hardware: Intel P4 1.7GHz mobile processor, Dell Laptop
OS: Red Hat Linux 9.0, Kernel 2.4.20-28.9
Compiler: gcc version 3.2.2 20030222 (Red Hat Linux 3.2.2-5)

Bug behaviour:
   The result of integration of the square of the radial component of 
the normalized hydrogenic bound state wave function (double 
gsl_sf_hydrogenicR_1(double Z, double r)  from gsl_sf_coulomb.h) over 
range of r from 0 to INFINITY should be 1.0

   Using gsl_integration_qagiu to perform this semi-infinite integral 
yields the  result 0.0 or core-dumps with an underflow. See file 
qagiu.txt (attached)

   Using gsl_integration_qag or gsl_integration_qags with a 
sufficiently large upper limit for r ( 10*a where a is calculated in the 
attached code) gives the expected result of 1.0. See files qag.txt and 
qags.txt (both attached).


Code snippet: integrl-zero2inf.cc , Makefile (both attached)
  Compiled using g++
  Make command: % make integrl-zero2inf
  Run command:  % ./integrl-zero2inf

   Important lines of code: lines 53 - 66 where you can choose one of 
qagiu, qag or qags to integrate the function with.

  Note: Makefile may have to be modified to tell the compiler where the 
GSL includes and libraries are located.

Script started on Thu 26 Feb 2004 03:45:59 PM EST
address@hidden gsl]$ make clean
rm -f gsl-try coulomb_wavefn integrn integrn2 integrn2D integrn-circle 
integrl-zero2inf *.o
address@hidden gsl]$ make integrl-zero2inf
g++ -c  -DHAVE_INLINE integrl-zero2inf.cc
g++ -o integrl-zero2inf integrl-zero2inf.o  -lgsl -lgslcblas -lm -DHAVE_INLINE
address@hidden gsl]$ integrl-zero2inf 
Probability of finding an electron in a 1s orbital around a hydrogen-like 
nucleus
Enter atomic number: 1
APPROACH 2: Using QAG with 61 points
Z = 1   a = 5.29465e-11 Zbya = 1.8887e+10
Cumulative probability = 1
Absolute error = 1.11022e-14
GSL error code = 0
address@hidden gsl]$ exit

Script done on Thu 26 Feb 2004 03:46:25 PM EST

Script started on Thu 26 Feb 2004 03:44:48 PM EST
address@hidden gsl]$ make clean
rm -f gsl-try coulomb_wavefn integrn integrn2 integrn2D integrn-circle 
integrl-zero2inf *.o
address@hidden gsl]$ make integrl-zero2inf
g++ -c  -DHAVE_INLINE integrl-zero2inf.cc
g++ -o integrl-zero2inf integrl-zero2inf.o  -lgsl -lgslcblas -lm -DHAVE_INLINE
address@hidden gsl]$ integrl-zero2inf
Probability of finding an electron in a 1s orbital around a hydrogen-like 
nucleus
Enter atomic number: 1
APPROACH 1: Using QAGIU!
gsl: coulomb_bound.c:66: ERROR: underflow
Default GSL error handler invoked.
Abort (core dumped)
address@hidden gsl]$ exit

Script done on Thu 26 Feb 2004 03:45:27 PM EST

Script started on Thu 26 Feb 2004 03:46:57 PM EST
address@hidden gsl]$ make clean
rm -f gsl-try coulomb_wavefn integrn integrn2 integrn2D integrn-circle 
integrl-zero2inf *.o
address@hidden gsl]$ make integrl-zero2inf
g++ -c  -DHAVE_INLINE integrl-zero2inf.c
g++ -o integrl-zero2inf integrl-zero2inf.o  -lgsl -lgslcblas -lm -DHAVE_INLINE
address@hidden gsl]$ integrl-zero2inf
Probability of finding an electron in a 1s orbital around a hydrogen-like 
nucleus
Enter atomic number: 1
APPROACH 3: Using QAGS
Z = 1   a = 5.29465e-11 Zbya = 1.8887e+10
Cumulative probability = 1
Absolute error = 4.21826e-11
GSL error code = 0
address@hidden gsl]$ exit

Script done on Thu 26 Feb 2004 03:47:20 PM EST

#include <iostream>
#include <cmath>
#include <gsl/gsl_const_mksa.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_sf_coulomb.h>

// Verifies that the cumulative probability of finding an electron in space 
around a
// hydrogen-like nucleus is unity


double R(double r, void *p)
{
  // Radial probability density of finding an electron around a H nucleus
  // Levine, pg 145
  double Zbya = *static_cast<double*>(p);
  double Rval = gsl_sf_hydrogenicR_1(Zbya, r);
  return Rval*Rval*r*r;
}


using namespace std;

int main(void)
{
  // allocate workspace
  gsl_integration_workspace *gsl_wk(NULL);
  gsl_wk = gsl_integration_workspace_alloc(1000);
  assert(gsl_wk != NULL);

  // Intro
  cout << "Probability of finding an electron in a 1s orbital around a 
hydrogen-like nucleus" << endl
       << "Enter atomic number: ";
  int Z;
  cin >> Z;

  const double hbar = GSL_CONST_MKSA_PLANCKS_CONSTANT_HBAR;
  const double eps0 = GSL_CONST_MKSA_VACUUM_PERMITTIVITY;
  const double m_e = GSL_CONST_MKSA_MASS_ELECTRON;
  const double m_p = GSL_CONST_MKSA_MASS_PROTON;
  const double mu = (m_e*m_p)/(m_e+m_p);
  const double q0 = GSL_CONST_MKSA_ELECTRON_CHARGE;
  const double a = (hbar*hbar * 4.0*M_PI*eps0)/ (mu*q0*q0);
  double Zbya = Z/a;

  // setup integrand
  gsl_function integrand;
  integrand.function = R;
  integrand.params = static_cast<void*>(&Zbya);
  double result(-1), abserr(-1);
  
  // APPROACH 1: **** This causes the error: either zero result or underflow
  cout << "APPROACH 1: Using QAGIU!" << endl;
  int errcode = gsl_integration_qagiu(&integrand, 0.0, 1.0e-7, 1.0e-7, 1000, 
gsl_wk,
                                      &result, &abserr);

//   // APPROACH 2
//   cout << "APPROACH 2: Using QAG with 61 points" << endl;
//   int errcode = gsl_integration_qag(&integrand, 0.0, 10*a, 1.0e-7, 1.0e-7, 
1000,
//                                  GSL_INTEG_GAUSS61, gsl_wk, &result, 
&abserr);
  
//   // APPROACH 3
//   cout << "APPROACH 3: Using QAGS" << endl;
//   int errcode = gsl_integration_qags(&integrand, 0.0, 10*a, 1.0e-7, 1.0e-7, 
1000,
//                                   gsl_wk, &result, &abserr);

  // outputs
  cout << "Z = " << Z << "\ta = " << a << "\tZbya = " << Zbya << endl;
  cout << "Cumulative probability = " << result << endl;
  cout << "Absolute error = " << abserr << endl;
  cerr << "GSL error code = " << errcode << endl;

  // deallocate workspace
  gsl_integration_workspace_free(gsl_wk);
  return 0;
}

CXX = g++
CC = gcc

CXXFLAGS = -DHAVE_INLINE

INCLUDEDIRS = 
LIBDIRS = 

LIBS = -lgsl -lgslcblas -lm

.cc.o:
        $(CXX) -c $(INCLUDEDIRS) $(CXXFLAGS) $<

TARGETS = gsl-try coulomb_wavefn integrn integrn2 integrn2D \
        integrn-circle integrl-zero2inf

all:    $(TARGETS)


gsl-try: gsl-try.o
        $(CXX) -o $@ $< $(LIBDIRS) $(LIBS) $(CXXFLAGS)

coulomb_wavefn: coulomb_wavefn.o
        $(CXX) -o $@ $< $(LIBDIRS) $(LIBS) $(CXXFLAGS)

integrn: integrn.o
        $(CXX) -o $@ $< $(LIBDIRS) $(LIBS) $(CXXFLAGS)

integrn2: integrn2.o
        $(CXX) -o $@ $< $(LIBDIRS) $(LIBS) $(CXXFLAGS)

integrn2D: integrn2D.o
        $(CXX) -o $@ $< $(LIBDIRS) $(LIBS) $(CXXFLAGS)

integrn-circle: integrn-circle.o
        $(CXX) -o $@ $< $(LIBDIRS) $(LIBS) $(CXXFLAGS)

integrl-zero2inf: integrl-zero2inf.o
        $(CXX) -o $@ $< $(LIBDIRS) $(LIBS) $(CXXFLAGS)


clean:
        rm -f $(TARGETS) *.o