Subsections

3 スプライン補間

3.1 区分多項式

ラグランジュの補間はデータ点数が増えてくると関数が振動し,補間の精度が悪くなるの は先に述べたとおりである.そこで,補間する領域をデータ間隔 $ [x_i,x_{i+1}]$ に区切 り,その近傍の値を使い低次の多項式で近似することを考える.区分的に近似関数を使う わけですが,上手に近似をしないと境界でその導関数が不連続になる.導関数が連続にな るように,上手に近似する方法がスプライン補間(spline interpolation)である.

ここでは,通常よくつかわれる3次のスプライン補間について説明する.補間する関数が3次関数 を使うため,そう呼ばれている.これ以降の説明は,文献[1]を参考にした.

データは先と同じように $ (x_0,y_0) ,(x_1,y_1), (x_2,y_2), \cdots, (x_N,y_N)$ と する.そして,区間 $ [x_j,x_{j+1}]$ で補間に使う関数を$ S_j(x)$ とする.この様子を図 5に示す.

図 5: スプライン補間の区分
\includegraphics[keepaspectratio,scale=0.7]{figure/Spline.eps}

3.2 区分多項式の条件と計算方法

3次のスプライン補間を考えるので,区分多項式は

$\displaystyle S_j(x)=a_j(x-x_j)^3+b_j(x-x_j)^2+c_j(x-x_j)+d_j \qquad(j=0,1,2,3,\cdots,N-1)$ (5)

となる.この $ a_j, b_j, c_j, d_j$ を求めることが,スプライン補間の関数をきめる問題 となる.

$ N+1$ 個のデータ数があるため,区分多項式は$ N$ 個ある.したがって,区分多項式の係数 をである未知数は$ 4N$ 個あることになる.これを求めるためには,$ 4N$ 個の方程式が必要 となる.3次のスプライン補間に以下の条件を課して,その係数を求めることにする.

[条件1]
全てのデータ点を通る.各々の$ S_j(x)$ に対して両端での値が決まる ため,2N個の方程式ができる.
[条件2]
各々の区分補間式は,境界点の1次導関数は連続とする.これにより, N-1個の方程式ができる
[条件3]
各々の区分補間式は,境界点の2次導関数も連続とする.これにより, N-1個の方程式ができる.

以上の3つの条件を課すと$ 4N-2$ 個の方程式で未知数である係数の関係を表現できる.実際,未 知数は$ 4N$ 個なので,2個方程式が不足している.この不足を補うために,いろいろな条 件が考えられるが,通常は両端$ x_0$$ x_N$ での2次導関数の値を0とする.すなわち, $ S_0^{\prime\prime}(x_0)=S_{N-1}^{\prime\prime}(x_N)=0$ である.これを自然スプラ イン(natural spline)という.自然スプライン以外には,両端の1次導関数の値を指定す るものもある.

これで全ての条件が決まった.あとは,この条件に満たす連立方程式を求めるだけである. このように,必要な条件が決まった場合,$ 4N$ 個の未知数 $ a_j, b_j, c_j, d_j$ を既知の $ x_j,y_j$ を使って連立方程式を作るのが普通である.これも可能であるが,少し手間 を省くために,

  1. これら$ 4N$ 個の未知数を$ x_j$$ y_j$ ,さらに$ x=x_j$ における2次導関数の値を $ u_j$ で表現する.

    $\displaystyle u_j$ $\displaystyle =S_k^{\prime\prime}(x_j)$ (6)
         ただし,  $\displaystyle j=0,1,2,\cdots,N$    
      $\displaystyle \qquad\hspace{15mm}k=j-1,j$    

  2. $ u_j$ が満たす連立方程式を作り,$ u_j$ を解く.
  3. 既知の$ x_j$$ y_j$ と連立方程式により求められた$ u_j$ を用いて,区分多項式の 係数 $ a_j, b_j, c_j, d_j$ を計算する.
というアプローチで問題を解く.ここで,本当に未知数 $ (a_j, b_j, c_j, d_j)$ $ (x_j,y_j,u_j)$ で表現することができるのか,という疑問が湧く者もいるだろう.これは,先の連 立方程式を作る条件を上手に使うことにより可能なのである.また,このような方法は, 問題解決の遠回りをしているように思えるが,以降の説明を見ると実際にはかなり簡潔に なることがわかるので我慢して読んで欲しい.

3.2.1 係数の表現

3.2.1.1 $ b_j$ の表現

これは, $ u_j=S_j^{\prime\prime}(x_j)$ から求めることができ る.式(5)から,

$\displaystyle S_j^{\prime\prime}(x)=6a_j(x-x_j)+2b_j$ (7)

となる.$ x=x_j$ の時,これは$ u_j$ となるので,

$\displaystyle b_j=\frac{u_j}{2}$ (8)

が直ちに導かれる.これで,$ b_j$$ u_j$ で表現できたことになる.$ b_j$ を表現するた めには,$ x_j$$ y_j$ を使っても良かったが,ここではたまたま必要なかったのである.

3.2.1.2 $ a_j$ の表現

これは,2次導関数が区分多項式の境界で等しいという条件から 導くことができる.先ほどは区分多項式の左端$ x_j$ を考えた.次に右端$ x=x_{j+1}$ を考 えることにする.右端の導関数は,

