【S02】3次スプライン補間(変数の正規化ver.)

# 【S02】3次スプライン補間(変数の正規化ver.)
　
この記事は[アドカレ](https://qiita.com/advent-calendar/2023/hachinohe_kosen)に参加しています。

# スプライン補間
　3次スプライン補間とは`n`番目と`n+1`番目の点の間を`n-1`番目から`n+2`番目の4つの点で三次関数として補間するものでした。
![a_care08.jpeg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3005402/1fad63a4-00d1-44a6-7c87-49d5398ce1af.jpeg)

　実は3次スプライン補間にも様々な種類が存在し、使用用途によって端の加速度を変えたりすることがあります。マイナー処理です。
 
変数を正規化したver.のスプライン補間もそんなマイナーの1つですが、行列が綺麗なので僕は好きです。ベジェ曲線との互換性を持とうとしたら、こちらの方が便利だと思います。ということで、変数を正規化したver.のスプライン補間の紹介です。

# 三次関数の式から係数を求める(変数の正規化ver.)
　まずは点の数を`m+1`個として、$p_{0}$～$p_{m}$という点があったとします。このとき、$p_{n}$と$p_{n+1}$を補間する三次関数$S_{n}(t)$は以下のように表すことができます。
```math
\begin{align}
y&=S_{n}(t)　　(0 \leqq t \leqq 1)\\
y&=a_{n}t^3+b_{n}t^2+c_{n}t+d_{n}
\end{align}
```
$t$というのは三次関数の引数で、$t=0$のとき$S_{n}(0)=p_{n}$であり、$t=1$のとき$S_{n}(1)=p_{n+1}$となります。

## 係数d
<details><summary>係数$d_{n}$の求め方</summary>

$t=0$のときを考えます。
```math
\begin{align}
y&=S_{n}(t) \\
p_{n}&=S_{n}(0)
\end{align}
```
という関係を用いて、
```math
\begin{align}
p_{n}&=a_{n}t^3+b_{n}t^2+c_{n}t+d_{n} \\
p_{n}&=a_{n} \times 0^3+b_{n} \times 0^2+c_{n} \times 0+d_{n} \\
p_{n}&=d_{n}
\end{align}
```
よって、$d_{n}=p_{n}$となります。
</details>

## 係数c
<details><summary>係数$c_{n}$の求め方</summary>

関数$S_{n}(t)$の終点($p_{n+1}$)と$S_{n+1}(t)$の始点($p_{n+1}$)の値は同じです。

$t_{n}=1$、$t_{n+1}=0$として
```math
\begin{align}
S_{n}(t_{n})&=S_{n+1}(t_{n+1}) \\
a_{n}t_{n}^3+b_{n}t_{n}^2+c_{n}t_{n}+d_{n}&=a_{n+1}t_{n+1}^3+b_{n+1}t_{n+1}^2+c_{n+1}t_{n+1}+d_{n+1} \\
a_{n} \times 1+b_{n} \times  1+c_{n} \times 1+d_{n}&=a_{n+1} \times 0+b_{n+1} \times 0+c_{n+1} \times 0+d_{n+1} \\
a_{n}+b_{n}+c_{n}+d_{n}&=d_{n+1} \\
c_{n}&=d_{n+1}-a_{n}-b_{n}-d_{n}
\end{align}
```
よって、$c_{n}=d_{n+1}-a_{n}-b_{n}-d_{n}$となります。
</details>

## 係数a
<details><summary>係数$a_{n}$の求め方</summary>

$S_{n}(t)$の二次微分$\frac{d^2S_{n}(t)}{dt^2}$を考えます。
```math
\begin{align}
S_{n}(t)&=a_{n}t^3+b_{n}t^2+c_{n}t+d_{n} \\
\frac{dS_{n}(t)}{dt}&=3a_{n}t^2+2b_{n}t+c_{n} \\
\frac{d^2S_{n}(t)}{dt^2}&=6a_{n}t+2b_{n}
\end{align}
```

$S_{n}(t)$の始点と$S_{n+1}(t)$の終点でなめらかに曲線が繋がっているためには、一次微分と二次微分が接合点で等しくなっている必要があります。そこで、$t_{n}=1$、$t_{n+1}=0$として、
```math
\begin{align}
\frac{d^2S_{n}(t_{n})}{dt^2}&=\frac{d^2S_{n+1}(t_{n+1})}{dt^2} \\
6a_{n}t_{n}+2b_{n}&=6a_{n+1}t_{n+1}+2b_{n+1} \\
6a_{n} \times 1+2b_{n}&=6a_{n+1} \times 0+2b_{n+1} \\
6a_{n}+2b_{n}&=2b_{n+1} \\
6a_{n}&=2b_{n+1}-2b_{n} \\
a_{n}&=\frac{b_{n+1}-b_{n}}{3}
\end{align}
```
よって、$a_{n}=\frac{b_{n+1}-b_{n}}{3}$となります。
</details>

## 係数b
<details><summary>係数$b_{n}$の求め方</summary>

$S_{n}(t)$の一次微分$\frac{dS_{n}(t)}{dt}$を考えます。
```math
\begin{align}
S_{n}(t)&=a_{n}t^3+b_{n}t^2+c_{n}t+d_{n} \\
\frac{dS_{n}(t)}{dt}&=3a_{n}t^2+2b_{n}t+c_{n}
\end{align}
```

$S_{n}(t)$の始点と$S_{n+1}(t)$の終点でなめらかに曲線が繋がっているためには、一次微分と二次微分が接合点で等しくなっている必要があります。そこで、$t_{n}=1$、$t_{n+1}=0$として、
```math
\begin{align}
\frac{dS_{n}(t_{n})}{dt}&=\frac{dS_{n+1}(t_{n+1})}{dt} \\
3a_{n}t_{n}^2+2b_{n}t_{n}+c_{n}&=3a_{n+1}t_{n+1}^2+2b_{n+1}t_{n+1}+c_{n+1} \\
3a_{n} \times 1+2b_{n}\times 1+c_{n}&=3a_{n+1}\times 0+2b_{n+1}\times 0+c_{n+1} \\
3a_{n}+2b_{n}+c_{n}&=c_{n+1} \\
\end{align}
```
ここで、$c_{n}=d_{n+1}-a_{n}-b_{n}-d_{n}$、$c_{n+1}=d_{n+2}-a_{n+1}-b_{n+1}-d_{n+1}$となるので、
```math
\begin{align}
3a_{n}+2b_{n}+c_{n}&=c_{n+1} \\
3a_{n}+2b_{n}+(d_{n+1}-a_{n}-b_{n}-d_{n})&=(d_{n+2}-a_{n+1}-b_{n+1}-d_{n+1}) \\
2a_{n}+b_{n}+d_{n+1}-d_{n}&=d_{n+2}-a_{n+1}-b_{n+1}-d_{n+1}
\end{align}
```
項を整理します。
```math
2a_{n}+a_{n+1}+b_{n}+b_{n+1}=d_{n}-2d_{n+1}+d_{n+2}
```

更に$a_{n}=\frac{b_{n+1}-b_{n}}{3}$、$a_{n+1}=\frac{b_{n+2}-b_{n+1}}{3}$として、
```math
\begin{align}
2a_{n}+a_{n+1}+b_{n}+b_{n+1}&=d_{n}-2d_{n+1}+d_{n+2} \\
2(\frac{b_{n+1}-b_{n}}{3})+(\frac{b_{n+2}-b_{n+1}}{3})+b_{n}+b_{n+1}&=d_{n}-2d_{n+1}+d_{n+2}
\end{align}
```
分母が邪魔なので、両辺を3倍します。
```math
\begin{align}
(2b_{n+1}-2b_{n})+(b_{n+2}-b_{n+1})+3b_{n}+3b_{n+1}&=3d_{n}-6d_{n+1}+3d_{n+2} \\
b_{n}+4b_{n+1}+b_{n+2}&=3d_{n}-6d_{n+1}+3d_{n+2}
\end{align}
```

さて、$b_{n}+4b_{n+1}+b_{n+2}=3d_{n+2}-6d_{n+1}+3d_{n}$という関係式は、あるものを彷彿とさせます。見えてきませんか？視えてきますよね？笑

$b_{n}$は二次微分なので、加速度と捉えることができます。よって、$p_{0}$と$p_{m}$の加速度を0とします。($b_{0}=b_{m}=0$)

$b_{n}$について以下の連立方程式が成り立ちます。
```math
\begin{align}
b_{0}&=0 \\
b_{0}+4b_{1}+b_{2}&=3d_{0}-6d_{1}+3d_{2} \\
b_{1}+4b_{2}+b_{3}&=3d_{1}-6d_{2}+3d_{3} \\
b_{2}+4b_{3}+b_{4}&=3d_{2}-6d_{3}+3d_{4} \\
\vdots \\
b_{n-1}+4b_{n}+b_{n+1}&=3d_{n-1}-6d_{n}+3d_{n+1} \\
\vdots \\
b_{m-3}+4b_{m-2}+b_{m-1}&=3d_{m-3}-6d_{m-2}+3d_{m-1} \\
b_{m-2}+4b_{m-1}+b_{m}&=3d_{m-2}-6d_{m-1}+3d_{m} \\
b_{m}&=0 \\
\end{align}
```

はい、視えました。
```math
\begin{bmatrix}
1 & 0 & \cdots & \cdots & \cdots & \cdots & \cdots & 0 \\
1 & 4 & 1 & 0 & \cdots & \cdots & \cdots & 0\\
0 & 1 & 4 & 1 & 0 & \cdots & \cdots & 0\\
\vdots & 0 & 1 & 4 & 1 & 0 & \cdots & 0\\
\vdots && \ddots & \ddots & \ddots & \ddots & \ddots \\0 & \cdots & \cdots & 0 & 1 & 4 & 1 & 0\\
0 & \cdots & \cdots & \cdots & 0 & 1 & 4 & 1 \\
0 & \cdots & \cdots & \cdots & \cdots & \cdots & 0 & 1
\end{bmatrix}

\begin{bmatrix}
b_{0} \\ b_{1} \\ b_{2} \\ b_{3} \\
\vdots \\ b_{m-2} \\ b_{m-1} \\ b_{m} 
\end{bmatrix}
=
\begin{bmatrix}
0 \\
3d_{0}-6d_{1}+3d_{2} \\
3d_{1}-6d_{2}+3d_{3} \\
3d_{2}-6d_{3}+3d_{4} \\
\vdots \\
3d_{m-3}-6d_{m-2}+3d_{m-1} \\
3d_{m-2}-6d_{m-1}+3d_{m} \\
0 
\end{bmatrix}
```
係数$b_{n}$はこの行列式を解くことで求めることができます。
</details>

## 係数の式一覧
それぞれの係数を求める式をまとめると、以下のようになります。
```math
\begin{align}
d_{n}&=p_{n} \\
c_{n}&=d_{n+1}-a_{n}-b_{n}-d_{n} \\
b_{n-1}+4b_{n}+b_{n+1}&=3d_{n-1}-6d_{n}+3d_{n+1} \\
a_{n}&=\frac{b_{n+1}-b_{n}}{3}
\end{align}
```


# プログラム例
　luaで申し訳ないですが、プログラム例です。luaの配列は`1`から始まるので注意が必要ですね。
```lua
--------------------------------------------------
--a,b,c,dを求める関数
--------------------------------------------------
local function spline_M(y)

local max_p=#y

local a={}
local b={}
local c={}
local d=y

b[1]=0
b[max_p]=0
for n=2,max_p-1 do
 b[n]=3*(d[n-1]-2*d[n]+d[n+1])
end

--行列の下三角を0にする
local w={0}
for i=2,max_p do
 local tmp=4-w[i-1]
 b[i]=(b[i]-b[i-1])/tmp
 w[i]=1/tmp
end

--行列の上三角を0にして、単位行列に
local i=max_p-1
while i>1 do
 b[i]=b[i]-b[i+1]*w[i]
 i=i-1
end

for n=1,max_p-1 do
 a[n]=(b[n+1]-b[n])/3
 c[n]=d[n+1]-a[n]-b[n]-d[n]
end

local sol={}
sol[1]=a
sol[2]=b
sol[3]=c
sol[4]=d

return sol
end
--------------------------------------------------
--------------------------------------------------
```

# 参考文献
・[3次スプライン補間で軌跡生成：３次多項式と補間条件](https://tajimarobotics.com/cubic-spline-interpolation/)
・[3次スプライン補間で軌跡生成：３次多項式のパラメータを求める（その１）](https://tajimarobotics.com/cubic-spline-interpolation-2/)
・[3次スプライン補間で軌跡生成：３次多項式のパラメータを求める（その２）](https://tajimarobotics.com/cubic-spline-interpolation-3/)
・[3次スプライン補間で軌跡生成：連続的な軌跡を生成する（その１）](https://tajimarobotics.com/cubic-spline-interpolation-result/)
・[3次スプライン補間で軌跡生成：連続的な軌跡を生成する（その２）](https://tajimarobotics.com/cubic-spline-interpolation-result-2/)
・[曲線描画](https://fussy.web.fc2.com/algo/curve1_polynomial.htm)
・[C++を用いてスプライン補間を行う](https://qiita.com/jagaimon/items/dac038be100e65371600)
・[Thomas法による数値解析](https://qiita.com/keito_takaishi/items/18f6943eb857b3360472)
・[第3-2回　TDMA法（Thomas’algorithm）[python]](https://teru-hide.com/cfd-python-tdma/)
・[３次スプライン補間](https://linx.hatenadiary.jp/entry/cubic_spline)


# むすび
　特殊な3次スプライン補間の紹介でした。