正文

牛顿插值算法与实现2006-09-15 00:09:00

【评论】 【打印】 【字体: 】 本文链接:http://blog.pfan.cn/rickone/18540.html

分享到:

牛顿真是牛,拉格朗日插值法只能算是数学意义上的插值,从插值基函数的巧妙选取,已经构造性的证明了插值法的存在性和惟一性,但是从实现的角度看并不很好,而牛顿很好的解决了这个问题。

牛顿插值是基于下面这些的公式:

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[i])<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[i])+f[0][i];
 }
 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;
}

rickone 2006/09/15

阅读(18612) | 评论(3)


版权声明:编程爱好者网站为此博客服务提供商,如本文牵涉到版权问题,编程爱好者网站不承担相关责任,如有版权问题请直接与本文作者联系解决。谢谢!

评论

loading...
您需要登录后才能评论,请 登录 或者 注册