1 プログラム方法

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

• 各点での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
著者: 山本昌志