3次スプライン曲線の導出と計算（１）
はじめに
データのグラフに対して補助的に直線、曲線を書き込んで傾向を調べる場合、最小二乗法で近似直線や2次曲線のあてはめをするのが統計的に意味があるので一般的なやり方です。ただ、それなりにきれいな線を描きたいという見た目が欲しい時もあります。そういう場合は、補間曲線がよいので、今回は3次スプライン曲線について、曲線の式の導出とプログラムの作成を行いました（Excel VBA）。
3次スプライン曲線の定義
3次スプライン曲線は、与えられた点を必ず通り、点間をなめらかな曲線（条件を満たす3次曲線）でつなぐものです。物理的な意味などは各種のテキストを見ていただくとして、具体的には次のように定義されます(Fig. 1)。
点$${(x_i, y_i) i=0 ～ n}$$ の点に対して、
区間$${(x_i, y_i)～(x_{i+1}, y_{i+1})}$$ の曲線を $${y = g_i(x)}$$ とする。点$${(x_i, y_i), (x_{i+1}, y_{i+1})}$$ は曲線 $${y = g_i(x)}$$ を通る。
つまり $${y_i = g_i(x_i) , y_{i+1} = g_i(x_{i+1}), i = 0 ～ n-2}$$ 。点$${(x_{i+1}, y_{i+1})}$$ において、前後の区間の曲線の1階の微分は等しい。
つまり$${g'_i(x_{i+1}) = g'_{i+1}(x_{i+1}), i = 0 ～ n-2}$$ 。点$${(x_{i+1}, y_{i+1})}$$ において、前後の区間の曲線の２階の微分は等しい。
つまり $${g''_i(x_{i+1}) = g''_{i+1}(x_{i+1})}$$ 。終点$${(x_0, y_0)}$$ と端点$${(x_n, y_n)}$$ での２階の微分は０とする。
つまり $${g''_0(x_0) = g''_{n-1}(x_n) = g''_n(x_n) = 0}$$ 。
曲線 $${y = g_n(x_n) }$$ は定義されていませんが、あるとしてこれも２階の微分は０とします。
3次スプライン曲線の導出
3次曲線$${g_i(x)}$$を次のように定義します。
$$
y = g_i(x) = a_{3,i} (x - x_i)^3 + a_{2,i} (x - x_i)^2 + a_{1,i} (x - x_i)+ a_{0,i} …(1)
$$
$$
y' = g'_i(x) = 3a_{3,i} (x - x_i)^2 + 2 a_{2,i} (x - x_i) + a_{1,i} …(2)
$$
$$
y'' = g''_i(x) = 6a_{3,i} (x - x_i) + 2 a_{2,i} …(3)
$$
定義２と(1)より
$$
y_i = g_i(x_i)=a_{0,i} …(4) \\y_{i+1} = g_i(x_{i+1}) = h_i^3a_{3,i} +h_i^2 a_{2,i} + h_i a_{1,i} + a_{0,i} …(5)\\ここで h_i = x_{i+1} - x_i
$$
定義３と(2)より
$$
g'_i(x_{i+1}) = g'_{i+1}(x_{i+1}) \\ =3h_i^2 a_{3,i} + 2h_i a_{2,i} + a_{1,i} = a_{1,i+1} …(6)
$$
定義４と(3)より
$$
g''_i(x_{i+1}) = g''_{i+1}(x_{i+1}) \\ = 6h_i a_{3,i} + 2a_{2,i} = 2a_{2,i+1} …(7)
$$
ここからは係数$${a_{k,i}}$$をデータ$${(x_i,y_i)}$$で表すことを考えたいが見通しが悪いので、いったん$${g''}$$も使って表すことを考える。
(7)は
$$
a_{3,i} = \frac{2a_{2,i+1} - 2a_{2,i}}{6h_i} = \frac{g''_{i+1}-g''_i}{6h_i} …(8)
$$
ここで、$${g''_i(x_i) を g''_i}$$ と表記する。
つぎに(5)は、(6)、(7)、(8)を使って
$$
y_{i+1} = h_i^3a_{3,i} +h_i^2 a_{2,i} + h_i a_{1,i} + a_{0,i}\\ =h_i^3 \frac{g''_{i+1}-g''_i}{6h_i} + h_i^2 \frac{g''_i}{2} + h_ia_{1,i} + y_i …(9)
$$
よって
$$
a_{1,i} = \frac{1}{h_i}\Big(y_{i+1} - y_i - h_i^2 \frac{g''_{i+1} - g''_i}{6} - h_i^2 \frac{g''_i}{2}\Big)\\ = \frac{y_{i+1} - y_i}{h_i} - \frac{h_i}{6}\Big(2g''_i + g''_{i+1}\Big) …(10)
$$
さらに(6)は(7)、(8)、(10)より
$$
3h_i^2 a_{3,i} + 2h_i a_{2,i} + a_{1,i} = a_{1,i+1} \\3\frac{h_i^2}{6h_i}(g''_{i+1}-g''_i) + h_ig''_i + \frac{g''_{i+1}-g''_i}{h_i} - \frac{h_i}{6}(2g''_i + g''_{i+1}) = \frac{y_{i+2} - y_{i+1}}{h_{i+1}} - \frac{h_{i+1}}{6}(2g''_{i+1} + g''_{i+2})\\
$$
整理すると
$$
h_ig''_i + 2(h_{i+1}+ h_i) g''_{i+1} + h_{i+1} g''_{i+2} = K_{i+1} …(11)
$$
$$
ここで\\ K_i= 6\frac{y_{i+1}-y_i}{h_i}- 6\frac{y_i-y_{i-1}}{h_{i-1}}、 i=1～n-1\\ とおいた。
$$
$${i=0、および i=n-2}$$ の時、(11)は定義５より
$$
2(h_1+ h_0) g''_1 + h_1 g''_2 = K_1\\h_{n-2}g''_{n-2} + 2(h_{n-1}+ h_{n-2}) g''_{n-1} = K_{n-1}
$$
なので、$${g''_i}$$について整理して行列に書くと
$$
\begin{bmatrix}H_0 & h_1 &0… \\h_1 & H_1 & h_2 & 0…\\0& h_2 & H_2 & h_3 & 0…\\……&&&&&……\\……&&&&&……\\…&0& h_{n-4} & H_{n-4} & h_{n-3} &0\\…&0&0 & h_{n-3} & H_{n-3} & h_{n-2} \\…&0&0&0 & h_{n-2} & H_{n-2} \end{bmatrix}\begin{bmatrix}g''_1\\g''_2\\g_3\\…\\…\\g''_{n-3}\\g''_{n-2}\\g''_{n-1}\end{bmatrix}=\begin{bmatrix}K_1\\K_2\\K_3\\…\\…\\K_{n-3}\\K_{n-2} \\K_{n-1}\end{bmatrix} …(12)
$$
ここで、$${H_i = 2(h_i + h_{i+1})}$$ とおいた。
$${h_i、H_i、K_i}$$ はすべて $${x_i、y_i}$$ から求まるので、(12)を解けば$${g''_i}$$ が求まる。
これより、$${a_{3,i}は(8)、a_{2,i}は(7)、a_{1,i}は(10)、a_{0,i}は(4)}$$ によって求められる。
次回に上記をプログラムにしたものを載せます（ExcelVBA）。
以上
参考文献
１．スプライン補間、http://www.yamamo10.jp/yamamoto/lecture/2006/5E/interpolation/interpolation_html/node3.html?utm_source=pocket_mylist、2022.06.16閲覧
２．3次スプライン補完、https://qiita.com/YudaiSadakuni/items/e886ad237a1171e069a1?utm_source=pocket_mylist、2022.06.16閲覧
３．3次スプライン補間で軌跡生成：３次多項式のパラメータを求める（その１）および関連ページ、https://tajimarobotics.com/cubic-spline-interpolation-2/?utm_source=pocket_mylist、2022.06.16閲覧
４．3次スプライン補間を実装する、https://swdrsker.hatenablog.com/entry/2018/08/14/182915?utm_source=pocket_mylist、2022.06.16閲覧
５．伊理正雄、藤野和建. 数値計算の常識. 共立出版. 1985.
６．奥村晴彦. C言語によるアルゴリズム辞典. 技術評論社. 1991.
７．ニューメリカルレシピ・イン・シー日本語版. W.H.Pressら. 技術評論社. 1993.
#３次スプライン補間 , #Excel . #VBA
応援してやろうということで、お気持ちをいただければ嬉しいです。もっと勉強したり、調べたりする糧にしたいと思います。