Integrasi Numerik - menangani NaNs (C / Fortran)

12

Saya berurusan dengan integral rumit yang menunjukkan NaN pada nilai-nilai tertentu mendekati nol dan saat ini saya berurusan dengan mereka secara kasar menggunakan pernyataan ISNAN yang menetapkan integrand ke nol ketika ini terjadi. Saya telah mencoba ini dengan perpustakaan NMS di FORTRAN (rutin q1da - q1dax tidak berbeda) dan dengan perpustakaan GSL di C (menggunakan rutin QAGS).

Saya melihat ke CQUAD (bagian dari pustaka GSL untuk C) yang dirancang khusus untuk menangani NaNs dan INF di integrand, tetapi ada sedikit info yang berguna dalam referensi dan tidak ada contoh program online yang bisa saya temukan. Adakah yang tahu rutin integrasi numerik lainnya untuk C atau FORTRAN yang dapat melakukan pekerjaan itu?

Josh
sumber
^ Saya telah menghapus posting itu.
Josh

Jawaban:

10

Saya penulis CQUADdi GSL. Antarmuka hampir identik dengan itu QAGS, jadi jika Anda telah menggunakan yang terakhir, seharusnya tidak sulit sama sekali untuk mencoba yang pertama. Ingatlah untuk tidak mengonversikan angka NaNs dan Infnol menjadi angka di integand - kode akan menangani ini sendiri.

Rutin ini juga tersedia di Octave as quadcc, dan di Matlab di sini .

Bisakah Anda memberikan contoh integrand yang Anda hadapi?

Memperbarui

Berikut adalah contoh penggunaan CQUADuntuk mengintegrasikan fungsi dengan singularitas di salah satu titik akhir:

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

    }

yang saya susun gcc -g -Wall cquad_test.c -lgsl -lcblas. Outputnya adalah

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

Yang, mengingat hasil yang dihitung dalam Maple menjadi 20 digit, , benar menjadi 14 digit.0.94608307036718301494

Perhatikan bahwa tidak ada yang istimewa di sini, tidak ada yang tahu di CQUADmana singularitas itu, atau perlakuan khusus apa pun di dalam integrand itu sendiri. Saya biarkan saja kembali NaN, dan integrator merawatnya secara otomatis.

Perhatikan juga bahwa ada bug dalam versi GSL terbaru 1.15 yang dapat mempengaruhi perawatan singularitas. Sudah diperbaiki, tetapi belum sampai ke distribusi resmi. Saya menggunakan sumber terbaru, diunduh dengan bzr branch http://bzr.savannah.gnu.org/r/gsl/trunk/.

Pedro
sumber
Bagus, terima kasih atas jawabannya. Saya menggunakan integrator untuk menemukan Fungsi Green, dan integrand saya melibatkan eksponensial dan beberapa sinus / cosinus. Saya kemudian mengintegrasikan ini lagi wrt variabel yang berbeda dan di situlah saya mendapatkan NaN bermunculan. Apakah Anda mengetahui contoh program yang menggunakan CQUAD? Saya bingung tentang bagaimana dan di mana harus meletakkan fungsi workspace. Saya harus menyebutkan bahwa saya cukup pemula dalam hal semacam ini!
Josh
@Josh: Poin bagus, saya kira seseorang harus menjadi yang pertama menggunakannya, jadi saya telah menambahkan contoh minimal bagaimana itu bisa dipanggil.
Pedro
3

Anda juga bisa melihat rumus kuadratur eksponensial ganda. Mereka melakukan perubahan variabel (implisit), memastikan bahwa mereka "melonggarkan" singularitas batas. Implementasi yang sangat bagus (Fortran77 dan C) dapat ditemukan di situs web Ooura .

GertVdE
sumber