Go言語で関数の近似(ネヴィル補間)

をネヴィル補間で近似する

与えられた５個の関数値 f(x₀), f(x₁), f(x₂), f(x₃), f(x₄) を通る４次式は、次のようにして求めることができる.
(1)

(2)

(3)

(4)

(5)

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

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 まで, ５点だけ値をセット
    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 = neville(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
}

// Neville (ネヴィル) 補間
func neville(d float64, x []float64, y []float64) float64 {
    var w[N][N] float64

    for i := 0; i < N; i++ {
        w[i][i] = y[i]
    }

    for j := 1; j < N; j++ {
        for i := j; i < N; i++ {
            w[i][i-j] = w[i][i-j+1] + (w[i][i-j+1] - w[i-1][i-j]) * (d - x[i]) / (x[i] - x[i-j])
        }
    }
    
    return w[N-1][0]
}

Z:\>8g GO0703.go

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

Z:\>GO0703
-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

参考文献