３次スプラインのバリエーションをまとめました。

準備
- ３次スプラインの定義
- 解くべき式の導出
３次スプラインのバリエーションと計算法
各バリエーションの比較テスト
- 自然スプライン、固定スプライン、ノットなしスプライン
- 周期的スプライン
まとめ

準備

はじめに、３次スプラインの基本を整理します。

３次スプラインの定義

データ点 $(x_0,y_0)...(x_N,y_N)$ を、次の区分多項式 $f_j(x)$ で補間します。 $(j=1...N)$

$\displaystyle f_j(x)=a_j(x-x_j)^3+b_j(x-x_j)^2+c_j(x-x_j)+d_j, \;\;\;\;\; (x_{j-1} \le x \le x_j)$

ただし、 $f_j(x)$ は以下の条件を満たします。

条件1. 区間の右側の端点でデータ点に一致。 $(j=1...N)$

$\displaystyle f_j(x_j)=y_j$

条件2. 区間の左側の端点でデータ点に一致。 $(j=1...N)$

$\displaystyle f_j(x_{j-1})=y_{j-1}$

条件3. 隣接する区間の境界点で1階微分が一致。 $(j=1...N-1)$

$\displaystyle f'_j(x_j)=f'_{j+1}(x_j)$

条件4. 隣接する区間の境界点で2階微分が一致。 $(j=1...N-1)$

$\displaystyle f''_j(x_j)=f''_{j+1}(x_j)$

図１にイメージを示します。

解くべき式の導出

上の4つの条件式を係数の組 $(a_j, b_j, c_j, d_j)$ について解けば良いですが、これらを直接 $(x_j, y_j)$ の式で表すことは困難です。そこで先ず $y''_j=f''_j(x_j)$ を使って、 $(a_j, b_j, c_j, d_j)$ を $(x_j, y_j, y''_j)$ の式で表します。その次に、 $y''_j$ が満たすべき式を導きます。

[ $d_j$ の表式 ]

条件1.より

$\displaystyle f_j(x_j)=d_j=y_j$

よって

$\displaystyle d_j=y_j \tag{1}$

[ $b_j$ の表式 ]

$f_j(x)$ の定義式より

$\displaystyle \left.f''_j(x_j)=6a_j(x-x_j)+2b_j\,\right|_{x=x_j}=2b_j=y''_j$

よって

$\displaystyle b_j=\frac{y''_j}{2} \tag{2}$

[ $a_j$ の表式 ]

条件4.より

$\displaystyle \left.f''_{j+1}(x_j)=6a_{j+1}(x-x_{j+1})+2b_{j+1}\,\right|_{x=x_j}=6a_{j+1}(-\Delta x_{j+1})+y''_{j+1}=y''_j$

よって

$\displaystyle a_j=\frac{1}{6}\frac{y''_j-y''_{j-1}}{\Delta x_j} \tag{3}$

ただし、 $\Delta$ は１つ前の項との差分を表します（ $\Delta x_j=x_j-x_{j-1}$ ）。

[ $c_j$ の表式 ]

条件2.より

$\displaystyle f_j(x_{j-1})=a_j(x_{j-1}-x_j)^3+b_j(x_{j-1}-x_j)^2+c_j(x_{j-1}-x_j)+d_j=y_{j-1}$

すでに求めた係数 $(a_j, b_j, d_j)$ の表式を代入すると

$\displaystyle \frac{1}{6}\frac{y''_j-y''_{j-1}}{\Delta x_j}(-\Delta x_j)^3+\frac{y''_j}{2}(-\Delta x_j)^2+c_j(-\Delta x_j)+y_j=y_{j-1}$

よって

$\displaystyle c_j=s_{j-1}+\frac{1}{6}\Delta x_j(2y''_j+y''_{j-1}) \tag{4}$

ただし、

$\displaystyle s_j=\frac{\Delta y_{j+1}}{\Delta x_{j+1}}$

としました。添字は、後の式の見栄えを良くするため、あえて１つずらしました。

[ $y''_j$ の表式 ]

続けて $y''_j$ が満たすべき式を導きます。条件3.より

$\displaystyle f'_j(x_j)=\left.3a_j(x-x_j)^2+2b_j(x-x_j)+c_j\right|_{x=x_j}\\ =\left.3a_{j+1}(x-x_{j+1})^2+2b_{j+1}(x-x_{j+1})+c_{j+1}\right|_{x=x_j}=f'_{j+1}(x_j)$

式を整理して

$\displaystyle 3a_{j+1}(\Delta x_{j+1})^2-2b_{j+1}\Delta x_{j+1}+\Delta c_{j+1}=0$

係数 $(a_j,b_j,c_j)$ の表式を代入すると

$\displaystyle \frac{1}{2}\frac{y''_{j+1}-y''_j}{\Delta x_{j+1}}(\Delta x_{j+1})^2-y''_{j+1}\Delta x_{j+1}+\Delta \left(s_j+\frac{1}{6}\Delta x_{j+1}(2y''_{j+1}+y''_j)\right)=0$

最終的に、 $y''_j$ の隣接３項間漸化式の形に整理できます。 $(j=1...N-1)$

$\displaystyle \Delta x_j y''_{j-1}+2(\Delta x_j+\Delta x_{j+1})y''_j+\Delta x_{j+1} y''_{j+1}= 6\Delta s_j$

行列の形で書けば

