Integração Numérica - manipulação de NaNs (C / Fortran)

Estou lidando com uma integral complicada que exibe NaNs em certos valores próximos de zero e, no momento, estou lidando com eles de maneira bastante grosseira, usando uma instrução ISNAN que define o integrando como zero quando isso ocorre. Eu tentei isso com a biblioteca NMS em FORTRAN (a rotina q1da - q1dax não é diferente) e com a biblioteca GSL em C (usando a rotina QAGS).

Eu olhei para o CQUAD (parte da biblioteca GSL para C), que é projetado especificamente para lidar com NaNs e INF no integrando, mas há muito pouca informação útil na referência e nenhum exemplo de programas online que eu poderia encontrar. Alguém conhece outra rotina de integração numérica para C ou FORTRAN que possa fazer o trabalho?

Eu sou o autor da CQUADGSL. A interface é quase idêntica à de QAGS, portanto, se você usou o último, não deve ser difícil tentar o primeiro. Lembre-se de não converter se NaNe Infz em zeros no integrando - o código lidará com eles.

A rotina também está disponível no Octave as quadcce no Matlab aqui .

Você poderia fornecer um exemplo dos integrandos com os quais está lidando?

Atualizar

Aqui está um exemplo de uso CQUADpara integrar uma função com uma singularidade em um dos pontos de extremidade:

#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;
    gsl_integration_cquad_workspace *ws = NULL;
    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. */
    gsl_integration_cquad_workspace_free( ws );

    /* Bye. */
    return 0;

    }

com o qual eu compilei gcc -g -Wall cquad_test.c -lgsl -lcblas. A saída é

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

$0.94608307036718301494$

Observe que não há nada de especial aqui, nem para dizer CQUADonde está a singularidade ou qualquer tratamento especial dentro do próprio integrando. Apenas deixo que retorne NaNs, e o integrador cuida deles automaticamente.

Observe também que há um erro na versão 1.15 da GSL mais recente que pode afetar o tratamento de singularidades. Foi corrigido, mas não chegou à distribuição oficial. Eu usei a fonte mais recente, baixada com bzr branch http://bzr.savannah.gnu.org/r/gsl/trunk/.

Você também pode verificar as fórmulas de quadratura com dupla exponencial. Eles fazem uma alteração (implícita) de variáveis, certificando-se de que "facilitam" as singularidades dos limites. Uma implementação muito boa (Fortran77 e C) pode ser encontrada no site da Ooura .