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: スプライン補間の区分

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

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

$$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$で表現する．

   $$u_j =S_k^{\prime\prime}(x_j)$$ (6)

   $$j=0,1,2,\cdots,N$$

   $$\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)から，

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

となる．$x=x_j$の時，これは$u_j$となるので，

$$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}$を考 えることにする．右端の導関数は，

$$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{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}$$

と導くことができる．これで，$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)に代入すると

$$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)に代入すると

$$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{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}$$

となる．これで，$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次導関数が境界点で等しいという条件から決める．次の式を 使うのである．

$$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)から，

$$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$の連立方程式にするだけである．最終的に式は，

$$(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{aligned}&\left( \begin{array}{@{\,}ccccccc@{\,}} 2(h_0+h... ...\\ \vdots \\ v_j \\ \vdots \\ v_{N-1} \end{pmatrix} \end{aligned}$$

となる．ただし，$h_j$と$v_j$は以下のとおり．

$$h_j =x_{j+1}-x_j (j=0,1,2,\cdots,N-1)$$ (18)

$$v_j =6\left[\frac{y_{j+1}-y_j}{h_j} -\frac{y_j-y_{j-1}}{h_{j-1}}\right] (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
著者: 山本昌志