LoginSignup
0
0
八戸高専Advent Calendar 2023Day 4
@metaphysical_bard

【S01】3次スプライン補間

Posted at

【S01】3次スプライン補間

　
この記事はアドカレに参加しています。

スプライン補間とは

　$p_{0}～p_{m}$というm+1個の点があったとします。その点すべてを通る曲線は、m次の関数で表わすことができます。

2個の点があれば、一次関数です。
3個の点があれば、二次関数です。
11個の点があれば、十次関数です。
101個の点があれば、百次関数です。

この理論でいくと、1001個の点すべてを通る曲線は千次関数とかいうとんでもなものになってしまいます。千次関数なんて、係数を求めるのもグラフを描くのも大変です。

そこでn番目とn+1番目の点の間をn-1番目からn+2番目の4つの点で補間しようというのが、3次スプライン補間の考え方です。

3次スプライン補間は千次関数ほどの正確さはありませんが、少ない演算でそれっぽいグラフを作成できるので多くの分野で使用されています。

以下の図を見てイメージしてみましょう。
a_care08.jpeg
3次スプライン補間はn番目とn+1番目の点の間をn-1番目からn+2番目の4つの点で補間します。

$p_{3}$と$p_{4}$の間を補間する三次関数を求めるには、$p_{2}$～$p_{5}$の4点から三次関数を求めます。
$p_{4}$と$p_{5}$の間を補間する三次関数を求めるには、$p_{3}$～$p_{6}$の4点から三次関数を求めます。

$p_{n}$と$p_{n+1}$の間を補間する三次関数を求めるには、$p_{n-1}$～$p_{n+2}$の4点から三次関数を求めます。

三次関数の式から係数を求める

　三次関数の式から係数を求めてみましょう。

まずは点の数をm+1個として、$p_{0}$～$p_{m}$という点があったとします。このとき、$p_{n}$と$p_{n+1}$を補間する三次関数$S_{n}(x)$は以下のように表すことができます。

\begin{align}
y&=S_{n}(x)　　(x_{n} \leqq x \leqq x_{n+1})\\
y&=a_{n}(x-x_{n})^3+b_{n}(x-x_{n})^2+c_{n}(x-x_{n})+d_{n}
\end{align}

$S_{n}(x)$は、$p_{n}$～$p_{n+1}$を補間する関数なので、$x=x_{n}$のとき$S_{n}(x_{n})=y_{n}$であり、$x=x_{n+1}$のとき$S_{n}(x_{n+1})=y_{n+1}$となります。

変位h

変位$h_{n}$の求め方

この先の導出で頻出するので、変位$h_{n}$というものを定義しておきます。変位$h_{n}$は$x_{n+1}$と$x_{n}$の差を表すものです。

\begin{align}
h_{n}&=x_{n+1}-x_{n} \\
\end{align}

係数d

係数$d_{n}$の求め方

$x=x_{n}$のときを考えます。

\begin{align}
y&=S_{n}(x) \\
y_{n}&=S_{n}(x_{n})
\end{align}

という関係を用いて、

\begin{align}
y_{n}&=a_{n}(x-x_{n})^3+b_{n}(x-x_{n})^2+c_{n}(x-x_{n})+d_{n} \\
y_{n}&=a_{n}(x_{n}-x_{n})^3+b_{n}(x_{n}-x_{n})^2+c_{n}(x_{n}-x_{n})+d_{n} \\
y_{n}&=a_{n}(0)^3+b_{n}(0)^2+c_{n}(0)+d_{n} \\
y_{n}&=d_{n}
\end{align}

よって、$d_{n}=y_{n}$となります。

係数a

係数$a_{n}$の求め方

$S_{n}(x)$の二次微分$\frac{d^2S_{n}(x)}{dx^2}$を考えます。

\begin{align}
S_{n}(x)&=a_{n}(x-x_{n})^3+b_{n}(x-x_{n})^2+c_{n}(x-x_{n})+d_{n} \\
\frac{dS_{n}(x)}{dx}&=3a_{n}(x-x_{n})^2+2b_{n}(x-x_{n})+c_{n} \\
\frac{d^2S_{n}(x)}{dx^2}&=6a_{n}(x-x_{n})+2b_{n}
\end{align}

