牛頓真是牛，拉格朗日插值法只能算是數(shù)學(xué)意義上的插值，從插值基函數(shù)的巧妙選取，已經(jīng)構(gòu)造性的證明了插值法的存在性和惟一性，但是從實(shí)現(xiàn)的角度看并不很好，而牛頓很好的解決了這個(gè)問題。

牛頓插值是基于下面這些的公式：

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)

前兩個(gè)是均差的遞推關(guān)系式，而后一個(gè)就是牛頓插值公式，其中N(x)=f(x)-Rn(x)，即目標(biāo)多項(xiàng)式，Rn(x)是n階插值余項(xiàng)，我們就是用N(x)去近似f(x)。

可以構(gòu)造這樣一個(gè)均方差表：

xk?? f(xk)?? 一階均差?? 二階均差 ...
x0?? f(x0)
x1?? f(x1)???? f[x0,x1]
x2?? f(x2)???? f[x1,x2]???? f[x0,x1,x2]
...

如 果有n個(gè)點(diǎn)插值，表會(huì)有(n*n)/2+n個(gè)表項(xiàng)，如果直接編程會(huì)有O(n*n)的空間復(fù)雜度，編程時(shí)做個(gè)簡(jiǎn)單的改進(jìn)，不難發(fā)現(xiàn)在這個(gè)表中只有部分?jǐn)?shù)據(jù)有 用，對(duì)角線（斜行）它們是目標(biāo)值，用來表示多項(xiàng)式的，左邊的兩縱行（實(shí)際上只需要x一行）以及最底下的一行，表示當(dāng)前插值的狀態(tài)。經(jīng)過改進(jìn)后只需要O (n)的空間復(fù)雜度。

兩個(gè)過程：
1，新增加一個(gè)點(diǎn)時(shí)的更新。只須更新最底下一行數(shù)據(jù)，其遞推關(guān)系由均差公式給出，最后算出高一隊(duì)的均差值，需時(shí)O(n)
2，插入點(diǎn)完成后如何計(jì)算多項(xiàng)式在另外給定點(diǎn)的值N(x)。
由牛頓插值公式，最終的表達(dá)式為：
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)
如果直接將它展開，再算實(shí)在麻煩，實(shí)際上大可不必這樣做，還記得多項(xiàng)式求值的秦九韶算法嗎？將多項(xiàng)式‘疊’起來，從內(nèi)層括號(hào)往外一層層撥開，n次多項(xiàng)多的計(jì)算，只需要做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]}...}
就可以同樣簡(jiǎn)單的計(jì)算了，時(shí)間復(fù)雜度O(n)

綜合起來的性能：對(duì)于n個(gè)點(diǎn)的插值，產(chǎn)生多項(xiàng)式的時(shí)間復(fù)雜度是O(n*n)，最終進(jìn)行一個(gè)點(diǎn)的計(jì)算的時(shí)間復(fù)雜度是O(n)。

C++代碼實(shí)現(xiàn)
// 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 為最大插值點(diǎn)數(shù) 可任意設(shè)定
?~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);
?//重復(fù)點(diǎn)檢查
?for(i=0;i<n;++i)
??if(fabs(X-x[i])<1e-5)
???return;
?//如果確保不會(huì)有重復(fù)點(diǎn)可刪去上面語句
?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<<"輸入插入點(diǎn)個(gè)數(shù)(n<=20)\nn=";
?cin>>n;
?for(int i=1;i<=n;++i)
?{
??cout<<"輸入第"<<i<<"個(gè)點(diǎn)\nx=";
??cin>>x;
??cout<<"y=";
??cin>>y;
??nt.InsertPoint(x,y);
?}
?while(1)
?{
??cout<<"計(jì)算N(x)\nx=";
??cin>>x;
??cout<<"N("<<x<<")=\n"<<nt.GetValue(x)<<endl;
??if(x==0.0)
???break;
?}
?return 0;
}