Subsections

1 プログラム方法

理論から分かるように,スプライン補間を行うためには以下の2つの計算が必要である. それぞれ,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次導関数の値が入る.

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

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

\begin{equation*}\begin{aligned}&\left( \begin{array}{@{ }ccccccc@{ }} 2(h_0+h...
...\ \vdots \ v_j \ \vdots \ v_{N-1} \end{pmatrix} \end{aligned}\end{equation*}

ただし,$ h_j$$ v_j$ は以下のとおり.

$\displaystyle h_j$ $\displaystyle =x_{j+1}-x_j$   $\displaystyle (j=0,1,2,\cdots,N-1)$ (2)
$\displaystyle v_j$ $\displaystyle =6\left[\frac{y_{j+1}-y_j}{h_j} -\frac{y_j-y_{j-1}}{h_{j-1}}\right]$   $\displaystyle (j=1,2,\cdots,N-1)$ (3)

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

通常は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$ を求める関数.
\includegraphics[keepaspectratio,scale=0.7]{figure/spline_cal_u.eps}

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次関数の係数を求めることである. 2次導関数と標本点の値から,それらの係数は

$\displaystyle a$ $\displaystyle =\frac{u_{j+1}-u_j}{6(x_{j+1}-x_j)}$ (4)
$\displaystyle b$ $\displaystyle =\frac{u_j}{2}$ (5)
$\displaystyle c$ $\displaystyle =\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)
$\displaystyle d$ $\displaystyle =y_j$ (7)

と計算できる.これらの係数から,補間の値$ y$

$\displaystyle 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: 補間された値を求める関数.
\includegraphics[keepaspectratio,scale=0.7]{figure/spline.eps}

ホームページ: Yamamoto's laboratory
著者: 山本昌志
Yamamoto Masashi
2005-12-21


no counter