[Help-gsl] integration w/ qagp, bessel_Jnu

[Jorge Talamantes]

Dear all,
 
I am trying to integrate the following function:
 
I = \int_0^{x1} M (D, alpha, x, n) dx,
 
where D, alpha and n are parameters to be passed to M, and
 
M = x^(D-alpha-1) * [ j(x,n+0.5) ]^2.
 
Here, j is the Bessel function of order (n + 0.5).
 
For some combinations of D and alpha, the integrand M diverges at the
origin. So, I am trying to use gsl_integration_qagp -- adaptive
integration with known singular points.
 
The problem I am having is that, for a given x, there is a maximum n for
which I can compute j(x,n+0.5) --  increasing n leads to an underflow
error from gsl_sf_bessel_Jnu.
 
Naturally, qagp chooses where to evaluate the integrand M. As it turns
out, as I increase the upper limit of integration (x1), I also need to
increase the value of n. Now, at some point, qagp chooses to evaluate M
for some small value of x -- but since n is set to some large value, I get
the underflow error.
 
I am appending a sample program below to illustrate my problem. Notice
that it works for x1 = 110.39, n = 107. But it dies for x=106.40, n = 107.
Evidently, gsl_sf_bessel_Jnu can compute j(110.39, 107.5), but not
j(106.40, 107.5).
 
Any insight would be deeply appreciated.
 
Jorge
 
 
#include <stdio.h>
#include <iostream>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_sf_bessel.h>
 
using namespace std;
 
// structure definition to pass parameters to gsl_integration_XXX
 
struct M{
  double D;
  int alpha;
  int n;
};
 
// function to be integrated is defined
 
double I_x1 (double x, void * p) {
 
  struct M * params = (struct M *)p;
  double D = (params->D);
  int alpha = (params->alpha);
  int n = (params->n);
 
  double N = 0.5 + static_cast<double>(n);
 
  double j = gsl_sf_bessel_Jnu (N,x);
  double f = pow(x, D-static_cast<double>(alpha)-1)*j*j;
  return f;
}
 
/* implementation of gsl's adaptive integration with known singular
points*/
 
double QAGP (double x1, double D, int alpha, int n) {
 
  gsl_integration_workspace * w
    = gsl_integration_workspace_alloc (1000);
   
  double result, error;
  double max_rel_error = 1.0e-12;
 
  struct M VM;
  VM.D     = D;
  VM.alpha = alpha;
  VM.n     = n;
 
  gsl_function F;
  F.function = &I_x1;
  F.params = &VM;
 
  double pts[2];
  pts[0]=0.;
  pts[1]=x1;
 
  size_t npts = 2;
 
  gsl_integration_qagp (&F, pts, npts, 0, max_rel_error, 1000,
                        w, &result, &error);
  return result;
}
 
// driver for testing purposes
 
int main(void){
 
  double x1    = 110.39;
  double D     = 1.8;
  int    alpha = 3;
  int    n     = 107;
 
  double I_n_alpha_D = QAGP (x1, D, alpha, n);
  printf("Integral = % .18f\n", I_n_alpha_D) ;
 
  x1 = 106.4;
 
  I_n_alpha_D = QAGP (x1, D, alpha, n);
  printf("Integral = % .18f\n", I_n_alpha_D) ;
}


---
Jorge Talamantes, Physics and Geology, California State University,
Bakersfield