3次Spline補間

 


理論

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

FilenameDescription
f90_SPL.txtProgram for Cubic Spline Interpolation
inp_data.txtInput data sample
out_data.txtOutput data sample
gpl_spline.txtScript 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)
...  ...


解析事例

fig_png_spline.png


pic
inserted by FC2 system