C#で関数の近似(ラグランジュ補間)

をラグランジュ補間で近似する

n+1個の点 (x₀, y₀), (x₁, y₁) … (x_n, y_n)
が与えられているとき, これらすべての点を通る n次式は次のように表すことができる.

この式を使って, 与えられた点以外の点の値を求める.

using System;

public class CS0701
{
    private const int N = 5; // データ点の数

    public static void Main()
    {
        double[] x = new double[N];
        double[] y = new double[N];

        // 0.25刻みで -0.5〜0.5 まで, ５点だけ値をセット
        for (int i = 0; i < N; i++)
        {
            double d = i * 0.25 - 0.5;
            x[i] = d;
            y[i] = f(d);
        }

        // 0.05刻みで 与えられていない値を補間 
        for (int i = 0; i <= 20; i++)
        {
            double d  = i * 0.05 - 0.5;
            double d1 = lagrange(d, x, y);
            double d2 = f(d);

            // 元の関数と比較
            Console.WriteLine(string.Format("{0,5:F2}: {1,8:F5} - {2,8:F5} = {3,8:F5}", d, d1, d2, d1 - d2));
        }
    }

    // 元の関数
    private static double f(double x)
    {
        return 4 * Math.Pow(x,4) + 3 * Math.Pow(x,3) - 2 * Math.Pow(x,2) - x + 1;
    }

    // Lagrange (ラグランジュ) 補間
    private static double lagrange(double d, double[] x, double[] y)
    {
        double sum = 0;
        for (int i = 0; i < N; i++)
        {
            double prod = y[i];
            for (int j = 0; j < N; j++)
            {
                if (j != i)
                    prod *= (d - x[j]) / (x[i] - x[j]);
            }
            sum += prod;
        }
        return sum;
    }
}

Z:\>csc CS0701.cs -nologo

Z:\>CS0701
-0.50:  0.87500 -  0.87500 =  0.00000
-0.45:  0.93565 -  0.93565 =  0.00000
-0.40:  0.99040 -  0.99040 =  0.00000
-0.35:  1.03640 -  1.03640 =  0.00000
-0.30:  1.07140 -  1.07140 =  0.00000
-0.25:  1.09375 -  1.09375 =  0.00000
-0.20:  1.10240 -  1.10240 =  0.00000
-0.15:  1.09690 -  1.09690 =  0.00000
-0.10:  1.07740 -  1.07740 =  0.00000
-0.05:  1.04465 -  1.04465 =  0.00000
 0.00:  1.00000 -  1.00000 =  0.00000
 0.05:  0.94540 -  0.94540 =  0.00000
 0.10:  0.88340 -  0.88340 =  0.00000
 0.15:  0.81715 -  0.81715 =  0.00000
 0.20:  0.75040 -  0.75040 =  0.00000
 0.25:  0.68750 -  0.68750 =  0.00000
 0.30:  0.63340 -  0.63340 =  0.00000
 0.35:  0.59365 -  0.59365 =  0.00000
 0.40:  0.57440 -  0.57440 =  0.00000
 0.45:  0.58240 -  0.58240 =  0.00000
 0.50:  0.62500 -  0.62500 =  0.00000

参考文献