$S_{n}(x)$の始点と$S_{n+1}(x)$の終点でなめらかに曲線が繋がっているためには、一次微分と二次微分が接合点で等しくなっている必要があります。そこで、

\begin{align}
\frac{d^2S_{n}(x_{n+1})}{dx^2}&=\frac{d^2S_{n+1}(x_{n+1})}{dx^2} \\
6a_{n}(x_{n+1}-x_{n})+2b_{n}&=6a_{n+1}(x_{n+1}-x_{n+1})+2b_{n+1} \\
6a_{n}h_{n}+2b_{n}&=6a_{n+1}(0)+2b_{n+1} \\
6a_{n}h_{n}+2b_{n}&=2b_{n+1} \\
6a_{n}h_{n}&=2b_{n+1}-2b_{n} \\
a_{n}&= \frac{b_{n+1}-b_{n}}{3h_{n}}
\end{align}

よって、$a_{n}=\frac{b_{n+1}-b_{n}}{3h_{n}}$となります。

係数c

係数$c_{n}$の求め方

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

重なっている点の位置は$x_{n+1}$なので、

\begin{align}
S_{n}(x_{n+1})&=S_{n+1}(x_{n+1}) \\
a_{n}(x_{n+1}-x_{n})^3+b_{n}(x_{n+1}-x_{n})^2+c_{n}(x_{n+1}-x_{n})+d_{n}&=a_{n+1}(x_{n+1}-x_{n+1})^3+b_{n+1}(x_{n+1}-x_{n+1})^2+c_{n+1}(x_{n+1}-x_{n+1})+d_{n+1} \\
a_{n}h_{n}^3+b_{n}h_{n}^2+c_{n}h_{n}+d_{n}&=a_{n+1}(0)^3+b_{n+1}(0)^2+c_{n+1}(0)+d_{n+1} \\
a_{n}h_{n}^3+b_{n}h_{n}^2+c_{n}h_{n}+d_{n}&=d_{n+1} \\
c_{n}h_{n}&=d_{n+1}-d_{n}-a_{n}h_{n}^3-b_{n}h_{n}^2 \\
c_{n}&=\frac{d_{n+1}-d_{n}}{h_{n}}-a_{n}h_{n}^2-b_{n}h_{n} \\
\end{align}

ここで、$a_{n}=\frac{b_{n+1}-b_{n}}{3h_{n}}$として、

\begin{align}
c_{n}&=\frac{d_{n+1}-d_{n}}{h_{n}}-a_{n}h_{n}^2-b_{n}h_{n} \\
c_{n}&=\frac{d_{n+1}-d_{n}}{h_{n}}-(\frac{b_{n+1}-b_{n}}{3h_{n}})h_{n}^2-b_{n}h_{n} \\
c_{n}&=\frac{d_{n+1}-d_{n}}{h_{n}}-(\frac{b_{n+1}-b_{n}}{3})h_{n}-b_{n}h_{n} \\
c_{n}&=\frac{d_{n+1}-d_{n}}{h_{n}}-(\frac{b_{n+1}-b_{n}}{3}+b_{n})h_{n} \\
c_{n}&=\frac{d_{n+1}-d_{n}}{h_{n}}-(\frac{b_{n+1}-b_{n}}{3}+\frac{3b_{n}}{3})h_{n} \\
c_{n}&=\frac{d_{n+1}-d_{n}}{h_{n}}-\frac{b_{n+1}+2b_{n}}{3}h_{n} \\
\end{align}

よって、$c_{n}=\frac{d_{n+1}-d_{n}}{h_{n}}-\frac{b_{n+1}+2b_{n}}{3}h_{n}$となります。

係数b

係数$b_{n}$の求め方

$S_{n}(x)$の一次微分$\frac{dS_{n}(x)}{dx}$を考えます。

\begin{align}
S_{n}(x)&=a_{n}(x-x_{n})^3+b_{n}(x-x_{n})^2+c_{n}(x-x_{n})+d_{n} \\
\frac{dS_{n}(x)}{dx}&=3a_{n}(x-x_{n})^2+2b_{n}(x-x_{n})+c_{n}
\end{align}

$S_{n}(x)$の始点と$S_{n+1}(x)$の終点でなめらかに曲線が繋がっているためには、一次微分と二次微分が接合点で等しくなっている必要があります。そこで、

