Neville 演算法解多項式插值

原文連結：http://blog.csdn.net/zhangxaochen/article/details/8100012

原理圖如下：

 

Numerical Recipes 隨書附帶的代碼：（xa, ya 是n個樣本點的座標值數組， x 是待求點的橫座標， 輸出值為 y， dy， 其中dy 表示誤差）

void NR::polint(Vec_I_DP &xa, Vec_I_DP &ya, const DP x, DP &y, DP &dy){ int i,m,ns=0; DP den,dif,dift,ho,hp,w; int n=xa.size(); Vec_DP c(n),d(n); dif=fabs(x-xa[0]); for (i=0;i<n;i++) {  if ((dift=fabs(x-xa[i])) < dif) {   ns=i;   dif=dift;  }  c[i]=ya[i];  d[i]=ya[i]; } y=ya[ns--]; for (m=1;m<n;m++) {  for (i=0;i<n-m;i++) {   ho=xa[i]-x;   hp=xa[i+m]-x;   w=c[i+1]-d[i];   if ((den=ho-hp) == 0.0) nrerror("Error in routine polint");   den=w/den;   d[i]=hp*den;   c[i]=ho*den;  }  y += (dy=(2*(ns+1) < (n-m) ? c[ns+1] : d[ns--])); }}

其中有兩個數組 c[n], d[n], 分別表示後一列相對前一列的增量（因為有兩個分支），如上面第一張圖所示，事實上我一直不明白弄這個增量有什麼深意，好像沒有降低時間複雜度， 反而增加了編程的複雜度。 如果有前輩恰好瞭解這個深意， 盼賜教哇。
個人實現的代碼就比較簡單：

void MyPolint(const double* xa, const double* ya, const int n, const double x, double& y){ double* p=new double[n]; for(int i=0; i<n; i++){  p[i]=ya[i]; } for(int k=1; k<n; k++){  for(int i=0; i<n-k; i++){   double factor1=x-xa[i+k],    factor2=xa[i]-x;   p[i]=(factor1*p[i]+factor2*p[i+1])/(factor2+factor1);  } } y=p[0];}

測試代碼：

#include "nr.h"using namespace std;using namespace NR;void MyPolint(const double* xa, const double* ya, const int n, const double x, double& y);void main(){ const int n=4; DP xarr[n]={-1, 0, 1, 2}; //DP yarr[n]={-2, 0, 2, 10};//曲線為： y=x^3+x DP yarr[n]={1, 0, -1, 4}; //曲線為： y=x^3-2x Vec_I_DP xa(xarr, n); Vec_I_DP ya(yarr, n); DP x=1.5; DP y=INT_MAX,  dy=INT_MAX; polint(xa, ya, x, y, dy); printf("y, dy: %lf, %lf\n", y, dy); // 4.875， 1.875， √ //-------------我的函數 MyPolint(xarr, yarr, n, x, y); printf("y: %lf\n", y);}

輸出結果:

y, dy: 0.375000, 1.875000
y: 0.375000

因為個人也沒大理解那個 dy憑什麼就是誤差，所以自己實現的函數裡就沒算這個 dy

 

原文連結：http://blog.csdn.net/zhangxaochen/article/details/8100012

{{OVER}}