$\displaystyle \begin{pmatrix} \bigcirc & \cdots & \cdots & \cdots & \cdots & \bigcirc \\ \hline \Delta x_1 & 2(\Delta x_1+\Delta x_2) & \Delta x_2 & & \huge{0} & \\ &\Delta x_2 & \ddots & \ddots & & \\ & & \ddots & \ddots & \Delta x_{N-1} & \\ & \huge{0} & & \Delta x_{N-1} & 2(\Delta x_{N-1}+\Delta x_N) & \Delta x_{N} \\ \hline \bigcirc & \cdots & \cdots & \cdots & \cdots & \bigcirc \\ \end{pmatrix} \begin{pmatrix} y''_0 \\ \vdots \\ y''_N \\ \end{pmatrix} = \begin{pmatrix} \bigcirc \\ \hline 6\Delta s_1 \\ \vdots \\ 6\Delta s_{N-1} \\ \hline \bigcirc \\ \end{pmatrix} \tag{5}$

です（横線はただの仕切り線です）。

式 (5) に示した上下２つの空白 ( $\bigcirc$ ) の行に適当な $y''_j$ の条件式を入れると、すべての $y''_j$ が定まり、式 (1)-(4) から係数の組 $(a_j,b_j,c_j,d_j)$ が求まります。このとき加える $y''_j$ の条件式により、３次スプラインには複数のバリエーションが存在します。

３次スプラインのバリエーションと計算法

３次スプラインのバリエーションを挙げます。なお、周期的スプライン以外は、左端と右端とで異なるバリエーションを組み合わせることも可能です。

自然スプライン (natural spline)

最端点 $x=x_0, x_{N}$ で２階微分がゼロとなるように $y''_j$ を決めます。自然スプライン曲線は、区間 $[x_0, x_N$ ] の外で直線（=２階微分がゼロ）に滑らかにつながります。

[ $y''_j$ の条件式]

$y''_0=y''_N=0$

[ 解くべき式 ]

$\displaystyle \begin{pmatrix} 1 & 0 & \cdots & \cdots & \cdots & 0 \\ \hline \Delta x_1 & 2(\Delta x_1+\Delta x_2) & \Delta x_2 & & \huge{0} & \\ &\Delta x_2 & \ddots & \ddots & & \\ & & \ddots & \ddots & \Delta x_{N-1} & \\ & \huge{0} & & \Delta x_{N-1} & 2(\Delta x_{N-1}+\Delta x_N) & \Delta x_{N} \\ \hline 0 & \cdots & \cdots & \cdots & 0 & 1 \\ \end{pmatrix} \begin{pmatrix} y''_0 \\ \vdots \\ y''_N \\ \end{pmatrix} = \begin{pmatrix} 0 \\ \hline 6\Delta s_1 \\ \vdots \\ 6\Delta s_{N-1} \\ \hline 0 \\ \end{pmatrix}$

このような、係数行列が３重対角行列で与えられる線形連立方程式は、３重対角行列アルゴリズム（Tri-Diagonal Matrix Algorithm, TDMA）により効率的に解かれます。プログラムの実装は極めて簡単です。

固定スプライン (clamped spline)

日本語訳が見当たらなかったので、ここでは仮に固定スプラインと呼びます。固定スプラインでは、最端点 $x=x_0, x_N$ での１階微分の値を与えます。

$\displaystyle f'_1 ( x = x_0 ) = a_1 ( - \Delta x_1 )^2+b_1( - \Delta x_1 ) + c_1 = \alpha \\ f'_N(x=x_N)=c_N=\beta$

$\alpha, \beta$ の値は任意です。 $y''_j$ の条件式は、式(2)-(4)を代入して直ちに求まります。

[ $y''_j$ の条件式]

$2\Delta x_1 y''_0+\Delta x_1 y''_1=6(s_0-\alpha)$
$2\Delta x_N y''_N+\Delta x_N y''_{N-1}=6(\beta-s_{N-1})$

[ 解くべき式 ]

$\displaystyle \begin{pmatrix} 2\Delta x_1 & \Delta x_1 & 0 & \cdots & \cdots & 0 \\ \hline \Delta x_1 & 2(\Delta x_1+\Delta x_2) & \Delta x_2 & & \huge{0} & \\ &\Delta x_2 & \ddots & \ddots & & \\ & & \ddots & \ddots & \Delta x_{N-1} & \\ & \huge{0} & & \Delta x_{N-1} & 2(\Delta x_{N-1}+\Delta x_N) & \Delta x_N \\ \hline 0 & \cdots & \cdots & 0 & \Delta x_N & 2\Delta x_N \\ \end{pmatrix} \begin{pmatrix} y''_0 \\ \vdots \\ y''_N \\ \end{pmatrix} = \begin{pmatrix} 6(s_0-\alpha) \\ \hline 6\Delta s_1 \\ \vdots \\ 6\Delta s_{N-1} \\ \hline 6(\beta-s_{N-1}) \\ \end{pmatrix}$

ここでも係数行列が３重対角行列になっており、３重対角行列アルゴリズムが有効です。

ノットなしスプライン (not-a-knot spline)

ノット（節点）とは、隣接する区間の曲線が互いに結合する点（データ点）のことです。ノットなしスプラインでは、 $x=x_1,x_{N-1}$ で３階微分が連続になるように $y''_j$ を決めます。自然スプラインや固定スプラインと異なり、最端点 $x=x_0, x_N$ での１階、２階微分の値を指定しません。

$\displaystyle f'''_1 ( x = x_1 ) = 6a_1 = 6a_2 = f'''_2 ( x = x_1 ) \\ f'''_{N-1}(x=x_{N-1}) = 6a_{N-1} = 6a_N = f'''_N(x=x_{N-1})$

