Welcome to ErranLi's Blog!

牛頓真是牛，拉格朗日插值法只能算是數學意義上的插值，從插值基函數的巧妙選取，已經構造性的證明了插值法的存在性和惟一性，但是從實現的角度看并不很好，而牛頓很好的解決了這個問題。

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

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;
}