两种多项式插值方法

泰勒公式表明许多复杂的三角函数、指数甚至对数方程都可以转换为多项式的形式，或在某一区间用多项式方程代替，多项式不但因为结构清晰简单，而且其各级导数连续、平滑且易于计算，所以在回归、曲线拟合、插值等领域普遍使用，工业机器人的路径规划、关节参数计算就是典型应用。

插值(Interpolation)是在有限的离散数据之间补充一些数据，使这组离散数据符合某个连续函数，是数学计算中最基本和常用的手段，插值的方法很多，如埃尔米特(Hermite)、埃特金(Aitken)、阿克玛(Akima)等等，以前介绍的多项式最小二乘回归也可以视作多项式(Polynomial)的插值形式，其实各类样条曲线(Spline)也是插值的应用，这里主要介绍一些常用插值的程序实现。

一、拉格朗日(Lagrange) 多项式插值

拉格朗日多项式函数表示为：





K为多项式的阶次；

举例：由函数y = x3 生成数据集x：[1，3，5，7] y：[1，27，125，343]

K = 2时

X:[1，3，5]区间：



=9x2-23x+15

X:[3，5，7]区间：



=15x2-71x+105

K = 3时

/// <summary>
/// Lagrange polynomial
/// </summary>
/// <param name="order">多项式阶次</param>
/// <param name="x">x数据集</param>
/// <param name="y">y数据集</param>
/// <param name="xp">插值x</param>
/// <param name="dy">一阶导数</param>
/// <returns>插值y</returns>
public static double Lagrange(int order, double[] x, double[] y, double xp,out double dy)
{
//            x = new double[] { 1, 3, 4, 7, 9, 10, 12, 15, 17 };
//            y = new double[] { 0, 2, 15, 12, 5, 3, 1, 0, -1 };
int len = x.Length;
int k = 0;
if (order < len)
{
while ((x[k] < xp) && k < len - 1)
k++;
k -= order / 2;
if (k < 0)
k = 0;
if (k + order > len - 1)
k = len - order - 1;
}
else
{
order = len - 1;
k = 0;
}
double nominator;
double denominator;
double yp = 0.0;
dy = 0;
for (int i = k; i <= k + order; i++)
{
nominator = 1.0;
denominator = 1.0;
double td = 0;
for (int j = k; j <= k + order; j++)
{
if (j != i)
{
nominator *= (xp - x[j]);
denominator *= (x[i] - x[j]);
double tdx = 1;
for (int m = k; m <= k + order; m++)
{
if (m != j&&m!=i)
{
tdx *= (xp - x[m]);
}
}
td += tdx;
}
}
if (denominator != 0.0)
{
yp += y[i] * nominator / denominator;
dy += y[i] * td / denominator;
}
else
{
yp = 0.0;
dy = 0;
}
}
return yp;
}

一、牛顿(Newton)多项式插值

牛顿多项式函数表示为：



b0 = y0

n0(x) = 1

j>0:



n1 = (x-x0)

n2 = (x-x0)(x-x1)

n3 = (x-x0)(x-x1)(x-x2)





K为多项式的阶次；

举例：由函数y = x3 生成数据集x：[1，3，5，7] y：[1，27，125，343]

K = 2时

X:[1，3，5]区间：





X:[3，5，7]区间：





K = 3时

/// <summary>
/// Newton Polynomial
/// </summary>
/// <param name="order">多项式阶次</param>
/// <param name="x">x数据集</param>
/// <param name="y">y数据集</param>
/// <param name="xp">插值x</param>
/// <param name="dy">一阶导数</param>
/// <returns>插值y</returns>
public static double NewTon(int order, double[] x, double[] y, double xp,out double dy)
{
//            x = new double[] { 1, 3, 4, 7, 9, 10, 12,15,17 };
//            y = new double[] { 0, 2, 15, 12, 5, 3, 1,0 ,-1};

int len = x.Length;
int k = 0;
if (order < len)
{
while ((x[k] < xp) && k < len - 1)
k++;
k -= order / 2;
if (k < 0)
k = 0;
if (k + order > len - 1)
k = len - order - 1;
}
else
{
order = len - 1;
k = 0;
}

//求多阶导数
double[] f = new double[order + 1];
f[0] = y[k];
for (int i = 1; i <= order; i++)
{
double temY = 0;
for (int j = 0; j <= i; j++)
{
double tmpD = 0;
{
tmpD = y[j + k];
for (int n = 0; n <= i; n++)
{
if (n != j)
tmpD /= x[j + k] - x[n + k];
}
}
temY += tmpD;
}
f[i] = temY;
}
dy = 0;
//  求Polynomial插值及一阶导数
double tempYp = 0;
double yp = f[0];
for (int i = 1; i <= order; i++)
{
tempYp = f[i];
double td = 0;
for (int j = 0; j < i; j++)
{
tempYp *= (xp - x[j + k]);
{
double tdx = 1;
for (int m = 0; m < i; m++)
{
if (m != j)
tdx *= (xp - x[m + k]);
}
td += tdx;
}
}
yp += tempYp;
dy += f[i] * td;
}
return yp;
}

拉格朗日与牛顿多项式插值易于在临近边界的区域出现数据振荡，特别是阶次较高时，也称为龙格(Runge's Phenomenon)现象，但中间的数据平滑，并且计算量小，在实际使用中最好先进行图形化验证，根据插值数据动态选取计算数据，保证较好的拟合性能。