多項式補間
定義 1
互いに異なる $x_{0} , \cdots , x_{n}$ のデータ $(x_{0}, y_{0}) , \cdots , (x_{n} , y_{n})$ に対して $p (x_{i} ) = y_{i}$ と $\deg p \le n$ を満たす 補間 である 多項式 $p$ を ポリノミアル・インターポレーションpolynomial Interpolationと呼ぶ。
定理
存在性と一意性
- [1]: 与えられたデータに対して $p$ は一意に存在する。
ラグランジュの公式
- [2]: $$p_{n} (x) = \sum_{i=0}^{n} y_{i} l_{i} (X)$$
ニュートンの差商公式
- [3]: $$p_{n} (x) = f(x_{0}) + \sum_{i=1}^{n} f [ x_{0} , \cdots , x_{i} ] \prod_{j=0}^{i-1} (x - x_{j} )$$
誤差解析
- [4]: $(n+1)$ 回微分可能な $f : \mathbb{R} \to \mathbb{R}$ とある $\xi \in \mathscr{H} \left\{ x_{0} , \cdots , x_{n} \right\}$ に対して $f$ のポリノミアル・インターポレーション $p_{n}$ はある $t \in \mathbb{R}$ に対して次を満たす。 $$\displaystyle f(t) - p_{n} (t) = {{ (t - x_{0}) \cdots (t - x_{n}) } \over { (n+1)! }} f^{(n+1)} ( \xi )$$
- $\mathscr{H} \left\{ a,b,c, \cdots \right\}$ は $a,b,c, \cdots$ を含む最小の区間を表す。
説明
ポリノミアル・インターポレーションは韓国語で多項式補間法と訳される。
条件 $\deg p \le n$
条件 $\deg p \le n$ とは、必ずしも $\deg p = n$ を満たす $p$ が常に存在するという保証がないということである。簡単な例として $n=2$ のとき上図のように三点が一直線上にあると、$(n+1)=3$ 個の点を通る $p_{2} (x) = a_{2} x^{2} + a_{1} x + a_{0}$ は存在しないが、それより次数が低い $p_{1} (x) = a_{1} x + a_{0}$ は存在する。これは実際には $p_{2} (x) = a_{2} x^{2} + a_{1} x + a_{0}$ を求めたが $a_{2} = 0$ である場合を意味する。
ラグランジュの公式とニュートンの差商公式は同じである
公式 [2] と [3] は見た目は異なっても、実は [1] によって互いに等しいことが分かる。本質的に両公式は文字通り $p_{n}$ をどう表すかの違いに過ぎず、用途上の違いはニュートンの差商公式が新しいデータが追加されたときに更新が容易であるという利点がある程度である。
実際の関数との誤差
定理 [4] はある関数 $f$ を補間する $p_{n}$ が $f$ とどれだけ差があるかを示す。通常、$(n+1)!$ の発散速度は非常に速いため、データが多ければ多いほど補間 $p_{n}$ の精度は上がる。ただしこの公式は厳密に収束性を論じるものではない点に注意する必要がある。簡単な例として $t$ が $\mathscr{H} \left\{ x_{0} , \cdots , x_{n} \right\}$ で本当に非常に遠くにある場合を考えられる。また $f$ があまり良くなく、微分するたびに値が大きくなり、それが $(n+1)!$ の発散速度より速ければ、いくらでもおかしくなる可能性がある。
証明
[1]
戦略: $(n+1)$ 個の連立方程式を行列表示にし、その逆行列の存在性と一意性を示す。
すべての $i = 0, \cdots , n$ に対して $y_{i} = p_{n} (x_{i}) = a_{0} + a_{1} x_{i} + \cdots + a_{n} x_{i}^{n}$ を満たす $p_{n}$ の係数 $a_{0} , a_{1} , \cdots , a_{n}$ が一意であることを示せばよい。
$$
\mathbf{y} := \begin{bmatrix}
y_{0}
\\ y_{1}
\\ \vdots
\\ y_{n}
\end{bmatrix}, X := \begin{bmatrix}
1 & x_{0} & \cdots & x_{0}^{n}
\\ 1 & x_{1} & \cdots & x_{1}^{n}
\\ \vdots & \vdots & \ddots & \vdots
\\ 1 & x_{n} & \cdots & x_{n}^{n}
\end{bmatrix} , \mathbf{a} := \begin{bmatrix}
a_{0}
\\ a_{1}
\\ \vdots
\\ a_{n}
\end{bmatrix}
$$ を定義し連立方程式を行列で表すと
$$
\begin{bmatrix}
y_{0}
\\ y_{1}
\\ \vdots
\\ y_{n}
\end{bmatrix} = \begin{bmatrix}
1 & x_{0} & \cdots & x_{0}^{n}
\\ 1 & x_{1} & \cdots & x_{1}^{n}
\\ \vdots & \vdots & \ddots & \vdots
\\ 1 & x_{n} & \cdots & x_{n}^{n}
\end{bmatrix} \begin{bmatrix}
a_{0}
\\ a_{1}
\\ \vdots
\\ a_{n}
\end{bmatrix}
$$
これで $\mathbf{y} = X \mathbf{a}$ を満たす解 $\mathbf{a}$ を求める問題に変わった。
ヴァンデルモンド行列の行列式: $X$ の行列式は $\displaystyle \det X = \prod_{1 \le i < j \le n } (x_{j} - x_{i})$
仮定より $x_{0} , \cdots , x_{n}$ は互いに異なるので $\det X \ne 0$ であり、 $X^{-1}$ が存在する。したがって $\mathbf{a} = X^{-1} \mathbf{y}$ も存在する。一方、与えられた行列の逆行列は一意なので、 $\mathbf{a}$ もまた一意である。
■
[2]
■
[3]
■
[4]
戦略: ダミー関数を新たに定義し、それらの微分可能性を利用して直接的な計算を回避する。設定が複雑なので実際には後ろから逆に理解する方が楽である。
主張: $E (x) := f(x) - p_{n} (X)$ について次が成り立つ。
$$
E(t) = {{ (t - x_{0}) \cdots (t - x_{n}) } \over { (n+1)! }} f^{(n+1)} ( \xi )
$$
Part 1.
まず $t$ がノード点 $x_{0} , \cdots , x_{n}$ と等しいなら自明に成り立つので、$t$ はこれらのノード点とは異なると仮定する。
$$
\begin{align*}
\Psi (x) :=& (x - x_{0} ) (x - x_{1}) \cdots (x - x_{n})
\\ G (x) :=& E(x) - {{ \Psi (x) } \over { \Psi (t) }} E(t)
\end{align*}
$$
上記のように新しい関数を定義すると $f$、$p_{n}$、$\Psi$ は $(n+1)$ 回微分可能だから $G$ も $(n+1)$ 回微分可能である。
Part 2.
$G$ に $x = x_{i}$ を代入すると $E(x_{i}) = f(x_{i}) - p_{n} (x_{i}) = 0$ であり $\Psi (x_{i} ) = 0$ なので
$$
G(x_{i}) = E(x_{i} ) - {{ \Psi (x_{i}) } \over { \Psi (t) }} E(t) = 0
$$
$G$ に $x = t$ を代入すると $\displaystyle {{ \Psi (t) } \over { \Psi (t) }} = 1$ なので
$$
G(t) = E(t) - E(t) =0
$$
したがって $G$ は $(n+2)$ 個の互いに異なる零点 $x_{0} , \cdots , x_{n} , t$ を持つ。
Part 3.
便宜上 $x_{n+1} :=t$ とする。
ロールの定理: 関数 $f(x)$ が $[a,b]$ で連続かつ $(a,b)$ で微分可能であり $f(a)=f(b)$ ならば $f ' (c)=0$ を満たす $c$ が $(a,b)$ に少なくとも一つ存在する。
すべての $i=0, \cdots , n$ に対して
$$
G(x_{i}) = 0 = G(x_{i+1})
$$
なのでロールの定理により $g ' ( x’_{i}) = 0$ を満たす $x'_{i} \in \mathscr{H} \left\{ x_{i} , x_{i+1} \right\}$ が存在し、同様にすべての $i=0, \cdots , (n-1)$ に対して
$$
g ' (x’_{i}) = 0 = g ' (x’_{i+1})
$$
なのでロールの定理により $G''( x’’_{i}) = 0$ を満たす $x''_{i} \in \mathscr{H} \left\{ x’_{i} , x’_{i+1} \right\}$ が存在する。このようにして帰納的にロールの定理を $(n+1)$ 回繰り返し適用すると
$$
G^{(n+1)} ( \xi ) = 0
$$
を満たす $\xi \in \mathscr{H} \left\{ x_{0} , \cdots , x_{n+1} \right\}$ が存在することが保証される。
一方 $p_{n}$ は $n$ 次の多項式なので
$$
E^{(n+1)} (x) = f^{(n+1)} ( x)
$$
$\Psi$ の最高次項は $x^{n+1}$ なので
$$
\Psi^{(n+1)} (x) = (n+1)!
$$
$\displaystyle G (x) = E(x) - {{ \Psi (x) } \over { \Psi (t) }} E(t)$ の両辺を $x$ について $(n+1)$ 回微分すると
$$
G^{(n+1)} (x) = f^{(n+1)} ( x) - {{ (n+1)! } \over { \Psi (t) }} E(t)
$$
$G^{(n+1)}$ と $f^{(n+1)}$ に $x=\xi$ を代入すると
$$
0 = f^{(n+1)} ( \xi ) - {{ (n+1)! } \over { \Psi (t) }} E(t)
$$
$E (t) = f(t) - p_{n} (t)$ なので次を得る。
$$
f(t) - \sum_{j=0}^{n} f( x_{j} ) l_{j} (t) = {{ (t - x_{0}) \cdots (t - x_{n}) } \over { (n+1)! }} f^{(n+1)} ( \xi )
$$
■
コード
以下はポストで式的に示した方法と同様に多項式近似 $p : \mathbb{R} \to \mathbb{R}$ を Julia で実装したものである。
using Plots, Polynomials
function interpolate(xs, ys)
V = [x^(i-1) for x in xs, i in 1:length(xs)]
coeffs = V \ ys
return Polynomial(coeffs)
end
foo = sort(rand(5))
bar = randn(5)
p = interpolate(foo, bar)
scatter(foo, bar, legend = :none)
plot!(0:0.01:1, p.(0:0.01:1))
Atkinson. (1989). An Introduction to Numerical Analysis(2nd Edition): p131. ↩︎