数值计算算法-多项式插值算法的实现与分析

数值计算是指在数值分析领域中的算法。数值分析是专门研究和数字以及近似值相关的数据问题，数值计算在数值分析的研究中发挥了特别重要的作用。

多项式插值是计算函数近似值的一种方法。其中函数值仅在几个点上已知。

该算法的基础是建立级数小于等于n的一个插值多项式p_n(z)，其中n+1是已知函数值的点的个数。

多项式插值法

许多问题都可以按照函数的方式来描述。然而，通常这个函数是未知的，我们只能通过少量的已知点来推断函数的大致模型。为了实现这个目的，在已知点之间做插值处理。如图所示，关于函数f(x)，已知的点为x₀...x₈，在图中以黑色的圆点表示。通过插值法的帮助，我们能够获取函数在z0、z1、z2处的值，在图中以白色小方块表示。本节主要讨论多项式插值法。

多项式插值法的根本点就是要建立一个特殊形式的多项式，称为插值多项式。

为了深入理解插值多项式的意义，我们先来看看多项式的一些基本法则：

首先，多项式是具有如下形式的函数：

p(x) = a₀ + a₁x + a₂x² + ... + a_nxⁿ

这里的a₀，...，a_n是系数。当a_n为非零整数时，这种形式的多项式称为n阶多项式。这是多项式的指数形式，在数学问题中尤为常见。但是，在某些特定的环境中其他形式的多项式则更为简便。比如，在有关多项式插值问题中，牛顿插值多项式就是一个很好的例子：

p(x) = a₀ + a₁(x - c₁) + a₂(x - c₁)(x - c₂) + ... + a_n(x - c₁)(x - c₂)...(x - c_n)

这里a₀，...，a_n是系数，而c₀，...，c_n是中值。注意到，当c₀，...，c_n全为0时，牛顿插值多项式就退化为前面定义的n阶多项式。

构建插值多项式

下面我们来看看如何对函数f(x)构建一个插值多项式。

为了对函数f(x)进行插值，要构建一个阶数小于等于n的多项式p_n(z)，而这又需要用到函数f(x)的n+1个已知点：x₀，...，x_n。这些已知点x₀，...，x_n就称为插值点。通过插值多项式p_n(z)可以计算出函数f(x)在x=z处的近似值。插值法需要满足点z在[x₀,x_n]内。可以采用如下公式来构建插值多项式p_n(z)。

p_n(z) = f[x₀] + f[x₀,x₁](z-x₀) + f[x₀,x₁,x₂](z-x₀)(z-x₁) + ... + f[x₀,...,x_n](z-x₀)(z-x₁)...(z-x_n-1)

其中函数f(x)在点x₀,...,x_n处的值已知，而f[x₀],...,f[x₀,...,x_n]则称为差商。

差商可以通过点x₀,...,x_n以及函数f(x)在这些点处的值来计算得出。这就是牛顿插值多项式的计算公式。注意到该公式同牛顿插值多项式的相同点。差商的计算公式为：

f[xi,...,xj] = f(xi) 如果i=j

f[x_i,...,x_j] = ( f[x_i+1,...x_j] - f[x_i,...x_j-1] ) / (x_j - x_i) 如果i < j

通过这个公式不难看出，当i < j时，必须预先计算出其他的差商值。例如，要计算f[x₀,x₁,x₂,x₃]，就需要先计算出f[x₁,x₂,x₃]和f[x₀,x₁,x₂]的值。幸运的是，可以通过一个差商表来帮助我们以一种系统的方式来计算差商值。如下图。

差商表由多行组成。最顶行保存已知点x₀,...,x_n 的值。第二行保存f[x₀],...,f[x_n]的值。要计算出表中其他的差商值，从每个待求的差商值处画一条对角线，使其回到f[xi]和f[xj]（如下图中差商f[x1,x2,x3]处的虚线）。要得到分母中的xi和xj，直接通过xi和xj求得。分子中的两个差商就是前一阶段计算出来的结果。

当完成了整个差商表的计算后，插值多项式的系数就是从第二行开始，每行最左边的那一项。

计算插值多项式

一旦确定了插值多项式的系数，对于函数f，如果我们想知道某个点处的函数值，只需要对多项式求值即可。

比如，已知函数f在4个点处的函数值：x0=-3.0，f（x0）=-5.0；x1=-2.0，f（x1）=-1.1；x2=2.0，f（x2）=1.9；x3=3.0，f（x3）=4.8；现在要求出点z0=-2.5，z1=0.0，z2=1.0，z3=2.5处的函数值。由于已经知道函数f在4个点处的值，因此插值多项式为3阶。下图是3阶插值多项式p3(z)的差商表。

一旦从差商表中得到了系数，就可以采用前面介绍过的牛顿公式来构建插值多项式p3(z)：

p3(z)=-5.0 + 3.9(z+3.0) + (-0.63)(z+3.0)(z+2.0) + 0.1767(z+3.0)(z+2.0)(z-2.0)

下一步可以通过该多项式计算出每个点z处的函数值。比如，在点z=-2.5处，经过如下计算得到：

p3(z)=-5.0 + 3.9(-2.5+3.0) + (-0.63)(-2.5+3.0)(-2.5+2.0) + 0.1767(-2.5+3.0)(-2.5+2.0)(-2.5-2.0) = -2.694

