Go言語で関数の近似(ニュートン補間)

をニュートン補間で近似する

与えられた４個の関数値 f(x₀), f(x₁), f(x₂), f(x₃) を通る３次式を求める場合, まず次のような表を作る.

このとき,

と定義する.
これを x₀とx₁ の第１差分商, または１階差分商という.
同様に

と定義し, これを x₀, x₁, x₂ の第２差分商, または２階差分商という.
第ｎ差分商を

で表すと, 与えられた n+1点を通る 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 まで, ５点だけ値をセット
    for i := 0; i < N; i++ {
        var d float64 = float64(i) * 0.25 - 0.5
        x[i] = d
        y[i] = f(d)
    }

    // 差分商の表を作る
    var d[7][7] float64
    for j := 0; j < N; j++ {
        d[0][j] = y[j]
    }
    for i := 1; i < N; i++ {
        for j := 0; j < N - i; j++ {
            d[i][j] = (d[i-1][j+1] - d[i-1][j]) / (x[j+i] - x[j])
        }
    }

    // ｎ階差分商
    var a [N]float64
    for j := 0; j < N; j++ {
        a[j] = d[j][0]
    }

    // 0.05刻みで 与えられていない値を補間 
    for i := 0; i <= 20; i++ {
        var d  float64 = float64(i) * 0.05 - 0.5
        var d1 float64 = newton(d, x[:], a[:])
        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
}

// Newton (ニュートン) 補間
func newton(d float64, x []float64, a []float64) float64 {
    var sum float64 = a[0]
    for i := 1; i < N; i++ {
        var prod float64 = 1
        for j := 0; j < i; j++ {
            prod *= (d - x[j])
        }
        sum += a[i] * prod
    }
    return sum
}

Z:\>8g GO0702.go

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

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

参考文献