最端点での１階、２階微分の値が不明な場合、ノットなしスプラインが有用です。

ノットなしスプラインでは、 $x=x_1,x_{N-1}$ で０～３階微分がすべて連続になるため、 $f_1(x)$ と $f_2(x)$ 、 $f_{N-1}(x)$ と $f_N(x)$ はそれぞれ全く同じ関数になります。つまり、見かけ上 $x=x_1,x_{N-1}$ のノットが消失します。 $y''_j$ の条件式は、式(3)を代入して直ちに求まります。

[ $y''_j$ の条件式]

$\Delta x_2 y''_0 - (\Delta x_1 + \Delta x_2) y''_1 + \Delta x_1 y''_2 = 0$
$\Delta x_N y''_{N-2} - (\Delta x_{N-1} + \Delta x_N) y''_{N-1} + \Delta x_{N-1} y''_N = 0$

[ 解くべき式 ]

$\displaystyle \begin{pmatrix} \Delta x_2 & -(\Delta x_1+\Delta x_2) & [\Delta x_1] & 0 & \cdots & 0 \\ \hline \Delta x_1 & 2(\Delta x_1+\Delta x_2) & \Delta x_2 & & \huge{0} & \\ &\Delta x_2 & \ddots & \ddots & & \\ & & \ddots & \ddots & \Delta x_{N-1} & \\ & \huge{0} & & \Delta x_{N-1} & 2(\Delta x_{N-1}+\Delta x_N) & \Delta x_N \\ \hline 0 & \cdots & 0 & [\Delta x_N] & -(\Delta x_{N-1}+\Delta x_N) & \Delta x_{N-1} \\ \end{pmatrix} \begin{pmatrix} y''_0 \\ \vdots \\ y''_N \\ \end{pmatrix} = \begin{pmatrix} 0 \\ \hline 6\Delta s_1 \\ \vdots \\ 6\Delta s_{N-1} \\ \hline 0 \\ \end{pmatrix}$

[ ] で囲った要素により、係数行列が３重対角行列になっていません。第２行の $\Delta x_1 / \Delta x_2$ 倍を第１行から引き、第 $N$ 行の $\Delta x_N / \Delta x_{N-1}$ 倍を第 $(N+1)$ 行から引くことで、[ ] で囲った要素を消去します。適当に定数倍して式を整えると、

$\displaystyle \begin{pmatrix} \Delta x_1-\Delta x_2 & 2\Delta x_1+\Delta x_2 & 0 & \cdots & \cdots & 0 \\ \hline \Delta x_1 & 2(\Delta x_1+\Delta x_2) & \Delta x_2 & & \huge{0} & \\ &\Delta x_2 & \ddots & \ddots & & \\ & & \ddots & \ddots & \Delta x_{N-1} & \\ & \huge{0} & & \Delta x_{N-1} & 2(\Delta x_{N-1}+\Delta x_N) & \Delta x_N \\ \hline 0 & \cdots & \cdots & 0 & 2\Delta x_N+\Delta x_{N-1} & \Delta x_N-\Delta x_{N-1} \\ \end{pmatrix} \begin{pmatrix} y''_0 \\ \vdots \\ y''_N \\ \end{pmatrix} = \begin{pmatrix} 6\frac{\Delta s_1}{1+\Delta x_2 / \Delta x_1} \\ \hline 6\Delta s_1 \\ \vdots \\ 6\Delta s_{N-1} \\ \hline 6\frac{\Delta s_{N-1}}{1+\Delta x_{N-1} / \Delta x_N} \\ \end{pmatrix}$

を得ます。これで、３重対角行列アルゴリズムを使える形になりました。

＜注意＞
$\Delta x_1-\Delta x_2 = 0$ または $\Delta x_N-\Delta x_{N-1} = 0$ の場合、３重対角行列アルゴリズムをそのまま適用するとゼロ除算が起きます。この場合は、 $y''_1, y''_{N-1}$ についてそれぞれ独立に解けますので、これらを消去してから適用します。

$\displaystyle \begin{pmatrix} \hline \Delta x_1 & \Delta x_2 & & & & \\ 0 & 2(\Delta x_2+\Delta x_3) & \Delta x_3 & & \huge{0} & \\ &\Delta x_3 & \ddots & \ddots & & \\ & & \ddots & \ddots & \Delta x_{N-2} & \\ & \huge{0} & & \Delta x_{N-2} & 2(\Delta x_{N-2}+\Delta x_{N-1}) & 0 \\ & & & & \Delta x_{N-1} & \Delta x_N \\ \hline \end{pmatrix} \begin{pmatrix} y''_0 \\ y''_2 \\ \vdots \\ y''_{N-2} \\ y''_N \\ \end{pmatrix} = \begin{pmatrix} \hline 6\Delta s_1 -2(\Delta x_1+\Delta x_2)y''_1 \\ 6\Delta s_2 -\Delta x_2y''_1 \\ 6\Delta s_3 \\ \vdots \\ 6\Delta s_{N-3} \\ 6\Delta s_{N-2} -\Delta x_{N-1}y''_{N-1}\\ 6\Delta s_{N-1} -2(\Delta x_{N-1}+\Delta x_N)y''_{N-1}\\ \hline \end{pmatrix}$

$\displaystyle y''\_1 = 6\frac{\Delta s_1}{1+\Delta x_2 / \Delta x_1}\frac{1}{2\Delta x_1+\Delta x_2}$

