5. 補間

5.1 概要

補間法はしばしば最小二乗法と混同されます. 最小二乗法ではデータ点に誤差が含まれていることを前提とし, 近似式はぴったりデータ点を通るわけではありません. しかし, 補間法は誤差が含まれていないデータ点を滑らかにつなぐ方法で, 補間式はぴったり全データ点を通ります. データ点に誤差が含まれていると不自然に曲がった補間曲線が得られてしまいます.

いろいろな関数値を表にした数表というのがあります. これを使う際に, 表にのっていない点の値を求めるために補間 (内挿) が使われていました. 現在では関数電卓やコンピュータがあるためその需要はなくなりましたが, 例えば, 1 点を計算するのに多大な労力(計算時間)が必要な関数の値を素早く求めたい場合などに現在でも使われています.

なお, 見えないところで数値解法の基盤として補間法が使われています. 例えば, 数値積分では分点の補間式を求めて補間関数を積分した関数の値を積分値として採用する方法が使われています.

5.2 補間法の解法

5.2.1 多項式補間 (ラグランジュ補間)

\(n\) 次多項式 \(p_n(x)\) を次のように定義します.
\[
p_n(x) = a_0x^n + a_1x^{n-1} + \dots + a_{n-1}x + a_n
\] この \(p_n(x)\) で関数 \(f(x)\) を補間するものとします.

\(n + 1\) 個の相異なる点 \(x_0, x_1, \dots, x_n\) における関数値 \(f(x_0), f(x_1), \dots, f(x_n)\) が与えられていれば, \((n + 1)\)元連立一次方程式を解くことにより \(n + 1\) 個の係数 \(a_0, a_1, \dots, a_n\) を定めることができます.

多項式補間を行うには, このように連立一次方程式を解いて多項式の係数を直接求めるのが直観的な方法ですが, 実務的には方程式の条件数が大きくなりやすいなどの理由により, 等価な多項式であるラグランジュ補間公式が用いられます.

ラグランジュ補間公式は次式で表されます.
\[
f(x) = \sum_{k=0}^{n}L_k(x) f(x_k)
\] ただし
\[
L_k(x) = \prod_{i=0}^{n (i \neq k)}(x – x_i)/(x_k – x_i)
\] \(L_k(x_i)\) は \(i = k\) のとき 1, \(i \neq k\) のとき 0 になります. したがって, ラグランジュ補間公式は各標本点において与えられたデータに一致します. また, 標本点の並んでいる順序に関係なく成り立ち, 標本点は等間隔な必要はありません.

ラグランジュ補間公式は \(n\) 次多項式であり, バラバラにして整理するとその係数は連立一次方程式を解いて求めた多項式の係数と一致します. しかし, 補間においては関数値が求まればよくて多項式の係数は計算に直接は必要ないので, この形のまま計算に使われます.

5.2.2 3次スプライン補間

全体を小さい区間に分割しそれぞれの区間で多項式補間を行う方法を区分的多項式補間といいます. その中で最もよく使われるのは 3 次スプライン補間です. 連続するデータ点の組を 3 次多項式で近似し, つなぎ目が滑らかになるように 2 次までの導関数が連続になるように構成されたものです.

データ \((x_1, y_1), (x_2, y_2), \dots, (x_n, y_n)\) が与えられているものとします. これらを節点といいます. 3 次スプラインは, 区間 \([x_1, x_n]\) において, 隣り合う節点の間では 3 次多項式であり, 隣の小区間の多項式とはその間の節点で 1 次と 2 次微分が連続になります.

このような \(n – 1\) 個の 3 次多項式を定めるためにはそれぞれ \(4\) 個のパラメータがあるので, 全部で \(4 \times (n – 1)\) 個のパラメータを決める必要があります. これに対して, \(n\) 個の節点で \(y_i\) の値が決められているから \(n\) 個, そして内側の \(n – 2\) 個の節点において関数値, 1 次微分, 2 次微分が連続になるから \(3 \times (n – 2)\) 個, 合わせて \(4n – 6\) 個のパラメータが決まり, あと \(2\) 個の自由度があることになります.

