C#版本拉格朗日插值演算法

        實際中有時需要對資料進行分析, 最近就遇到了這樣一個情況: 有一系列橫座標是時間, 縱座標是記錄值的一些資料, 但橫座標卻不是等距記錄的. 就是說在第一分鐘記錄一次, 第二分鐘記錄一次, 第四分鐘記錄一次...不等距. 需求是根據現在的這些資料進行計算, 獲得等距時間所對應的記錄值, 也就是說第三分鐘沒有記錄, 得通過計算來獲得第三分鐘的值.

        原本想把這些資料顯示在Chart控制項中, 然後再從Chart控制項中擷取各個X座標所對應的Y值, 但後來發現這樣不可行, 只能得到有記錄的資料點的值, 新點的值得不到. 只能使用插值來計算了, 下面是拉格朗日插值的演算法(C#版本):

/// <summary>/// 根據離散點進行二次的拉格朗日插值/// http://www.cnblogs.com/technology/// </summary>class Lagrange{/// <summary>/// X各點座標組成的數組/// </summary>public int[] x { get; set; }/// <summary>/// X各點對應的Y座標值組成的數組/// </summary>public double[] y { get; set; }/// <summary>/// x數組或者y數組中元素的個數, 注意兩個數組中的元素個數需要一樣/// </summary>public int itemNum { get; set; }/// <summary>/// 初始化拉格朗日插值/// </summary>/// <param name="x">X各點座標組成的數組</param>/// <param name="y">X各點對應的Y座標值組成的數組</param>public Lagrange(int[] x, double[] y){    this.x = x; this.y = y;    this.itemNum = x.Length;}/// <summary>/// 獲得某個橫座標組應的Y座標值/// </summary>/// <param name="xValue">x座標值</param>/// <returns></returns>public double GetValue(int xValue){    //用於累乘數組始末下標    int start, end;    //傳回值    double value = 0.0;    //如果初始的離散點為空白, 返回0    if (itemNum < 1) { return value; }    //如果初始的離散點只有1個, 返回該點對應的Y值    if (itemNum == 1) { value = y[0]; return value; }    //如果初始的離散點只有2個, 進行線性插值並返回插值    if (itemNum == 2)    {        value = (y[0] * (xValue - x[1]) - y[1] * (xValue - x[0])) / (x[0] - x[1]);        return value;    }    //如果插值點小於第一個點X座標, 取數組前3個點做插值    if (xValue <= x[1]) { start = 0; end = 2; }    //如果插值點大於等於最後一個點X座標, 取數組最後3個點做插值    else if (xValue >= x[itemNum - 2]) { start = itemNum - 3; end = itemNum - 1; }    //除了上述的一些特殊情況, 通常情況如下    else    {        start = 1; end = itemNum;        int temp;        //使用二分法決定選擇哪三個點做插值        while ((end - start) != 1)        {            temp = (start + end) / 2;            if (xValue < x[temp - 1])                end = temp;            else                start = temp;        }        start--; end--;        //看插值點跟哪個點比較靠近        if (Math.Abs(xValue - x[start]) < Math.Abs(xValue - x[end]))            start--;        else            end++;    }    //這時已經確定了取哪三個點做插值, 第一個點為x[start]    double valueTemp;    //注意是二次的插值公式    for (int i = start; i <= end; i++)    {        valueTemp = 1.0;        for (int j = start; j <= end; j++)            if (j != i)                valueTemp *= (double)(xValue - x[j]) / (double)(x[i] - x[j]);        value += valueTemp * y[i];    }    return value;}

        等時間距擷取值的問題就這樣解決了, 但要注意這是二次插值, 也就是擬合的是二次及二次方程以下的函數.