三次拉格朗日插值多项式

本文属于代数插值的内容，对于很多数据处理，比如物理实验之类的，得到了很多个数据点，但是并不能找到一个解析的公式来描述，但是我们确实知道该函数是连续的，这样，我们可以通过一系列的点，来插值，找出某个点的值，并且我们使用多项式来逼近，这是多项式的插值。
本次算法使用多项式的插值公式，并且是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;
}

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