[Top][All Lists]

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

[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.
#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
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];
  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,

reply via email to

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