$\displaystyle u_{j+1}=S_j^{\prime\prime}(x_{j+1})= 6a_j(x_{j+1}-x_j)+2b_j\qquad(j=0,1,2,\cdots,N-1)$ (9)

となる.これから,$ a_j$

\begin{equation*}\begin{aligned}a_j&=\frac{u_{j+1}-2b_j}{6(x_{j+1}-x_j)} &=\fr...
...j+1}-u_j}{6(x_{j+1}-x_j)}\qquad(j=0,1,2,\cdots,N-1) \end{aligned}\end{equation*}

と導くことができる.これで,$ a_j$$ x_j$$ y_j$$ u_j$ で表現できたことになる.

3.2.1.3 $ d_j$ の表現

これは,区分多項式は全てのデータ点上を通過するという条件 から導くことができる.まずは,多項式$ S_j(x)$ の左端$ x_j$ を考える.ここでは, $ S_j(x_j)=y_j$ なので,式(5)に代入すると

$\displaystyle d_j=y_j$ (11)

が直ちに導ける.これで,$ d_j$$ y_j$ で表現できたことになる.$ d_j$ の表現には, $ x_j$$ u_j$ はたまたま不要であった.

3.2.1.4 $ c_j$ の表現

これもまた,区分多項式は全てのデータ点上を通過するという条件 から導くことができる.今度は,多項式$ S_j(x)$ の右端$ x_{j+1}$ である.ここでは, $ S_j(x_{j+1})=y_{j+1}$ なので,式(5)に代入すると

$\displaystyle 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}$ (12)

が導ける.式 (8),(10),(11)を用いる と,

\begin{equation*}\begin{aligned}c_j&=\frac{1}{x_{j+1}-x_j}\left[ y_{j+1}-a_j(x_{...
..._{j+1}-x_j} -\frac{1}{6}(x_{j+1}-x_j)(2u_j+u_{j+1}) \end{aligned}\end{equation*}

となる.これで,$ a_j$$ x_j$$ y_j$$ u_j$ で表現できたことになる.

以上で,$ a_j$$ b_j$$ c_j$$ d_j$$ x_j$$ y_j$$ u_j$ で表せたことになる.$ x_j$$ y_j$ はデータ点なので,既知である.したがって,$ u_j$ が分かれば,補間に必要な係 数が全て分かるのである.また,連立方程式の[条件1][条件3]も満 たしている.従って,[条件2]を満たすように$ u_j$ を決めれば良いことになる. すると全ての区分多項式の係数が分かるのである.

3.2.2 連立方程式

先に述べたように$ u_j$ は,1次導関数が境界点で等しいという条件から決める.次の式を 使うのである.

$\displaystyle S^\prime(x_{j+1})=S_j^\prime(x_{j+1})=S_{j+1}^\prime(x_{j+1})\qquad(j=0,1,2,\cdots,N-2)$ (14)

これと式(5)から,

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

となる.あとは,この式の$ a_j$$ b_j$$ c_j$$ x_j$$ y_j$$ u_j$ で 表して,$ u_j$ の連立方程式にするだけである.最終的に式は,

$\displaystyle (x_{j+1}-x_j)u_j+2(x_{j+2}-x_j)u_{j+1}+(x_{j+2}-x_{j+1})u_{j+2}= ...
...\frac{y_{j+2}-y_{j+1}}{x_{j+2}-x_{j+1}} -\frac{y_{j+1}-y_j}{x_{j+1}-x_j}\right]$ (16)

となる.この方程式は, $ j=0,1,2,\cdots,N-2$ で成り立つ.したがって,式の数は, $ N-1$ 個である.$ u_j$ の数は$ N+1$ 個であるが,$ u_0=u_N=0$ としたので,未知の$ u_j$$ N-1$ 個である.式(16)を解くことにより,全ての$ u_j$ が 決定できる.これが決まれば,$ a_j$$ b_j$ ,そして$ c_j$ が計算できる.

$ u_0=u_N=0$ を代入した連立1次方程式は,

\begin{equation*}\begin{aligned}&\left( \begin{array}{@{ }ccccccc@{ }} 2(h_0+h...
... \vdots  v_j  \vdots  v_{N-1} \end{pmatrix} \end{aligned}\end{equation*}

となる.ただし,$ h_j$$ v_j$ は以下のとおり.

$\displaystyle h_j$ $\displaystyle =x_{j+1}-x_j$   $\displaystyle (j=0,1,2,\cdots,N-1)$ (18)
$\displaystyle v_j$ $\displaystyle =6\left[\frac{y_{j+1}-y_j}{h_j} -\frac{y_j-y_{j-1}}{h_{j-1}}\right]$   $\displaystyle (j=1,2,\cdots,N-1)$ (19)

3.2.3 区分多項式の決定

いうまでもないと思うが,スプライン補間を行う区分多項式(5)は, 以下の手順で求める.
  1. 連立方程式(17)を計算して,各点での二次の導関数の値 $ u_j$ を求める.
  2. 区分多項式の係数$ a_j$$ b_j$$ c_j$$ d_j$ を式(10), (8),(13),(11)を 計算する.
区分多項式が分かれば,補間の値計算できる.
ホームページ: Yamamoto's laboratory
著者: 山本昌志
Yamamoto Masashi
2005-12-21


no counter