牛顿插值法及其代码

简介

Lagrange 插值方法虽然易算，但若要增加一个节点时，全部基函数 φ(x) 都需重新算过。本节介绍另外一种插值方法-牛顿插值法。牛顿插值法是一种重要的插值方法，与拉格朗日插值法在同阶时产生的多项式在化简以后是一样的，余项也是一样的，因此两者只是写法形式不同而已。在针对变化的数据时，其各自独特的写法形式给其适合的应用场景带来了区别。牛顿插值法适用于被插的点不断增多，而且可以很好的复用以前的计算结果。

来源

我们知道 Lagrange 插值多项式的插值基函数为：
$${\phi }_i\left(x\right):=\prod _{j=0,j\ne i}^n\frac{x-x_j}{x_i-x_j}(i=0,1,2,3,\cdot \cdot \cdot ,n)$$
上述基函数形式上比较复杂、计算量大，并且重复计算很多，由线性代数可知任何一个 n 次多项式都可以表示成：
$$1,(x-x_0),(x-x_0)(x-x_1),\cdot \cdot \cdot ,(x-x_0)(x-x_1)\cdot \cdot \cdot (x-x_{n-1})$$
共 n + 1 个多项式的线性组合。那么，是否可以将这 n + 1 个多项式作为插值基函数呢？答案是肯定的。多项式组$1,(x-x_0),(x-x_0)(x-x_1),\cdot \cdot \cdot ,(x-x_0)(x-x_1)\cdot \cdot \cdot (x-x_{n-1})$线性无关，因此可以作为插值基函数。

设已知 n + 1 个不同节点 $x_0,x_1,\cdot \cdot \cdot ,x_n$ 上的函数值$f(x_0)=y_0,f(x_1)=y_1,\cdot \cdot \cdot ,f(x_n)=y_n$构造满足插值条件 $p_n(x_i)=f(x_i)=y_i(i=0,1,\cdot \cdot \cdot ,n)$ 的次数不超过 n 次多项式。

$$\begin{gathered} p_n(x)=a_0+a_1(x-x_0)+a_2(x-x_0)(x-x_1)+a_3(x-x_0)(x-x_1)(x-x_2)+\cdot \cdot \cdot + \\ {a}_n(x-x_0)(x-x_1)(x-x_2)\cdot \cdot \cdot (x-x_{n-1}) \end{gathered}$$

为了便于计算系数值$a_0,a_1,a_2,\cdot \cdot \cdot ,a_n$，我们引入差商的概念。

我们将差商的值作为牛顿插值公式的系数
利用差商的递推定义，可以用递推来计算差商。

如此，我们便得到了牛顿插值公式的基函数以及它的系数，如此我们就可以写出来牛顿插值公式的函数，其 n 次插值多项式为：
$$\begin{gathered} p_n(x)=f(x_0)+f(x_0,x_1)(x-x_0)+f(x_0,x_1,x_2)(x-x_0)(x-x_1)+\cdot \cdot \cdot + \\ f(x_0,x_1,\cdot \cdot \cdot ,x_n)(x-x_0)(x-x_1)\cdot \cdot \cdot (x-x_{n-1}) \end{gathered}$$

实例

例：已知节点 0, 2, 3, 5 对应的函数值为 1, 3, 2, 5，构造三次Newton 插值多项式，当 x = 2.5 时，计算 Newton 多项式插值结果。
解：由前面计算的差商值，我们很容易写出三次 Newton 插值多项式如下：
$$\begin{gathered} p3(x)=f(x0)+f(x0,x1)(x-x0)+f(x0,x1,x2)(x-x0)(x-x1)+ \\ f(x0,x1,x2,x3)(x-x0)(x-x1)(x-x2) \\ =1+1\times (x-0)-\frac{2}{3}\times (x-0)(x-2)+\frac{3}{10}\times (x-0)(x-2)(x-3) \\ =\frac{3}{10}\times x^3-\frac{13}{6}\times x^2+\frac{62}{15}\times x+1 \end{gathered}$$
当 x = 2.5 时，Newton 多项式插值结果为：
$$p3(2.5)=\frac{3}{10}\times 2.5^3-\frac{13}{6}\times 2.5^2+\frac{62}{15}\times 2.5+1=2.479167$$

代码

#include<iostream>
using namespace std;
double n,x[100],y[100],l[100][100] ;
void input()
{
	for(int i=0;i<n;i++) 
		cin>>x[i]>>y[i]; 
} 

double newton(double test)
{
	double re = 0;
	for(int i=0;i<100;i++) l[i][0]=y[i]; 
	for(int j =1;j<n;j++)
	{ 
		for(int i=j;i<n;i++)
		{
			
			l[i][j] = (l[i][j-1]-l[i-1][j-1])/(x[i]-x[i-j]);
			cout<<i<<" "<<j<<" "<<l[i][j]<<endl; 
		}	
	} 
	double g[100];
	for(int i=0;i<100;i++) g[i]=1.0;
	for(int i=0;i<n;i++)
	{ 
		for(int j=0;j<i;j++)
		{
			g[i] *= (test -x[j]);
		} 
		re += g[i]*l[i][i];
	}
	return re;
}
 

int main()
{
//	for(int i=0;i<100;i++) cout<<l[i]<<endl;
	cout<<"输入多少个点"<<endl;
	cin>>n;
	cout<<"输入插值点："<<endl;
	input(); 
	double test;
	while(1)
	{
		cout<<"输入要计算的值或输入 -1 退出"<<endl; 
		cin>>test;
		if(test == -1) break;
		double re = newton(test);
		cout<<"result: "<<re <<endl;
	}
	return 0;
}

结果：