$\displaystyle y''_{N-1} = 6\frac{\Delta s_{N-1}}{1+\Delta x_{N-1} / \Delta x_N}\frac{1}{2\Delta x_N+\Delta x_{N-1}}$

上記の場合以外で、 $\Delta x_1-\Delta x_2$ または $\Delta x_N-\Delta x_{N-1}$ が $0$ に近い場合は、計算が不安定になると思われます。この場合は、ガウスの掃き出し法など、別のアルゴリズムを使った方が良いかもしれません。

周期的スプライン (periodic spline)

周期的スプラインでは、両最端点 $x=x_0, x_N$ での１階、２階微分の値が互いに一致するように $y''_j$ を決めます。 $y_0=y_N$ であれば、スプライン曲線は区間 $[x_0,x_N$ ] の外に周期的に繋がります。

$\displaystyle f'_1(x=x_0) = 3a_1(-\Delta x_1)^2+2b_1(-\Delta x_1)+c_1 = c_N = f'_N(x=x_N) \\ f''_1(x=x_0) = 6a_1(-\Delta x_1)+2b_1 = 2b_N = f''_N(x=x_N)$

$y''_j$ の条件式は、式(2)-(4)を代入して直ちに求まります。

[ $y''_j$ の条件式]

$2\Delta x_1 y''_0+\Delta x_1 y''_1+\Delta x_N y''_{N-1}+2\Delta x_N y''_N=6(s_0-s_{N-1})$
$y''_0-y''_N=0$

[ 解くべき式 ]

$\displaystyle \begin{pmatrix} 2\Delta x_1 & \Delta x_1 & 0 & \cdots & 0 & \Delta x_N & 2\Delta x_N \\ \hline \Delta x_1 & 2(\Delta x_1+\Delta x_2) & \Delta x_2 & & & \huge{0} & \\ &\Delta x_2 & \ddots & & \ddots & & \\ & & \ddots & & \ddots & \Delta x_{N-1} & \\ & \huge{0} & & & \Delta x_{N-1} & 2(\Delta x_{N-1}+\Delta x_N) & \Delta x_N \\ \hline 1 & 0 & \cdots & \cdots & \cdots & 0 & -1 \\ \end{pmatrix} \begin{pmatrix} y''_0 \\ \vdots \\ y''_N \\ \end{pmatrix} = \begin{pmatrix} 6(s_0-s_{N-1}) \\ \hline 6\Delta s_1 \\ \vdots \\ 6\Delta s_{N-1} \\ \hline 0 \\ \end{pmatrix}$

３重対角行列アルゴリズムを適用するために変形を行います。第 $(N+1)$ 行の $2\Delta x_N$ 倍を第 $1$ 行、 $\Delta x_N$ 倍を第 $N$ 行に加えます。このとき、第 $(N+1)$ 行の式 $y''_0 - y''_N=0$ を別に書くと、解くべき式は次の通りになります。

$\displaystyle \begin{pmatrix} 2(\Delta x_N+\Delta x_1) & & \Delta x_1 & 0 & \cdots & 0 & & [\Delta x_N] \\ \hline \Delta x_1 & & 2(\Delta x_1+\Delta x_2) & \Delta x_2 & & & \huge{0} & \\ & & \Delta x_2 & \ddots & & \ddots & & \\ & \huge{0} & & \ddots & & \ddots & & \Delta x_{N-1} \\ [\Delta x_N] & & & & & \Delta x_{N-1} & & 2(\Delta x_{N-1}+\Delta x_N) \\ \hline \end{pmatrix} \begin{pmatrix} y''_0 \\ \vdots \\ y''_{N-1} \\ \end{pmatrix} = \begin{pmatrix} 6(s_0-s_{N-1}) \\ \hline 6\Delta s_1 \\ \vdots \\ 6\Delta s_{N-1} \\ \hline \end{pmatrix} \\$

$\displaystyle y''\_0 - y''\_N=0$

右上と左下の角に非ゼロの要素（[ ] で囲った要素）がありますが、これも３重対角行列アルゴリズムを適用できる形です。例えば、次のWikipediaのページ ("Variants") を参照。

en.wikipedia.org

各バリエーションの比較テスト

自然スプライン、固定スプライン、ノットなしスプライン

FORTRANで即席のプログラムを作成しました。SPLINE_INT() を呼び出して係数を決定したのち、SPLINE() を呼び出して区間 $[x_0, x_N$ ] の任意の点を計算します。また、SPLINE_INT() では、左端と右端とで別々に３次スプラインの種類を指定します。

▶クリックで展開 ("splines.f")

      PROGRAM MAIN
      IMPLICIT REAL*8(A-H, O-Z)
      ALLOCATABLE XJ(:), YJ(:),
     $            X(:), YN(:), YC(:), YK(:)

C ファイルからデータを読み込む
      NDATA = 0  ! データ点の数
      OPEN(UNIT=10, FILE="input1.txt", STATUS="OLD", ACTION="READ")
      DO
         READ(UNIT=10, FMT=*, END=10) DUMMY
         NDATA = NDATA + 1
      ENDDO
10    N = NDATA - 1
      ALLOCATE(XJ(0:N), YJ(0:N))
      REWIND(UNIT=10)
      DO J=0, N
         READ(UNIT=10, FMT=*) XJ(J), YJ(J)
      ENDDO
      CLOSE(UNIT=10)
      
