[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.
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
- [Help-gsl] integration w/ qagp, bessel_Jnu,
Jorge Talamantes <=