如何使用gsl_interp获得多项式插值系数

How do I get the Polynomial Interpolation coefficients using gsl_interp?

本文关键字：多项式 插值 何使用 gsl interp      更新时间：2023-10-16

所以我有下面的代码。它完美地计算多项式的所有 y 点（并打印它们以使用 gnuplot 绘制），但是我如何获得结果多项式（在本例中为 1-x²）？

void twoDegreePoly() {
    int n = 3;
    double x[n],y[n];
    printf ("#m=0,S=16n");
    for (int i=0; i<n ;i++) {
        x[i] = ((double)2*i)/2 -1;
        y[i] = f(x[i]);
        printf ("%g %gn", x[i], y[i]);
    }
    printf ("#m=1,S=0n");
    gsl_interp_accel *acc = gsl_interp_accel_alloc ();
    const gsl_interp_type *t = gsl_interp_polynomial;
    gsl_interp* poly = gsl_interp_alloc(t,n);
    gsl_interp_init (poly, x, y,n);
    for (double xi=x[0]; xi<x[n-1]; xi+= 0.01) {
        double yi = gsl_interp_eval (poly, x, y, xi, acc);
        printf ("%g %gn", xi, yi);
    }
}

快速浏览文档后，GSL 中似乎没有这样的功能。这可能是由两个原因引起的：首先，获得多项式协效是特殊的，因为这种插值方法不适合一般设计（可以处理任意函数）。第二，引用数值配方：

但是，请确保系数是您需要的。通常，插值多项式的系数的确定精度远低于其在所需横坐标处的值。因此，确定系数仅用于计算插值不是一个好主意。这样计算的值不会完全通过表格点，例如，...

这样做的原因是，原则上，计算系数涉及求解具有范德蒙德矩阵的线性系统，该线性系统是高度不合格的。

尽管如此，数值配方给出了一个例程polcoe，您可以通过该获得插值多项式。您可以在第 3.5 章中找到它。在免费的第二版中。

我对Akima的插值做了类似的事情。首先，将状态定义为 GSL 执行的操作：

typedef struct
{
  double *b;
  double *c;
  double *d;
  double *_m;
}akima_state_t;

然后，创建插值器

spline = gsl_spline_alloc (gsl_interp_akima, M_size);
gsl_spline_init (spline, x, y, M_size);

之后，您可以执行以下操作：

const akima_state_t *state = (const akima_state_t *) ( spline -> interp -> state);
  double _b,_c,_d;
  for (int i = 0; i < M_size; i++)
  {
    _b = state->b[i];
    _c = state->c[i];
    _d = state->d[i];
    std::cout << "(x>"<<x[i]<<")*(x<"<<x[i+1]<<")*("<<y[i]<< "+ (x-"<< x[i]<<")*("<<_b<<"+(x-"<< x[i]<<")*("<<_c<<"+"<<_d<<"*(x-"<<x[i]<<")))) + ";
  }

我没有尝试过多项式插值，但这里的多项式的状态结构，它应该是一个很好的起点。

typedef struct
{
  double *d;
  double *coeff;
  double *work;
}
polynomial_state_t;
相关文章：