c++ - 如何使用 gsl_interp 获得多项式插值系数？

Question

所以我有下面的代码。它完美地计算了多项式的所有 y 点(并打印它们以使用 gnuplot 绘制)，但我如何获得结果多项式(在这种情况下为 1-x²)？

void twoDegreePoly() {
    int n = 3;
    double x[n],y[n];
    printf ("#m=0,S=16\n");
    for (int i=0; i<n ;i++) {
        x[i] = ((double)2*i)/2 -1;
        y[i] = f(x[i]);
        printf ("%g %g\n", x[i], y[i]);
    }
    printf ("#m=1,S=0\n");

    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 %g\n", xi, yi);
    }
}

Answer 1

快速扫描 documentation 之后, GSL 中似乎没有这样的功能。这可能是由两个原因引起的:首先，获取多项式系数对于这种插值方法来说是特殊的，不适合一般设计(可以处理任意函数)。二、引用数值食谱:

Please be certain, however, that the coefficients are what you need. Generally, the coefficients of the interpolating polynomial can be determined much less accurately than its value at a desired abscissa. Therefire, it is not a good idea to determine the coefficients only for use in calculating interpolating values. Values thus calculated will not pass exactly through the tabulated points, for example, ...

原因是原则上，计算系数涉及求解具有高度病态的 Vandermonde 矩阵的线性系统。

不过，Numerical Recipes 提供了一个例程 polcoe，您可以通过它获得插值多项式。您可以在第 3.5 章中找到它。在free second edition .