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

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

拉格朗日系数多项式

假设$n$个点$(x_0,y_0),(x_1,y_1),\cdots ,(x_{n-1},y_{n-1})$，并且有$x_0<x_1<\cdots <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})}$.
很显然$L_k(x)$具有下面的性质：
$$\begin{array}{r}{L}_k(x_i)=\{\begin{array}{r}1,i=k\\ 0,i\ne k\end{array}\end{array}$$
而且$L_k(x)$只跟横坐标有关，与纵坐标无关，有时也叫做拉格朗日基函数。
根据拉格朗日系数多项式，可以快速求出过点$(x_0,y_0),(x_1,y_1),\cdots ,(x_{n-1},y_{n-1})$的多项式$p(x)$.
$$p(x)=y_0{L}_0(x)+y_1{L}_1(x)+\cdots +y_{n-1}{L}_{n-1}(x).$$
根据$L_k(x)$的性质，$p(x)$显然过点$(x_0,y_0),(x_1,y_1),\cdots ,(x_{n-1},y_{n-1})$。

c++演示代码

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

#include<iostream>
#include<vector>
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;
	}
}