拉格朗日插值

拉格朗日插值

——------------------------------------------------------------------------------------------------

定义

对于一个多项式函数，以知给定的k+1个点((x_{0},y_{0}),....,(x_{k},y_{k}))
设任意的两个(x_j)都不相同，那么应用拉格朗日插值公式所得到的多项式为：
(L(x)=sum_{j=0}^{k}y_{i}ell_{i}(x))
其中每个ell {j}(x)为拉格朗日基本多项式，其表达式为：
(ell _{j}(x):=prod _{{i=0,\,i eq j}}^{{k}}{frac {x-x_{i}}{x_{j}-x_{i}}})
拉格朗日基本多项式ell {j}(x)的特点是在x{j}上取值为1，在其它的点x{i},,i eq j上取值为0。

例子

假设有某个二次多项式函数f，已知它在三个点上的取值为：

(f(4)=10)
(f(5)=5.25)
(f(6)=1)
要求(f(18))的值。

首先写出每个拉格朗日基本多项式：

(ell _{0}(x)={frac {(x-5)(x-6)}{(4-5)(4-6)}})
(ell _{1}(x)={frac {(x-4)(x-6)}{(5-4)(5-6)}})
(ell _{2}(x)={frac {(x-4)(x-5)}{(6-4)(6-5)}})
然后应用拉格朗日插值法，就可以得到p的表达式（p为函数f的插值函数）：

(p(x)=f(4)ell _{0}(x)+f(5)ell _{1}(x)+f(6)ell _{2}(x))
(.\,\,\,\,\,\,\,\,\,\,=10cdot {frac {(x-5)(x-6)}{(4-5)(4-6)}}+5.25cdot {frac {(x-4)(x-6)}{(5-4)(5-6)}}+1cdot {frac {(x-4)(x-5)}{(6-4)(6-5)}})
(.\,\,\,\,\,\,\,\,\,\,={frac {1}{4}}(x^{2}-28x+136))
此时代入数值(18)就可以求出所需之值：( f(18)=p(18)=-11)

证明

存在性：

对于给定的k+1个点：((x_{0},y_{0}),ldots ,(x_{k},y_{k}))
拉格朗日插值法的思路是找到一个在一点(x_{j})取值为1，而在其他点取值都是0的多项式(ell _{j}(x))
这样，多项式(y_{j}ell _{j}(x))在点(x_{j})取值为(y_{j})，而在其他点取值都是0
而多项式(L(x):=sum _{{j=0}}^{{k}}y_{j}ell _{j}(x))就可以满足
(L(x_{j})=sum _{{i=0}}^{{k}}y_{i}ell _{i}(x_{j})=0+0+cdots +y_{j}+cdots +0=y_{j})
在其它点取值为0的多项式容易找到，例如：
((x-x_{0})cdots (x-x_{{j-1}})(x-x_{{j+1}})cdots (x-x_{{k}}))
在点(x_{j})取值为：
((x_{j}-x_{0})cdots (x_{j}-x_{{j-1}})(x_{j}-x_{{j+1}})cdots (x_{j}-x_{{k}}))。
因为假定(x_{i})两两互不相同，因此上面的取值不等于0。
将多项式除以这个取值，得到一个满足“在(x_{j})取值为1，而在其他点取值都是0的多项式”：
(ell _{j}(x):=prod _{{i=0,\,i eq j}}^{{k}}{frac {x-x_{i}}{x_{j}-x_{i}}}={frac {(x-x_{0})}{(x_{j}-x_{0})}}cdots {frac {(x-x_{{j-1}})}{(x_{j}-x_{{j-1}})}}{frac {(x-x_{{j+1}})}{(x_{j}-x_{{j+1}})}}cdots {frac {(x-x_{{k}})}{(x_{j}-x_{{k}})}})
完

唯一性

次数不超过(k)的拉格朗日多项式至多只有一个:
对任意两个次数不超过k的拉格朗日多项式：(P_{1}和P_{2})，它们的差
(P_{1}-P_{2})在所有(k+1)个点上取值都是(0)，因此必然是多项式((x-x_{0})(x-x_{{1}})cdots (x-x_{{k}}))的倍数。
因此，如果这个差(P_{1}-P_{2})不等于0，次数就一定不小于(k+1)。
但是(P_{1}-P_{2})是两个次数不超过(k)的多项式之差，它的次数也不超过(k)。
所以(P_{1}-P_{2}=0)，也就是说(P_{1}=P_{2})。
得证

还剩代码啦

#include<iostream>  
#include<string>  
#include<vector>  
using namespace std;
double Lagrange(int n,vector<double>&X,vector<double>&Y,double x){  
  double ret=0;  
    for(int i=1;i<=n;i++) {  
        double tmp=Y[i];  
        for(int j=1;j<=n;j++){  
            if(i!=j) {  
                tmp*=(x-X[j]);  
                tmp/=(X[i]-X[j]);  
            }  
        }  
        ret+=tmp;  
    }  
    return ret;  
};  
int main(){  
    int n;
    std::cin>>n;
    vector<double>X(n,1);  
    vector<double>Y(n,1);
    for(int i=1;i<=n;i++){  
        cin>>X[i]>>Y[i];  
    }
    double x;  
    std::cin>>x;  
    std::cout<<Lagrange(n,X,Y,x)<<endl;  
    return 0;  
} 

注，以上部分转自http://www.cnblogs.com/ECJTUACM-873284962