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

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

これらの未知数は、4N個あります。従って、4N個の方程式が必要になります。 そのために、3次のスプライン補間に以下の条件を課します。

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

以上の条件を課すと方程式は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$とします。すなわち、

$$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次導関数の条件は満足されたことになります。この式から、

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

となります。これから、

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

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

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

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

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

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

$$d_j=y_j$$ (10)

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

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

3.3 連立方程式

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

これと式(4)から、

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

$$(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{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)$$ (17)

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

ホームページ: Yamamoto's laboratory
著者: 山本昌志