スプライン補間のプログラム

山本昌志^{∗}

2005 年 1 月 18 日

1 プログラム方法

理論から分かるように、スプライン補間を行うためには以下の 2 つの計算が必要である。

• 各点での 2 次導関数の値 u_{i} の計算

• 係数 a, b, c, d を求めて、補間の値の計算

それぞれ、 C 言語の関数を作成することにする。それぞれの関数の作成方法について、以下に述べる。

1.1 2 次導関数の値を計算する関数

まず最初の関数で、標本点 x_{i} と y_{i} から、 x_{i} での 2 次導関数の値 u_{i} を計算するプログラムである。この関数は、最初に 1 回呼び出すだけでよく、結果は補間値を計算する次の関数に渡す。この 2 次導関数の値を計算するプログラムは、独立した関数としてモジュール化するのが常套手段である。そのほうが、プログラムが分かりやすくなり、メンテナンスも容易である。

C 言語の関数でこれを実現するための関数のプロトタイプ宣言は、

     void spline_cal_u(int n, double x[], double y[], double u[]);

のようにすればよい。 n は標本点の最大の番号である。標本点の番号は、 0 から始まり n で終わる。配列 x[] と y[] に標本の x と y 座標の値を入れる。配列 u[] は計算された 2 次導関数の値が入る。

この関数での処理の内容は、

• 標本点での 2 次導関数の値 u_{i} が満たす連立方程式を作成する。

• その連立方程式を解く。である。要するに連立方程式を作り、そして解いているのである。解くべき連立方程式は、次の通りである。

∗独立行政法人 秋田工業高等専門学校 電気工学科

( \begin{matrix}2(h_{0} + h_{1}) & h_{1} \\ h_{1} & 2(h_{1} + h_{2}) & h_{2} & & & 0 \\ & h_{2} & 2(h_{2} + h_{3}) & h_{3} \\ ^{.} & ^{.} & . \\ & & & h_{j−1} & 2(h_{j−1} + h_{j}) & h_{j} \\ & 0 & & & & ^{.}^{.}. \\ & & & & & h_{N−2} & 2(h_{N−2} + h_{N−1}) \end{matrix} ) ( u_{1} ; u_{2} ; u_{3} ; ⋮ ; u_{j} ; ⋮ ; u_{N−1})_{= ( v_{1} ; v_{2} ; v_{3} ; ⋮ ; v_{j} ; ⋮ ; v_{N−1})}

(1)

ただし、 h_{j} と v_{j} は以下のとおり。

h_{j} = x_{j+1} − x_{j}  (j = 0, 1, 2, ⋯, N − 1)  (2)
v_{j} = 6 [\frac{y_{j+1} − y_{j}}{h_{j}} − \frac{y_{j} − y_{j−1}}{h_{j−1}} ]  (j = 1, 2, ⋯, N − 1)  (3)

この関数のフローチャートを図 1 に示す。この関数では、連立方程式を解く必要がある。この連立方程式は、以下の特徴がある。

• 係数は、 3 重対角行列である。

• 対角成分は、ゼロにならない。

通常は LU 分解、あるいは反復法で解くことになる。しかし、そのプログラムを書くのも面倒なので、以前作成したピボット選択のないガウス・ジョルダンを使えばよい。連立方程式により解かれた 2 次導関数の値は、配列 u[] に入れられて呼び出し側に戻る。

連立方程式により計算される 2 次導関数の値は、 u_{1}, u_{2}, u_{3}, ⋯, u_{n−1} である。両端は自然境界条件で、 u_{0} = 0 u_{n} = 0 とする。ここの関数では最後にその条件を配列を u[0]=0, u[n]=0 とすることで実現している。

1.2 補間の値を計算する関数

2 次導関数の値が分かったので、残りの処理は、その値を利用して任意の x の補関値を求めることである。スプライン補間は、 3 次の区分多項式で標本点の間を補完するのは今まで述べたとおりである。そのためには、標本点の座標と 2 次導関数の値、即ち x_{i}, y_{i}, u_{i} がわかれば計算できる。 2 次導関数の値 u_{i} は予め、

先の関数 spline cal u(n, x, y, u) で計算しておくものとする。

スプライン補間の計算に必要な x_{i}, y_{i}, u_{i} から、補間点 x の値を求める関数を作ろう。 C 言語の関数でこれを実現するための関数のプロトタイプ宣言は、

     double spline(int n, double x[], double y[], double u[], double xx);

のようにすればよい。この関数の実引数 xx に補関値を求めたい x を入れるのである。そして、この関数の戻り値が x での補間値 y となる。

この関数の処理の内容は、

• 補間値を求める x が、どの標本点の間にあるか探す。

• 見つかった標本点の場所から、 3 次関数の係数を計算する。

• 3 次関数から x の時の値、 y を計算する。

である。ここで重要なことは、 x がある区間を探し、 3 次関数の係数を求めることである。 2 次導関数と標本点の値から、それらの係数は

a = \frac{u_{j+1} − u_{j}}{6(x_{j+1} − x_{j})}  (4)
b = \frac{u_{j}}{2}  (5)
c = \frac{y_{j+1} − y_{i}}{x_{j+1} − x_{j}} − \frac{1}{6}(x_{j+1} − x_{j})(2u_{j} + u_{j+1})  (6)
d = y_{j}  (7)

と計算できる。これらの係数から、補間の値 y は

y = ax^{3} + bx^{2} + cx + d  (8)

となる。

この関数のフローチャートを図 2 に示す。まず初めに、 x が存在する区間を 2 分探索により探している。これがわかれば、 x を挟む場所の x, y, u の値がわかる。そうすると、式 (4)〜(7) を用いて、 3 次関数の補間式の係数 a, b, c, d の値が計算できる。そして、式 (8) から補間値 y がわかる。補間したい点 x が変わる毎に

spline(n, x, y, u, xx) を呼び出せばよいわけである。

図 1: 2 次導関数の u_{i} を求める関数。 4 図 2: 補間された値を求める関数。

5