2

我正在研究依赖于计算多项式插值的内维尔算法。有关它的更多信息，您可以在http://en.wikipedia.org/wiki/Neville%27s_algorithm找到。在某些时候计算多项式对我来说不是问题。有很多关于它的来源。我的问题依赖于我不想在某些时候计算多项式，我想以这种形式获得多项式a_0 + a_1x + ... + a_nx^n：我不知道如何开始。你能给我一些建议吗？

4

2 回答 2

2

Neville 算法的一种变体允许计算一些常数（不是多项式系数），以便与另一个函数一起使用，该函数可以在任何点评估插值多项式。后一种功能很容易区分。下面的 C 代码是我使用的，我相信它可以正常工作。但是，我怀疑您可能会对多项式插值的效果感到失望，除非您提供的数据确实来自多项式。如果您可以轻松地在许多点对数据进行采样，那么最好通过对数据进行最小二乘拟合来找到多项式（表示为切比雪夫多项式之和）。

// fill C (allocated if null) with params for interpolating polynomial
// use params with interp_poly_eval
// !! these are NOT polynomial coefficients. 

double* interp_poly( Int deg, const double* x, const double* y, double* restrict C)
{
double* c = C ?: calloc( deg+1, sizeof *c);
Int i, j;
    memcpy( c, y, (deg+1)*sizeof *y);
    for (i=1; i<=deg; i++) 
    {   for (j=deg; j>=i; j--)
        {   c[j] = (c[j]-c[j-1]) / (x[j]-x[j-i]);
        }
    }
    return c;
}


double  interp_poly_eval( Int deg, const double* c, const double* x, double X)
{
double  p = c[deg];
Int i = deg;
    while( --i >= 0)
    {   p = c[i] + (X-x[i])*p;
    }
    return p;
}

// as above but also returns derivative of the polynomial through pdp
double  interp_poly_eval_d( Int deg, const double* c, const double* x, double X, double* pdp)
{
double  p = c[deg];
double  dp = 0.0;
Int i = deg;
    while( --i >= 0)
    {   dp = (X-x[i])*dp + p;
        p = c[i] + (X-x[i])*p;
    }   
    *pdp = dp;
    return p;
}
于 2014-10-31T12:32:09.527 回答
1

您需要有一个可以表示单变量（超过一个变量）多项式的类/数据结构。然后很容易计算实际的多项式（即它的系数），而不仅仅是使用 Neville 算法本身的单个点，即

  P[0,0] = y[0]   ; constant
  P[1,1] = y[1]   ; constant

然后递归

  P[a,b] = P[a,b-1] * x[b]/C + P[a,b-1] * -X/C + P[a+1,b] * -x[a]/C + P[a+1,b] * X/C

其中 X 是多项式“x”（即一阶单项式），C 是常数 X[b]-X[a]。

当您像往常一样遵循算法中的递归时，这将为您提供实际的多项式。

请注意，以上所有算术都是针对多项式的，即 P[a,b-1]*x[b]/C 表示多项式 P[a,b-1] 乘以常数 x[b]/C（即第 b 个点的 x 坐标除以 C)。

如果你想得到一个精确的结果，请使用任意精度的有理算术包（例如 C/C++ 的 GMP）。或者，使用浮点数进行计算，但最终可能会出现影响您的舍入错误。

于 2014-12-13T19:37:53.397 回答