6.1 多項式補間

補間式として多項式を使用するものを多項式補間という. ここでは, (誤差がなければ) 同じ結果を与えるが計算手順が異なるラグランジュ補間公式とニュートン補間公式も広い意味では多項式補間ととらえておくことにする.

ラグランジュ補間公式はよく使われる代表的な補間公式である. ニュートン補間公式はすでに求めた補間式に新たに標本点を追加するときに適している.

6.1.1 多項式補間

\(n\) 次多項式 \(p_n(x)\) を次のように定義する.
\[
p_n(x) = a_0x^n + a_1x^{n-1} + \dots + a_{n-1}x + a_n
\] \(n + 1\) 個の相異なる点 \(x_0, x_1, \dots, x_n\) における関数 \(f(x)\) の値 \(f(x_0), f(x_1), \dots, f(x_n)\) が が与えられているとき, 次の連立一次方程式を解くことにより \(f(x)\) を補間する \(n\) 次多項式 \(p_n(x)\) の係数 \(a_0, a_1, \dots, a_n\) を求めることができる.
\[
\boldsymbol{Va} = \boldsymbol{f}
\] ただし,
\[
\boldsymbol{V} =
\begin{pmatrix}
x_0^n & x_0^{n-1} & \dots & 1 \\
x_1^n & x_1^{n-1} & \dots & 1 \\
\vdots & \vdots & \dots & \vdots \\
x_n^n & x_n^{n-1} & \dots & 1 \\
\end{pmatrix}
,
\boldsymbol{a} =
\begin{pmatrix}
a_0 \\
a_1 \\
\vdots \\
a_n \\
\end{pmatrix}
,
\boldsymbol{f} =
\begin{pmatrix}
f_0 \\
f_1 \\
\vdots \\
f_n \\
\end{pmatrix}
\] である.

ここで, \(\boldsymbol{V}\) は連立一次方程式の係数行列, \(\boldsymbol{f}\) は連立一次方程式の右辺ベクトル (各点における関数値 \(f_i = f(x_i)\) を表す), \(\boldsymbol{a}\) は多項式の係数で, 連立一次方程式の解ベクトルになる.

この連立一次方程式の係数行列 \(\boldsymbol{V}\) はファンデルモンド行列と呼ばれる形をしているが, ファンデルモンド行列の行列式は次式で求められることが知られている.
\[
det(\boldsymbol{V}) = \prod_{0 \le i < j \le n}(x_j - x_i) \] したがって, \(x_0, x_1, \dots, x_n\) が相異なる点であれば, \(det(\boldsymbol{V}) \ne 0\) であり, 上の連立一次方程式は解を持つ. そして, 解は \(n\) 次多項式 \(p_n(x)\) の係数である. すなわち, 上の連立一次方程式を解くことにより補間多項式が求められる.

6.1.2 数値実験 (1)

区間 \([-1, 1]\) において次の指数関数を多項式補間する.
\[
f(x) = 2e^{x-1} – 1
\] \(n = 2 \sim 5\) として (すなわち, \(2 \sim 5\) 次多項式を用いて) 補間を行う. ただし, 区間 \([-1, 1]\) を \(n\) 等分して, 両端を含め等間隔に並んだ \(n + 1\) 点における関数値をデータとして使用することにする. 多項式補間としては各標本点が相異なっていればよく等間隔である必要はないが, 数表などを想定すればこのように標本点を選ぶのは自然である.

上の連立一次方程式を解いて多項式の係数を求め, 区間内における多項式の値をプロットする. 横軸は \(x\), 縦軸は関数値である.

縦軸を誤差 (補間値 – 厳密値) としてプロットすると次のようになる.

多項式の次数が上がると精度も上がっているのがわかる.

6.1.3 数値実験 (2)

多項式補間ではどんな関数についても精度を上げたいときには次数を上げていけばよいというわけではないことを示す.

区間 \([-1, 1]\) において次の関数 (ルンゲの関数という) を多項式補間する.
\[
f(x) = 1/(1 + 25x^2)
\] \(n = 4, 8, 12, 20\) として (すなわち, \(4, 8, 12, 20\) 次多項式を用いて) 補間を行う. ただし, 区間 \([-1, 1]\) を \(n\) 等分して, 両端を含め等間隔に並んだ \(n + 1\) 点における関数値をデータとして使用する.