C スプライン補間を実行
C 元のデータ点数の100倍の数の等間隔点で補間
      KMAX = NDATA * 100 - 1
      DX = (XJ(N) - XJ(0)) / KMAX
      ALLOCATE(X(0:KMAX), YN(0:KMAX), YC(0:KMAX), YK(0:KMAX))

C 固定スプラインの両端での傾き
      ALPHA = -1.D0
      BETA  =  1.D0

C 自然スプライン
C 'N'は自然(NATURAL)
C ALPHA, BETA はダミー
      CALL SPLINE_INI('N', ALPHA, 'N', BETA, XJ, YJ, N)
      XTMP = XJ(0) - DX
      DO K=0, KMAX
         XTMP = XTMP + DX
         CALL SPLINE(XJ, N, XTMP, Y)
         X(K) = XTMP
         YN(K) = Y
      ENDDO

C 固定スプライン
C 'C'は固定(CLAMPED)
      CALL SPLINE_INI('C', ALPHA, 'C', BETA, XJ, YJ, N)
      DO K=0, KMAX
         CALL SPLINE(XJ, N, X(K), Y)
         YC(K) = Y
      ENDDO

C ノットなしスプライン
C 'K'はノット(KNOT)
C ALPHA, BETA はダミー
      CALL SPLINE_INI('K', ALPHA, 'K', BETA, XJ, YJ, N)
      DO K=0, KMAX
         CALL SPLINE(XJ, N, X(K), Y)
         YK(K) = Y
      ENDDO

      OPEN(UNIT=20, FILE="spline.txt", STATUS="REPLACE", ACTION="WRITE")
      DO K=0, KMAX
         WRITE(UNIT=20, FMT=*)X(K), YN(K), YC(K), YK(K)
      ENDDO
      CLOSE(UNIT=20)
      
      DEALLOCATE(XJ, YJ, X, YN, YC, YK)
      
      STOP
      END
      
      
C 係数を決定、保存
      SUBROUTINE SPLINE_INI(KEY0, ALPHA, KEY1, BETA, XJ, YJ, N)
      IMPLICIT REAL*8(A-H, O-Z)

      INTENT(IN) KEY0, KEY1, XJ, YJ, N
      
      CHARACTER*1 KEY0, KEY1
C スプラインの種類 (N/C/K)
C     KEY0: 左端
C     KEY1: 右端

      DIMENSION XJ(0:N), YJ(0:N),
     $          UJ(0:N), SJ(0:N-1), DXJ(1:N), DSJ(1:N-1),
     $          TDMA_A(1:N), TDMA_B(0:N), TDMA_C(0:N-1), TDMA_D(0:N)
      ALLOCATABLE COEFF(:, :)
      SAVE COEFF

      IF(.NOT.(KEY0.EQ.'N'.OR.KEY0.EQ.'C'.OR.KEY0.EQ.'K').OR.
     $   .NOT.(KEY1.EQ.'N'.OR.KEY1.EQ.'C'.OR.KEY1.EQ.'K'))
     $   STOP "! KEY ERROR"
     
      IF(ALLOCATED(COEFF)) DEALLOCATE(COEFF)

      DXJ(1) = XJ(1) - XJ(0)
      SJ(0) = (YJ(1) - YJ(0)) / DXJ(1)
      DO J=1, N-1
         DXJ(J+1) = XJ(J+1) - XJ(J)
         SJ(J) = (YJ(J+1) - YJ(J)) / DXJ(J+1)
         DSJ(J) = SJ(J) - SJ(J-1)
      ENDDO

C 式(5)を作る。
C TDMA は TriDiagonal Matrix Algorithm (３重対角行列アルゴリズム)。
C 
C [ TDMA_B(0) TDMA_C(0)                                 ] [UJ(0)]   [TDMA_D(0)]
C [ TDMA_A(1) TDMA_B(1) TDMA_C(1)            0          ] [     ]   [         ]
C [           TDMA_A(2) TDMA_B(2)    ...                ] [ ... ] = [   ...   ]
C [                        ...       ...                ] [     ]   [         ]
C [        0                                TDMA_C(N-1) ] [     ]   [         ]
C [                               TDMA_A(N) TDMA_B(N)   ] [UJ(N)]   [TDMA_D(N)]

      DO J=1, N-1
         TDMA_A(J) = DXJ(J)
         TDMA_B(J) = DXJ(J) + DXJ(J+1)
         TDMA_C(J) = DXJ(J+1)
         TDMA_D(J) = DSJ(J)
      ENDDO
      TDMA_B = 2.D0 * TDMA_B
      TDMA_D = 6.D0 * TDMA_D

C 左端
      SELECT CASE (KEY0)
         CASE ('N')
            TDMA_B(0) = 1.D0
            TDMA_C(0) = 0.D0
            TDMA_D(0) = 0.D0
         CASE ('C')
            TDMA_B(0) = 2.D0 * DXJ(1)
            TDMA_C(0) = DXJ(1)
            TDMA_D(0) = 6.D0 * (SJ(0) - ALPHA)
         CASE ('K')
            TDMA_B(0) = DXJ(1) - DXJ(2)
            TDMA_C(0) = 2.D0 * DXJ(1) + DXJ(2)
            TDMA_D(0) = 6.D0 * DSJ(1) / (1.D0 + DXJ(2) / DXJ(1))
      END SELECT
      
