理論
3次スプライン補間は,与えられた点列に対し,隣りあう点間を3次曲線で近似し補間する手法である.
今,n個の(x,y)座標が与えられた点を3次スプライン補間することを考える.
補間曲線の方程式
以下の3次式およびその微分を考える.
\begin{align}
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) =&b_i + 2 c_i (x-x_i) + 3 d_i (x-x_i)^2 \\
S''_i(x)=&2 c_i + 6 d_i (x-x_i)
\end{align}
補間式作成の条件
\begin{align}
&\text{条件1:与えられた点の値が一致} & S_i(x_i)=y_i \\
&\text{条件2:隣りあう曲線の値が連続する} & S_{i}(x_{i+1})=S_{i+1}(x_{i+1}) \\
&\text{条件3:隣りあう曲線の傾きが連続する} & S'_{i}(x_{i+1})=S'_{i+1}(x_{i+1}) \\
&\text{条件4:隣りあう曲線の2階微分値が連続する}& S''_{i}(x_{i+1})=S''_{i+1}(x_{i+1}) \\
\end{align}
ただし,条件4においては,端点 $x_1$ および $x_n$ において,$c_1=0$,$c_n=0$ である.
連立方程式
各条件を満足するために必要な関係式は以下のとおり.
\begin{align}
&\text{条件1より} &a_i =&y_i \\
&\text{条件2より} &a_{i+1}=&a_i + b_i h_i + c_i h_i^2 + d_i h_i^3 \qquad ( h_i = x_{i+1} - x_i ) \\
&\text{条件3より} &b_{i+1}=&b_i + 2 c_i h_i + 3 d_i h_i^2 \\
&\text{条件4より} &d_i =&\frac{c_{i+1}-c_i}{3 h_i}
\end{align}
上記関係式から,$b_i$と$d_i$を消去することにより,下式が得られる.
\begin{equation}
h_{i-1} c_{i-1} + 2(h_{i-1} + h_i) c_i + h_i c_{i+1} = \frac{3(a_{i+1}-a_i)}{h_i}-\frac{3(a_i-a_{i-1})}{h_{i-1}} \qquad (i=2,\dots, n-1)
\end{equation}
上式を,$c_1=0$,$c_n=0$の関係を考慮して,マトリックス形式で表現すれば以下のとおりとなる.
\begin{equation}
\begin{bmatrix}
2(h_1+h_2) & h_2 & 0 & 0 & \dots & 0 & 0 & 0 \\
h_2 & 2(h_2+h_3) & h_3 & 0 & \dots & 0 & 0 & 0 \\
0 & h_3 & 2(h_3+h_4) & h_4 & \dots & 0 & 0 & 0 \\
\dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots \\
0 & 0 & 0 & 0 & \dots & h_{n-3} & 2(h_{n-3}+h_{n-2}) & h_{n-2} \\
0 & 0 & 0 & 0 & \dots & 0 & h_{n-2} & 2(h_{n-2}+h_{n-1}) \\
\end{bmatrix}
\begin{Bmatrix}
c_2 \\ c_3 \\ c_4 \\ \dots \\ c_{n-2} \\ c_{n-1}
\end{Bmatrix}
=
\begin{Bmatrix}
u_2 \\ u_3 \\ u_4 \\ \dots \\ u_{n-2} \\ u_{n-1}
\end{Bmatrix}
\end{equation}
\begin{equation}
u_i = \frac{3(a_{i+1}-a_i)}{h_i}-\frac{3(a_i-a_{i-1})}{h_{i-1}} \qquad (i=2,\dots, n-1)
\end{equation}
上記連立方程式は,Cholesky法などにより,効率的に解くことができる.
未知数 $c_i$ が定まれば,$a_i$ は既知であるので,係数 $d_i$ および $b_i$ を求め,曲線を確定することができる.
プログラム
Programs
Filename | Description |
---|---|
f90_SPL.txt | Program for Cubic Spline Interpolation |
inp_data.txt | Input data sample |
out_data.txt | Output data sample |
gpl_spline.txt | Script for gnuplot drawing |
Fortranプログラム実行用コマンド
gfortran -o f90_SPL f90_SPL.f90 ./f90_SPL inp_data.txt 5 > out_data.txt gnuplot gpl_spline.txt
./f90_SPL fnameR num > fnameW
f90_SPL | 実行プログラム |
fnameR | データ入力ファイル名 |
num | 出力時の点間分割数 |
fnameW | 出力ファイル名 |
データ入力書式
x(1) y(1) x(2) y(2) x(3) y(3) ... ...