牛顿插值法

差商定义：

f(x)在  点处的0阶差商  

f(x)在  点处的1阶差商  

f(x)在  点处的1阶差商  

f(x)在  点处的2阶差商  

......

f(x)在  点处的n阶差商  

 

差商表：

		1阶差商	2阶差商	3阶差商	4阶差商	......
......	......	......	......	......	......	......

 

牛顿插值函数:

过 n + 1 个插值节点的n次牛顿插值函数的一般形式为： 

 

所以利用牛顿插值法首先需要计算差商，再求近似值。在代入插值公式时，可以采用秦九韶算法提取公因数进行算法优化。另外牛顿插值函数比拉格朗日插值函数更具有优势。因为当增加新插值节点时，牛顿插值函数能避免重复计算。

 

代码：

#include <stdio.h>

const int maxn = 50;

// 计算差商 
void DifferenceQuotient(int n, double *x, double **f){
	int k;
	for(int j = 1; j < n; j++){                    // 第 j 列
		for(int i = j; i < n; i++){            // 第 i 行
			f[i][j] = (f[i][j-1] - f[i-1][j])/(x[i]-x[i-j]);
		}
	}
}

// 秦九韶算法 
double Newton(int n, double *x, double **f, double _x){
	double _y = f[n-1][n-1];
	for(int i = n-2; i >= 0; i--){
		_y = _y * (_x - x[i]) + f[i][i];
	}
	return _y;
}

int main(){
	double x[maxn], _x;
	double **f = new double *[maxn];
	for(int i = 0; i < maxn; i++) f[i] = new double[maxn];
	int n;
	printf("请输入插值结点的个数：");
	scanf("%d",&n);
	for(int i = 0; i < n; i++){
		printf("请输入插值结点x%d,y%d：",i,i);
		scanf("%lf %lf",&x[i],&f[i][0]);
	}
	DifferenceQuotient(n,x,f);
	printf("请输入_x:"); 
	scanf("%lf",&_x);
	printf("N(%lf)的近似值:%lf",_x,Newton(n,x,f,_x));
	delete f;
	return 0;
}

空间优化：牛顿插值公式，只用了差商表的各列的最上面一个差商，因此可以只存储各列的最上面的差商，使空间复杂度从O(n^2)降为O(n)。但增加新节点时，不能重复利用之前的差商，需要重新计算差商表。

 

思路：用一个一维数组f[n]存放差商表，f[i] 表示第 i 阶差商。初始化时，0阶差商等于函数本身。即，f[i] = f(xi)。因计算新的差商时，需要用到旧的差商。所以采取从下往上计算的方式，保证在计算新差商时，旧差商没有被覆盖。由于计算过程会保证 i >= j，所以新的差商f[i+1],不会覆盖旧差商 f[i]。当计算结束后，所以阶的差商也就被计算出来了。

 

代码：

// 优化
for(int j = 1; j < n; j++){
    for(int i = n-1; i >= j; i--){                // 自下往上计算
        f[i] = (f[i] - f[i-1])/(x[i] - x[i-j]);
    }  
}

 

测试用例:

3
2.56 1.6
2.89 1.7
3.24 1.8
3
近似值：1.732

 

测试结果：