C 右端
      SELECT CASE (KEY1)
         CASE ('N')
            TDMA_B(N) = 1.D0
            TDMA_A(N) = 0.D0
            TDMA_D(N) = 0.D0
         CASE ('C')
            TDMA_B(N) = 2.D0 * DXJ(N)
            TDMA_A(N) = DXJ(N)
            TDMA_D(N) = 6.D0 * (BETA - SJ(N-1))
         CASE ('K')
            TDMA_B(N) = DXJ(N) - DXJ(N-1)
            TDMA_A(N) = 2.D0 * DXJ(N) + DXJ(N-1)
            TDMA_D(N) = 6.D0 * DSJ(N-1) / (1.D0 + DXJ(N-1) / DXJ(N))
      END SELECT
      

      NMIN=0
      NMAX=N

C ゼロ除算の回避
      IF(TDMA_B(0).EQ.0.D0)THEN
         UJ(1) = TDMA_D(0) / TDMA_C(0)
         TDMA_D(1) = TDMA_D(1) - TDMA_B(1) * UJ(1)
         TDMA_D(2) = TDMA_D(2) - TDMA_A(2) * UJ(1)
         TDMA_B(1) = TDMA_A(1)
         TDMA_A(2) = 0.D0
         NMIN=1
      ELSEIF(TDMA_B(N).EQ.0.D0)THEN
         UJ(N-1) = TDMA_D(N) / TDMA_A(N)
         TDMA_D(N-1) = TDMA_D(N-1) - TDMA_B(N-1) * UJ(N-1)
         TDMA_D(N-2) = TDMA_D(N-2) - TDMA_C(N-2) * UJ(N-1)
         TDMA_B(N-1) = TDMA_C(N-1)
         TDMA_A(N-2) = 0.D0
         NMAX=N-1
      ENDIF

C TDMA を実行
      DO J=NMIN+1, NMAX
         W = TDMA_A(J) / TDMA_B(J-1)
         TDMA_B(J) = TDMA_B(J) - W * TDMA_C(J-1)
         TDMA_D(J) = TDMA_D(J) - W * TDMA_D(J-1)
      ENDDO
      
      UJ(NMAX) = TDMA_D(NMAX) / TDMA_B(NMAX)
      DO J=NMAX-1, NMIN, -1
         UJ(J) = (TDMA_D(J) - TDMA_C(J) * UJ(J+1)) / TDMA_B(J)
      ENDDO
      
      IF(NMIN.EQ.1)THEN
         TMP = UJ(0)
         UJ(0) = UJ(1)
         UJ(1) = TMP
      ENDIF

      IF(NMAX.EQ.N-1)THEN
         TMP = UJ(N)
         UJ(N) = UJ(N-1)
         UJ(N-1) = TMP
      ENDIF
      
C 係数を計算
      ALLOCATE(COEFF(4,N))
C     COEFF(1, J): A_J
C     COEFF(2, J): B_J
C     COEFF(3, J): C_J
C     COEFF(4, J): D_J
      DO J=1, N
         COEFF(1, J) = (UJ(J)-UJ(J-1)) / DXJ(J) / 6.D0
         COEFF(2, J) = .5D0 * UJ(J)
         COEFF(3, J) = SJ(J-1)
     $               + DXJ(J) * (2.D0 * UJ(J) + UJ(J-1)) / 6.D0
         COEFF(4, J) = YJ(J)
      ENDDO
      
      RETURN
      
      
C スプライン補間の計算
      ENTRY SPLINE(XJ, N, X, Y)

C X がどの区間にあるか判断
      JMIN = 0
      JMAX = N
      DO
         IF(JMAX.EQ.JMIN+1) EXIT
         JAVE = (JMIN + JMAX) / 2
         IF(XJ(JAVE).LT.X)THEN
            JMIN = JAVE
         ELSE
            JMAX = JAVE
         ENDIF
      ENDDO
      J = JMAX
      
C 計算を実行
      XDIF = X - XJ(J)
      Y = COEFF(1, J)
      DO NORD=2, 4
         Y = Y * XDIF + COEFF(NORD, J)
      ENDDO
      
      RETURN
      
      END

"input1.txt" からデータ点を読み取ってスプライン補間し、結果を "spline.txt" に出力します。ここでは次の適当なデータ点で試してみます。

"input1.txt"

0.2   0.4392
0.7   0.8638
1.6   0.5449
2.3   0.2019
3.0   0.0190
4.0  -0.0374

補間の結果を、下の図２に示します。

固定スプラインの結果は、最端点 $x=0.2, 4.0$ でやや無理な傾きを指定しているため $(\alpha=-1, \beta=1)$ 、自然スプラインとノットなしスプラインの結果に比べて大きくたわんでいます。一方、自然スプラインとノットなしスプラインの結果の差は、図の左側で大きくなっているように見えます。最左端 $x=0.2$ の付近を拡大してみます。

図３では、 $x=0.2$ での接線を点線で描き加えました（係数が分かっているので、当然接線を求めることができます）。 $x$ が増加したときの接線との乖離の程度を見比べると、明らかに自然スプラインの方が直線に近い形をしていることが分かります。

図４では、 $f_1(x)$ を点線で描き加えました。自然スプラインは $x\gt x_1$ で乖離していますが、一方でノットなしスプラインは $x\gt x_2$ で乖離しています。見かけ上 $x=x_1$ のノットが消失していることが分かります。

周期的スプライン

FORTRANで即席のプログラムを作成しました。SPLINE_INT() を呼び出して係数を決定したのち、SPLINE() を呼び出して区間 [x0,xN] の任意の点を計算します。

▶クリックで展開 ("spline_periodic.f")

      PROGRAM MAIN
      IMPLICIT REAL*8(A-H, O-Z)
      ALLOCATABLE XJ(:), YJ(:)

