データ分析

スプライン補間の具体的な計算

More than 1 year has passed since last update.

3次スプライン補間を具体的に計算していきます.

3次スプライン補間の条件

条件は5つある.詳細は後述する.

  1. $S_j(x_j) = y_j$
  2. $S_j(x_{j + 1}) = y_{j + 1}$
  3. $S_j'(x_{j + 1}) - S_{j + 1}'(x_{j + 1}) = 0$
  4. $S_j''(x_{j + 1}) - S_{j + 1}''(x_{j + 1}) = 0$
  5. $S_0''(x_0) = 0 , \ S_{n - 1}''(x_{n}) = 0$

具体的な計算の準備(方程式編)

3点$(x_0, y_0), (x_1, y_1), (x_2, y_2)$が与えられたと仮定する.
3点により2つの区間$S_0, S_1$が得られる.
3次スプライン補間を考えるので,
 $S_j = a_jx^3+b_jx^2+c_jx+d_j$
とする.
これの1階微分は 
 $S_j'(x)=3a_jx^2+2b_jx+c_j$
2階微分は
 $S_j''(x)=6a_jx+2b_j$ 
である.

条件1

この条件は求める区間の多項式は対応する与えられた点を通るという意味である.
区間$S_0$は$(x_0, y_0)$を通り,区間$S_1$は$(x_1, y_1)$を通る.
従って
 $S_0(x_0) = y_0$
 $S_1(x_1) = y_1$
を得る.
すなわち
 $x_0^3a_0+x_0^2b_0+x_0c_0+1d_0=y_0$
 $x_1^3a_1+x_1^2b_1+x_1c_1+1d_1=y_1$

条件2

この条件は求める区間の多項式は右端の与えられた点を通るという意味である.
区間$S_0$は$(x_1, y_1)$を通り,区間$S_1$は$(x_2, y_2)$を通る.
従って
 $S_0(x_1) = y_1$
 $S_1(x_2) = y_2$
を得る.

すなわち
 $x_1^3a_0+x_1^2b_0+x_1c_0+1d_0=y_1$
 $x_2^3a_1+x_2^2b_1+x_2c_1+1d_1=y_2$

条件3

この条件は1次導関数がすべての接点上で連続であるという意味である.
区間が2本のため,接点は$(x_1, y_1)$の1点だけである.
従って
 $S_0'(x_1) - S_1'(x_1) = 0$
を得る.

すなわち
 $3x_1^2a_0+2x_1b_0+1c_0+0-3x_1^2-2x_1b_1-1c_1-0=0$

条件4

この条件は2次導関数がすべての接点上で連続であるという意味である.
同様に区間が2本のため,接点は$(x_1, y_1)$の1点だけである.
従って
 $S_0''(x_1) - S_1''(x_1) = 0$
を得る.

すなわち
 $6x_1a_0+2b_0+0+0-6x_1a_1-2b_1-0-0=0$

条件5

この条件は両端の2次導関数を0とするという意味である.
(この条件は,以上の4つの条件だけでは多項式が求まらないので,仕方なく?付け足した条件である.)
両端は区間$S_0$の$(x_0, y_0)$と,区間$S_1$の$(x_2, y_2)$である.
従って
 $S_0''(x_0) = 0, \ S_1''(x_2) = 0$
を得る.

すなわち
 $6x_0a_0+2b_0+0+0=0, \ 6x_2a_1+2b_1+0+0=0$

具体的な計算の準備(行列編)

求めたいのは多項式の係数$(a_0, b_0, c_0, d_0, a_1, b_1, c_1, d_1)$である.
これを行列$coef$とおくと,
 $A*coef = b$
となるように行列$A$と$b$を定めたい.
連立方程式を解いているだけなので,行の順序はどうでもいいが,以下の追加の方が実装しやすいと思う.

条件1を追加する

  A = \left[
    \begin{array}{rrrrrrrr}
      x_0^3 & x_0^2 & x_0 & 1 & 0 & 0 & 0 & 0 \\
      * & * & * & * & * & * & * & * \\
      * & * & * & * & * & * & * & * \\
      * & * & * & * & * & * & * & * \\
      0 & 0 & 0 & 0 & x_1^3 & x_1^2 & x_1 & 1 \\
      * & * & * & * & * & * & * & * \\
      * & * & * & * & * & * & * & * \\
      * & * & * & * & * & * & * & * \\
    \end{array}
  \right]

