两种多项式插值方法

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

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

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

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

K为多项式的阶次；

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

K = 2时

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

=9x²-23x+15

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

=15x²-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)多项式插值

牛顿多项式函数表示为：

 b₀ = y₀

n₀(x) = 1

j>0: 

n₁ = (x-x₀₎

n₂ = (x-x₀)(x-x₁)

n₃ = (x-x₀)(x-x₁)(x-x₂)

K为多项式的阶次；

举例：由函数y = x³ 生成数据集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)现象，但中间的数据平滑，并且计算量小，在实际使用中最好先进行图形化验证，根据插值数据动态选取计算数据，保证较好的拟合性能。