C ファイルからデータを読み込む
      NDATA = 0  ! データ点の数
      OPEN(UNIT=10, FILE="input2.txt", STATUS="OLD", ACTION="READ")
      DO
         READ(UNIT=10, FMT=*, END=10) DUMMY
         NDATA = NDATA + 1
      ENDDO
10    N = NDATA - 1
      ALLOCATE(XJ(0:N), YJ(0:N))
      REWIND(UNIT=10)
      DO J=0, N
         READ(UNIT=10, FMT=*) XJ(J), YJ(J)
      ENDDO
      CLOSE(UNIT=10)
      
C スプライン補間を実行
C 元のデータ点数の100倍の数の等間隔点で補間
      KMAX = NDATA * 100 - 1
      DX = (XJ(N) - XJ(0)) / KMAX

C 周期的スプライン
      OPEN(UNIT=20, FILE="splineP.txt",
     $     STATUS="REPLACE", ACTION="WRITE")
      CALL SPLINE_INI(XJ, YJ, N)
      XX = XJ(0) - DX
      DO K=0, KMAX
         XX = XX + DX
         CALL SPLINE(XJ, N, XX, YY)
         WRITE(UNIT=20, FMT=*)XX, YY
      ENDDO
      CLOSE(UNIT=20)
      
      DEALLOCATE(XJ, YJ)
      
      STOP
      END
      
      
C 係数を決定、保存
      SUBROUTINE SPLINE_INI(XJ, YJ, N)
      IMPLICIT REAL*8(A-H, O-Z)

      INTENT(IN) XJ, YJ, N
      
      DIMENSION XJ(0:N), YJ(0:N),
     $          UJ(0:N), SJ(0:N-1), DXJ(1:N), DSJ(1:N-1),
     $          TDMA_A(0:N-1), TDMA_B(0:N-1),
     $          TDMA_C(0:N-1), TDMA_D(0:N-1),
     $          PJ(N-1), QJ(N-1), RJ(N-1)
      ALLOCATABLE COEFF(:, :)
      SAVE COEFF

     
      IF(ALLOCATED(COEFF)) DEALLOCATE(COEFF)

      DXJ(1) = XJ(1) - XJ(0)
      SJ(0) = (YJ(1) - YJ(0)) / DXJ(1)
      DO J=1, N-1
         DXJ(J+1) = XJ(J+1) - XJ(J)
         SJ(J) = (YJ(J+1) - YJ(J)) / DXJ(J+1)
         DSJ(J) = SJ(J) - SJ(J-1)
      ENDDO

C 解くべき式
C 
C [ TDMA_B(0)   TDMA_C(0)                       TDMA_A(0)   ] [UJ(0)  ]   [TDMA_D(0)  ]
C [ TDMA_A(1)   TDMA_B(1)   TDMA_C(1)        0              ] [       ]   [           ]
C [             TDMA_A(2)   TDMA_B(2)    ...                ] [  ...  ] = [    ...    ]
C [                            ...       ...                ] [       ]   [           ]
C [                 0                           TDMA_C(N-2) ] [       ]   [           ]
C [ TDMA_C(N-1)                     TDMA_A(N-1) TDMA_B(N-1) ] [UJ(N-1)]   [TDMA_D(N-1)]

      TDMA_A(0) = DXJ(N)
      TDMA_B(0) = DXJ(N) + DXJ(1)
      TDMA_C(0) = DXJ(1)
      TDMA_D(0) = SJ(0) - SJ(N-1)
      DO J=1, N-1
         TDMA_A(J) = DXJ(J)
         TDMA_B(J) = DXJ(J) + DXJ(J+1)
         TDMA_C(J) = DXJ(J+1)
         TDMA_D(J) = DSJ(J)
      ENDDO
      TDMA_B = 2.D0 * TDMA_B
      TDMA_D = 6.D0 * TDMA_D


C 下記２つの線形連立方程式に TDMA (３重対角行列アルゴリズム) を適用。
C
C [ TDMA_B(1) TDMA_C(1)                                   ] [PJ(1)  ]   [TDMA_D(1)  ]
C [ TDMA_A(2) TDMA_B(2) TDMA_C(2)              0          ] [       ]   [           ]
C [           TDMA_A(3) TDMA_B(3)    ...                  ] [  ...  ] = [    ...    ]
C [        0               ...       ...      TDMA_C(N-2) ] [       ]   [           ]
C [                               TDMA_A(N-1) TDMA_B(N-1) ] [PJ(N-1)]   [TDMA_D(N-1)]
C

      CALL TDMA(N-1, TDMA_A(1), TDMA_B(1), TDMA_C(1), TDMA_D(1), PJ)

C
C [ TDMA_B(1) TDMA_C(1)                                   ] [QJ(1)  ]   [-TDMA_A(1)  ]
C [ TDMA_A(2) TDMA_B(2) TDMA_C(2)              0          ] [       ]   [     0      ]
C [           TDMA_A(3) TDMA_B(3)    ...                  ] [  ...  ] = [     ...    ]
C [        0               ...       ...      TDMA_C(N-1) ] [       ]   [     0      ]
C [                               TDMA_A(N-1) TDMA_B(N-1) ] [QJ(N-1)]   [-TDMA_C(N-1)]

      RJ = 0.D0
      RJ(1) = -TDMA_A(1)
      RJ(N-1) = -TDMA_C(N-1)
      
      CALL TDMA(N-1, TDMA_A(1), TDMA_B(1), TDMA_C(1), RJ, QJ)

