ここ最近、個人的に"Numerical Recipies in C"¹ という、C言語の数値計算の本の一部をRustで書き直すことにハマっています。久しぶりに積分のコードを実装し、その際に補間・補外の章を読むことになりました。この記事では、その際に理解に詰まったNevilleのアルゴリズムについて自分なりにまとめています。

Nevilleのアルゴリズムとは？

互いに異なる$N$個の点が与えられているとき、それらの点を通るLagrange補間多項式を計算し、任意の$x$における補間多項式の値を返すアルゴリズムです。

文字や式の説明は後に回しますが、次の漸化式がNevilleのアルゴリズムの鍵になります：

$$\begin{align} P_{i(i+1)\cdots(i+m)} = \frac{ (x - x_{i + m}) P_{i(i + 1) \cdots (i + m - 1)} - (x_i - x) P_{i(i + 1) \cdots (i + m)} }{ x_i - x_{i + m} } \end{align}$$

Lagrange補間多項式

$N$個の点$\{(x_1, y_1), \ldots, (x_N, y_N)\}$が与えられているとき、Lagrange補間多項式は

$$\begin{align} P(x) &= \frac{ (x - x_2) (x - x_3) \cdots (x - x_N) }{ (x_1 - x_2) (x_1 - x_3) \cdots (x_1 - x_N) } y_1 + \frac{ (x - x_1) (x - x_3) \cdots (x - x_N) }{ (x_2 - x_1) (x_2 - x_3) \cdots (x_2 - x_N) } y_2 \notag \\ &+ \cdots + \frac{ (x - x_1) (x - x_2) \cdots (x - x_{N-1}) }{ (x_N - x_1) (x_N - x_2) \cdots (x_N - x_{N-1}) } y_N \\ \end{align}$$

で与えられます。これは一意的です。

$L_i(x)$を

$$L_i(x) = \prod_{j = 1 (j \neq i)} \frac{ (x - x_j)}{(x_i - x_j)}$$

のように表現しておけば、これは$L_i(x_j) = \delta_{ij}$（$\delta_{ij}$はKroneckerのデルタ）を満たします。 よって、$P(x) = \sum_{i} L_i y_i$のように書くこともできますね。

例：4点の場合でイメージをつかむ

残念ながら私は数式だけであまりイメージができなかったので、まずは簡単な場合を考えて意味ます。 フリーレンも言っていますね。

魔法はイメージの世界だ

帰納的な図（Nevilleの三角表）

$$\begin{array}{cccccccc} x_1:& y_1 = P_1 & & & \\ & & P_{12} & & \\ x_2:& y_2 = P_2 & & P_{123} & \\ & & P_{23} & & P_{1234}\\ x_3:& y_3 = P_3 & & P_{234} & \\ & & P_{34} & & \\ x_4:& y_4 = P_4 & & & \\ \end{array}$$

今、$\{ (x_1, y_1), (x_2, y_2), (x_3, y_3), (x_4, y_4) \}$の4点について考えていみましょう。まずは0次の補間は定数で、$P_1 = y_1, P_2 = y_2, P_3 = y_3, P_4 = y_4$となります。さてそれぞれの隣り合う点同士の1次の補間関数は例えば、

$$\begin{align} P_{12}(x) = \frac{(x - x_2) }{(x_1 - x_2)} y_1 + \frac{(x - x_1) }{(x_2 - x_1)} y_2 \end{align}$$

となります。明らかに1次関数であり、$P_{12}(x_1) = y_1$と$P_{12}(x_2) = y_2$を満たすので、$(x_1, y_1), (x_2, y_2)$を通る直線を想像できるでしょう。

このセンスで、$P_{12}(x_1)$と$P_{23}(x_3)$を通る2次曲線も構築できるでしょう。つまり、

$$\begin{align} P_{123}(x) = \frac{(x - x_3) }{(x_1 - x_3)} P_{12}(x) + \frac{(x - x_1) }{(x_3 - x_1)} P_{23}(x) \end{align}$$

ということです。新しくできた2次関数（1次関数にさらに$x$の項をかけているので）は明らかに$P_{123}(x_1) = y_1, P_{123}(x_2) = y_3$です。 もちろん$P_{123}(x_2)$のときも

$$\begin{align*} P_{123}(x_2) &= \frac{(x_2 - x_3) }{(x_1 - x_3)} P_{12}(x_2) + \frac{(x_2 - x_1) }{(x_3 - x_1)} P_{23}(x_2) \\ &= \frac{(x_2 - x_3)y_2 - (x_2 - x_1)y_2}{(x_1 - x_3)} \\ &= y_2 \end{align*}$$

です。

このようにして、1次補間を何度も繰り返すことで、4次の補完関数まで求めることができます。

---

Footnotes

1. William H. Press et al （丹慶勝市 et al.）, "Numerical Recipes in C", 1988 (日本語版 1993) ↩