\begin{align}
\frac{dS_{n}(x_{n+1})}{dx}&=\frac{dS_{n+1}(x_{n+1})}{dx} \\
3a_{n}(x_{n+1}-x_{n})^2+2b_{n}(x_{n+1}-x_{n})+c_{n}&=3a_{n+1}(x_{n+1}-x_{n+1})^2+2b_{n+1}(x_{n+1}-x_{n+1})+c_{n+1} \\
3a_{n}h_{n}^2+2b_{n}h_{n}+c_{n}&=3a_{n+1}(0)^2+2b_{n+1}(0)+c_{n+1} \\
3a_{n}h_{n}^2+2b_{n}h_{n}+c_{n}&=c_{n+1} \\
3a_{n}h_{n}^2+2b_{n}h_{n}&=-c_{n}+c_{n+1}
\end{align}

ここで、$a_{n}=\frac{b_{n+1}-b_{n}}{3h_{n}}$、$c_{n}=\frac{d_{n+1}-d_{n}}{h_{n}}-\frac{b_{n+1}+2b_{n}}{3}h_{n}$、$c_{n+1}=\frac{d_{n+2}-d_{n+1}}{h_{n+1}}-\frac{b_{n+2}+2b_{n+1}}{3}h_{n+1}$となるので、

\begin{align}
3a_{n}h_{n}^2+2b_{n}h_{n}&=-c_{n}+c_{n+1} \\
3\Big(\frac{b_{n+1}-b_{n}}{3h_{n}}\Big)
h_{n}^2+2b_{n}h_{n}&=
-\Big( \frac{d_{n+1}-d_{n}}{h_{n}}-\frac{b_{n+1}+2b_{n}}{3}h_{n}\Big)+\Big(\frac{d_{n+2}-d_{n+1}}{h_{n+1}}-\frac{b_{n+2}+2b_{n+1}}{3}h_{n+1}\Big)
\end{align}

分母が邪魔なので両辺を3倍します。

\begin{align}
3\Big(\frac{b_{n+1}-b_{n}}{h_{n}}\Big)
h_{n}^2+6b_{n}h_{n}&=
-\Big( 3\frac{d_{n+1}-d_{n}}{h_{n}}-(b_{n+1}+2b_{n})h_{n}\Big)+\Big(3\frac{d_{n+2}-d_{n+1}}{h_{n+1}}-(b_{n+2}+2b_{n+1})h_{n+1}\Big) \\
3(b_{n+1}-b_{n})h_{n}+6b_{n}h_{n}&=
-3\frac{d_{n+1}-d_{n}}{h_{n}}+3(b_{n+1}+2b_{n})h_{n}+3\frac{d_{n+2}-d_{n+1}}{h_{n+1}}-3(b_{n+2}+2b_{n+1})h_{n+1}
\end{align}

$b$を左辺に、$d$を右辺に持っていき、項を整理します。

\begin{align}
3(b_{n+1}-b_{n})h_{n}+6b_{n}h_{n}
-3(b_{n+1}+2b_{n})h_{n}+3(b_{n+2}+2b_{n+1})h_{n+1}&=
-3\frac{d_{n+1}-d_{n}}{h_{n}}+3\frac{d_{n+2}-d_{n+1}}{h_{n+1}} \\
(b_{n}+2b_{n+1})(h_{n})+(b_{n+2}+2b_{n+1})h_{n+1}&=
3\frac{d_{n+2}-d_{n+1}}{h_{n+1}}-3\frac{d_{n+1}-d_{n}}{h_{n}} \\
b_{n}h_{n}+2b_{n+1}(h_{n}+h_{n+1})+b_{n+2}h_{n+1}&=
3 \big(\frac{d_{n+2}-d_{n+1}}{h_{n+1}}-\frac{d_{n+1}-d_{n}}{h_{n}} \big)
\end{align}

$b_{n}$についての関係式$b_{n}h_{n}+2b_{n+1}(h_{n}+h_{n+1})+b_{n+2}h_{n+1}=3 \big(\frac{d_{n+2}-d_{n+1}}{h_{n+1}}-\frac{d_{n+1}-d_{n}}{h_{n}} \big)$が成り立ちました。

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

$b_{n}$について以下の連立方程式が成り立ちます。

