プログラミングと数値計算

スプライン補間（1, 2次元）

7件のコメント

3次スプライン補間に関する簡単な説明とfortran90による計算コードです。

概要

・与えられたデータ点すべてを通るように補間
・与えられたデータ点間を3次関数で近似する
・与えられたデータ点間の隣り合う補間関数（スプライン関数）の一階微分、二階微分が連続
・与えられたデータ点の両端は二階微分がゼロと仮定

1次元の場合

区間\([x_i, x_{i+1}]\)の間を3次関数
\(
S_i(x)=a_i(x-x_i)^3+b_i(x-x_i)^2+c_i(x-x_i)+d_i,~(i=0,1,\cdots N-1)
\)
で補間します。補間は与えられたデータ点間の隣り合う区間の補間関数の一階微分と二階微分が一致するように決めます。
すなわち
\(
\displaystyle \frac{d S_i(x_i)}{dx}=\frac{d S_{i-1}(x_i)}{dx} \\
\displaystyle \frac{d^2 S_i(x_i)}{dx^2}=\frac{d^2 S_{i-1}(x_i)}{dx^2} \\
\)
という式が成り立っているとします。与えられたデータ点の端、点\(x_0\)と\(x_N\)では自然境界条件と呼ばれる境界条件
\(
\displaystyle \frac{d^2 S_0(x_0)}{dx^2}=\frac{d^2 S_{N-1}(x_N)}{dx^2}=0
\)
を課します…この境界条件を課す理由は良く分かりません。[1]に端ではない領域で補間が滑らかになるようにこう決める、とあります。良く分かりません。きっと解く都合上こう決めると簡単に解けるからでしょう。実際そうですし…

詳細は[1]が詳しいのでそちらを参照してください。
係数は以下の通り与えられます。
\(
\begin{align}
a_i&=\frac{v_{i+1}-v_i}{6(x_{i+1}-x_i)} \\
b_i&=\frac{v_i}{2}\\
c_i&=\frac{y_{i+1}-y_i}{x_{i+1}-x_i}-\frac{1}{6}(x_{i+1}-x_i)(2v_i+v_{i+1}) \\
d_i&=y_i
\end{align}
\)
ここで、\(v_i\)は以下の方程式の解です。
\(
\scriptsize
\begin{eqnarray}
\left(
\begin{array}{ccccc}
2(h_0+h_1)& h_1 & 0 & \cdots & 0 \\
h_1 & 2(h_1+h_2) & h_2 & \cdots & 0\\
0 & h_2 & 2(h_2+h_3) & \cdots & 0\\
\ldots &&&& \\
0& \cdots & 0 & h_{N-2} & 2(h_{N-2}+h_{N-1})
\end{array}
\right)
\left(
\begin{array}{c}
v_1 \\
v_2 \\
v_3 \\
\ldots \\
v_{N-1}
\end{array}
\right)=\left(
\begin{array}{c}
w_1 \\
w_2 \\
w_3 \\
\ldots \\
w_{N-1}
\end{array}
\right)
\end{eqnarray}
\normalsize
\)
ここで、
\(
\begin{align}
h_i&=x_{i+1}-x_i,~ (i=0,1,\cdots,N-1)\\
w_i&=6\left(\frac{y_{i+1}-y_{i}}{h_i}-\frac{y_i-y_{i-1}}{h_{i-1}}\right),~ (i=1,2,\cdots,N-1)
\end{align}
\)
です。

計算コード

Fortranコードを以下に置いておきます。
注意しておきますが、[1]を元にしたコードではありません。[2]を元にしたコードです。

・配列 fdata(0:N),xdata(0:N) は与えられたデータで、既知の値です。
・ c3spline1p は任意の点 x の補間値 f 、一階微分 df 、二階微分 df2 を fdata(0:N),xdata(0:N) を元に計算します。
・ c3spline は任意の点配列 xa(0:M) の補間値 fa(0:M) を fdata(0:N),xdata(0:N) を元に計算します。
・ c3spline_integral はスプライン補間された関数の積分 s を計算します。

c3spline1p を用いて複数の点を計算することもできますが、c3spline に比べて遅いです。
計算時間は  c3spline1p : c3spline = 1 : 0.03です。

2次元の補間

2次元の補間を考えます…が、少し横着します。
「横着」は2次元のスプラインを調べるのが面倒で、思いついてしまったのでこれでいいだろう、という意味です。

格子状に与えられたデータ列\(f(x_i,y_j), i=0,\cdots,N_x, j=0,\cdots,N_y\)があり、点\((x,y)\)での補間された値が欲しいとします。
点\((x,y)\)は、それぞれ
\(
\begin{align}
& x_i\le x \lt x_{i+1} \\
& y_j\le y \lt y_{j+1}
\end{align}
\)
に存在するとします。
この時、考え方は以下の通りです。

まず、\(j\)を固定して\(x\)上の補間値を\(N_y\)個(点\(f(x,y_0),f(x,y_1),f(x,y_2),\cdots,f(x,y_{N_y})\))を得ます。
その後、補間された\(N_y\)個の補間値に対して補間を行い、点\((x,y)\)に補間値\(f(x,y)\)を得ることが出来ます。

ここで問題なのが、\(x\)方向、\(y\)方向のどちらを先に補間するのか?
そしてその2つの結果は変わるのか?ということです。

直感的に考えて良さそうなのは２つの平均を取ってしまう事でしょう。しかし、単純に考えて計算時間が2倍になります。
直ぐには分からなそうなので、プログラムではそれを指定できる形にする、にとどめましょう。

偏微分を求められます。最後に補間した方向の偏微分を得ることが出来ます。