条件2を追加する

  A = \left[
    \begin{array}{rrrrrrrr}
      x_0^3 & x_0^2 & x_0 & 1 & 0 & 0 & 0 & 0 \\
      x_1^3 & x_1^2 & x_1 & 1 & 0 & 0 & 0 & 0 \\
      * & * & * & * & * & * & * & * \\
      * & * & * & * & * & * & * & * \\
      0 & 0 & 0 & 0 & x_1^3 & x_1^2 & x_1 & 1 \\
      0 & 0 & 0 & 0 & x_2^3 & x_2^2 & x_2 & 1 \\
      * & * & * & * & * & * & * & * \\
      * & * & * & * & * & * & * & * \\
    \end{array}
  \right]

条件3を追加する

  A = \left[
    \begin{array}{rrrrrrrr}
      x_0^3 & x_0^2 & x_0 & 1 & 0 & 0 & 0 & 0 \\
      x_1^3 & x_1^2 & x_1 & 1 & 0 & 0 & 0 & 0 \\
      * & * & * & * & * & * & * & * \\
      * & * & * & * & * & * & * & * \\
      0 & 0 & 0 & 0 & x_1^3 & x_1^2 & x_1 & 1 \\
      0 & 0 & 0 & 0 & x_2^3 & x_2^2 & x_2 & 1 \\
      3x_1^2 & 2x_1 & 1 & 0 & -3x_1^2 & -2x_1 & -1 & 0 \\
      * & * & * & * & * & * & * & * \\
    \end{array}
  \right]

条件4を追加する

  A = \left[
    \begin{array}{rrrrrrrr}
      x_0^3 & x_0^2 & x_0 & 1 & 0 & 0 & 0 & 0 \\
      x_1^3 & x_1^2 & x_1 & 1 & 0 & 0 & 0 & 0 \\
      * & * & * & * & * & * & * & * \\
      * & * & * & * & * & * & * & * \\
      0 & 0 & 0 & 0 & x_1^3 & x_1^2 & x_1 & 1 \\
      0 & 0 & 0 & 0 & x_2^3 & x_2^2 & x_2 & 1 \\
      3x_1^2 & 2x_1 & 1 & 0 & -3x_1^2 & -2x_1 & -1 & 0 \\
      6x_1 & 2 & 0 & 0 & -6x_1 & -2 & 0 & 0 \\
    \end{array}
  \right]

条件5を追加する

  A = \left[
    \begin{array}{rrrrrrrr}
      x_0^3 & x_0^2 & x_0 & 1 & 0 & 0 & 0 & 0 \\
      x_1^3 & x_1^2 & x_1 & 1 & 0 & 0 & 0 & 0 \\
      6x_0 & 2 & 0 & 0 & 0 & 0 & 0 & 0 \\
      0 & 0 & 0 & 0 & 6x_2 & 2 & 0 & 0 \\
      0 & 0 & 0 & 0 & x_1^3 & x_1^2 & x_1 & 1 \\
      0 & 0 & 0 & 0 & x_2^3 & x_2^2 & x_2 & 1 \\
      3x_1^2 & 2x_1 & 1 & 0 & -3x_1^2 & -2x_1 & -1 & 0 \\
      6x_1 & 2 & 0 & 0 & -6x_1 & -2 & 0 & 0 \\
    \end{array}
  \right]

行列bを定める

  b = \left[
    \begin{array}{rrrrrrrr}
      y_0 & y_1 & 0 & 0 & y_1 & y_2 & 0 & 0 \\
    \end{array}
  \right]

最後に

以上でスプライン補間の具体的な計算を終わります.
具体的な数値を代入しようと思いましたが,ここまですればもういいやってなりました.

実装しやすい行の追加は上記の通りだと思うのですが,その辺の説明をどうするか,余力があれば追記したいと思います.