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

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

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

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

object Scala0701 {
    // データ点の数 - 1
    val N = 6 

    def main(args: Array[String]) {
        // 1.5刻みで -4.5〜4.5 まで, ７点だけ値をセット
        val x = (0 to N).map(_ * 1.5 - 4.5)
        val y = x.map(f)

        // 0.5刻みで 与えられていない値を補間
        val d1 = (0 to 18).map(_ * 0.5 - 4.5)
        val d2 = d1.map(f)
        val d3 = d1.map(lagrange(_, x, y))

        (d1 zip d2 zip d3).foreach { 
            case ((d1, d2), d3) => 
                println("%5.2f\t%8.5f\t%8.5f\t%8.5f".format(d1, d2, d3, d2 - d3))
        }
    }

    // 元の関数
    def f(x:Double) = {
        x - Math.pow(x,3) / (3 * 2) + Math.pow(x,5) / (5 * 4 * 3 * 2)
    }

    // Lagrange (ラグランジュ) 補間
    def lagrange(d:Double, x:IndexedSeq[Double], y:IndexedSeq[Double]) = {
        var sum_list = List[Double](0)
        for (i <- 0 to N) {
            var prod_list = List(y(i))
            for (j <- 0 to N) {
                if (i != j) {
                    prod_list = ((d - x(j)) / (x(i) - x(j)))::prod_list
                }
            }
            sum_list = prod_list.product::sum_list
        }
        sum_list.sum
    }
}

Z:\>scala Scala0701.scala
-4.50   -4.68984        -4.68984         0.00000
-4.00   -1.86667        -1.86667        -0.00000
-3.50   -0.73099        -0.73099        -0.00000
-3.00   -0.52500        -0.52500         0.00000
-2.50   -0.70964        -0.70964        -0.00000
-2.00   -0.93333        -0.93333        -0.00000
-1.50   -1.00078        -1.00078         0.00000
-1.00   -0.84167        -0.84167        -0.00000
-0.50   -0.47943        -0.47943        -0.00000
 0.00    0.00000         0.00000         0.00000
 0.50    0.47943         0.47943         0.00000
 1.00    0.84167         0.84167         0.00000
 1.50    1.00078         1.00078         0.00000
 2.00    0.93333         0.93333         0.00000
 2.50    0.70964         0.70964         0.00000
 3.00    0.52500         0.52500         0.00000
 3.50    0.73099         0.73099         0.00000
 4.00    1.86667         1.86667         0.00000
 4.50    4.68984         4.68984         0.00000

参考文献