ONLY DO WHAT ONLY YOU CAN DO

こけたら立ちなはれ 立ったら歩きなはれ

Go言語で関数の近似(ラグランジュ補間)



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

n+1個の点 (x0, y0), (x1, y1) … (xn, yn)
が与えられているとき, これらすべての点を通る n次式は次のように表すことができる.

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

package main

import "fmt"
import "math"

const N = 5 // データ点の数

func main() {
    var x [N]float64
    var y [N]float64

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

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

        // 元の関数と比較
        fmt.Printf("%5.2f: %8.5f - %8.5f = %8.5f\n", d, d1, d2, d1 - d2)
    }
}

// 元の関数
func f(x float64) float64 {
    return 4 * math.Pow(x,4) + 3 * math.Pow(x,3) - 2 * math.Pow(x,2) - x + 1
}

// Lagrange (ラグランジュ) 補間
func lagrange(d float64, x []float64, y []float64) float64 {
    var sum float64 = 0
    for i := 0; i < N; i++ {
        var prod float64 = y[i]
        for j := 0; j < N; j++ {
            if (j != i) {
                prod *= (d - x[j]) / (x[i] - x[j])
            }
        }
        sum += prod
    }
    return sum
}
Z:\>8g GO0701.go

Z:\>8l -o GO0701.exe GO0701.8

Z:\>GO0701
-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
参考文献