\begin{align}
b_{0}&=0 \\
b_{0}h_{0}+2b_{1}(h_{0}+h_{1})+b_{2}h_{1}&=
3 \big(\frac{d_{2}-d_{1}}{h_{1}}-\frac{d_{1}-d_{0}}{h_{0}} \big) \\
b_{1}h_{1}+2b_{2}(h_{1}+h_{2})+b_{3}h_{2}&=
3 \big(\frac{d_{3}-d_{2}}{h_{2}}-\frac{d_{2}-d_{1}}{h_{1}} \big) \\
b_{2}h_{2}+2b_{3}(h_{2}+h_{3})+b_{4}h_{3}&=
3 \big(\frac{d_{4}-d_{3}}{h_{3}}-\frac{d_{3}-d_{2}}{h_{2}} \big) \\
\vdots \\
b_{n-1}h_{n-1}+2b_{n}(h_{n-1}+h_{n})+b_{n+1}h_{n}&=
3 \big(\frac{d_{n+1}-d_{n}}{h_{n}}-\frac{d_{n}-d_{n-1}}{h_{n-1}} \big) \\
\vdots \\
b_{m-3}h_{m-3}+2b_{m-2}(h_{m-3}+h_{m-2})+b_{m-1}h_{m-2}&=
3 \big(\frac{d_{m-1}-d_{m-2}}{h_{m-2}}-\frac{d_{m-2}-d_{m-3}}{h_{m-3}} \big) \\
b_{m-2}h_{m-2}+2b_{m-1}(h_{m-2}+h_{m-1})+b_{m}h_{m-1}&=
3 \big(\frac{d_{m}-d_{m-1}}{h_{m-1}}-\frac{d_{m-1}-d_{m-2}}{h_{m-2}} \big) \\
b_{m}&=0 \\
\end{align}

行列にすると以下のようになります。

