help-gsl
[Top][All Lists]

## [Help-gsl] integration w/ qagp, bessel_Jnu

 From: Jorge Talamantes Subject: [Help-gsl] integration w/ qagp, bessel_Jnu Date: Thu, 10 Nov 2005 12:02:58 -0800

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;
pts=0.;
pts=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