使用 GSL 库制作样条曲线并使用它们进行集成

Using the GSL Libraries to Make Splines and Using them for Integration

本文关键字：集成 GSL 样条曲线 使用      更新时间：2023-10-16

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

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

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

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

由于样条是平滑的，我从眼睛判断，我希望这种方法比评估样条和使用梯形等算法更有效

这是一个谬论。您假设数据通过某种阶数的样条比通过阶跃函数更好地近似，但您没有任何支持。你唯一拥有的是一堆(x, f(x))对。使用中点积分是在这里近似积分的一种完全值得尊敬的方法。一大优点：您可以自己轻松实现它。

我想出了如何使这种情况比梯形规则更快、更准确地发生。关键是使用样条对象作为结构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 nn", time_spent);
printf ("Result for Spline Integration     =  %.18fn", result);
printf ("Estimated error                   =  %.18fn", error);
printf ("Intervals                         =  %zunn", w->size);
printf("Result for Trapezoidal Integration = %f n", b);
gsl_integration_workspace_free (w);
free(x); free(f);
return 0;
}
相关文章：