\begin{bmatrix}
1 & 0 & \cdots & \cdots & \cdots & \cdots & \cdots & 0 \\
h_{0} & 2(h_{0}+h_{1}) & h_{1} & 0 & \cdots & \cdots & \cdots & 0\\
0 & h_{1} & 2(h_{1}+h_{2}) & h_{2} & 0 & \cdots & \cdots & 0\\
\vdots & 0 & h_{2} & 2(h_{2}+h_{3}) & h_{3} & 0 & \cdots & 0\\
\vdots && \ddots & \ddots & \ddots & \ddots & \ddots \\
0 & \cdots & \cdots & 0 & h_{m-3} & 2(h_{m-3}+h_{m-2}) & h_{m-2} & 0\\
0 & \cdots & \cdots & \cdots & 0 & h_{m-2} & 2(h_{m-2}+h_{m-1}) & h_{m-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 \\
3 \big(\frac{d_{2}-d_{1}}{h_{1}}-\frac{d_{1}-d_{0}}{h_{0}} \big) \\
3 \big(\frac{d_{3}-d_{2}}{h_{2}}-\frac{d_{2}-d_{1}}{h_{1}} \big) \\
3 \big(\frac{d_{4}-d_{3}}{h_{3}}-\frac{d_{3}-d_{2}}{h_{2}} \big) \\
\vdots \\
3 \big(\frac{d_{m-1}-d_{m-2}}{h_{m-2}}-\frac{d_{m-2}-d_{m-3}}{h_{m-3}} \big) \\
3 \big(\frac{d_{m}-d_{m-1}}{h_{m-1}}-\frac{d_{m-1}-d_{m-2}}{h_{m-2}} \big) \\
0 
\end{bmatrix}

係数$b_{n}$はこの行列式を解くことで求めることができます。

係数の式一覧

それぞれの係数を求める式をまとめると、以下のようになります。

\begin{align}
h_{n}&=x_{n+1}-x_{n} \\
d_{n}&=y_{n} \\
c_{n}&=\frac{d_{n+1}-d_{n}}{h_{n}}-\frac{b_{n+1}+2b_{n}}{3}h_{n} \\
b_{n-1}h_{n-1}+2b_{n}(h_{n-1}+h_{n})+b_{n+1}h_{n}&=
3 \big(\frac{d_{n+1}-d_{n}}{h_{n}}-\frac{d_{n}-d_{n-1}}{h_{n-1}} \big) \\
a_{n}&= \frac{b_{n+1}-b_{n}}{3h_{n}}
\end{align}

三重対角行列を解くアルゴリズム

　係数$b_{n}$は行列を解くことで求めることができます。普通にガウスの消去法などを用いて行列を解いてもいいのですが、$b_{n}$の行列の形は三重対角行列という特殊な形なのでThomas法を使うことで効率的に行列を解くことができます。

Thomas法では以下のような行列を解くことができます。

\begin{bmatrix}
b_{0} & c_{0} & \cdots & \cdots & \cdots & \cdots & \cdots & 0 \\
a_{1} & b_{1} & c_{1} & 0 & \cdots & \cdots & \cdots & 0\\
0 & a_{2} & b_{2} & c_{2} & 0 & \cdots & \cdots & 0\\
\vdots & 0 & a_{3} & b_{3} & c_{3} & 0 & \cdots & 0\\
\vdots && \ddots & \ddots & \ddots & \ddots & \ddots \\
0 & \cdots & \cdots & 0 & a_{m-2} & b_{m-2} & c_{m-2} & 0\\
0 & \cdots & \cdots & \cdots & 0 & a_{m-1} & b_{m-1} & c_{m-1} \\
0 & \cdots & \cdots & \cdots & \cdots & \cdots & a_{m} & b_{m}
\end{bmatrix}

\begin{bmatrix}
x_{0} \\ x_{1} \\ x_{2} \\ x_{3} \\
\vdots \\ x_{m-2} \\ x_{m-1} \\ x_{m} 
\end{bmatrix}
=
\begin{bmatrix}
d_{0} \\ d_{1} \\ d_{2} \\ d_{3} \\
\vdots \\ d_{m-2} \\ d_{m-1} \\ d_{m}
\end{bmatrix}

$0～m$のすべての自然数$n$で、$a_{n}=c_{n}=0$かつ$b_{n}=1$であれば$x_{n}=d_{n}$となるので、行列を解いたことになります。

Thomas法での行列の解き方

まずは$a_{n}=0$、$b_{n}=1$となることを目指します。

一行目は$b_{0}$で割るだけでいいですね。演算後の値には$'$を付けています。

\begin{align}
c_{0}'&=c_{0}/b_{0} \\
d_{0}'&=d_{0}/b_{0}
\end{align}

一行目以降の$n$行目に対しては、上の行から順番に$a_{n}=0$となるようにします。さらに、$b_{n}=1$になるようにします。(上の行に$a_{n}$を掛けて引いた後に、$b_{n}$で割るという操作になる。)

\begin{align}
a_{n}'&=\frac{a_{n}-b_{n-1}'\times a_{n}}{b_{n}-c_{n-1}'\times a_{n}}
=\frac{a_{n}-1\times a_{n}}{b_{n}-c_{n-1}'\times a_{n}}=\frac{0}{b_{n}-c_{n-1}'\times a_{n}}=0\\
b_{n}'&=\frac{b_{n}-c_{n-1}'\times a_{n}}{b_{n}-c_{n-1}'\times a_{n}}=1 \\
c_{n}'&=\frac{c_{n}-0 \times a_{n}}{b_{n}-c_{n-1}'\times a_{n}}
=\frac{c_{n}}{b_{n}-c_{n-1}'\times a_{n}}\\
d_{n}'&=\frac{d_{n}-d_{n-1}'\times a_{n}}{b_{n}-c_{n-1}'\times a_{n}}
\end{align}

これで上三角行列となりました。

\begin{align}
\begin{bmatrix}
1 & c_{0}' & \cdots & \cdots & \cdots & \cdots & \cdots & 0 \\
0 & 1 & c_{1}' & 0 & \cdots & \cdots & \cdots & 0\\
0 & 0 & 1 & c_{2}' & 0 & \cdots & \cdots & 0\\
\vdots & 0 & 0 & 1 & c_{3}' & 0 & \cdots & 0\\
\vdots && \ddots & \ddots & \ddots & \ddots & \ddots \\
0 & \cdots & \cdots & 0 & 0 & 1 & c_{m-2}' & 0\\
0 & \cdots & \cdots & \cdots & 0 & 0 & 1 & c_{m-1}' \\
0 & \cdots & \cdots & \cdots & \cdots & \cdots & 0 & 1
\end{bmatrix}
\begin{bmatrix}x_{0} \\ x_{1} \\ x_{2} \\ x_{3} \\ \vdots \\ x_{m-2} \\ x_{m-1} \\ x_{m} \end{bmatrix} &= 
\begin{bmatrix}d_{0}' \\ d_{1}' \\ d_{2}' \\ d_{3}' \\ \vdots \\ d_{m-2}' \\ d_{m-1}' \\ d_{m}' \end{bmatrix}
\end{align}

今度は単位行列を目指すために$c_{n}=0$となるようにしていきます。最後の行ではすでに$x_{m}=d_{m}''=d_{m}'$となっているので、下の行から処理していきます。(下の行に$c_{n}'$を掛けて引く操作。)