点z1、z2、z3处的函数值可以通过相似的方法计算得出。最终结果以表格和函数图像的方式表达。如下图。

和任何其他近似算法一样，通常会有一些与插值多项式相关的误差出现，理解这一点很重要。定性的来讲，如果要使误差降至最小，构建的插值多项式必须要在函数f(x)上获取足够多的已知点才行。并且点与点之前的距离要适当，这样最终得到的多项式才能精确地表示出函数的特性。

多项式插值的接口定义

interpol

int interpol ( const double *x, const double *fx, int n, double *z, double *pz, int m );

返回值：如果插值操作成功，返回0；否则返回-1；

描述：采用多项式插值法来求得函数在某些特定点上的值。

由调用者在参数x处指定函数值已知的点集。每个已知点所对应的函数值都在fx中指定。对应待求的点由参数z来指定，而z所对应的函数值将在pz中返回。x和f（x）中的元素个数由参数n来表示。z中的待求点的个数（以及pz中返回值个数）由参数m来表示。由调用者管理x、fx、z以及pz所关联的存储空间。

复杂度：O（mn²），这里m代表待求值的个数，而n代表已知点的个数。

多项式插值的实现与分析

多项式插值法主要基于对一系列期望点的插值多项式的确定。要得到这个多项式，首先必须通过计算差商来确定该多项式的系数。

首先，为差商以及待确定的系数分配存储空间。注意，由于差商表中每一行的每一项都仅依赖于其前一行的计算结果，因此，并不需要一次性保留所有的表项。所以，只为最占用空间的行分配空间即可，该行将有n个条目。

接下来，用fx中的值来初始化差商表的第一行。这是为计算差商表中的第三行做准备。（前两行不需要计算，因为这两行中的条目都已经保存在x和fx中）。

初始化的最后一步是在coeff[0]中保存fx[0]的值，因为这是插值多项式的第一个系数。

计算差商的过程涉及一个嵌套循环，我们在循环中根据前面介绍过的公式来计算差商。在外层循环中，k用来统计正在计算的是哪一行（排除x和fx所代表的行）。在内层循环中，i表示在当前行中正在计算的是哪一个条目。一旦计算完一行的条目，table[0]中的值就成为插值多项式的下一个系数。因此，保存该值到coeff[k]中。一旦得到插值多项式的所有系数，就可以计算出z中每个目标点的值，最后将这些值保存在pz中。

我们命名这个函数为interpol，它的时间复杂度为O（mn^₂），这里m代表z中的元素个数（也是pz中值的个数），n代表x中（也是fx中）的元素个数。复杂度因子n²是这样得到的，变更j控制循环中的每次迭代，当前迭代中需要乘以的因子比上一轮要多一个。也就是说，当j=1时，term需要做1次乘法，当j=2时，term需要做2次乘法，持续这个过程直到j=n-1时，term需要做n-1次乘法。实际上，这就成了对1~n-1的整数求和，得到的计算时间为T（n）=（n（n+1）/2）-n，再乘以某段固定的时间。（这是由计算等差数列的著名公式得到的）。在大O记法中可以简化为O（n²）。O（mn²）中的因子m来源于针对z中的每个点计算多项式插值的过程。在第一个嵌套循环中，计算出所有的差商，其复杂度为O（n²）。因此，最终的复杂度有一个额外的因子m，该因子对实际的复杂度没有多大的影响。

示例：多项式插值的实现

/*interpol.c*/
#include <stdlib.h>
#include <string.h>#include "nummeths.h"/*interpol  */
int interpol(const double *x, const double *fx, int n, double *z, double *pz, int m)
{double term, *table, *coeff;int i,j,k;/*为差商和待确定的系数分配空间*/if((table = (double *)malloc(sizeof(double)*n)) == NULL)return -1;if((coeff = (double *)malloc(sizeof(double)*n)) == NULL){free(table);return -1;}/*初始化差商表*/memcpy(table,fx,sizeof(double)*n);/*重点：确定差商表的系数*/coeff[0] = table[0];for(k=1; k<n; k++)/*外层循环k用来统计正在计算的是哪一行*/{for(i=0; i<n-k; i++)/*内层循环i表示在当前行中正在计算的是哪一个条目（随之行数的增加，条目减少）*/{j=i+k;/*当前行的每一项中分子的差商就是其前一阶段计算的结果*/table[i] = (table[i+1] - table[i]) / (x[j] - x[i]);}/*当前行的首个条目计算结果即是多项式的下一个系数*/coeff[k]=table[0];}free(table);/*在指定点上对插值多项式进行求值(循环构造插值多项式)*/for(k=0; k<m; k++)          /*最外层：遍历z点数组*/{/*插值多项式的第一个因子*/pz[k] = coeff[0];       for(j=1; j<n; j++)      /*嵌套：构造多项式（新算式等于上一步的结果加上新因子）*/{term = coeff[j];    /*因子构成：以多项式系数为基础*/for(i=0; i<j; i++)  /*嵌套：新因子以上一步的结果乘以（z[k] - x[i]）*/term=term*(z[k] - x[i]);pz[k]=pz[k] + term;}}free(coeff);return 0;
}