数值计算算法-多项式插值算法的实现与分析
Posted DreamGo
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了数值计算算法-多项式插值算法的实现与分析相关的知识，希望对你有一定的参考价值。
数值计算是指在数值分析领域中的算法。数值分析是专门研究和数字以及近似值相关的数据问题，数值计算在数值分析的研究中发挥了特别重要的作用。
多项式插值是计算函数近似值的一种方法。其中函数值仅在几个点上已知。
该算法的基础是建立级数小于等于n的一个插值多项式pn(z)，其中n+1是已知函数值的点的个数。
多项式插值法
许多问题都可以按照函数的方式来描述。然而，通常这个函数是未知的，我们只能通过少量的已知点来推断函数的大致模型。为了实现这个目的，在已知点之间做插值处理。如图所示，关于函数f(x)，已知的点为x0...x8，在图中以黑色的圆点表示。通过插值法的帮助，我们能够获取函数在z0、z1、z2处的值，在图中以白色小方块表示。本节主要讨论多项式插值法。
多项式插值法的根本点就是要建立一个特殊形式的多项式，称为插值多项式。
为了深入理解插值多项式的意义，我们先来看看多项式的一些基本法则：
首先，多项式是具有如下形式的函数：
p(x) = a0 + a1x + a2x2 + ... + anxn
这里的a0，...，an是系数。当an为非零整数时，这种形式的多项式称为n阶多项式。这是多项式的指数形式，在数学问题中尤为常见。但是，在某些特定的环境中其他形式的多项式则更为简便。比如，在有关多项式插值问题中，牛顿插值多项式就是一个很好的例子：
p(x) = a0 + a1 (x - c1) + a2 (x - c1)(x - c2) + ... + an(x - c1)(x - c2)...(x - cn)
这里a0，...，an是系数，而c0，...，cn是中值。注意到，当c0，...，cn全为0时，牛顿插值多项式就退化为前面定义的n阶多项式。
构建插值多项式
下面我们来看看如何对函数f(x)构建一个插值多项式。
为了对函数f(x)进行插值，要构建一个阶数小于等于n的多项式pn(z)，而这又需要用到函数f(x)的n+1个已知点：x0，...，xn。这些已知点x0，...，xn就称为插值点。通过插值多项式pn(z)可以计算出函数f(x)在x=z处的近似值。插值法需要满足点z在[x0,xn]内。可以采用如下公式来构建插值多项式pn(z)。
pn(z) = f[x0] + f[x0,x1](z-x0) + f[x0,x1,x2](z-x0)(z-x1) + ... + f[x0,...,xn](z-x0)(z-x1)...(z-xn-1)
其中函数f(x)在点x0,...,xn处的值已知，而f[x0],...,f[x0,...,xn]则称为差商。
差商可以通过点x0,...,xn以及函数f(x)在这些点处的值来计算得出。这就是牛顿插值多项式的计算公式。注意到该公式同牛顿插值多项式的相同点。差商的计算公式为：
f[xi,...,xj] = f(xi) 如果i=j
f[xi,...,xj] = ( f[xi+1,...xj] - f[xi,...xj-1] ) / (xj - xi) 如果i < j
通过这个公式不难看出，当i < j时，必须预先计算出其他的差商值。例如，要计算f[x0,x1,x2,x3]，就需要先计算出f[x1,x2,x3]和f[x0,x1,x2]的值。幸运的是，可以通过一个差商表来帮助我们以一种系统的方式来计算差商值。如下图。
差商表由多行组成。最顶行保存已知点x0,...,xn 的值。第二行保存f[x0],...,f[xn]的值。要计算出表中其他的差商值，从每个待求的差商值处画一条对角线，使其回到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（mn2），这里m代表待求值的个数，而n代表已知点的个数。
多项式插值的实现与分析
多项式插值法主要基于对一系列期望点的插值多项式的确定。要得到这个多项式，首先必须通过计算差商来确定该多项式的系数。
首先，为差商以及待确定的系数分配存储空间。注意，由于差商表中每一行的每一项都仅依赖于其前一行的计算结果，因此，并不需要一次性保留所有的表项。所以，只为最占用空间的行分配空间即可，该行将有n个条目。
接下来，用fx中的值来初始化差商表的第一行。这是为计算差商表中的第三行做准备。（前两行不需要计算，因为这两行中的条目都已经保存在x和fx中）。
初始化的最后一步是在coeff[0]中保存fx[0]的值，因为这是插值多项式的第一个系数。
计算差商的过程涉及一个嵌套循环，我们在循环中根据前面介绍过的公式来计算差商。在外层循环中，k用来统计正在计算的是哪一行（排除x和fx所代表的行）。在内层循环中，i表示在当前行中正在计算的是哪一个条目。一旦计算完一行的条目，table[0]中的值就成为插值多项式的下一个系数。因此，保存该值到coeff[k]中。一旦得到插值多项式的所有系数，就可以计算出z中每个目标点的值，最后将这些值保存在pz中。
我们命名这个函数为interpol，它的时间复杂度为O（mn2），这里m代表z中的元素个数（也是pz中值的个数），n代表x中（也是fx中）的元素个数。复杂度因子n2是这样得到的，变更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（n2）。O（mn2）中的因子m来源于针对z中的每个点计算多项式插值的过程。在第一个嵌套循环中，计算出所有的差商，其复杂度为O（n2）。因此，最终的复杂度有一个额外的因子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; }
以上是关于数值计算算法-多项式插值算法的实现与分析的主要内容，如果未能解决你的问题，请参考以下文章