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)$ (4)

となります。この $ a_j, b_j, c_j, d_j$を決めなくてはなりません。

これらの未知数は、4N個あります。従って、4N個の方程式が必要になります。 そのために、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次導関数の値を指定するものもあります。

これで全ての条件が決まりました。あとは、この条件に満たす連立方程式を求 めるだけです。まずは、2次導関数が区分関数の境界で等しいという条件から です。$ x=x_j$における2次導関数の値を$ u_j$とします。すなわち、

$\displaystyle u_j=S^{\prime\prime}(x_j)\qquad(j=0,1,2,\cdots,N)$ (5)

です。 $ u_j=S_{j-1}^{\prime\prime}(x_j)=S_j^{\prime\prime}(x_j)$としま すので、2次導関数の条件は満足されたことになります。この式から、

$\displaystyle u_j=S_j^{\prime\prime}(x_j)=2b_j\qquad(j=0,1,2,\cdots,N-1)$ (6)

となります。これから、

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

が、直ちに導けます。ここで、スプライン補間の係数、すなわち計算で求める べき$ b_j$$ u_j$で表した理由があります。以降の議論を見ると分かるように、 $ u_j$を連立方程式で計算することにより、他の係数を求めることができます。 そのようなわけで、できるだけ$ u_j$で表現するようにします。

さらに2次導関数の計算から、

$\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)$ (8)

が導けます。この式から、$ 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*}

となります。これで、2次導関数の条件は終わりです。

つぎに、全てのデータ点上を通過する(最初の条件)という条件を考えます。ま ずは、区分の左端の点から、

$\displaystyle d_j=y_j$ (10)

が直ちに導けます。つぎに、区分の右端の点から

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

が導けます。式 (7),(9),(10)を用いると、

\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$$ b_j$$ c_j$$ d_j$$ x_j$$ y_j$$ u_j$で表せました。 $ x_j$$ y_j$はデータ点なので、値はわかっています。したがって、$ u_j$が 分かれば、補間に必要な係数が全て分かります。

3.3 連立方程式

それでは、$ 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)$ (13)

これと式(4)から、

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

となります。あとは、この式の$ 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]$ (15)

となります。この方程式は、 $ j=0,1,2,\cdots,N-2$で成り立ちます。従って、 式の数は、N-1個です。$ u_j$の数はN+1個ですが、$ u_0=u_N=0$ですので、未知 の$ u_j$はN-1個となります。式(15)を解くことに より、全ての$ 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)$ (17)
$\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)$ (18)


ホームページ: Yamamoto's laboratory
著者: 山本昌志
Yamamoto Masashi
平成17年1月14日


no counter