Go言語で関数の近似(エルミート補間)

をエルミート補間で近似する

２個の関数値 f(x₀), f(x₁) と,それぞれの微分係数f'(x₀), f'(x₁) とを
与えられたとき, 与えられた２個の関数値を通る５３次式を求めるには, まず次のような表を作る.

このとき,
.
あとは, ニュートン補間と同じ


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

package main

import "fmt"
import "math"

// データ点の数
const N  = 3
const N2 = 6

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

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

    // 差分商の表を作る
    var z[N2]     float64
    var d[N2][N2] float64
    for i := 0; i < N2; i++ {
        j := i / 2
        z[i]    = x[j]
        d[0][i] = y[j]
    }

    for i := 1; i < N2; i++ {
        for j := 0; j < N2 - i; j++ {
            if (i == 1 && j % 2 == 0) {
                d[i][j] = yy[j / 2]
            } else {
                d[i][j] = (d[i-1][j+1] - d[i-1][j]) / (z[j+i] - z[j])
			}
        }
    }

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

    // 0.1刻みで 与えられていない値を補間 
    for i := 0; i <= 80; i++ {
        var d  float64 = float64(i) * 0.1 - 4.0
        var d1 float64 = hermite(d, z[:], 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 x - math.Pow(x,3)/6 + math.Pow(x,5)/120
}
// 導関数
func fd(x float64) float64 {
    return (24 - 12 * math.Pow(x,2) + math.Pow(x,4)) / 24
}

// Hermite (エルミート) 補間
func hermite(d float64, z []float64, a []float64) float64 {
    var sum float64 = a[0]
    for i := 1; i < N2; i++ {
        var prod float64 = 1
        for j := 0; j < i; j++ {
            prod *= (d - z[j])
        }
        sum += a[i] * prod
    }

    return sum
}

Z:\>8g GO0704.go

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

Z:\>GO0704
-4.00: -1.86667 - -1.86667 = -0.00000
-3.90: -1.53218 - -1.53218 = -0.00000
-3.80: -1.25760 - -1.25760 = -0.00000
-3.70: -1.03650 - -1.03650 = -0.00000
-3.60: -0.86285 - -0.86285 = -0.00000
-3.50: -0.73099 - -0.73099 = -0.00000
-3.40: -0.63562 - -0.63562 = -0.00000
-3.30: -0.57178 - -0.57178 = -0.00000
-3.20: -0.53487 - -0.53487 = -0.00000
-3.10: -0.52060 - -0.52060 = -0.00000
-3.00: -0.52500 - -0.52500 = -0.00000
-2.90: -0.54443 - -0.54443 = -0.00000
-2.80: -0.57553 - -0.57553 = -0.00000
-2.70: -0.61524 - -0.61524 = -0.00000
-2.60: -0.66078 - -0.66078 = -0.00000
-2.50: -0.70964 - -0.70964 = -0.00000
-2.40: -0.75955 - -0.75955 = -0.00000
-2.30: -0.80853 - -0.80853 = -0.00000
-2.20: -0.85480 - -0.85480 = -0.00000
-2.10: -0.89684 - -0.89684 = -0.00000
-2.00: -0.93333 - -0.93333 = -0.00000
以下略

エルミート補間に３個の関数値と,それぞれの微分係数を与えた結果

ちなみに、ニュートン補間に６点を与えた場合...

参考文献