# neco-spline -- 数理詳細

この文書では neco-spline が使う自然三次スプライン補間の導出をまとめる。中心になるのは、補間点をただ折れ線で結ぶのではなく、滑らかにつなぎながら各点を通る曲線を作ろうとすると、なぜ三重対角連立方程式が現れるかである。使い方と API は [README-ja.md](README-ja.md) を参照。

ここで扱う前提は、補間点 $(x_i, y_i)$ の $x_i$ が狭義単調増加であること、各区間で三次多項式を使うこと、両端で二次導関数を 0 にする自然境界条件を採ることである。実装はこの前提のもとで係数を一度求め、各区間の多項式を評価する。

## 区分的三次多項式

補間で最初に欲しいのは、与えられた点を通ることと、点の間でも不自然に折れないことである。一次式では区間ごとに折れ曲がり、二次式 1 本では多くの点を柔軟に通しにくい。そこで各区間ごとに三次多項式を置き、それらを滑らかにつなぐ。

各区間 $[x_i, x_{i+1}]$ でスプラインを

$$
S_i(x) = a_i + b_i(x - x_i) + c_i(x - x_i)^2 + d_i(x - x_i)^3
$$

と書く。

補間条件 $S_i(x_i) = y_i$ を課すと

$$
a_i = y_i
$$

が直ちに決まる。したがって未知数は実質的に $b_i, c_i, d_i$ である。

## 導関数

区間を滑らかにつなぐには、隣接区間で傾きや曲がり方がどうつながるかを調べる必要がある。そのため一次導関数と二次導関数を使う。

一次導関数と二次導関数は

$$
S_i'(x) = b_i + 2c_i(x - x_i) + 3d_i(x - x_i)^2
$$

$$
S_i''(x) = 2c_i + 6d_i(x - x_i)
$$

となる。左端 $x = x_i$ では

- $S_i'(x_i) = b_i$
- $S_i''(x_i) = 2c_i$

である。

## 自然境界条件

両端の条件を決めないと、内部でどれだけ滑らかさを課しても解は一意に決まらない。自然三次スプラインでは、端で余計な曲げモーメントを持たない最も素直な条件として、両端の二次導関数を 0 にする。

したがって両端で

$$
S''(x_0) = 0,
\qquad
S''(x_n) = 0
$$

を課す。

左端では $S_0''(x_0) = 2c_0$ なので

$$
c_0 = 0
$$

が従う。同様に右端から

$$
c_n = 0
$$

が得られる。したがって未知数は内部節点に対応する

$$
c_1, c_2, \ldots, c_{n-1}
$$

へ縮約される。

## 連続性条件から三重対角系を導く

各区間で三次式を使うだけでは、区間のつなぎ目で傾きや曲率が飛ぶ可能性がある。そこで隣接区間の接続点で、値、傾き、曲がり方が一致するように条件を課す。

スプラインには隣接区間で

1. $C^0$ 連続性: 値が一致する
2. $C^1$ 連続性: 一次導関数が一致する
3. $C^2$ 連続性: 二次導関数が一致する

を課す。

以下では

$$
h_i = x_{i+1} - x_i,
\qquad
f_i = \frac{y_{i+1} - y_i}{h_i}
$$

とおく。

### ステップ 1: $C^0$ 連続性

まず値がつながることを使って、区間係数の一部を消去する。$x_{i+1}$ で補間条件を使うと

$$
a_i + b_i h_i + c_i h_i^2 + d_i h_i^3 = y_{i+1}
$$

である。$a_i = y_i$ を代入すると

$$
d_i = \frac{y_{i+1} - y_i - b_i h_i - c_i h_i^2}{h_i^3}
$$

を得る。

### ステップ 2: $C^2$ 連続性

次に、曲がり方がつながる条件を使う。$S_i''(x_{i+1}) = S_{i+1}''(x_{i+1})$ から

$$
2c_i + 6d_i h_i = 2c_{i+1}
$$

