Web教材一覧数値解析
(BOK大区分:1 基礎理論、中区分:2 アルゴリズムとプログラミング、中区分:2 アルゴリズム)

補間(スプライン)


スプライン補間とは

補間に関する基本概念は「補間」を参照してください。
 そこで紹介したラグランジュの補間公式やニュートンの補間公式は、与えられたn+1の点を通る唯一のn次式を求めることがポイントでした。そのため、次数が高くなると、極端な補間値になってしまうことがあります。

スプライン補間は、n+1個の点 (x0,y0), (x1,y1), … , (xn,yn) を、2つの点の区間 [xj~xj+1]ごとに別々の式(一般には3次式)Sj(x) を設定して緩やかな曲線にし、極端な補間値にならないようにする補間法です。計算はかなり複雑になります。

定式化

区間 [xj~xj+1] における3次式を
  Sj=aj(x-xj)3 + bj(x-xj)2 + cj(x-xj) + dj   (j=0~n-1)
とします。
 これらの3次式が滑らかにつながっている必要があります。それで、xj におけるSj-1とSjの1次導関数(傾き)S(1) と2次導関数(滑らかさ)S(2) の値が等しいという条件を加えます。さらに、始点 x0 と終点 xn での2次導関数の値が0であるとの仮定をします。

決定すべき件数 aj, bj, cj, dj の個数は4n個です。
 それに対して、上の条件を満たす方程式は、次の4n個になりますので、連立方程式を解くことにより、各係数の値を決定することができます。

連立方程式

上の関係に、xj, yj の値を代入して連立方程式にすればよいのでが、複雑なので、次のように変形して解く方法がとられます。この変形過程はここでの目的ではないので、結果だけを示します。
  uj=2bj  (j=1~n-1)
    (u0=0, un=0 とする)
  hj=xj+1 - xj  (j=0~n-1)
  vj=6( (yj+1-yj)/hj - (yj-yj-1)/hj-1)  (j=1~n-1)
と置換すると、uj を求める連立方程式は、次のようになります。変数側の要素がすべてxの差の1次式になっていること、対角線とその両側にだけ値を持つ(それ以外は0)である特徴があります。

  u1   u2   u3   … uj-1   uj   uj+1 … un-2   un-1
2(h0+h1)  h1                             = v1
  h1  2(h1+h2)  h2                         = v2
      h2  2(h2+h3) h3                      = v3
             :::                       :
               hj-1 2(hj-1+hj) hj            = vj
                         :::           :
                           hn-2 2(hn-2+hn-1) = vn-1

各係数への復元

条件Aより、dj=yj は直ちに求められます。
 連立方程式を解いて、ujが求まれば、
 bj=(1/2)*uj で計算できます。
 cjは、かなりの式操作がありますが、
  cj=(yj+1-yj)/(xj+1-xj) - (1/6)*(2uj+uj+1)*(xj+1-xj)
から計算できます。
 aj は、条件Dから
  6aj(xj+1-xj) + 2bj=2bj+1
  ∴ aj=(1/6)*(uj+1-uj)/(xj+1-xj)
から計算できます。

数値例

次の6点(n=5)(このページの先頭にあるグラフの赤点)を入力とする補間式を考えます。
  j    0   1   2   3   4   5
  xj   0   2   3   5   7   8
  yj   4   2   6   8   6   8

 hj=xj+1-xj, vj=6( (yj+1-yj)/hj - (yj-yj-1)/hj-1) を計算して付け加えます。
  hj   2   1   2   2   1
  vj      30  -18 -12 18

これより、左の連立方程式が作られ、右の解が得られます。
  u1  u2  u3  u4          解
  6  1       = 30    u1= 5.544
  2  8  2    =-18    u2=-3.262
     2  6  1 =-12    u3=-1.496
        2  6 = 18    u4= 3.499 

次の逆変換により、係数 aj, bj, cj, dj が得られます。
  aj=(1/6)*(uj+1-uj)/(xj+1-xj)
  bj=(1/2)*uj
  cj=(yj+1-yj)/(xj+1-xj) - (1/6)*(2uj+uj+1)*(xj+1-xj)
  dj=yj

その結果、区間ごとの補間式は次のようになります。
 j 区間   aj      bj     cj     dj
 0 0~2: 0.4620(x-0)3 +0  (x-0)2 -2.848(x-0) + 4
 1 2~3: -1.4680(x-2)3 +2.772(x-2)2 +2.696(x-2) + 2
 2 3~5: 0.1472(x-3)3 -1.631(x-3)2 +3.673(x-3) + 6
 3 5~7: 0.4162(x-5)3 -0.748(x-5)2 -1.169(x-5) + 8
 4 7~8: -0.5831(x-7)3 +1.749(x-7)2 +0.834(x-7) + 6

これらをつなぐと上のグラフの曲線になり、赤点近傍で滑らかに接続しています。

この補間式で、x=1, 4, 6 のときの y を計算すると次のようになり、グラフの◇の点になります。
  x=1: 0.4620(1-0)3 +0  (1-0)2 -2.848(1-0) + 4= 1.614
  x=4: 0.1472(4-3)3 -1.631(4-3)2 +3.673(4-3) + 6= 8.189
  x=6: 0.4162(6-5)3 -0.748(6-5)2 -1.169(6-5) + 8= 6.499


計算プログラム

区間数n(データの組数は0~nのn+1組) 補間するxの個数
x[j]
y[j]
補間値x[k]

数値解析へ