牛顿插值算法与实现.doc

牛顿插值算法与实现牛顿真是牛，拉格朗日插值法只能算是数学意义上的插值，从插值基函数的巧妙选取，已经构造性的证明了插值法的存在性和惟一性，但是从实现的角度看并不很好，而牛顿很好的解决了这个问题。牛顿插值是基于下面这些的公式：fx0,x1,.xk=(fx1,.xk-fx0,.xk-1)/(xk-x0)fx=f(x)f(x)=fx0+fx0,x1(x-x0)+fx0,x1,x2(x-x0)(x-x1)+.fx0,.xn(x-x0).(x-xn-1)+Rn(x)前两个是均差的递推关系式，而后一个就是牛顿插值公式，其中N(x)=f(x)-Rn(x)，即目标多项式，Rn(x)是n阶插值余项，我们就是用N(x)去近似f(x)。可以构造这样一个均方差表：xk f(xk) 一阶均差 二阶均差 .x0 f(x0)x1 f(x1) fx0,x1x2 f(x2) fx1,x2 fx0,x1,x2.如果有n个点插值，表会有(n*n)/2+n个表项，如果直接编程会有O(n*n)的空间复杂度，编程时做个简单的改进，不难发现在这个表中只有部分数据有用，对角线（斜行）它们是目标值，用来表示多项式的，左边的两纵行（实际上只需要x一行）以及最底下的一行，表示当前插值的状态。经过改进后只需要O(n)的空间复杂度。两个过程：1，新增加一个点时的更新。只须更新最底下一行数据，其递推关系由均差公式给出，最后算出高一队的均差值，需时O(n)2，插入点完成后如何计算多项式在另外给定点的值N(x)。由牛顿插值公式，最终的表达式为：N(x)=fx0+fx0,x1(x-x0)+fx0,x1,x2(x-x0)(x-x1)+.fx0,.xn(x-x0).(x-xn-1)如果直接将它展开，再算实在麻烦，实际上大可不必这样做，还记得多项式求值的秦九韶算法吗？将多项式叠起来，从内层括号往外一层层拨开，n次多项多的计算，只需要做n次乘法，同样的思想，将上式改写成：N(x)=fx0+(x-x0)fx0,x1+(x-x1)fx0,x1,x2+(x-x2).fx0,.xn-1+(x-xn-1)fx0,.xn.就可以同样简单的计算了，时间复杂度O(n)综合起来的性能：对于n个点的插值，产生多项式的时间复杂度是O(n*n)，最终进行一个点的计算的时间复杂度是O(n)。C+代码实现/ file: newton.h#ifndef NEWTON_DEF_#define NEWTON_DEF_class CNewtondouble *f2;double *x;int max;int n;public:CNewton(int MaxN);/MaxN 为最大插值点数 可任意设定CNewton();void InsertPoint(double X,double Y);double GetValue(double X);#endif/ file: newton.cpp#include newton.h#include assert.h#include math.h#ifndef NULL#define NULL 0#endifCNewton:CNewton(int MaxN)max=MaxN+1;n=0;x=new doublemax;f0=new doublemax;f1=new doublemax;assert(x!=NULL);assert(f0!=NULL);assert(f1!=NULL);CNewton:CNewton()if(x)deletex;if(f0)deletef0;if(f1)deletef1;void CNewton:InsertPoint(double X,double Y)int i;double fw;assert(nmax);/重复点检查for(i=0;in;+i)if(fabs(X-xi)1e-5)return;/如果确保不会有重复点可删去上面语句xn=X;fw=Y;for(i=1;i=0;-i)s=s*(X-xi)+f0i;return s;/ file: test cpp#include newton.h#include iostream.hint main(void)int n;double x,y;CNewton nt(20);cout输入插入点个数(nn;for(int i=1;i=n;+i)cout输入第ix;couty;nt.Insert