c3spline2d1pは元となるデータ配列を入れると、点(x,y)の補間された値を返します。
c3spline2dは元となるデータ配列を入れると、配列で指定された格子状の点を返します。

上記コードを動かすと、以下の画像を作ることが出来ます。

赤い点は補間に用いたオリジナルのデータ列を表し、黒の線は補間された関数を表します。また、底面のカラーマップは本当の関数との相対誤差を表しています。xが大きいところの端で補間した結果の関数が本当の関数とあまり一致していません。

これは、元にするデータ点の両端において、元の関数の二階微分がゼロという仮定が満たされていないためです。

上の図を再現するgnuplotのコードは以下のものです。

set palette defined ( -1 '#ffffff', 0 '#000090',1 '#000fff',2 '#0090ff',3 '#0fffee',4 '#90ff70',5 '#ffee00',6 '#ff7000',7 '#ee0000',8 '#7f0000')

set ticslevel 0
set view 46,40,1,1

set xl "x"
set yl "y"
set cbr[0:1]
unset key
splot "interpolate_data1.d" u 1:2:6 w pm3d at b
replot "interpolate_data1.d" u 1:2:3 w l lw 1
replot "original_data.d" u 1:2:3 w p pt 7 ps 0.5

参考文献

[1]https://www.mk-mode.com/rails/docs/INTERPOLATION_SPLINE.pdf

[2]http://150.19.250.16/MULTIMEDIA/numeanal/node16.html

アイキャッチ画像のフォント：キルゴhttp://getsuren.com/killgoU.html

関連記事

このエントリーをはてなブックマークに追加

「スプライン補間（1, 2次元）」への7件のフィードバック

  1. 2次元スプライン補間の導出を参考にさせていただきました．
    記事でデータ点の端の二階微分値を0にする理由がわからないとおっしゃっていましたが，
    これは材料力学の梁(はり*)のたわみ計算に近いものと理解しました．
    梁のたわみ計算にもスプライン曲線と似た考え方が登場するのですが，
    このとき，二階微分値は梁にかかるモーメントに相当します．
    そのため，梁の両端にモーメントをかけると不自然に曲がってしまうように，
    スプライン曲線もデータ点の端の二階微分値を0にしないと
    不自然な形になってしまうのではないかと考えます．
    *あるいは宙に浮かんだ棒を考えてください．

    ご参考になれば幸いです．

    返信

    1. コメントありがとうございます。
      たしかに、物理的な視点を加えればそのような解釈が可能ですね。
      スプライン補間をある特定の物理的な現象（二階微分がゼロである境界条件を満たす現象）に施す場合ならば、二階微分がゼロで補間する方法は妥当な境界条件であり、その通りだと思いました。

      しかし、
      　１．スプライン補間の考えそれ自体に物理現象は関係ない点
      　２．例えば周期境界条件を満たすような物理現象を記述する関数を補間する際には適切ではない点
      が説明できません。
      にもかかわらず、よく見るスプライン補間の説明では、データ点の端の二階微分値を0にするというのが通例となっています。
      その意味でなぜ二階微分がゼロの境界条件が採用されているのか分からない、と記述しました。
      二階微分値を0に置く物理現象が多いので採用しているのかな、と思っています。

      返信

  2. 大変分かりやすく素晴らしい記事でした、ありがとうございます。

    境界条件に関してですが、この記事で紹介されているものの他にも、様々なパターンがあるようです。両端で2階微分が0という条件を課したものは特にnatural splineやsimple splineと呼ばれているようです。他にも両端の1階微分の係数を与えるclamped splineや、周期境界条件を課すようなperiodic splineなど色々あるようです。

    同じデータでも、この境界条件の取り方でスプライン曲線は変わりますので、状況に応じて使い分ける必要がありますね。色々比較してみると興味深いかと思います。

    返信

    1. コメントありがとうございます。

      なるほど！キーワードで検索したところ、たくさん出てきました。私の調べが足りなかったようで、ご指摘ありがとうございます。
      固定端だと数値計算的に三重対角行列になるので、最も簡単に解けるから教育的にも良いのかな？と考えていました。
      特に端が重要になる場合は使い分ける必要がありますね。

      返信

  3. プログラム勉強中に出会い、素晴らしいサイトに感動しております。
    subroutine spline_abcd中のl, mu, zの部分はどのような仕組み（アルゴリズム）なのでしょうか？大変恐縮なのですが、お教え頂ければ幸いです。

    返信

    1. この部分は、三重対角行列の解法で『TDMA(Tri-Diagonal Matrix Algorithm)』または『Thomasのアルゴリズム』と呼ばれています。
      googleの検索で、『三重対角行列　TDMA』と調べれば出てくるかと思います。

      この記事を書いたのが少し前なのであまり覚えていませんが、私もどこかの擬コードを元に作成したか、
      今回の場合に特化させている段階で途中の値をl,mu,zという変数で置いていたか、どちらかと思います。

      私のページでは、クランク・ニコルソン法（https://slpr.sakura.ne.jp/qp/tdse-crank-nicolson/）で三重対角行列を使っています。
      その中に以下の似たコードがありました。

        u=b(1)
        x(1)=r(1)/u
        do j=2,n
           w(j)=c(j-1)/u
           u=b(j)-a(j)*w(j)
           x(j)=(r(j)-a(j)*x(j-1))/u
        enddo

      spline_ancdの中のzが上のコードのxに対応しているような感じです。

      返信

  4. 早速のお返事ありがとうございます。TDMAを調べ簡単なプログラムを組んで勉強させていただきました。Sikinote様のwebsiteは非常に勉強になります。このようなサイトを作ってくださりありがとうございます。

    返信

コメントを残す

メールアドレスが公開されることはありません。 ※ が付いている欄は必須項目です