Numerische Integration - Umgang mit NaNs (C / Fortran)

12

Ich habe es mit einem trickreichen Integral zu tun, das NaNs bei bestimmten Werten nahe Null aufweist, und ich habe es im Moment ziemlich grob mit einer ISNAN-Anweisung zu tun, die den Integranden auf Null setzt, wenn dies auftritt. Ich habe dies mit der NMS-Bibliothek in FORTRAN (die Q1DA-Routine - Q1dax ist nicht anders) und mit der GSL-Bibliothek in C (mit der QAGS-Routine) versucht.

Ich habe mich mit CQUAD (Teil der GSL-Bibliothek für C) befasst, das speziell für den Umgang mit NaNs und INF im Integranden entwickelt wurde, aber in der Referenz gibt es nur sehr wenige nützliche Informationen und keine Online-Beispielprogramme, die ich finden konnte. Kennt jemand eine andere numerische Integrationsroutine für C oder FORTRAN, die diese Aufgabe übernehmen könnte?

Josh
quelle
^ Ich habe diesen Beitrag gelöscht.
Josh

Antworten:

10

Ich bin der Autor CQUADin der GSL. Die Benutzeroberfläche ist fast identisch mit der von QAGS. Wenn Sie also die letztere verwendet haben, sollte es nicht schwierig sein, die erstere zu testen. Denken Sie daran, Ihr NaNs und Infs im Integranden nicht in Nullen umzuwandeln - der Code wird sich selbst darum kümmern.

Die Routine ist auch in Octave als quadccund in Matlab hier .

Können Sie ein Beispiel für die Integranden nennen, mit denen Sie sich befassen?

Aktualisieren

Im Folgenden finden Sie ein Beispiel für die CQUADIntegration einer Funktion mit einer Singularität an einem der Endpunkte:

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

    }

mit denen ich kompiliert habe gcc -g -Wall cquad_test.c -lgsl -lcblas. Die Ausgabe ist

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

0.94608307036718301494

Beachten Sie, dass es hier nichts Besonderes gibt, weder um zu sagen, CQUADwo sich die Singularität befindet, noch um eine spezielle Behandlung innerhalb des Integranden. Ich lasse es einfach NaNs zurückgeben und der Integrator kümmert sich automatisch um sie.

Beachten Sie auch, dass die neueste GSL-Version 1.15 einen Fehler enthält, der sich auf die Behandlung von Singularitäten auswirken kann. Es wurde behoben, hat es aber nicht in die offizielle Distribution geschafft. Ich habe die neueste Quelle verwendet, die mit heruntergeladen wurde bzr branch http://bzr.savannah.gnu.org/r/gsl/trunk/.

Pedro
quelle
Super, danke für die Antwort. Ich benutze den Integrator, um die Funktionen von Green zu finden, und mein Integrand enthält Exponentiale und einige Sinus / Cosinus. Ich integriere diese dann wieder in eine andere Variable und bekomme dort die NaNs hoch. Kennen Sie Beispielprogramme mit CQUAD? Ich bin verwirrt darüber, wie und wo ich die Funktionen des Arbeitsbereichs einfügen soll. Ich sollte erwähnen, dass ich so ziemlich ein Anfänger in dieser Art von Dingen bin!
Josh
@Josh: Guter Punkt, ich denke, jemand muss der Erste sein, der es benutzt, also habe ich ein minimales Beispiel hinzugefügt, wie es aufgerufen werden kann.
Pedro
3

Sie können auch die doppelt exponentiellen Quadraturformeln überprüfen. Sie führen eine (implizite) Änderung von Variablen durch, um sicherzustellen, dass sie Rand-Singularitäten "beseitigen". Eine sehr schöne Implementierung (Fortran77 und C) finden Sie auf der Website von Ooura .

GertVdE
quelle