Численная интеграция - обработка NaN (C / Fortran)

Я имею дело с хитрым интегралом, который демонстрирует NaN при определенных значениях, близких к нулю, и в настоящий момент я имею дело с ними довольно грубо, используя выражение ISNAN, которое устанавливает подынтегральное выражение в ноль, когда это происходит. Я пробовал это с библиотекой NMS в FORTRAN (процедура q1da - q1dax не отличается) и с библиотекой GSL в C (используя процедуру QAGS).

Я посмотрел на CQUAD (часть библиотеки GSL для C), которая специально разработана для обработки NaN и INF в подынтегральном выражении, но в справочнике очень мало полезной информации и нет примеров программ в Интернете, которые я мог бы найти. Кто-нибудь знает какие-либо другие процедуры числовой интеграции для C или FORTRAN, которые могли бы сделать эту работу?

Я автор CQUADв GSL. Интерфейс почти идентичен интерфейсу QAGS, поэтому, если вы использовали последний, вам не составит труда попробовать первый. Только помните, что не нужно преобразовывать ваши NaNs и Infs в нули в подынтегральном выражении - код сам справится с ними.

Процедура также доступна в Octave as quadccи Matlab здесь .

Не могли бы вы привести пример с интегратами, с которыми вы имеете дело?

Обновить

Вот пример использования CQUADдля интеграции функции с особенностью в одной из конечных точек:

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

    }

который я собрал gcc -g -Wall cquad_test.c -lgsl -lcblas. Выход

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

$0.94608307036718301494$

Обратите внимание, что здесь нет ничего особенного, ни для того, чтобы сказать, CQUADгде находится особенность, ни для какой-либо особой обработки внутри самой подынтегральной функции. Я просто позволю ему вернуть NaNs, и интегратор позаботится о них автоматически.

Также обратите внимание, что в последней версии GSL 1.15 есть ошибка, которая может повлиять на обработку сингулярностей. Это было исправлено, но не дошло до официального распространения. Я использовал самый последний источник, скачанный с bzr branch http://bzr.savannah.gnu.org/r/gsl/trunk/.

Вы также можете проверить квадратурные формулы с двумя экспонентами. Они делают (неявное) изменение переменных, удостоверяясь, что они «ослабляют» граничные особенности. Очень хорошая (Fortran77 и C) реализация может быть найдена на веб-сайте Ooura .