上の連立一次方程式を解いて多項式の係数を求め, 区間内における多項式の値をプロットする. 横軸は \(x\), 縦軸は関数値である.

図のように両端で発散して (ルンゲの現象という) 補間式としては使えない. 特に \(n\) が大きくなるほど両端での発散が大きくなった. 最大誤差だけを見るとすべての \(n\) の中で \(n = 4\) のときが最小になる. ただし, 図のように中央付近では次数が上がるほどよい近似をしていることがわかるので, この領域でだけ補間式として使うという方法はあるかもしれない.

この問題は標本点を等間隔ではなくチェビシェフ多項式のゼロ点にとると解決することが知られている. 「6.3 直交多項式補間」を参照せよ.

なお, 一般に区間内で関数値の変化が十分ゆるやかであればこのような振動現象を起こすことはなく, 通常は等間隔であっても標本点を増やすと補間の精度は向上する.

6.1.4 ラグランジュ補間

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

ラグランジュ補間公式は次式で表される.
\[
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 \ne k\) のとき \(0\) である. したがって, ラグランジュ補間公式は各標本点において与えられたデータに一致する. また, 標本点の並んでいる順序に関係なく成り立ち, 標本点は等間隔でなくてもよい.

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

ラグランジュ補間公式の計算は次のようにして行うとよい.

まず, 与えられたデータについて
\[
w_k = f(x_k)/\prod_{i=0}^{n (i \ne k)}(x_k – x_i)
\] を \(k = 0 \sim n\) について計算しておく.

次に, 求めたい点 \(x\) それぞれについて次のようにして補間値を計算する.
\[
\begin{align}
& w(x) = \prod_{k=0}^{n}(x – x_k) \\
& p_n(x) = w(x)\sum_{k=0}^{n}w_k/(x – x_k) \\
\end{align}
\] ただし, \(x = x_k (k = 0, \dots, n)\) であれば計算せずに \(p_n(x) = f(x_k)\) とすればよい.

6.1.2 および 6.1.3 の数値実験についてラグランジュ補間公式により計算すると多項式補間のときと誤差の範囲で同じ結果が得られる.

6.1.5 差分商

