補間と微分

# 微分可能な補間

![interpolation.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/706455/d94b9a18-d219-132f-b440-b5f307e87a11.png)

補間とは，$t$-$x$空間で$N$個の点列$\{(t_{i},x_{i})\}~(i=0\sim N-1)$が与えられた時に，それらをつなぐ関数$x(t)~~(t_{0}\leq t\leq t_{N-1})$を求めることです．点列が何らかの連続プロセスによる出力であると考えれば，出力から元の連続プロセスを逆推定する単純な方法，とも解釈できます．

以下，$t_{0}<\cdots<t_{N-1}$を仮定します．例えば点列を次のような折れ線（区分線形関数）で結ぶのは，最も簡単な補間です．
```math
x(t)=\frac{x_{i+1}-x_{i}}{t_{i+1}-t_{i}}(t-t_{i})+x_{i}
\quad\mbox{(ただし $t_{i}\leq t<t_{i+1}$のとき)}
```

しかしこの関数は，$t=t_{i}$においては微分不可能です．元のプロセスが滑らかさを伴うものであると分かっているならば，この推定方法はあまり素性が良くありません．

本記事は，滑らかな補間を行う代表的方法として
 - Lagrange補間
 - 3次スプライン補間
 - 区分3次エルミート補間

を採り上げ概説します．いずれも有名な方法なので既に優れた解説がウェブ上に多く存在しますが，筆者自身が納得できる説明を試みます．

# Lagrange補間

例えば2点与えられたらそれらを通る直線（1次関数）を，3点与えられたらそれらを通る放物線（2次関数）をそれぞれ一意に引くことができます．一般的に，$N$個の点列が与えられたら，それらを全て通る$N-1$次関数を一意に定めることができます．これがLagrange補間です．具体的には$x(t)$は次式で表せます．
```math
x(t)=\sum_{i=0}^{N-1}x_{i}f_{i}(t)
```
ただし，
```math
f_{i}(t)=
\frac{\displaystyle\prod_{j=0,j\neq i}^{N-1}(t-t_{j})}
{\displaystyle\prod_{j=0,j\neq i}^{N-1}(t_{i}-t_{j})}
=
\frac{(t-t_{1})\cdots(t-t_{i-1})(t-t_{i+1})\cdots(t-t_{N-1})}
{(t_{i}-t_{1})\cdots(t_{i}-t_{i-1})(t_{i}-t_{i+1})\cdots(t_{i}-t_{N-1})}
```
とおきました．$f_{i}(t)$は$N-1$次関数なので，それらの和をとった$x(t)$も$N-1$次関数です．また，$f_{i}(t)$は$t=t_{i}$のとき$1$，それ以外の$t=t_{j} (j\neq i)$のとき$0$となることから，上式の正しさが分かります．これは次のようにも書けます．
```math
x(t)=\sum_{i=0}^{N-1}z_{i}p_{i}(t)
```
ただし，
```math
z_{i}=
\frac{x_{i}}{\displaystyle\prod_{j=0,j\neq i}^{N-1}(t_{i}-t_{j})}
,\quad
p_{i}(t)=\prod_{j=0,j\neq i}^{N-1}(t-t_{j})
```
とそれぞれおきました．$z_{i} (i=1,\cdots,N-1)$は定数であり，予め計算しておけることに注意して下さい．

Lagrange補間の1階微分は
```math
\dot{x}(t)=\sum_{i=0}^{N-1}z_{i}\dot{p}_{i}(t)
```
であり，関数の積の微分より
```math
\dot{p}_{i}(t)
=\sum_{k=0}^{N-1}\prod_{j=0,j\neq i,j\neq k}^{N-1}(t-t_{j})
```
です．これより，
```math
\dot{x}(t)
=\sum_{i=0}^{N-1}z_{i}
\sum_{k=0}^{N-1}\prod_{j=0,j\neq i,j\neq k}^{N-1}(t-t_{j})
=\sum_{i=0}^{N-2}\sum_{j=i+1}^{N-1}
(z_{i}+z_{j})\prod_{k=0,k\neq i,k\neq j}^{N-1}(t-t_{k})
```
ですので，定義式通りに計算する場合に比べて計算量を半分程度にすることができます．

