已知 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_0<x_1<\cdots <x_{n-1} x0<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)=y0L0(x)+y1L1(x)+⋯+yn−1Ln−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<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;
}
}