C UJが求まる
      UJ(0) = (TDMA_D(0) - TDMA_A(0) * PJ(N-1) - TDMA_C(0) * PJ(1))
     $      / (TDMA_B(0) + TDMA_A(0) * QJ(N-1) + TDMA_C(0) * QJ(1))
      DO J=1, N-1
         UJ(J) = PJ(J) + UJ(0) * QJ(J)
      ENDDO
      UJ(N) = UJ(0)

C 係数を計算
      ALLOCATE(COEFF(4,N))
C     COEFF(1, J): A_J
C     COEFF(2, J): B_J
C     COEFF(3, J): C_J
C     COEFF(4, J): D_J
      DO J=1, N
         COEFF(1, J) = (UJ(J)-UJ(J-1)) / DXJ(J) / 6.D0
         COEFF(2, J) = .5D0 * UJ(J)
         COEFF(3, J) = SJ(J-1)
     $               + DXJ(J) * (2.D0 * UJ(J) + UJ(J-1)) / 6.D0
         COEFF(4, J) = YJ(J)
      ENDDO

      RETURN
      
      
C スプライン補間の計算
      ENTRY SPLINE(XJ, N, X, Y)

C X がどの区間にあるか判断
      JMIN = 0
      JMAX = N
      DO
         IF(JMAX.EQ.JMIN+1) EXIT
         JAVE = (JMIN + JMAX) / 2
         IF(XJ(JAVE).LT.X)THEN
            JMIN = JAVE
         ELSE
            JMAX = JAVE
         ENDIF
      ENDDO
      J = JMAX
      
C 計算を実行
      XDIF = X - XJ(J)
      Y = COEFF(1, J)
      DO NORD=2, 4
         Y = Y * XDIF + COEFF(NORD, J)
      ENDDO
            
      RETURN
      
      END
      
      
C TDMA
      SUBROUTINE TDMA(M, A, B, C, D, X)
      IMPLICIT REAL*8(A-H, O-Z)
      DIMENSION A(M), B(M), C(M), D(M), X(M)
      
      DO J=2, M
         W = A(J) / B(J-1)
         B(J) = B(J) - W * C(J-1)
         D(J) = D(J) - W * D(J-1)
      ENDDO
      
      X(M) = D(M) / B(M)
      DO J=M-1, 1, -1
         X(J) = (D(J) - C(J) * X(J+1)) / B(J)
      ENDDO
      
      RETURN

      END

"input2.txt" からデータ点を読み取ってスプライン補間し、結果を "splineP.txt" に出力します。ここでは次の適当なデータ点 $(y=\sin x)$ で試してみます。

"input2.txt"

0.7854   0.7071
1.0472   0.8660
1.5708   1.0000
2.0944   0.8660
2.3562   0.7071
2.6180   0.5000
3.1416  -0.0004
3.6652  -0.5000
3.9270  -0.7071
4.1888  -0.8660
4.7124  -1.0000
5.2360  -0.8660
5.4978  -0.7071
5.7596  -0.5000
6.2832   0.0000
6.8068   0.5000
7.0686   0.7071

補間の結果を、下の図５に示します。

補間を行った区間は、中央の縦線 ( $x=2.25\pi$ ) より左側です。右側の曲線は、左側の補間の結果を１周期ずらしてプロットしたものです。１階微分を破線、２階微分を点線であわせて示しました。いずれも、縦線のところで滑らかにつながっています。

比較のため、自然スプライン、固定スプライン、ノットなしスプラインによる補間の結果を示します。

図６は自然スプラインの結果です。実線の結果は悪くないように見えますが、１階微分は明らかに不連続です。２階微分は、最端点でゼロになる条件に無理があり、縦線の付近で不自然な振る舞いをしています。

図７は固定スプラインの結果です。両最端点での１階微分は、どちらも $\sin'(2.25\pi) =1/\sqrt{2}$ としました。１階、２階微分とも滑らかにつながっているように見えます（実際に滑らかにつながっている訳ではありません）。最端点での１階微分の値が予め分かっているため、それらしい結果が得られました。

図８はノットなしスプラインの結果です。１階微分まで滑らかにつながっているように見えます（実際に滑らかにつながっている訳ではありません）。２階微分は明らかに不連続ですが、自然スプラインの結果と比べると、こちらの結果の方が無理がありません。

まとめ

３次スプラインの４つのバリエーションをまとめました。また、即席のプログラムで実際にスプライン補間を行い、それぞれの特徴を見ました（表１）。

表１　３次スプラインのバリエーション
周期的スプライン	両最端点の１階、２階微分の値を互いに一致させる。周期的なデータの補間に適する。
自然スプライン	最端点での２階微分の値をゼロとする。直線で無理なく補外できるデータの補間に適する。
固定スプライン	最端点での１階微分の値を指定する。最端点での傾きが分かっているデータの補間に適する。
ノットなしスプライン	最端点での条件を指定しない。一般的なデータの補間に用いられる。

Doc's log

ドクログ

３次スプライン補間 (cubic spline interpolation)

準備

３次スプラインの定義

解くべき式の導出

３次スプラインのバリエーションと計算法

自然スプライン (natural spline)

固定スプライン (clamped spline)

ノットなしスプライン (not-a-knot spline)

周期的スプライン (periodic spline)

各バリエーションの比較テスト

自然スプライン、固定スプライン、ノットなしスプライン

周期的スプライン

まとめ