拉格朗日多项式插值及其c++演示

已知 n n n个点，寻找一条穿过所有点的曲线的过程叫做插值。如果这 n n n个点在某一条多项式的曲线上，那么，寻找这个多项式的解析式的过程就是多项式插值。
拉格朗日（Lagrange）插值是一种快速求解穿过 n n n个点的多项式的方法。这 n n n个点必须满足其横坐标两两不相同，否则不存在过这 n n n个点的多项式。而且，满足条件的 n − 1 n-1 n−1次多项式只有一个。
最简单的例子就是两点确定一条直线，三个点确定一条抛物线。

拉格朗日系数多项式

假设 n n n个点 ( x 0 , y 0 ) , ( x 1 , y 1 ) , ⋯   , ( x n − 1 , y n − 1 ) (x_0,y_0),(x_1,y_1),\cdots,(x_{n-1},y_{n-1}) (x0​,y0​),(x1​,y1​),⋯,(xn−1​,yn−1​)，并且有 x 0 < x 1 < ⋯ < x n − 1 x_0x0​<x1​<⋯<xn−1​.
定义函数 L k ( x ) = ( x − x 0 ) ( x − x 1 ) ⋯ ( x − x k − 1 ) ( x − x k + 1 ) ⋯ ( x − x n − 1 ) ( x k − x 0 ) ( x k − x 1 ) ⋯ ( x k − x k − 1 ) ( x k − x k + 1 ) ⋯ ( x k − x n − 1 ) L_k(x)=\frac{(x-x_0)(x-x_1)\cdots(x-x_{k-1})(x-x_{k+1})\cdots(x-x_{n-1})}{(x_k-x_0)(x_k-x_1)\cdots(x_k-x_{k-1})(x_k-x_{k+1})\cdots(x_k-x_{n-1})} Lk​(x)=(xk​−x0​)(xk​−x1​)⋯(xk​−xk−1​)(xk​−xk+1​)⋯(xk​−xn−1​)(x−x0​)(x−x1​)⋯(x−xk−1​)(x−xk+1​)⋯(x−xn−1​)​.
很显然 L k ( x ) L_k(x) Lk​(x)具有下面的性质：
L k ( x i ) = { 1 , i = k 0 , i ≠ k \begin{aligned} L_k(x_i)=\left \{ \begin{aligned} 1,i=k\ 0,i\neq k \end{aligned} \right. \end{aligned} Lk​(xi​)={1,i=k0,i​=k​​
而且 L k ( x ) L_k(x) Lk​(x)只跟横坐标有关，与纵坐标无关，有时也叫做拉格朗日基函数。
根据拉格朗日系数多项式，可以快速求出过点 ( x 0 , y 0 ) , ( x 1 , y 1 ) , ⋯   , ( x n − 1 , y n − 1 ) (x_0,y_0),(x_1,y_1),\cdots,(x_{n-1},y_{n-1}) (x0​,y0​),(x1​,y1​),⋯,(xn−1​,yn−1​)的多项式 p ( x ) p(x) p(x).
p ( x ) = y 0 L 0 ( x ) + y 1 L 1 ( x ) + ⋯ + y n − 1 L n − 1 ( x ) . p(x)=y_0L_0(x)+y_1L_1(x)+\cdots+y_{n-1}L_{n-1}(x). p(x)=y0​L0​(x)+y1​L1​(x)+⋯+yn−1​Ln−1​(x).
根据 L k ( x ) L_k(x) Lk​(x)的性质， p ( x ) p(x) p(x)显然过点 ( x 0 , y 0 ) , ( x 1 , y 1 ) , ⋯   , ( x n − 1 , y n − 1 ) (x_0,y_0),(x_1,y_1),\cdots,(x_{n-1},y_{n-1}) (x0​,y0​),(x1​,y1​),⋯,(xn−1​,yn−1​)。

c++演示代码

代码中使用系数表示法表示一个多项式，数组下标表示次数，对应的值是该次数的系数。如 1 + 2 x + 3 x 2 + 4 x 3 1+2x+3x^2+4x^3 1+2x+3x2+4x3表示为 [ 1 , 2 , 3 , 4 ] [1,2,3,4] [1,2,3,4].

#include
#include
using namespace std;
class Lagrange{
	vector<vector<double>> points;// a point is (x,y)
	public:
		void setPoints(vector<vector<double>> points){
			this->points=points;
		}
		vector<double> interpolate(){
			vector<vector<double>> ls(points.size());
			for(int i=0;i<ls.size();i++){
				ls[i]=Lk(i);
				cout<<"The L"<<i<<"(x) is "<<endl; 
				for(auto it:ls[i]){
					cout<<it<<' ';
				}
				cout<<endl;
			}
			vector<double> res(points.size(),0);
			for(int i=0;i<ls.size();i++){
				for(int j=0;j<res.size();j++){
					res[j]+=points[i][1]*ls[i][j];
				}
			}
			return res;
		}
		vector<double> Lk(int k){
			double denominator=1;
			for(int i=0;i<points.size();i++){
				if(i==k){
					continue;
				}
				denominator*=points[k][0]-points[i][0];
			}
			vector<double> res(points.size(),0);
			Lk_help(k,0,0,1,res);
			for(int i=0;i<res.size();i++){
				res[i]=res[i]/denominator;
			}
			return res;
		}
		void Lk_help(int k,int i,int exp,double val,vector<double>& l){
			if(i>=points.size()){
				l[exp]+=val;
				return;
			}
			if(i==k){
				Lk_help(k,i+1,exp,val,l);
				return;	
			}
			Lk_help(k,i+1,exp,-val*points[i][0],l);
			Lk_help(k,i+1,exp+1,val,l);
		}
};
int main(){
	Lagrange ll;
	vector<vector<double>> p={{0,1},{1.2,0.362358},{2.6,34},{3.2,7.8}};
	ll.setPoints(p);
	vector<double> res=ll.interpolate();
	cout<<"The result is"<<endl;
	for(auto i:res){
		cout<<i<<' ';
	}
	cout<<endl<<endl;;
	for(int i=0;i<p.size();i++){
		double y=0;
		for(int k=res.size()-1;k>=0;k--){
			y=y*p[i][0]+res[k];
		}
		cout<<"The original y is "<<p[i][1]<<", computed y is "<<y<<endl;
	}
}