1. 程式人生 > >Neville 演算法解多項式插值

Neville 演算法解多項式插值

原理圖如下:

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

{{OVER}}