したがって

$$
d_i = \frac{c_{i+1} - c_i}{3h_i}
$$

である。

### ステップ 3: $b_i$ を $c_i$ で表す

ここまでで $d_i$ が消えたので、さらに $b_i$ も $c_i$ で表す。ステップ 1 にステップ 2 を代入すると

$$
y_i + b_i h_i + c_i h_i^2 + \frac{c_{i+1} - c_i}{3h_i} h_i^3 = y_{i+1}
$$

これを整理して

$$
b_i h_i + \frac{h_i^2}{3}(2c_i + c_{i+1}) = y_{i+1} - y_i
$$

さらに $f_i$ を使うと

$$
b_i = f_i - \frac{h_i}{3}(2c_i + c_{i+1})
$$

を得る。

### ステップ 4: $C^1$ 連続性

最後に傾きの連続性を使うと、未知数を $c_i$ だけに絞った連立方程式が現れる。$S_i'(x_{i+1}) = S_{i+1}'(x_{i+1})$ から

$$
b_i + 2c_i h_i + 3d_i h_i^2 = b_{i+1}
$$

である。ここへステップ 2 と 3 を代入すると

$$
h_i c_i + 2(h_i + h_{i+1}) c_{i+1} + h_{i+1} c_{i+2}
= 3(f_{i+1} - f_i)
$$

を得る。添え字を振り直すと、内部節点 $i = 1, \ldots, n-1$ に対して

$$
h_{i-1} c_{i-1} + 2(h_{i-1} + h_i)c_i + h_i c_{i+1}
= 3(f_i - f_{i-1})
$$

となる。各節点の係数 $c_i$ が隣接する $c_{i-1}, c_{i+1}$ としか結び付かないので、ここで三重対角構造が現れる。

## 三重対角行列

$c_0 = c_n = 0$ を使うと、未知数は $c_1, \ldots, c_{n-1}$ の $n-1$ 個になる。行列形式では

$$
\begin{pmatrix}
2(h_0+h_1) & h_1 & & \\
h_1 & 2(h_1+h_2) & h_2 & \\
& \ddots & \ddots & \ddots \\
& & h_{n-2} & 2(h_{n-2}+h_{n-1})
\end{pmatrix}
\begin{pmatrix}
c_1 \\ c_2 \\ \vdots \\ c_{n-1}
\end{pmatrix}
=
\begin{pmatrix}
3(f_1 - f_0) \\ 3(f_2 - f_1) \\ \vdots \\ 3(f_{n-1} - f_{n-2})
\end{pmatrix}
$$

となる。

この行列は対角成分が正で、自然三次スプラインの標準設定では解が一意に定まる。

## Thomas アルゴリズム

三重対角行列の利点は、一般の密行列として解く必要がないことである。実装はこの構造を使って Thomas アルゴリズムで $O(n)$ で解く。副対角、対角、上対角、右辺をそれぞれ $l_i, m_i, u_i, r_i$ と書く。

### 前進消去

$$
w = \frac{l_i}{m_{i-1}}
$$

$$
m_i \leftarrow m_i - w u_{i-1}
$$

$$
r_i \leftarrow r_i - w r_{i-1}
$$

を順に適用して、副対角成分を消す。

### 後退代入

最後の未知数を

$$
c_{n-1} = \frac{r_{n-1}}{m_{n-1}}
$$

で求め、残りを

$$
c_i = \frac{r_i - u_i c_{i+1}}{m_i}
$$

で後ろから復元する。

## 係数の復元

$c_i$ が求まれば、各区間の係数は

$$
a_i = y_i
$$

$$
b_i = f_i - \frac{h_i}{3}(2c_i + c_{i+1})
$$

$$
d_i = \frac{c_{i+1} - c_i}{3h_i}
$$

で決まる。したがって実装は、構築時にこの係数列を保存し、評価時には該当区間の三次多項式をそのまま計算すればよい。
