2 変数関数のスプライン関数による取り扱いについて
中野武雄 nakano@apm.seikei.ac.jp

2010-1-27
概 要
スパッタリング法のモンテカルロシミュレーションでは、粒子の散乱角をエネルギーと衝突パラメータの関数として求める。本文書は、この取り扱いをスプライン関数で行うプログラムルーチン作成に当たってのメモである。

1B スプライン
B スプラインとは、スプライン関数の基底関数 (basis) を構成するものである。 n 次の B スプラインは n + 2 個のノード間でのみ値を持ち、その他の領域では 0 になるような関数である。また、有限領域の両端では n − 1 次の微係数が 0 という性質を持つ。
この B スプラインは、切断べき級数の差分商によって与えられる。 M 次の切断べき級数は{0(x < qi )(1)(x − qi )M+ =M(x − qi ) (x ≥ qi )と定義される関数である。
一方差分商は以下のように定義される。いま関数 f (x) が異なる座標 qi , qi+1 , . . . , qi+N に対して定義されているとする。

1. 関数 f (x) に関する 1 階差分商 f (qi , qi+1 ) を以下のように定義する。
f (qi , qi+1 ) =f (qi+1 ) − f (qi )qi+1 − qi(2)

すなわちこれは yi = f (qi ) としたときのデータ点 (qi , yi ) と (qi+1 , yi+1 ) を結ぶ直線の傾きである。

2. r + 1 階の差分商 f (qi , qi+1 , . . . , qi+r+1 ) は r 階の差分商から以下の漸化式で求められる。
f (qi , qi+1 , . . . , qi+r+1 ) =f (qi+1 , qi+2 , . . . , qi + r + 1) − f (qi , qi+1 , . . . , qi + r)
qi+r+1 − qi(3)

これらを用いて、K − 1 次の B スプライン Bi,k (x) を以下のように定義する。

Bi,k (x) = (qi+K − qi )MK (x : qi , qi+1 , . . . , qi+K )(4)

ここで MK (x : qi , qi+1 , . . . , qi+K ) はMK (x : q) = (q − x)K−1+1(5)に対して qi , qi+1 , . . . , qi+K の接点で K 階の差分商を取ったものである。このとき各接点は qi ≤qi+1 ≤ · · · ≤ qi+K を満たしているものとする。
すると 0 次のスプライン関数 Bi,1 (x) は{Bi,1 (x) =(qi ≤ x < qi+1 )(x < qi , qi+1 ≤ x)

10(6)となり、 K − 1 次の B スプラインは漸化式

Bi,K (x) =x − qiqi+K − x
Bi,K−1 (x) +Bi+1,K−1 (x)qi+K−1 − qi
qi+K − qi+1(7)と表すことができる。

21 次元スプライン関数の導出
前節で定義された B スプライン関数を用いて、 K-1 次のスプライン関数を構成する。スプライン関数 s(x) はs(x) =N−1∑

αi Bi,K (x)

(8)

i=0

で定義され、与えられた N 個の点 (x0 , y0 ), (x1 , y1 ), . . . , (xN −1 , yN −1 ) を通るものとする。このとき

Bi,K が qi , . . . , qi+K の K 個のデータ点から構成されるから、式 8 を構成するためには N + K 個の接点が必要となる。さらに N 個の各データ点を s(x) が通るようにするには、 B スプラインはBi,K (xi ) 6= 0,(i = 0, 1, . . . , N − 1)(9)を満たす必要がある。各 B スプラインを構成する接点から式 9 を書き換えると

q0 < x0 < q0+K
q1 < x1 < q1+K..
.

(10)

qN −1 < xN −1 < qN −1+Kとなる。この条件は Shoenberg-Whittny の条件と呼ばれる。通常 K − 1 次のスプライン関数の接点は、簡単のために以下のように取ることが多い。

q0qi+KqN= q1 = · · · = qK−1 = x0(i = 0, 1, . . . , N − 1 − K)

= (xi + xi+K )/2= qN +1 = · · · = xN +K−1 = xN −1(11)接点が決定できたら、連立 N 次方程式
N−1∑

αi Bi,K (xj ) = yj (j = 0, 1, . . . , N − 1)i=0を解いて αi を求めれば良い。

2(12)2 次 (K=3) のスプライン関数における実作業

2.1

前提として、xi (i = 0 . . . N − 1) は等間隔であり、各々に対する yi = f (xi ) は既知とする。
まず qi (i = 0 . . . N + 2) を以下のように定める。

q0= q1 = q2 = x0

qi+3

= (xi + xi+3 )/2qN= qN +1 = qN +2 = xN −1(i = 0 . . . N − 4)

N × N の行列 B を用意し、以下のように Bij を決定する。
# 1st kindfor(j=0;j<N;j++){for(i=0;i<N;i++){if (x[j].ge.q[i] .and. x[j].le.q[i+1]){
B[i][j] = 1}else{
B[i][j] = 0}}}# 2nd kindfor(j=0;j<N;j++){for(i=0;i<N;i++){if (q[i+1].eq.q[i]){BP = B[i][j]}else{BP = (x[j]-q[i])/(q[i+1]-q[i])*B[i][j]}if (i+1.ge.N) {BL = 1}else{
BL = B[i+1][j]}if (q[i+2].eq.q[i+1]){
BL = BL}else{
BL = (q[i+2]-x[j]/(q[i+2]-q[i+1])*BL

3(13)

}
B[i][j] = BP + BL}}# 2nd kindfor(j=0;j<N;j++){for(i=0;i<N;i++){if (q[i+2].eq.q[i]){BP = B[i][j]}else{BP = (x[j]-q[i])/(q[i+2]-q[i])*B[i][j]}if (i+1.ge.N) {BL = 1}else{
BL = B[i+1][j]}if (q[i+3].eq.q[i+1]){
BL = BL}else{
BL = (q[i+3]-x[j]/(q[i+3]-q[i+1])*BL}
B[i][j] = BP + BL}}
このようにして Bij が求まったら、−1
αi = Bjiyj(14)

