使用GSL进行光滑单峰函数的数值积分？

Question

有人可以发布一个在有限区间内对一个光滑单峰函数进行数值积分的GSL简单示例吗？

- a06e

---

你所说的“smooth”是什么意思？是Lipschitz吗？ - Steve Cox

@SteveCox 这意味着它具有连续的导数，直到某个高阶，比如说10阶导数。 - a06e

分析工具？还是只有n次可微性？ - Steve Cox

@SteveCox 我无法对解析性做出任何评论。你可以假设只有连续的导数。 - a06e

你不希望导数也是单峰的，对吧？ - Steve Cox

Answer 1

以下是一个例子，对 1/(t^2 + 1) 在 [0,1000] 区间上进行积分。由于没有奇点存在，因此使用最简单的自适应积分规则。

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>

double f (double x, void * params) {
    double alpha = *(double *) params;
    double f = alpha / (x * x + 1);
    return f;
}

    int
main (void)
{
    gsl_integration_workspace * w 
        = gsl_integration_workspace_alloc (1000);

    double result, error;
    double alpha = 1.0;


    gsl_function F;
    F.function = &f;
    F.params = &alpha;

    gsl_integration_qag (&F,
                         0.0, 1000.0,
                         0.0, 1e-7, 1000,
                         GSL_INTEG_GAUSS15,
                         w,
                         &result, &error);

    printf ("result          = % .18f\n", result);
    printf ("estimated error = % .18f\n", error);

    gsl_integration_workspace_free (w);

    return 0;
} 

结果如下：

result          =  1.569796327128230029
estimated error =  0.000000000092546021

这是有道理的，因为积分应该约为 pi/2。

- Steve Cox

---

你应该使用宏“GSL_INTEG_GAUSS15”来明确第七个参数。 - Vivian Miranda

@ViniciusMiranda 好的呼叫 - Steve Cox