\(2\) 個の条件として, 例えば \(x_1\) および \(x_n\) における微分値がわかればそれを使うとよいことになります. そうでなければ \(x_1\) および \(x_n\) における 2 次微分値を 0 にする方法があり, これを自然スプラインといいます. また, \(x_1\) および \(x_n\) における 3 次微分を連続にする not-a-knot 条件と呼ばれる方法もあります.

以上のパラメータを求めることは, 係数が正定値対称な 3 重対角行列の連立一次方程式を解くことに帰着されます.

5.3 XLPack を使った補間法の解き方

3 次スプライン補間係数を求めるために Pchse を, 補間値を求めるために Pchfe を使用することができます. ワークシート関数では, それぞれ WPchse と WPchfe を使用することができます.

例題: 対数表

下表は自然対数表の一部を抜き出したものです.

5.3.1 ワークシート関数を使用した解き方

ワークシートの適当な場所にデータ \((x_i, y_i = ln(x_i))\) を入力します(左のオレンジ色のセル). ここで, \(x_i\) は昇順になっていなければなりません. スプライン補間に必要な係数 \(d_i\) を求めるためにデータ数個のセルを選択してワークシート関数 WPchse を入力します(左の緑色のセル).

WPchse の必要なパラメータは N, X, Y です. N はデータ数(この例では 8), X と Y は補間するデータ \((x_i, y_i)\) のセル範囲です.

入力が終了したら Ctrl + Shift + Enter を押します. これで, スプライン補間係数が求められました.

次に, 別の場所に求めたい点のデータ(この場合, \(n = 1.55 \sim 2.15\))を入力しておき(右のオレンジ色のセル), そこにおける関数値を WPchfe 関数を使って求めます(右の緑色のセル).

WPchfe の必要なパラメータは Ne, Xe, N, X, Y, D です. Ne は求めたい点のデータ数(この例では 7), Xe は求めたい点のデータのセル範囲です. N, X, Y は補間係数を求めたときの WPchse の入力データと同じです. D は求められた補間係数です.

Ctrl + Shift + Enter を押すと, 表に載っていない点における自然対数の値が求められます. ここでは, 表に載っていない点を丸で, 数表の値を直線で表すグラフも作成してみました.

この例では, データ点は等間隔でしたが, 必ずしも等間隔である必要はありません. ただし, 間隔が大きく空いている区間の補間値は精度が悪くなります. また, データ点の範囲外での値(外挿値)を計算することもできますが精度は落ちます.

5.3.2 VBAプログラムを使用した解き方

上と同じ例を VBA プログラムにより解きます. Pchse および Pchfe を使ったプログラム例を示します.

Sub Start()
    Const NMax = 10
    Dim N As Long, Ne As Long, I As Long, Info As Long
    Dim X(NMax) As Double, Y(NMax) As Double, D(NMax) As Double
    Dim Xe(NMax) As Double, Fe(NMax) As Double
    '--- Input data
    N = 8
    For I = 0 To N - 1
        X(I) = Cells(5 + I, 1)
        Y(I) = Cells(5 + I, 2)
    Next
    '--- Compute coefficients of a cubic spline
    Call Pchse(N, X(), Y(), D(), Info)
    If Info <> 0 Then
        MsgBox "** Error in Pchse **  Info = " + Str(Info)
        Exit Sub
    End If
    '--- Compute interporated values
    Ne = N - 1
    For I = 0 To Ne - 1
        Xe(I) = Cells(5 + I, 4)
    Next
    Call Pchfe(N, X(), Y(), D(), Ne, Xe(), Fe(), Info)
    If Info <> 0 Then
        MsgBox "** Error in Pchfe **  Info = " + Str(Info)
        Exit Sub
    End If
    '--- Output result
    For I = 0 To Ne - 1
        Cells(5 + I, 5) = Fe(I)
    Next
End Sub

所定の位置(下のオレンジ色のセル)にデータを入力し, マクロ Start としてプログラムを実行すると次の結果が得られます. ワークシート関数を使用したときと異なり, 補間データを入力しただけでは結果が得られず, VBA プログラムを実行してやる必要があります.