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

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

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

• 係数は，3重対角行列である．
• 対角成分は，ゼロにならない．

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

連立方程式により計算される2次導関数の値は， $u_1,\/u_2,\/u_3,\/\cdots \/,u_{n-1}$ である．両端は自然境界条件で， $u_0=0\quad u_n=0$ とする．こ この関数では最後にその条件を配列をu[0]=0, u[n]=0とすることで実 現している．

図 1: 2次導関数$u_i$ を求める関数．

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)を呼び出せばよいわけで ある．

図 2: 補間された値を求める関数．

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