2階微分も，同様に考えれば
```math
\dot{x}(t)
=\sum_{i=0}^{N-3}\sum_{j=i+1}^{N-2}\sum_{k=j+1}^{N-1}
2(z_{i}+z_{j}+z_{k})\prod_{l=0,l\neq i,l\neq j,l\neq k}^{N-1}(t-t_{l})
```
と計算できます．

# 3次スプライン補間

Lagrange補間は，与えられた全ての点をただ一つの多項式関数で補間しますので，点列の配置によっては蛇行した曲線が引かれてしまうこともあります．これは[Runge現象](https://ja.wikipedia.org/wiki/%E3%83%AB%E3%83%B3%E3%82%B2%E7%8F%BE%E8%B1%A1)としてよく知られています．点が増えれば次数が上がり，次数が上がると蛇行する可能性も上がってしまうので，なるべく次数の低い多項式関数を使いたいところです．

点列が機械システムに由来するものならば，実用上は2階微分（加速度）まで考慮すれば（つまり3次関数を用いれば）十分なことが多いでしょう．ただしこの場合，全ての点をただ一つの関数で補間するのではなく，区分ごとに異なる関数を用いることになります．これはあたかも雲形定規を使って点列を滑らかにつなぐ曲線を引く作業に似ています．スプライン補間はこのような発想によるものです．

今，仮に全ての$t=t_{i}$における関数の1階微分値が$v_{i}$と与えられているものとしましょう．すなわち区間$\mathrm{I}\_{i}=\left[t\_{i},t\_{i+1}\right]$は，次の境界条件を満たす3次関数で表されることになります．
```math
\displaylines{
x^{[i]}(t_{i})=x_{i}
\\
\dot{x}^{[i]}(t_{i})=v_{i}
\\
x^{[i]}(t_{i+1})=x_{i+1}
\\
\dot{x}^{[i]}(t_{i+1})=v_{i+1}
}
```
このような関数は次式で表せます．
```math
x^{[i]}(t)
=x_{i}\phi\left(\frac{t_{i+1}-t}{\varDelta t_{i}}\right)
+x_{i+1}\phi\left(\frac{t-t_{i}}{\varDelta t_{i}}\right)
-v_{i}\varDelta t_{i}\psi\left(\frac{t_{i+1}-t}{\varDelta t_{i}}\right)
+v_{i+1}\varDelta t_{i}\psi\left(\frac{t-t_{i}}{\varDelta t_{i}}\right)
```
ただし，$\varDelta t_{i}=t_{i+1}-t_{i}$とおきました．また，$\phi(\tau)$，$\psi(\tau)$は次の境界条件を満たす3次関数です．
```math
\begin{cases}
\phi(0)=0 \\
\phi(1)=1 \\
\psi(0)=0 \\
\psi(1)=0 \\
\end{cases}
,\quad
\begin{cases}
\dot{\phi}(0)=0 \\
\dot{\phi}(1)=0 \\
\dot{\psi}(0)=0 \\
\dot{\psi}(1)=1 \\
\end{cases}
```
このことは，
```math
\dot{x}^{[i]}(t)
=-\frac{x_{i}}{\varDelta t_{i}}\dot{\phi}\left(\frac{t_{i+1}-t}{\varDelta t_{i}}\right)
+\frac{x_{i+1}}{\varDelta t_{i}}\dot{\phi}\left(\frac{t-t_{i}}{\varDelta t_{i}}\right)
+v_{i}\dot{\psi}\left(\frac{t_{i+1}-t}{\varDelta t_{i}}\right)
+v_{i+1}\dot{\psi}\left(\frac{t-t_{i}}{\varDelta t_{i}}\right)
```
より
```math
\displaylines{
x^{[i]}(t_{i})=x_{i}\phi(1)+x_{i+1}\phi(0)-v_{i}\varDelta t_{i}\psi(1)+v_{i+1}\varDelta t_{i}\psi(0)
\\
x^{[i]}(t_{i+1})=x_{i}\phi(0)+x_{i+1}\phi(1)-v_{i}\varDelta t_{i}\psi(0)+v_{i+1}\varDelta t_{i}\psi(1)
\\
\dot{x}^{[i]}(t_{i})=-\frac{x_{i}}{\varDelta t_{i}}\dot{\phi}(1)+\frac{x_{i+1}}{\varDelta t_{i}}\dot{\phi}(0)+v_{i}\dot{\psi}(1)+v_{i+1}\dot{\psi}(0)
\\
\dot{x}^{[i]}(t_{i+1})=-\frac{x_{i}}{\varDelta t_{i}}\dot{\phi}(0)+\frac{x_{i+1}}{\varDelta t_{i}}\dot{\phi}(1)+v_{i}\dot{\psi}(0)+v_{i+1}\dot{\psi}(1)
}
```
に上記を代入すれば正しさが分かります．

$\phi(\tau)$，$\psi(\tau)$は比較的簡単に求まります．まず$\phi(\tau)$について，$\dot{\phi}(0)=\dot{\phi}(1)=0$より
```math
\dot{\phi}(\tau)=A\tau(\tau-1)
```
と表せます（$\phi(\tau)$が3次関数なので$\dot{\phi}(\tau)$は2次関数です）．ただし$A$は未知定数です．これを1階積分し，$\phi(0)=0$，$\phi(1)=1$に代入することで$A=-6$を得ます．つまり
```math
\phi(\tau)=3\tau^{2}-2\tau^{3}
```
です．$\psi(\tau)$についても同様で，$\psi(0)=\psi(1)=\dot{\psi}(0)=0$より
```math
\psi(\tau)=B\tau^{2}(\tau-1)
```
と表せます．これを1階微分し，$\dot{\psi}(1)=1$に代入することで$B=1$を得ます．つまり
```math
\psi(\tau)=\tau^{3}-\tau^{2}
```
です．

以上のように，両端の値と微分係数から決まる関数を区分3次エルミート関数といいます．また，これによって引かれる曲線はFerguson曲線と呼ばれます．

3次スプライン補間は，各区分の境界における2階微分が可能となる，つまり$x^{[i]}(t\_{i+1})=x^{[i+1]}(t\_{i+1})$を満たすよう$v_{i}$を求める方法です．$x^{[i]}(t)$を2階微分すると
```math
\ddot{x}^{[i]}(t)
=\frac{x_{i}}{\varDelta t_{i}^{2}}\ddot{\phi}\left(\frac{t_{i+1}-t}{\varDelta t_{i}}\right)
+\frac{x_{i+1}}{\varDelta t_{i}^{2}}\ddot{\phi}\left(\frac{t-t_{i}}{\varDelta t_{i}}\right)
-\frac{v_{i}}{\varDelta t_{i}}\ddot{\psi}\left(\frac{t_{i+1}-t}{\varDelta t_{i}}\right)
+\frac{v_{i+1}}{\varDelta t_{i}}\ddot{\psi}\left(\frac{t-t_{i}}{\varDelta t_{i}}\right)
```
で，
```math
\displaylines{
\ddot{\phi}(\tau)=6-12\tau
\quad\Rightarrow\quad
\ddot{\phi}(0)=6,\quad \ddot{\phi}(1)=-6
\\
\ddot{\psi}(\tau)=6\tau-2
\quad\Rightarrow\quad
\ddot{\psi}(0)=-2,\quad \ddot{\psi}(1)=4
}
```
ですから，
```math
\displaylines{
\ddot{x}^{[i]}(t_{i+1})
=\frac{6x_{i}}{\varDelta t_{i}^{2}}-\frac{6x_{i+1}}{\varDelta t_{i}^{2}}
+\frac{2v_{i}}{\varDelta t_{i}}+\frac{4v_{i+1}}{\varDelta t_{i}}
\\
\ddot{x}^{[i+1]}(t_{i+1})
=-\frac{6x_{i+1}}{\varDelta t_{i+1}^{2}}+\frac{6x_{i+2}}{\varDelta t_{i+1}^{2}}
-\frac{4v_{i+1}}{\varDelta t_{i+1}}-\frac{2v_{i+2}}{\varDelta t_{i+1}}
}
```
より，次式を得ます．
```math
\frac{v_{i+2}}{\varDelta t_{i+1}}
+2\left(\frac{1}{\varDelta t_{i+1}}+\frac{1}{\varDelta t_{i}}\right)v_{i+1}
+\frac{v_{i}}{\varDelta t_{i}}
=3\left(
\frac{x_{i+2}-x_{i+1}}{\varDelta t_{i+1}^{2}}+
\frac{x_{i+1}-x_{i}}{\varDelta t_{i}^{2}}
\right)
```
これは$\{v_{i}\}$に関する三項方程式となっています．未知数が$v_{0},\cdots,v_{N-1}$の$N$個であるのに対し，この方程式は両端を除く$N-2$個の点において成り立つので，このままでは式が2本足りません．そこでこれに，両端における適当な境界条件を加える必要があります．例えば点$0$における1階微分値が$v$と与えられているならば，次を加えます．
```math
v_{0}=v
```
あるいは点$0$における2階微分値が$a$と与えられているならば，次を加えます．
```math
\frac{v_{0}}{\varDelta t_{0}}+\frac{2v_{1}}{\varDelta t_{0}}
=\frac{1}{2}a
+3\left(\frac{x_{1}}{\varDelta t_{0}^{2}}-\frac{x_{0}}{\varDelta t_{0}^{2}}\right)
```
点$N-1$についても同様に考えます．なお、$v=0$の場合を**固定端条件**，$a=0$の場合を**自由端条件**と呼びます．両端の境界条件2つが決まれば方程式の数と未知数の数が一致しますので，全ての$v_{i}$が求まります．$x(t)$，$\dot{x}(t)$，$\dot{x}(t)$は次式となります．
```math
\begin{cases}
x(t)=x^{[i]}(t) \\
\dot{x}(t)=\dot{x}^{[i]}(t) \\
\ddot{x}(t)=\ddot{x}^{[i]}(t) \\
\end{cases}
,\quad t\in\mathrm{I}_{i}
```
$t$がどの区間に所属するかは，2分法により効率的に調べることができます．

# 区分3次エルミート補間（Piecewise Cubic Hermite Interpolating Polynomial, PCHIP）

PCHIPは，与えられた値列が単調に変化する，すなわち$x_{0}\leq\cdots\leq x_{N_1}$もしくは$x_{0}\geq\cdots\geq x_{N_1}$である場合に，補間曲線がその単調性を損ねないよう工夫したもので，次の文献が出典です．

F. N. Fritsch and R. E. Carlson, Monotone Piecewise Cubic Interpolation, SIAM Journal of Numerical Analysis, Vol. 17, No. 2, pp. pp. 238-246, 1980.

3次スプライン補間のときと同様に，仮に全ての$t=t_{i}$における関数の1階微分値が$v_{i}$と与えられているものとすれば，各区間は区分3次エルミート関数で表せるのでした．あとは，$x(t)$の単調性を保証するよう全ての$v_{i}$を求めれば良いわけです．

```math
\bar{v}_{i}\overset{\mathrm{def}}{=}\frac{x_{i+1}-x_{i}}{\varDelta t_{i}}
```
とおきましょう．
$x^{[i]}(t)$が単調であるためには，$v\_{i}$，$v\_{i+1}$，$\bar{v}\_{i}$が全て同符号であることが必要条件です．特に，$\bar{v}\_{i}=0$ならば$v\_{i}=v\_{i+1}=0$でなければなりません．以降は$\bar{v}\_{i}\neq 0$を仮定します．$t-t_{i}$を法として区分3次エルミート関数をべき級数展開すると，
```math
\displaylines{
x^{[i]}(t)
=\left(\frac{v_{i}+v_{i+1}-2\bar{v}_{i}}{\varDelta t_{i}^{2}}\right)(t-t_{i})^{3}
-\left(\frac{2v_{i}+v_{i+1}-3\bar{v}_{i}}{\varDelta t_{i}}\right)(t-t_{i})^{2}
+d_{i}(t-t_{i})+x_{i}
\\
\dot{x}^{[i]}(t)
=\left(\frac{3(v_{i}+v_{i+1}-2\bar{v}_{i})}{\varDelta t_{i}^{2}}\right)(t-t_{i})^{2}
-\left(\frac{2(2v_{i}+v_{i+1}-3\bar{v}_{i})}{\varDelta t_{i}}\right)(t-t_{i})
+v_{i}
\\
\ddot{x}^{[i]}(t)
=\left(\frac{6(v_{i}+v_{i+1}-2\bar{v}_{i})}{\varDelta t_{i}^{2}}\right)(t-t_{i})
-\left(\frac{2(2v_{i}+v_{i+1}-3\bar{v}_{i})}{\varDelta t_{i}}\right)
}
```

を得ます．これらより次が言えます．
 - $v\_{i}+v\_{i+1}-2\bar{v}\_{i}=0$ならば$x(t)$は2次関数となる．$\min\\{v\_{i},v\_{i+1}\\}\leq\dot{x}^{[i]}(t)\leq\max\\{v\_{i},v\_{i+1}\\}$なので，$v\_{i}$，$v\_{i+1}$，$\bar{v}\_{i}$が全て同符号であることは$x^{[i]}(t)$が単調であるための十分条件．
 - $v\_{i}+v\_{i+1}-2\bar{v}\_{i}\neq 0$ならば$\dot{x}^{[i]}(t)$は2次関数であり，$v\_{i}+v\_{i+1}-2\bar{v}\_{i}>0$ならば下に凸，$v\_{i}+v\_{i+1}-2\bar{v}\_{i}<0$ならば上に凸となる．$x_{i}< x_{i+1}$ （つまり$\bar{v}_{i}>0$） かつ$v_{i}+v_{i+1}-2\bar{v}_{i}<0$ならば，$0\leq\min\{v_{i},v_{i+1}\}\leq\dot{x}^{[i]}(t)$なので$x^{[i]}(t)$は単調増加．同様に$x_{i}>x_{i+1}$（つまり$\bar{v}_{i}<0$）かつ$v_{i}+v_{i+1}-2\bar{v}_{i}>0$ならば，$\dot{x}^{[i]}(t)\leq\max\left\\{v\_{i},v\_{i+1}\right\\}\leq 0$なので$x^{[i]}(t)$は単調減少．

ここで$\alpha_{i}\overset{\mathrm{def}}{=}v_{i}/\bar{v}\_{i}$，$\beta_{i}\overset{\mathrm{def}}{=}v_{i+1}/\bar{v}\_{i}$を定義すると，$v_{i}+v_{i+1}-2\bar{v}\_{i}=(\alpha_{i}+\beta_{i}-2)\bar{v}\_{i}$なので，$\alpha_{i}+\beta_{i}-2\leq 0$ならば$x^{[i]}(t)$は単調関数となることが分かります．

$\alpha_{i}+\beta_{i}-2>0$の場合を考えましょう．$\dot{x}^{[i]}(t)$は
```math
t^{*}=t_{i}+\frac{\varDelta t_{i}}{3}\left(\frac{2\alpha_{i}+\beta_{i}-3}{\alpha_{i}+\beta_{i}-2}\right)
```
において唯一の極値
```math
\dot{x}^{[i]}(t^{*})=
\left\{\alpha-\frac{1}{3}\frac{(2\alpha_{i}+\beta_{i}-3)^{2}}{\alpha_{i}+\beta_{i}-2}\right\}
\bar{v}_{i}
```
をとります．したがって，$x^{[i]}(t)$が区間$I_{i}$において単調であるための必要十分条件は，
 - $t^{*}\notin\mathrm{I}_{i}$
 - $t^{\*}\in\mathrm{I}\_{i}$かつ$\dot{x}^{[i]}(t^{*})$と$\bar{v}\_{i}$が同符号

のいずれかを満たすことです．これは
 - $2\alpha_{i}+\beta_{i}-3\leq 0$
 - $\alpha_{i}+2\beta_{i}-3\leq 0$
 - $\displaystyle\alpha_{i}-\frac{1}{3}\frac{(2\alpha_{i}+\beta_{i}-3)^{2}}{\alpha_{i}+\beta_{i}-2}\geq 0$

のいずれかを満たすことと等価に置き換えられます．このような$(\alpha_{i},\beta_{i})$の領域$\mathcal{M}$は下図のようになります．
![PCHIP_alpha_beta.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/706455/3d7751ff-deec-e060-2bf9-c7efd368a4d6.png)

したがって，全ての区間$I_{i}$（$i=1,\cdots,n-1$）について対応する$(\alpha_{i},\beta_{i})$が領域$\mathcal{M}$に収まるように$\\{v_{i}\\}$（$i=1,\cdots,n$）を配置すれば良いということになります．ただし，$I_{i}$と$I_{i-1}$の間には$\beta_{i-1}\bar{v}\_{i-1}=\alpha_{i}\bar{v}\_{i}=v_{i}$の関係があるので，$\alpha_{i}$を変えると$\beta_{i-1}$も変わってしまいます．そこで$\mathcal{M}$を，次の条件を満たす集合$\mathcal{I}$で置き換えることにします．
 - $(\alpha,\beta)\in\mathcal{I}$ならば，$0\leq\alpha^{\*}\leq\alpha$，$0\leq\beta^{\*}\leq\beta$を満たす任意の$(\alpha^{\*},\beta^{\*})$が$(\alpha^{\*},\beta^{\*})\in\mathcal{I}$を満たす．
 - $(\alpha,\beta)\in\mathcal{I}$ならば$(\beta,\alpha)\in\mathcal{I}$である．

このような$\mathcal{I}$の候補は無数にありますが，経験的に上図の太鎖線部のような4分の1円（半径3）の内側の値が比較的性質が良いと分かっているようです．この領域は次式で表せます．
```math
\mathcal{I}\overset{\mathrm{def}}{=}\left\{
(\alpha,\beta)\left|\alpha\geq 0,\beta\geq 0,\alpha^{2}+\beta^{2}\leq 3^{2}\right.\right\}
```

これを用いて，次のアルゴリズムで$v_{i} (i=0,\cdots,N-1)$を順番に決めていきます．
 1. $\\{v_{i}\\}$の適当な初期値を与える．
 1. $i=0$とする．
 1. $\bar{v}=(x_{i+1}-x_{i})/\varDelta t_{i}$
 1. $\alpha=v_{i}/\bar{v}$
 1. $\beta=v_{i+1}/\bar{v}$
 1. $l=\sqrt{\alpha^{2}+\beta^{2}}$
 1. もし$l\geq 3$ならば$v_{i}\leftarrow 3\alpha\bar{v}/l$，$v_{i+1}\leftarrow 3\beta\bar{v}/l$とする．
 1. $i=N-1$ならば終了．
 1. $i\leftarrow i+1$として3に戻る．

$\\{v_{i}\\}$の初期値は，例えば次のようなアルゴリズムで与えれば比較的滑らかな補間曲線が得られます．
 1. $v_{0}=x_{0}\cdot\frac{2t_{0}-t_{1}-t_{2}}{(t_{0}-t_{1})(t_{0}-t_{2})}+x_{1}\cdot\frac{t_{0}-t_{2}}{(t_{1}-t_{0})(t_{1}-t_{2})}+x_{2}\cdot\frac{t_{0}-t_{1}}{(t_{2}-t_{0})(t_{2}-t_{1})}$とする．
 1. $i=1$とする．
 1. $v_{i}=x_{i-1}\cdot\frac{t_{i}-t_{i+1}}{(t_{i-1}-t_{i})(t_{i-1}-t_{i+1})}+x_{i}  \cdot\frac{2t_{i}-t_{i-1}-t_{i+1}}{(t_{i}-t_{i-1})(t_{i}-t_{i+1})}+x_{i+1}\cdot\frac{t_{i}-t_{i-1}}{(t_{i+1}-t_{i-1})(t_{i+1}-t_{i})}$とする．
 1. $i\leftarrow i+1$とする．
 1. $i<N-1$ならば3に戻る．
 1. $v_{N-1}=x_{N-3}\cdot\frac{t_{N-1}-t_{N-2}}{(t_{N-3}-t_{N-2})(t_{N-3}-t_{N-1})}+x_{N-2}\cdot\frac{t_{N-1}-t_{N-3}}{(t_{N-2}-t_{N-3})(t_{N-2}-t_{N-1})}+x_{N-1}\cdot\frac{2t_{N-1}-t_{N-3}-t_{N-2}}{(t_{N-1}-t_{N-3})(t_{N-1}-t_{N-2})}$とする．

区分3次エルミート補間の構成要素は3次スプライン補間と同じですので，1階微分，2階微分も同様に得られます．ただし3次スプライン補間と異なり1階導関数が$t=t_{i}$において不連続となります．したがって，2階微分は右微分係数か左微分係数のどちらかを選ぶことになります．

# 比較

上述の補間法を定性的に比較するために，適当な点列を補間した例を示します．
![ip.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/706455/b7c72a26-6235-4624-b214-03dee5929010.png)
左からLagrange補間，3次スプライン補間，秋間スプライン補間，区分3次エルミート補間です．秋間スプライン補間というのはやはり区分3次エルミート関数を用いる方法で，人の手書きに近い性質を持つと言われています．理屈は幾何学的な議論を伴いますが、計算は簡単です。詳しくは次の原著論文をご参照下さい．

H. Akima, A new method of interpolation and smooth curve fitting based on local procedures, J. Assoc. Comput. Mach., Vol. 17, pp. 589--602, 1970.

また，上から順に補間曲線，1階微分，2階微分を示しています．元の点列は一番上の図に黒丸で記されています．5点あるのでLagrange補間は4次曲線となっています．このため2階微分まで滑らかです．この例だと3次スプライン補間と秋間スプライン補間の違いはあまりありませんが，前者は2階微分が連続になっています．点列が単調でないので，区分3次エルミート補間の見た目はLagrange補間とあまり変わらず，効果がよく分かりません．

別の例を示しましょう．Runge現象が現れる例としてよく知られているCauchy-Lorentz関数の補間です．
![ip.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/706455/61f67b3c-4f87-46b3-b9a6-a2216128907b.png)
Lagrange補間は外側に向かって大きく波打ってしまっています．残り三つの方法は比較的穏やかな変化を示しています．ただし，元の点列が左右対称に配置されているのに，これらの曲線はわずかながら対称性を損ねています．1階微分は三者三様といったところでしょうか．

最後に，区分3次エルミート補間の原著論文で示されていた例です．
![ip.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/706455/09db4bf8-cede-4bf7-a36e-949dbc797b8d.png)
最初7点は同じ値が続いているので，補間曲線はなるべく真っすぐになって欲しいところですが、Lagrange補間はまずだめ，3次スプライン補間も後段の点列に引っ張られて多少凹凸が出てしまっています．秋間スプライン補間と区分3次エルミート補間の結果は，好みが分かれるところかも知れません．

以上から分かるように，ここで示した補間法のいずれも万能ではありません．欠落した情報を無理矢理補うので，それぞれになにがしかの癖が出てしまうのは致し方ありません．元のプロセスの性質をなんとか理解した上で，それと親和性の高い方法を選択するのが現実的なやり方でしょう．

なお，以上の補間法は全て[ZM](https://github.com/mi-lib/zm)に実装されています．
 - Lagrange補間 [zm_ip_lagrange.c](https://github.com/mi-lib/zm/blob/main/src/zm_ip_lagrange.c)
 - 3次スプライン補間 [zm_ip_spline.c](https://github.com/mi-lib/zm/blob/main/src/zm_ip_spline.c)
 - 秋間スプライン補間 [zm_ip_akima.c](https://github.com/mi-lib/zm/blob/main/src/zm_ip_akima.c)
 - 区分3次エルミート補間 [zm_ip_pchip.c](https://github.com/mi-lib/zm/blob/main/src/zm_ip_pchip.c)