によって係数ベクトル α を求めれば良い。行列 B は疎 (sparse) な行列となるので、 Stanford 大学で開発された SPARSE ライブラリを用いるのがメモリの有効利用などの点で有利であろうと考えている（まだ実装はしていない）。

32 次元スプライン関数の導出

3.1

一般的な議論

モンテカルロシミュレーションでは、• 散乱角が衝突のエネルギー E と衝突パラメータ b の関数となる。
4

• 衝突ガス粒子の速度、角度が乱数 p とスパッタ粒子の速度 vg の関数となる。
これら二変数の適当な間隔ごとに関数値を求め、この (Ei , bj , θij )(0 ≤ i, j ≤ N −1) や (pi , vpi , vgi )を元にスプライン関数を構成することを考える。以降では変数は (x, y, z(x, y)) に置き換える。
まず x 方向の接点 qi と y 方向の接点 rj を式 11 に従って求める。スプライン関数 s(x,y) はs(x, y) =−1
N−1 N
∑
∑

αij Bi,K (x)Bj,K (y)(15)

i=0 j=0と与えられるから、これを前節と同様の方法で解けば良い。

3.22 次のスプライン関数 (K = 3) の実装

Berkley SPARSE ライブラリではなく、通常の行列解法を用いたバージョンを作成した。サブルーチン bd でマトリックス b を作り、これをベクトル z(nx*ny) に対して ba = z のように解いて係数ベクトル a を求める。関数値は spf によって得る。
これを 2 変数関数z = cos(x) sin(y)に対して使ってみた結果を添付図に示す。与えたデータ点は x, y ともに 0.5 刻みとした。

subroutine bd(mdim, m, ndimx, nx, ndimy, ny, ichix, ichiy,$ccx, y, qx, qy, b, bx, by, bxj, byj)

create the matrix of N x N dimension for the determination ofccc

parameters of spline function order ICHI using Deboor Cox argorism

cc(i.e. x[j] = ja)

ccc

input
mdimmassumption.. the data point is placed periodically.. declared row dimension of B in calling program.. order of the matrix (namely, nx * ny)

ccndimx .. declared row demension for xnx.. number of data points along x axis

cccndimy .. declared row demension for yny.. number of data points along y axis
ichix .. order of spline function for x axis

ccichiy ..
x(n) .. data points along x axis

cccy(n)cc.. data points along y axis

outputqx, qy .. node datab.. created matrix to be solved
5y(16)

ccworkccbx
by.. B spline for X axis..
Ycccbxj
byj.. B spline for X axis.. B spline for Y axis

integer mdim, m, ndimx, nx, ndimy, ny, ichix, ichiyreal b(mdim, m), x(nx), y(ny), qx(nx+ichix), qy(nx+ichiy)real bx(ndimx, nx), by(ndimy, ny)
real bxj(nx), byj(ny)===real function spf (nx, ny, ichix, ichiy, qx, qy, a, x, y)cc
cc

calculate the value of 2D spline function at (x, y).
input

ccnx, ny
ichix, ichiy.. number of data points.. order of spline functioncccqx, qy
ax, y.. node point.. parameter vector for spline function.. input datacinteger nx, ny, ichix, ichiyreal qx(nx+ichix), qy(ny+ichiy)real a(nx*ny), x, y

参考文献[1] 桜井 明 (監), 菅野敬祐, 吉村和美, 高山文雄「C によるスプライン関数」東京電機大学出版局(1993)

6