【数值分析】拉格朗日插值法与牛顿插值法的C++实现

数值分析——拉格朗日插值法与牛顿插值法的C++实现

---

---

一、插值法

1.1 插值法定义

设函数$y=f(x)$在区间$[a,b]$上有定义，且$a\le x_0<x_1<\cdots <x_n\le b$，已知在$x_0\dots x_n$点处的值分别为$y_0...y_n$，若存在一个简单的函数$p(x)$，使得$p(x_i)=y_i$，称$p(x)$为$f(x)$的插值函数，点$x_0...x_1$称为插值节点，$[a,b]$称为插值区间,$p(x_i)=y_i$称为插值条件,求插值函数$p(x)$的方法称为插值法。
若$p(x)$为代数不超过n的代数多项式（$p(x)=a_0+a_{1}x+a_2{x}^2+\cdots +a_n{x}^n$），其中$a_i$为实数，称$p(x)$为插值多项式；若 $p(x)$为分段多项式，就是分段插值。若$p(x)$为三角多项式，就称为三角插值。

本质上：根据给出的一系列点$(x_i,f(x_i))$，构造函数$p(x)$来对$f(x)$近似，其中要求$p(x)$必须经过这些插值点。

1.2 插值多项式唯一性定理

满足$p(x_i)=y_i,i=0,...,n$次数不超过 n的插值多项式是唯一存在的。

二、拉格朗日（Lagrange）插值法

2.1 拉格朗日多项式

求n次多项式$L_n(x)=y_0{l}_0(x)+y_1{l}_1(x)+...+y_n{l}_n(x)$使得$L_n(x_i)=y_i,i=0,1,...,n$——即插值多项式一定经过过插值节点$(x_i,y_i),i=0,1,...,n$。$l_i(x)$称为拉格朗日基函数（ Lagrange Basis）。
要求：无重合节点，即$i\ne j$时，$x_i不等于x_j$。
拉格朗日基函数的构造：
$l_i(x)=\prod _{i=0,i\ne j}^n\frac{(x-x_i)}{(x_i-x_j)}$

拉格朗日多项式：
$L_n(x)=l_i(x)y_i$

---

举个例子：
构造三点二次插值（抛物线）多项式：

三点：$(x_0,y_0)、(x_1,y_1)、(x_2,y_2)$

$l_0(x)=\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}$

$l_1(x)=\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}$

$l_2(x)=\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}$

三点二次拉格朗日多项式：$L_2(x)=l_0(x)y_0+l_1(x)y_1+l_2(x)y_2$

2.2 拉格朗日插值法的C++实现

拉格朗日多项式C++实现：

/*
	X:插值节点的x坐标
	Y:插值节点的y坐标
	len:使用插值节点数量（插值多项式阶数为len-1）
	x:待求节点处的x坐标
	输出:拉格朗日多项式在x处的值
	注意：X中坐标值不可重复（否则拉格朗日基函数中分母会出现0项）
*/
double lagrangePolynomial(vector<double> X, vector<double> Y, int len, double x) {
	double res = 0;
	for (int i = 0; i < len; i++) {
		double l_i = 1;//拉格朗日基函数li(x)
		for (int j = 0; j < len; j++) {
			if (j != i) {
				l_i *= (x - X[j]) / (X[i] - X[j]);
			}
		}
		res += l_i * Y[i];
	}
	return res;
}

测试代码：

使用四个点$(\frac{\pi }{6},\frac{1}{2})$、$(\frac{\pi }{4},\frac{1}{\sqrt{2}})$、$(\frac{\pi }{3},\frac{\sqrt{3}}{2})$、$(\frac{\pi }{2},1)$构造三次拉格朗日多项式$L_3(x)$来近似原函数$f(x)=sinx$，并近似求出$sin(\frac{\pi }{5})$的值

#define pi 3.14159265358979323846   //pi
int main() {

	//测试f(x) = sin(x)
	vector<double> X = {pi / 6 , pi / 4 , pi / 3, pi / 2};
	vector<double> Y = { 0.5, 1 / sqrt(2), sqrt(3) / 2, 1};
	int len = X.size();

	cout << "拉格朗日多项式sin(pi/5):" << lagrangePolynomial(X, Y, len, pi / 5) << endl;
	cout << "sin(pi/5)的准确值:" << sin(pi / 5) << endl;
	/*
		拉格朗日多项式sin(pi/5):0.587997
		sin(pi/5)的准确值:0.587785
	*/
}

三、牛顿插值法

拉格朗日插值法当增加一个节点时，所有的拉格朗日基函数$l_i(x)$需要重新进行计算；
使用牛顿插值法，增加一个节点时只需要在原来的牛顿插值函数的基础上加上一项即可。

3.1 差商

一阶差商：
$f[x_i,x_j]=\frac{f(x_i)-f(x_j)}{x_i-x_j}(i\ne j,x_i\ne x_j)$

二阶差商：
$f[x_i,x_j,x_k]=\frac{f[x_i,x_j]-f[x_j,x_k]}{x_i-x_k}(i\ne k)$

k+1阶差商：
$f[x_0,\dots ,x_n]=\frac{f[x_0,\dots ,x_k]-f[x_1,\dots ,x_k,x_{k+1}]}{x_0-x_{k+1}}=\frac{f[x_0,\dots ,x_{k-1},x_k]-f[x_0,\dots ,x_{k-1},x_{k+1}]}{x_k-x_{k+1}}$