\begin{align}
c_{n}''&=c_{n}'-b_{n+1}'\times c_{n}'=c_{n}'-1\times c_{n}'=0 \\
d_{n}''&=d_{n}'-d_{n+1}''\times c_{n}'
\end{align}

これで単位行列となったので、$0～m$の自然数$n$で$x_{n}=d_{n}''$となりました。

\begin{align}
\begin{bmatrix}
1 & 0 & \cdots & \cdots & \cdots & \cdots & \cdots & 0 \\
0 & 1 & 0 & 0 & \cdots & \cdots & \cdots & 0\\
0 & 0 & 1 & 0 & 0 & \cdots & \cdots & 0\\
\vdots & 0 & 0 & 1 & 0 & 0 & \cdots & 0\\
\vdots && \ddots & \ddots & \ddots & \ddots & \ddots \\
0 & \cdots & \cdots & 0 & 0 & 1 & 0 & 0\\
0 & \cdots & \cdots & \cdots & 0 & 0 & 1 & 0 \\
0 & \cdots & \cdots & \cdots & \cdots & \cdots & 0 & 1
\end{bmatrix}
\begin{bmatrix}x_{0} \\ x_{1} \\ x_{2} \\ x_{3} \\ \vdots \\ x_{m-2} \\ x_{m-1} \\ x_{m} \end{bmatrix} &= 
\begin{bmatrix}d_{0}'' \\ d_{1}'' \\ d_{2}'' \\ d_{3}'' \\ \vdots \\ d_{m-2}'' \\ d_{m-1}'' \\ d_{m}'' \end{bmatrix}
\end{align}
Thomas法まとめ

要点をまとめます。

まず、1行目の処理をします。

\begin{align}
c_{0}'&=c_{0}/b_{0} \\
d_{0}'&=d_{0}/b_{0}
\end{align}

次に、上から順に$n$行目$(0<n\leqq m)$に対して以下の処理をします。

