计算方法——牛顿插值方法的C++实现

牛顿插值多项式:

差商表：

牛顿插值多项式中的表示差商表中的项，具体里是0至n，就是差商表中的坐标项。

差商表的计算：

需要按列计算，即第一列、第二列按这个顺序计算。下面以已知插值节点是为例，计算差商表。

第0列：

	0
0	y[0]=0.41075
1	y[1]=0.57815
2	y[2]=0.69675
3	y[3]=0.88811
4	y[4]=1.02652

第1列：

x		0,a	1,b
x[0]=0.4	0	y[0]=0.41075
x[1]=0.55	1	y[1]=0.57815	b[1]=(y[1]-y[0])/(x[1]-x[0])
x[2]=0.65	2	y[2]=0.69675	b[2]=(y[2]-y[1])/(x[2]-x[1])
x[3]=0.8	3	y[3]=0.88811	b[3]=...
x[4]=0.9	4	y[4]=1.02652	b[4]=...

第2列：

x		0,a	0,b	0,c
x[0]=0.4	0	y[0]=0.41075
x[1]=0.55	1	y[1]=0.57815	b[1]=(y[1]-y[0])/(x[1]-x[0])
x[2]=0.65	2	y[2]=0.69675	b[2]=(y[2]-y[1])/(x[2]-x[1])	c[2]=(b[2]-b[1])/(x[2]-x[0])
x[3]=0.8	3	y[3]=0.88811	b[3]=...	c[3]=(b[3]-b[2])/(x[3]-x[1])
x[4]=0.9	4	y[4]=1.02652	b[4]=...	c[4]=...

第3，4列注意下标之间的规律后可自行计算

x		0,a	1,b	2,c	3,d	4,e
x[0]=0.4	0	y[0]=0.41075
x[1]=0.55	1	y[1]=0.57815	b[1]=(y[1]-y[0])/(x[1]-x[0])
x[2]=0.65	2	y[2]=0.69675	b[2]=(y[2]-y[1])/(x[2]-x[1])	c[2]=(b[2]-b[1])/(x[2]-x[0])
x[3]=0.8	3	y[3]=0.88811	b[3]=...	c[3]=(b[3]-b[2])/(x[3]-x[1])	d[3]=(c[3]-c[2])/(x[3]-x[0])
x[4]=0.9	4	y[4]=1.02652	b[4]=...	c[4]=...	d[4]=(c[4]-c[3])/(x[4]-x[1])	e[4]=(d[4]-d[3])/(x[4]-x[0])

        从上面的计算方法中我们可以看出，要想编程实现，可以用一个二维数组（这里设为z[n][n]）存储差商表的数值，将数组初始化第0列为y的值后，可以按表格规律z[j][i]=(z[j-1][i-1]-z[j][i-1])/(x[j]-x[i-j])，计算出后面其他列的值。

        最后把差商表中的值带入牛顿插值多项式公式，计算出x取值下插值多项式的值。

代码：

#include<iostream>
#define N 5
using namespace std;
int main(){
	int i,j;
	float varx =0.895,yx=0;
	float x[N]={0.4,0.55,0.65,0.8,0.9};
	float y[N]={0.41075,0.57815,0.69675,0.88811,1.02652};
	float z[N][N]={0};
//  计算差商表
//  初始化第一列
	for(i=0;i<N;i++){
		z[i][0]=y[i];
	}
//  计算后续列
	for(i=1;i<N;i++){
		for(j=i;j<N;j++){
			z[j][i]=(z[j][i-1]-z[j-1][i-1])/(x[j]-x[j-i]);
		}
	}
//  利用差商表中有用的值，计算插值多项式的值
	for(i=0;i<N;i++){
		float t=1;
		for(j=0;j<i;j++)
			t*=(varx-x[j]);	
		yx+=z[i][i]*t;
	}
	cout<<"Nn("<<varx<<")="<<yx<<endl;
}