差商表：

3.2 牛顿（Newton）插值法

$c_k(x)=f[x_0,x_1,\dots ,x_k]$

$N_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)+\cdots +f[x_0,\dots ,x_n](x-x_0)(x-x_1)\dots (x-x_{n-1})=c_0+c_1(x-x_0)+c_2(x-x_0)(x-x_1)+\cdots +c_n(x-x_0)(x-x_1)\dots (x-x_{n-1})$

且由唯一性定理可知$N_n(x)=L_n(x)$

---

举个例子：

四个点$(0.4,0.41075)$、$(0.55,0.57815)$、$(0.65,0.69675)$、$(0.80,0.88811)$

$f[x_0]=0.41075$、$f[x_1]=0.57815$、$f[x_2]=0.69675$、$f[x_3]=0.88811$

$f[x_0,x_1]=\frac{f[x_0]-f[x_1]}{x_0-x_1}=1.11600$、$f[x_1,x_2]=\frac{f[x_1]-f[x_2]}{x_1-x_2}=1.18600$、$f[x_2,x_3]=\frac{f[x_2]-f[x_3]}{x_2-x_3}=1.27573$

$f[x_0,x_1,x_2]=\frac{f[x_0,x_1]-f[x_1,x_2]}{x_0-x_2}=0.28000$、$f[x_1,x_2,x_3]=\frac{f[x_1,x_2]-f[x_2,x_3]}{x_1-x_3}=0.35893$

$f[x_0,x_1,x_2,x_3]=\frac{f[x_0,x_1,x_2]-f[x_1,x_2,x_3]}{x_0-x_3}=0.19733$

$N_3(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)+f[x_0,x_1,x_2,x_3](x-x_0)(x-x_1)(x-x_2)=0.41075+1.11600(x-0.40)+0.28000(x-0.40)(x-0.55)+0.19733(x-0.40)(x-0.55)(x-0.65)$

3.2 牛顿插值法的C++实现

牛顿插值法C++实现：

/*
	X:插值节点的x坐标
	Y:插值节点的y坐标
	n:使用插值节点数量
	x:待求节点处的x坐标

	输出:牛顿多项式在x处的值

*/
double newtonInterpolation(vector<double> X, vector<double> Y, int n, double x) {
	double res = 0;
	vector<vector<double>> table(X.size(), vector<double>(X.size(), 0));//构建一个完整的差商表[X.size(), X.size()]
	for (int i = 0; i < X.size(); i++) {
		table[0][i] = Y[i];
	}

	//一阶差商
	for (int i = 1; i < X.size(); i++) {
		table[1][i] = (table[0][i] - table[0][i - 1]) / (X[i] - X[i - 1]);//注意一阶差商下(X[i] - X[i - 1])与其它高阶差商(X[i] - X[0])有些区别
	}

	//二阶及以上差商
	for (int i = 2; i < X.size(); i++) {
		for (int j = i; j < X.size(); j++) {
			table[i][j] = (table[i - 1][j] - table[i - 1][j - 1]) / (X[i] - X[0]);
		}
	}

	//计算Newton多项式在x处的近似值(使用len个节点构造len-1次多项式)
	for (int i = 0; i < n; i++) {
		double temp = 1;
		for (int j = 0; j < i; j++) {
			temp *= x - X[j];//计算(x-x_0)...(x-x_i-1)
		}
		res += temp * table[i][i];
	}

	return res;
}

测试代码：

使用四个点$(\frac{\pi }{6},\frac{1}{2})$、$(\frac{\pi }{4},\frac{1}{\sqrt{2}})$、$(\frac{\pi }{3},\frac{\sqrt{3}}{2})$、$(\frac{\pi }{2},1)$构造差商表table，再使用len个节点构造len-1阶牛顿多项式来近似原函数$f(x)=sinx$，并近似求出$sin(\frac{\pi }{5})$的值

#define pi 3.14159265358979323846   //pi
int main() {

	//测试f(x) = sin(x)
	vector<double> X = {pi / 6 , pi / 4 , pi / 3, pi / 2};
	vector<double> Y = { 0.5, 1 / sqrt(2), sqrt(3) / 2, 1};
	int n = X.size();

	//使用前三个节点
	cout << "牛顿多项式sin(pi/5):" << newtonInterpolation(X, Y, n - 1, pi / 5) << endl;
	cout << "sin(pi/5)的准确值:" << sin(pi / 5) << endl;
	/*
		牛顿多项式sin(pi/5):0.588625
		sin(pi/5)的准确值:0.587785 
	*/

	//使用前四个节点
	cout << "牛顿多项式sin(pi/5):" << newtonInterpolation(X, Y, n, pi / 5) << endl;
	cout << "sin(pi/5)的准确值:" << sin(pi / 5) << endl;
	/*
		牛顿多项式sin(pi/5):0.586526
		sin(pi/5)的准确值:0.587785
	*/
	
}

总结

插值法就是利用一系列的点来构造插值函数对原函数进行近似，构造的插值函数必须满足插值节点的要求。拉格朗日插值法构造拉格朗日基函数较为简单，但是当增加或减少插值节点时拉格朗日基函数就需要重新计算；而牛顿插值法利用差商表（有点类似动态规划）来构建插值函数，当函数增加一个节点时只需要在原来的基础上增加一项即可。