\begin{align}
c_{n}'&=\frac{c_{n}}{b_{n}-c_{n-1}'\times a_{n}} \\
d_{n}'&=\frac{d_{n}-d_{n-1}'\times a_{n}}{b_{n}-c_{n-1}'\times a_{n}}
\end{align}

$m$行目の値は既に分かりました。

x_{m}=d_{m}''=d_{m}'

下から順に$n$行目$(m>n\geqq 1)$に対して以下の処理をします。

x_{n}=d_{n}''=d_{n}'-d_{n+1}''\times c_{n}'
Thomas法での行列の解き方(省メモリver.)

係数の分の配列を使用できるのであれば、こちらの方法の方が少ないメモリで処理することができます。

先ほどの方法との違いはいつ$b_{n}=1$とするかです。
先程の方法では、$a_{n}=0$、$b_{n}=1$、$c_{n}=0$の順で解いていきましたが、こちらの方法では$a_{n}=0$、$c_{n}=0$、$b_{n}=1$の順で解いていきます。

まずは一行目以外の$n$行目に対して、上の行から順番に$a_{n}=0$になるようにします。$w_{n}=\frac{a_{n}}{b_{n-1}'}$として、

\begin{align}
a_{n}'&=a_{n}-b_{n-1}'\times w_{n}=a_{n}-b_{n-1}' \times \frac{a_{n}}{b_{n-1}'}=0 \\
b_{n}'&=b_{n}-c_{n-1}\times w_{n} \\
c_{n}'&=c_{n}-0\times w_{n}=c_{n} \\
d_{n}'&=d_{n}-d_{n-1}'\times w_{n}
\end{align}

$a_{n}=0$となったので、上三角行列となりました。

\begin{align}
\begin{bmatrix}
b_{0} & c_{0} & \cdots & \cdots & \cdots & \cdots & \cdots & 0 \\
0 & b_{1}' & c_{1} & 0 & \cdots & \cdots & \cdots & 0\\
0 & 0 & b_{2}' & c_{2} & 0 & \cdots & \cdots & 0\\
\vdots & 0 & 0 & b_{3}' & c_{3} & 0 & \cdots & 0\\
\vdots && \ddots & \ddots & \ddots & \ddots & \ddots \\
0 & \cdots & \cdots & 0 & 0 & b_{m-2}' & c_{m-2} & 0\\
0 & \cdots & \cdots & \cdots & 0 & 0 & b_{m-1}' & c_{m-1} \\
0 & \cdots & \cdots & \cdots & \cdots & \cdots & 0 & b_{m}'
\end{bmatrix}
\begin{bmatrix}x_{0} \\ x_{1} \\ x_{2} \\ x_{3} \\ \vdots \\ x_{m-2} \\ x_{m-1} \\ x_{m} \end{bmatrix} &= 
\begin{bmatrix}d_{0} \\ d_{1}' \\ d_{2}' \\ d_{3}' \\ \vdots \\ d_{m-2}' \\ d_{m-1}' \\ d_{m}' \end{bmatrix}
\end{align}

このとき、最終行では$b_{m}'$で割ってあげると$x_{m}$が求まります。

x_{m}=d_{m}''=d_{m}'/b_{m}'

$x_{m}$が求まったので、下から順に$c_{n}=0$となるようにしていきます。

\begin{align}
b_{n}''&=b_{n}'-0\times \frac{c_{n}'}{b_{n+1}'}=b_{n}' \\
c_{n}''&=c_{n}-b_{n+1}''\times \frac{c_{n}}{b_{n+1}''}=0 \\
d_{n}''&=d_{n}'-d_{n+1}''\times \frac{c_{n}'}{b_{n+1}'}
\end{align}

同時に$b_{n}=1$となるように割ってあげます。

\begin{align}
b_{n}''&=b_{n}'/b_{n}'=1 \\
d_{n}''&=\frac{d_{n}'-d_{n+1}''\times \frac{c_{n}}{b_{n+1}'}}{b_{n}'}
\end{align}

一行目も同様に$c_{n}=0$、$b_{n}=1$となるようにします。

\begin{align}
b_{0}'&=b_{0}/b_{0}=1 \\
d_{0}'&=\frac{d_{0}-d_{1}''\times \frac{c_{0}}{b_{1}'}}{b_{0}}
\end{align}

これで単位行列となったので、$1～m$の自然数$n$で$x_{n}=d_{n}''$、一行目で$x_{0}=d_{n}'$となりました。

\begin{align}
\begin{bmatrix}
1 & 0 & \cdots & \cdots & \cdots & \cdots & \cdots & 0 \\
0 & 1 & 0 & 0 & \cdots & \cdots & \cdots & 0\\
0 & 0 & 1 & 0 & 0 & \cdots & \cdots & 0\\
\vdots & 0 & 0 & 1 & 0 & 0 & \cdots & 0\\
\vdots && \ddots & \ddots & \ddots & \ddots & \ddots \\
0 & \cdots & \cdots & 0 & 0 & 1 & 0 & 0\\
0 & \cdots & \cdots & \cdots & 0 & 0 & 1 & 0 \\
0 & \cdots & \cdots & \cdots & \cdots & \cdots & 0 & 1
\end{bmatrix}
\begin{bmatrix}x_{0} \\ x_{1} \\ x_{2} \\ x_{3} \\ \vdots \\ x_{m-2} \\ x_{m-1} \\ x_{m} \end{bmatrix} &= 
\begin{bmatrix}d_{0}' \\ d_{1}'' \\ d_{2}'' \\ d_{3}'' \\ \vdots \\ d_{m-2}'' \\ d_{m-1}'' \\ d_{m}'' \end{bmatrix}
\end{align}
Thomas法まとめ(省メモリver.)

要点をまとめます。

まず、上から順に$n$行目$(0<n\leqq m)$に対して以下の処理をします。

