三次拉格朗日插值多项式

本文属于代数插值的内容，对于很多数据处理，比如物理实验之类的，得到了很多个数据点，但是并不能找到一个解析的公式来描述，但是我们确实知道该函数是连续的，这样，我们可以通过一系列的点，来插值，找出某个点的值，并且我们使用多项式来逼近，这是多项式的插值。
本次算法使用多项式的插值公式，并且是3次插值，即，只用3个数据点进行插值。对于如果有很多歌数据点，那就找插值点最近的3个数据点进行插值。本算法会自动判断哪几个是插值节点。
下面贴代码：

//3点拉格朗日插值//采用抛物线插值，自动寻找插值点最近的3个点/* *输入是一个数组x[n]，该数组中存放数据点的横坐标 *还有一个输入数组y[n],该数组中存放数据点的纵坐标 *输入插值点t *返回值是double值 */#include<iostream>#include<cmath>#include<iomanip>using namespace std;double Lagrange_3(double* x, double* y, double t, int n){    int k, m,i;    if (n < 1)                                                      //如果少于一个数，抛出异常        throw("not have enough point !");                    if (n == 1)                                                     //只有一个数据点，直接返回那个点        return y[n - 1];    if (n == 2)                                                     //只有2个点，直接线性插值        return (y[0] * ((t - x[0]) / (x[0] - x[1])) + y[1] * (t - x[0]) / (x[1] - x[0]));    if(t<=x[1])                                                     //判断是否在头2个数据点之间    {        k = 0; m = 2;    }    else         if(t>x[n-2])                                                    //是否在末两个点之间        {            k = n - 3; m = n - 1;        }        else                                                          //不在就进行查找最近的2个点，二分查找        {            k = 1; m = n;            while (m - k > 1)                                      //循环条件            {                i = (k + m) / 2;                                   //二分                if (t < x[i - 1])                    m = i;                else                    k = i;            }            cout << "k " << k << "m " << m << endl;            k = k - 1;            m = m - 1;            if (fabs(t - x[k]) < fabs(t - x[m]))           //如果离k的那个点距离较近                k = k - 1;                                    //k在减一，则为3个点            else                m = m + 1;                                          }    double z = 0.0;      for(i=k;i<=m;i++)                                               //进行三次插值    {        double s = 1.0;        for (int j = k; j <= m; j++)            if (j != i)                s = s*(t - x[j]) / (x[i] - x[j]);        z = z + s*y[i];    }    return z;}int main(){    double aa = 1.682;    double bb = 1.813;    double c;    double x[5] = { 1.615,1.634,1.702,1.828,1.921 };    double y[5] = { 2.41450,2.46459,2.65271,3.03035,3.34066 };    c = Lagrange_3(x, y, aa, 5);    cout <<  setprecision(8)<<c<<endl;    double ff;    ff= Lagrange_3(x, y, bb, 5);    cout << ff<< endl;    return 0;}

如果对文章有什么见解，欢迎大家提出！