Numerical Integration - handling NaNs (C / Fortran)

I am dealing with a tricky integral that exhibits NaNs at certain values near zero and at the moment I am dealing with them quite crudely using an ISNAN statement which sets the integrand to zero when this occurs. I have tried this with the NMS library in FORTRAN (the q1da routine - q1dax is no different) and with the GSL library in C (using the QAGS routine).

I looked into CQUAD (part of the GSL library for C) which is specifically designed to handle NaNs and INF in the integrand, but there is very little useful info in the reference and no example programs online that I could find. Does anyone know any other numerical integration routine for either C or FORTRAN which could do the job?

-
Cross post on math.stackexchange.com/questions/242771/… –  user6 Nov 22 '12 at 19:43
^ I have deleted that post. –  Josh Nov 22 '12 at 19:52

I am the author of CQUAD in the GSL. The interface is almost identical to that of QAGS, so if you've used the latter, it should not be difficult at all to try the former. Just remember not to convert your NaNs and Infs to zeros in the integrand -- the code will deal with these itself.

The routine is also available in Octave as quadcc, and in Matlab here.

Could you provide an example of the integrands you are dealing with?

Update

Here's an example of using CQUAD to integrate a function with a singularity at one of the endpoints:

#include <stdio.h>
#include <gsl/gsl_integration.h>

/* Our test integrand. */
double thefunction ( double x , void *param ) {
return sin(x) / x;
}

/* Driver function. */
int main ( int argc , char *argv[] ) {

gsl_function f;
double res, abserr;
size_t neval;

/* Prepare the function. */
f.function = &thefunction;
f.params = NULL;

/* Initialize the workspace. */
if ( ( ws = gsl_integration_cquad_workspace_alloc( 200 ) ) == NULL ) {
printf( "main: call to gsl_integration_cquad_workspace_alloc failed.\n" );
abort();
}

/* Call the integrator. */
if ( gsl_integration_cquad( &f, 0.0 , 1.0 , 1.0e-10 , 1.0e-10 , ws , &res , &abserr , &neval ) != 0 ) {
printf( "main: call to gsl_integration_cquad failed.\n" );
abort();
}

/* Print the result. */
printf( "main: int of sin(x)/x in [0,1] is %.16e +/- %e (%i evals).\n" ,
res , abserr , neval );

/* Free the workspace. */

/* Bye. */
return 0;

}


which I compiled with gcc -g -Wall cquad_test.c -lgsl -lcblas. The output is

main: int of sin(x)/x in [0,1] is 9.4608307036718275e-01 +/- 4.263988e-13 (63 evals).


Which, given the result computed in Maple to 20 digits, $0.94608307036718301494$, is correct to 14 digits.

Note that there is nothing special here, neither to tell CQUAD where the singularity is, or any special treatment within the integrand itself. I just let it return NaNs, and the integrator takes care of them automatically.

Note also that there is a bug in the latest GSL version 1.15 which can affect the treatment of singularities. It has been fixed, but has not made it to the official distribution. I used the most recent source, downloaded with bzr branch http://bzr.savannah.gnu.org/r/gsl/trunk/.

-
Great, thanks for the reply. I'm using the integrator to find Green's Functions, and my integrand involves exponentials and some sine/cosines. I then integrate these again wrt a different variable and that is where I get the NaNs popping up. Do you know of any example programs using CQUAD? I am confused about how and where to put in the workspace functions. I should mention that I'm pretty much a beginner at this type of thing! –  Josh Nov 23 '12 at 10:47
@Josh: Good point, I guess somebody has to be the first to use it, so I've added a minimal example of how it can be called. –  Pedro Nov 26 '12 at 12:57

You could also check out the double-exponential quadrature formulas. They do an (implicit) change of variables, making sure that they "ease-away" boundary singularities. A very nice (Fortran77 and C) implementation can be found on Ooura's website.

-
Great, I'll give this a go. Thanks. –  Josh Nov 23 '12 at 10:49