\begin{align}
w_{n}&=\frac{a_{n}}{b_{n-1}'} \\
b_{n}'&=b_{n}-c_{n-1}\times w_{n} \\
d_{n}'&=d_{n}-d_{n-1}'\times w_{n}
\end{align}

最終行だけ$b_{m}'$で割ることで、$x_{m}$が求まります。

x_{m}=d_{m}''=d_{m}'/b_{m}'

下から順に$n$行目$(m>n\geqq 1)$に対して以下の処理をします。

x_{n}=d_{n}''=\frac{d_{n}'-d_{n+1}''\times c_{n}}{b_{n}'}

プログラム例

　プログラム例です。

/*

SplineCurve_M.hpp

*/

#pragma once
#include <vector>

class SplineCurve_M {

	int num;
	std::vector<double> a;
	std::vector<double> b;
	std::vector<double> c;
	std::vector<double> d;
	std::vector<double> px;

	void TDMA(int n, std::vector<double>& a, std::vector<double>& b, std::vector<double>& c, std::vector<double>& d) {
		for (int i = 1; i < n; i++) {
			double w = a[i] / b[i - 1];
			b[i] = b[i] - c[i - 1] * w;
			d[i] = d[i] - d[i - 1] * w;
		}

		d[n - 1] = d[n - 1] / b[n - 1];
		for (int i = n - 2; i >= 0; i--)
			d[i] = (d[i] - d[i + 1] * c[i]) / b[i];
	}

public:

	void set(int n, std::vector<double>& x, std::vector<double>& y) {
		num = n;
		a.resize(n); b.resize(n);
		c.resize(n); d.resize(n);
		px.resize(n);
		std::vector<double> h(n - 1);

		//変位h
		for (int i = 0; i < n - 1; i++) {
			h[i] = x[i + 1] - x[i];
		}

		//係数d、xのコピー
		for (int i = 0; i < n; i++) {
			d[i] = y[i];
			px[i] = x[i];
		}

		//係数b
		{
			std::vector<double> w(n);

			b[0] = 0; b[n - 1] = 0;//d
			a[0] = 0; a[n - 1] = 0;//a
			w[0] = 1; w[n - 1] = 1;//b
			c[0] = 0; c[n - 1] = 0;//c
			for (int i = 1; i < n - 1; i++) {
				b[i] = 3 * ((d[i + 1] - d[i]) / h[i] - (d[i] - d[i - 1]) / h[i - 1]);
				a[i] = h[i - 1];
				w[i] = 2 * (h[i - 1] + h[i]);
				c[i] = h[i];
			}

			TDMA(n, a, w, c, b);
		}

		//係数a、係数c
		for (int i = 0; i < n - 1; i++) {
			a[i] = (b[i + 1] - b[i]) / (3 * h[i]);
			c[i] = (d[i + 1] - d[i]) / h[i] - h[i] * (b[i + 1] + 2 * b[i]) / 3;
		}
	}

	double run(double x) {
		int i;
		for (i = 0; i < num - 1; i++)
			if (px[i] <= x && x <= px[i + 1]) break;
		double dx = x - px[i];
		return ((a[i] * dx + b[i]) * dx + c[i]) * dx + d[i];
	}

};

参考文献

・3次スプライン補間で軌跡生成：３次多項式と補間条件
・3次スプライン補間で軌跡生成：３次多項式のパラメータを求める（その１）
・3次スプライン補間で軌跡生成：３次多項式のパラメータを求める（その２）
・3次スプライン補間で軌跡生成：連続的な軌跡を生成する（その１）
・3次スプライン補間で軌跡生成：連続的な軌跡を生成する（その２）
・曲線描画
・C++を用いてスプライン補間を行う
・３次スプライン補間
・Thomas法による数値解析
・第3-2回　TDMA法（Thomas’algorithm）[python]
・Tridiagonal matrix algorithm

むすび

　コードとか参考文献は5月に準備していたのですが、その間に二つの参考文献のリンクが死んでいました。悲しい。ですが、無償で有料級の情報を提供してくださる方のおかげで、僕はプログラムを学べているので文句はいいません。智の開拓をしてくださった先人様に常に感謝です。

ちなみに、知識を共有してくれている方がここにもいますね。
僕に感謝してくれてもいいんですよ…？

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0