牛顿插值的java实现_转：牛顿插值算法与实现

原文：http://blog.csdn.net/nhczp/archive/2007/01/31/1498826.aspx

作者：

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

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

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

//重复点检查

for(i=0;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<

cin>>n;

for(int i=1;i<=n;++i)

{

cout<

cin>>x;

cout<

cin>>y;

nt.InsertPoint(x,y);

}

while(1)

{

cout<

cin>>x;

cout<

if(x==0.0)

break;

}

return 0;

}