Integracja numeryczna - obsługa NaNs (C / Fortran)

Mam do czynienia z trudną całką, która wykazuje wartości NaN przy pewnych wartościach zbliżonych do zera, a w tej chwili radzę sobie z nimi dość brutalnie, używając instrukcji ISNAN, która ustawia intandand na zero, gdy to nastąpi. Próbowałem tego z biblioteką NMS w FORTRAN (procedura q1da - q1dax nie różni się) i z biblioteką GSL w C (używając procedury QAGS).

Zajrzałem do CQUAD (część biblioteki GSL dla C), która jest specjalnie zaprojektowana do obsługi NaN i INF w integrandzie, ale w referencji jest bardzo mało przydatnych informacji i nie ma żadnych przykładowych programów online, które mogłem znaleźć. Czy ktoś zna jakąkolwiek inną procedurę integracji numerycznej dla C lub FORTRAN, która mogłaby wykonać to zadanie?

Jestem autorem CQUADw GSL. Interfejs jest prawie identyczny z interfejsem QAGS, więc jeśli używałeś tego drugiego, nie powinno być w ogóle trudne. Wystarczy pamiętać, aby nie konwertować NaNS i InfS do zera w podcałkowej - kod będzie radzić sobie z tymi siebie.

Procedura jest również dostępna w Octave as quadccoraz w Matlabie tutaj .

Czy możesz podać przykład integrandów, z którymi masz do czynienia?

Aktualizacja

Oto przykład użycia CQUADdo zintegrowania funkcji z osobliwością w jednym z punktów końcowych:

#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;

    }

z którą skompilowałem gcc -g -Wall cquad_test.c -lgsl -lcblas. Dane wyjściowe to

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

$0.94608307036718301494$

Zauważ, że nie ma tu nic specjalnego, ani powiedzieć, CQUADgdzie jest osobliwość, ani żadnego specjalnego traktowania w obrębie samego integranda. Po prostu pozwalam mu zwrócić NaNs, a integrator zajmie się nimi automatycznie.

Należy również zauważyć, że w najnowszej wersji GSL 1.15 występuje błąd, który może wpływać na leczenie osobliwości. Zostało to naprawione, ale nie dotarło do oficjalnej dystrybucji. Użyłem najnowszego źródła, pobranego z bzr branch http://bzr.savannah.gnu.org/r/gsl/trunk/.

Możesz także sprawdzić formuły kwadraturowe o podwójnej wykładniczej wartości. Dokonują (domyślnej) zmiany zmiennych, upewniając się, że „łagodzą” osobliwości graniczne. Bardzo ładną implementację (Fortran77 i C) można znaleźć na stronie internetowej Ooura .