牛顿插值是基于下面这些的公式：
f[x0,x1,...xk]=(f[x1,...xk]-f[x0,...xk-1])/(xk-x0)
f[x]=f(x)
f(x)=f[x0]+f[x0,x1](x-x0)+f[x0,x1,x2](x-x0)(x-x1)+...f[x0,...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) f[x0,x1]
x2 f(x2) f[x1,x2] f[x0,x1,x2]
...
如果有n个点插值，表会有(n*n)/2+n个表项，如果直接编程会有O(n*n)的空间复杂度，编程时做个简单的改进，不难发现在这个表中只有部分数据有用，对角线（斜行）它们是目标值，用来表示多项式的，左边的两纵行（实际上只需要x一行）以及最底下的一行，表示当前插值的状态。经过改进后只需要O(n)的空间复杂度。
两个过程：
1，新增加一个点时的更新。只须更新最底下一行数据，其递推关系由均差公式给出，最后算出高一队的均差值，需时O(n)
2，插入点完成后如何计算多项式在另外给定点的值N(x)。
由牛顿插值公式，最终的表达式为：
N(x)=f[x0]+f[x0,x1](x-x0)+f[x0,x1,x2](x-x0)(x-x1)+...f[x0,...xn](x-x0)...(x-xn-1)
如果直接将它展开，再算实在麻烦，实际上大可不必这样做，还记得多项式求值的秦九韶算法吗？将多项式‘叠’起来，从内层括号往外一层层拨开，n次多项多的计算，只需要做n次乘法，同样的思想，将上式改写成：
N(x)=f[x0]+(x-x0){f[x0,x1]+(x-x1){f[x0,x1,x2]+(x-x2){...{f[x0,...xn-1]+(x-xn-1)f[x0,...xn]}...}
就可以同样简单的计算了，时间复杂度O(n)
综合起来的性能：对于n个点的插值，产生多项式的时间复杂度是O(n*n)，最终进行一个点的计算的时间复杂度是O(n)。
C++代码实现
// file: newton.h
#ifndef NEWTON_DEF_
#define NEWTON_DEF_
class CNewton
{
double *f[2];
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
#endif
CNewton::CNewton(int MaxN)
{
max=MaxN+1;
n=0;
x=new double[max];
f[0]=new double[max];
f[1]=new double[max];
assert(x!=NULL);
assert(f[0]!=NULL);
assert(f[1]!=NULL);
}
CNewton::~CNewton()
{
if(x)
delete[]x;
if(f[0])
delete[]f[0];
if(f[1])
delete[]f[1];
}
void CNewton::InsertPoint(double X,double Y)
{
int i;
double fw;
assert(n<max);
//重复点检查
for(i=0;i<n;++i)
if(fabs(X-x)<1e-5)
return;
//如果确保不会有重复点可删去上面语句
x[n]=X;
fw=Y;
for(i=1;i<=n;++i)
{
double tmp=fw;
fw=(fw-f[1][i-1])/(x[n]-x[n-i]);
f[1][i-1]=tmp;
}
f[0][n]=f[1][n]=fw;
n++;
}
double CNewton::GetValue(double X)
{
if(n==0)
return 0.0;
double s=f[0][n-1];
for(int i=n-2;i>=0;--i)
{
s=s*(X-x)+f[0];
}
return s;
}
// file: test cpp
#include "newton.h"
#include "iostream.h"
int main(void)
{
int n;
double x,y;
CNewton nt(20);
cout<<"输入插入点个数(n<=20)\nn=";
cin>>n;
for(int i=1;i<=n;++i)
{
cout<<"输入第"<<i<<"个点\nx=";
cin>>x;
cout<<"y=";
cin>>y;
nt.InsertPoint(x,y);
}
while(1)
{
cout<<"计算N(x)\nx=";
cin>>x;
cout<<"N("<<x<<")=\n"<<nt.GetValue(x)<<endl;
if(x==0.0)
break;
}
return 0;
}
最新评论