我正在处理一个棘手的积分，它在接近零的某些值处显示 NaN，目前我正在使用 ISNAN 语句非常粗略地处理它们，该语句在发生这种情况时将被积函数设置为零。我已经用 FORTRAN 中的 NMS 库（q1da 例程 - q1dax 没有什么不同）和 C 中的 GSL 库（使用 QAGS 例程）尝试了这个。
我查看了专门设计用于处理被积函数中的 NaN 和 INF 的 CQUAD（C 的 GSL 库的一部分），但参考资料中几乎没有有用的信息，而且我找不到在线示例程序。有谁知道 C 或 FORTRAN 的任何其他数值积分例程可以完成这项工作？
我是CQUADGSL 的作者。界面与 的界面几乎一模一样QAGS，所以如果你用过后者，尝试前者应该一点都不难。请记住不要在被积函数中将您NaN的 s 和Infs 转换为零——代码将自己处理这些。
该例程在quadccOctave as和 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).
其中，给定在 Maple 中计算到 20 位的结果，, 正确到 14 位。
请注意，这里没有什么特别之处，既不告诉CQUAD奇点在哪里，也不告诉被积函数本身的任何特殊处理。我只是让它返回NaNs，积分器会自动处理它们。
另请注意，最新的 GSL 版本 1.15 中存在一个错误，可能会影响奇异点的处理。它已被修复，但尚未进入官方发行版。我使用了最新的源，用bzr branch http://bzr.savannah.gnu.org/r/gsl/trunk/.
您还可以查看双指数求积公式。他们对变量进行（隐式）更改，确保它们“缓解”边界奇点。在Ooura 的网站上可以找到一个非常好的（Fortran77 和 C）实现。