\(x = x_0, x_1\) に対する差分商を次のように定義する.
\[
f[x_0, x_1] = (f(x_1) – f(x_0))/(x_1 – x_0)
\] これは \(x_0\) と \(x_1\) が近ければ微分値 \(f'(x_0)\) を表す. あるいは, 区間 \([x_0, x_1]\) における微分値の平均値とみなすことができる.

同様に, \(x = x_0, x_1, x_2\) に対する 2 階の差分商は次のように定義できる.
\[
f[x_0, x_1, x_2] = (f[x_1, x_2] – f[x_0, x_1])/(x_2 – x_0)
\] すなわち, 2 階の差分商は 1 階の差分商の差分商である.

便宜上 \(f[x_i] = f(x_i)\) と表すと, 高階の差分商まで含めた一般形は n + 1 個の標本点 \(x_0, x_1, \dots, x_n\) について次のように表すことができる.
\[
f[x_i, \dots, x_{i+k}] = (f[x_{i+1}, \dots, x_{i+k}] – f[x_i, \dots, x_{i+k-1}])/(x_{i+k} – x_i) (k = 1 \sim n, i = 0 \sim n – k)
\]

ここで, 関数 \(f(x)\) を差分商を用いて表すことを考える.

\(x, x_0\) に対する 1 階の差分商
\[
f[x, x_0] = (f(x_0) – f(x))/(x_0 – x)
\] を変形すると次のようになる.
\[
f(x) = f(x_0) + (x – x_0)f[x, x_0] \] さらに, \(x, x_0, x_1\) に対する 2 階の差分商
\[
f[x, x_0, x_1] = (f[x_0, x_1] – f[x, x_0])/(x_1 – x)
\] を \(f[x, x_0]\) について解き, 代入すると次のようになる.
\[
f(x) = f(x_0) + (x – x_0)f[x_0, x_1] + (x – x_0)(x – x_1)f[x, x_0, x_1] \] 同様の代入を繰り返すと, n + 1 個の標本点 \(x_0, x_1, \dots, x_n\) について次式が得られる (\(f[x_i] = f(x_i)\) と表した).
\[
\begin{align}
f(x) & = f[x_0] \\
& + (x – x_0)f[x_0, x_1] \\
& + (x – x_0)(x – x_1)f[x_0, x_1, x_2] \\
& + \dots \\
& + (x – x_0)(x – x_1) \dots (x – x_{n-1})f[x_0, x_1, \dots, x_n] \\
& + (x – x_0)(x – x_1) \dots (x – x_{n-1})(x – x_n)f[x, x_0, x_1, \dots, x_n] \\
\end{align}
\] ここで, 最後の項は \(f[x, x_0, x_1, \dots, x_n]\) に変数 \(x\) を含むため, 通常は計算することが難しい.

6.1.6 ニュートン補間

関数 \(f(x)\) を差分商を用いて表した上式において最後の項を \(0\) とおいた次式をニュートン補間公式という.
\[
\begin{align}
f(x) & = f[x_0] \\
& + (x – x_0)f[x_0, x_1] \\
& + (x – x_0)(x – x_1)f[x_0, x_1, x_2] \\
& + \dots \\
& + (x – x_0)(x – x_1) \dots (x – x_{n-1})f[x_0, x_1, \dots, x_n] \\
\end{align}
\] この式は \(f(x)\) の近似式になっており, \(0\) とおいた最後の項は剰余項と呼ばれ近似誤差を表す.

ニュートン補間公式は, ラグランジュ補間公式を \(1, (x – x_0), (x – x_0)(x – x_1), \dots\) について展開した式に相当する. すでに求めた補間式に新たに標本点を追加するとき, ラグランジュ補間公式では全部計算し直す必要があるが, ニュートン補間公式では新たな項を追加するだけでよい.

ニュートン補間公式の計算は次のようにして行うとよい.

まず, 与えられたデータについて差分商を下図のように左から右の順に計算しておく. 実際に補間値の計算に必要なのは対角部分だけであるが, 標本点を追加するときに一番下の行が必要になる.
\[
\begin{align}
& f[x_0] \\
& & \ddots \\
& f[x_1] & \dots & \quad f[x_0, x_1] \\
& & \ddots & & \ddots \\
& f[x_2] & \dots & \quad f[x_1, x_2] & \dots & \quad f[x_0, x_1, x_2] \\
& & \ddots & & \ddots \\
& \vdots \\
& & \ddots & & \ddots & & & & \ddots \\
& f[x_n] & \dots & \quad f[x_{n-1}, x_n] & \dots & \quad f[x_{n-2}, x_{n-1}, x_n] & \dots & \quad & \dots & \quad f[x_0, \dots, x_n] \\
\end{align}
\] 次に, 求めたい点 \(x\) それぞれについて, 多項式を計算する要領でニュートン補間公式により補間値を計算する.

標本点を追加するときは差分商の一番下の行を追加計算し, ニュートン補間公式の最後の項を補間値に加える.

6.1.2 および 6.1.3 の数値実験についてニュートン補間公式により計算すると多項式補間およびラグランジュ補間公式と誤差の範囲で同じ結果が得られる.

6.1.7 反復補間法

補間により \(x\) における関数値を計算するときに, はじめは少ない標本点を使って計算し, 次第に標本点の数を増やして計算を繰り返し目標精度に達したら (標本点を増やす前後の差が許容範囲に入ったら) 収束とみなして終了する方法を反復補間法という.

計算は, ニュートン補間公式を使用して標本点を追加していくことにより行う. ニュートン補間公式は標本点の並び順に関係なく適用できるが, 反復補間法に適用する場合には, 標本点を \(x\) に近い順に追加するようにし, しかも標本点の数が \(x\) の両側でほぼ等しくなるようにするのが収束の速さの点で有利である.