c++ - 使用 GSL 库制作样条曲线并使用它们进行积分

Question

假设我有一组 N 个数据点。我可以使用 gsl 库 gsl_splines.h 例程来创建此数据的样条曲线。我想做的是使用这个样条和gsl集成库来找到这些数据的积分。我在这里用 C 语言工作。

在我的代码中，我生成了我要使用的样条线，并且由于样条线是平滑的，我用肉眼判断，我希望这种方法比评估样条线和使用梯形规则之类的算法更有效找到积分，但我无法想出一种将这两件事拼凑在一起的方法。

如果你能提供任何简单的例子，我将不胜感激！

如果 gsl 库不是您要使用的，我很高兴听到任何其他建议。

标签： c++cinterpolationgslnumerical-integration

Answer 1

我想出了如何使这种情况发生得更快，并且比梯形规则更准确。关键是使用样条对象作为结构 my_f_params 的成员，因为 GSL 集成需要一个 gsl_function 对象，

https://www.gnu.org/software/gsl/doc/html/integration.html#c.gsl_integration_qags。

这是一个将 1/x 从 1 积分到 1200 的示例代码，这不一定是漂亮的代码，我是物理学家而不是计算机科学家：

// Complie with 
// gcc -w -o Spline_Integration_Test Spline_Integration_Test.c -lgsl -lgslcblas -lm

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


double funk(double x, void *p);
struct my_f_params { gsl_interp_accel *facc; gsl_spline *fspline; };

double funk(double x, void *p)
{
    struct my_f_params * params = (struct my_f_params *)p;
    gsl_interp_accel *facc = (params->facc);
    gsl_spline *fspline = (params->fspline);

    double f = gsl_spline_eval (fspline, x, facc);
    return f;
}

int main()
{
    int i, N = 10000000;
    double *x; x = (double *)malloc((size_t)sizeof(double) * N); memset(x,0,(size_t) sizeof(double)*N);
    double *f; f = (double *)malloc((size_t)sizeof(double) * N); memset(f,0,(size_t) sizeof(double)*N);
    gsl_interp_accel *facc = gsl_interp_accel_alloc();
    gsl_spline *fspline = gsl_spline_alloc (gsl_interp_cspline, N);

    x[0] = 0.0;
    f[0] = x[0]*x[0];
    for(i=1; i<N; i++)
    {
        x[i] = x[i-1] + 0.001;
        f[i] = 1/x[i];
    }

    gsl_spline_init(fspline, x, f, N);

    clock_t begin = clock();
    gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);

    struct my_f_params params = { facc, fspline };

    double result, error;

    gsl_function F;
    F.function = &funk;
    F.params = &params;

    // Beginning GSL/spline integration part
    gsl_integration_qags (&F, 1, 1200, 0, 1e-7, 1000, w, &result, &error);

    clock_t end = clock();
    double time_spent = (double)(end - begin) / CLOCKS_PER_SEC;
    printf("Time for Spline Integration:         = %9f s \n", time_spent);


    // Begining trapezoidal integration part
    begin = clock();
    double a, b = 0.0;
    double delta;
    for(i=1000; i<1200*1000; i++)
    {
        delta = 0.001;
        if(i==1||i==N) a = 1.0/x[i];
        else a = 2.0/x[i];
        b += a*(delta)/2.0;
    }
    end = clock();
    time_spent = (double)(end - begin) / CLOCKS_PER_SEC;
    printf("Time for Trapezoidal Integration:    = %9f s \n\n", time_spent);

    printf ("Result for Spline Integration     =  %.18f\n", result);
    printf ("Estimated error                   =  %.18f\n", error);
    printf ("Intervals                         =  %zu\n\n", w->size);

    printf("Result for Trapezoidal Integration = %f \n", b);

    gsl_integration_workspace_free (w);
    free(x); free(f);

    return 0;
}