スプライン補間
入力に対する出力の値を決めておくと、その間を滑らかに補間します。これにより、すべての連続関数の近似値を求めることが可能です。
イベント命令の記述例と、相互変換モジュールを用いて貼り付けるためのコードを併記しています。実際に入力する際は、相互変換モジュールを用いてスイッチ・変数番号を任意の番号に書き換えて変換すると効率的です。
概要
y＝f(x)の関数があったとして、入力xに対する出力yがあらかじめデータとしていくつか与えられているとします。それぞれのデータを滑らかに補間し、連続的な関数にするのがこのスプライン補間です。データファイルをあらかじめ作成しておき（例：data.txt）、それをHSP構文で読み込ませてイベント命令を生成します。以下のイベント命令の生成には相互変換モジュールが必須になります。データファイルは1列目にx、2列目にyの値をそれぞれ入力し、それをデータの個数だけ続けます。入力に値を入れて、以下のプログラムにより生成されるイベント命令を実行すると、出力に補間後の値が返ります。
例として以下のdata.txtはy＝100e^(0.01x)のデータをいくつか列挙したものです（x＝25のときy＝128、など）。例えばx＝80を代入するとy＝222が返るので、e^0.8≒2.22であると求まります（真値は2.2255……）。
イベント命令
| データファイル例（data.txt） |
|---|
0 100
25 128
50 164
75 211
100 271
125 349
150 448
175 575
200 738
225 948
250 1218
275 1564
300 2008
325 2579
350 3311
375 4252
400 5459
425 7010
450 9001
475 11558
500 14841
525 19056
550 24469
575 31419
600 40342
625 51801
650 66514
675 85405
700 109663
725 140810
750 180804
775 232157
800 298095
825 382762
850 491476
875 631068
900 810308
921 999660
|
| 使用変数 |
|---|
|
| 相互変換モジュール用コード |
|---|
#include "rpgfunc.as"
/* スイッチ・変数番号を入力 */
v1 = 1 // 変数1:入力x
v2 = 2 // 変数2:出力y
v3 = 3 // 変数3:定数Ａ
v4 = 4 // 変数4:定数Ｂ
v5 = 5 // 変数5:定数Ｃ
v6 = 6 // 変数6:定数Ｄ
v7 = 7 // 変数7:Ａ指数部
v8 = 8 // 変数8:Ｂ指数部
v9 = 9 // 変数9:Ｃ指数部
v10 = 10 // 変数10:Ｄ指数部
fileName = "data.txt" // データファイル名
/* 以下、イベント命令の構文 */
goto *@f
#deffunc adj double di_i
mi = i
i = 0
di = di_i
repeat
if di = 0.0 : i = mi : break
if absf(di) < 100 : di *= 10.0 : i-- : continue
if absf(di) >= 1000 : di *= 0.1 : i++ : continue
break
loop
if int(10.0 * di) \ 10 >= 5 : di += 1.0
return
#deffunc shisuu int i_i
cy
co 1, v7, 1, i_i, 1
cybreak
coend
varsel v7 : var 1, 0, 1, 0
varsel v3 : var 0, 1, v2, 0
varsel v3 : var 5, 0, 10, 0
varsel v2 : var 4, 0, 10, 0
co 1, v3, 0, 5, 1
varsel v2 : var 1, 0, 1, 0
coend
cyend
return
#deffunc keta
cy
co 1, v2, 0, 1000, 4
cybreak
coend
varsel v7 : var 1, 0, 1, 0
varsel v3 : var 0, 1, v2, 0
varsel v3 : var 5, 0, 10, 0
varsel v2 : var 4, 0, 10, 0
co 1, v3, 0, 5, 1
varsel v2 : var 1, 0, 1, 0
coend
cyend
return
*@
data = ""
notesel data
exist fileName
if strsize = -1{
dialog "データファイルが存在しません。", 1
end
}
noteload fileName
n = notemax
ddim x, n + 1 : ddim y, n + 1
ddim v, n - 1 : ddim u, n + 1
// データの読み込み
repeat n
noteget k, cnt
x.cnt = double(strmid(k, 0, instr(k, 0, "\t")))
y.cnt = double(strmid(k, instr(k, 0, "\t") + 1, strlen(k)))
loop
// 連立方程式の準備
repeat n-1
v.(cnt + 1) = 6.0 * ((y.(cnt + 2) - y.(cnt + 1)) / (x.(cnt + 2) - x.(cnt + 1))-(y.(cnt + 1) - y.cnt) / (x.(cnt + 1) - x.cnt))
u.(cnt + 1) = 2.0 * (x.(cnt + 2) - x.cnt)
loop
// ガウスの消去法
repeat n - 2,2
di = 1.0 * (x.cnt - x.(cnt-1)) / u.(cnt-1)
u.cnt -= di * (x.cnt - x.(cnt-1))
v.(cnt) -= di * v.(cnt-1)
loop
// uを求める
repeat n - 1
if cnt > 0 : v.(n - cnt - 1) -= (x.(n - cnt) - x.(n - cnt - 1)) * v.(n - cnt)
v.(n - cnt - 1) /= u.(n - cnt - 1)
u.(n - cnt - 1) = v.(n - cnt - 1)
loop
// 係数の計算
i = 0
repeat n - 1
co 1, v1, 0, x.cnt, 1
co 1, v1, 0, x.(cnt + 1), 4 - 2 * (cnt = n - 2)
varsel v1 : var 2, 0, x.cnt, 0
adj (u.(cnt + 1) - u.cnt) / (6.0 * (x.(cnt + 1) - x.cnt))
varsel v3 : var 0, 0, int(di), 0
varsel v7 : var 0, 0, i, 0
adj u.cnt / 2.0
varsel v4 : var 0, 0, int(di), 0
varsel v8 : var 0, 0, i, 0
adj (y.(cnt + 1) - y.cnt) / (x.(cnt + 1) - x.cnt) - (x.(cnt + 1) - x.cnt) * (2.0 * u.cnt + u.(cnt + 1)) / 6.0
varsel v5 : var 0, 0, int(di), 0
varsel v9 : var 0, 0, i, 0
adj y.cnt
varsel v6 : var 0, 0, int(di), 0
varsel v10 : var 0, 0, i, 0
coend
coend
loop
varsel v2 : var 0, 1, v3, 0 ;Ａ
varsel v2 : var 3, 1, v1, 0 ;（ｘ－ｘj）を乗算
shisuu v8 ;指数部調整
varsel v2 : var 1, 1, v4, 0 ;＋Ｂ
keta ;桁合わせ
varsel v2 : var 3, 1, v1, 0 ;（ｘ－ｘj）を乗算
shisuu v9 ;指数部調整
varsel v2 : var 1, 1, v5, 0 ;＋Ｃ
keta ;桁合わせ
varsel v2 : var 3, 1, v1, 0 ;（ｘ－ｘj）を乗算
shisuu v10 ;指数部調整
varsel v2 : var 1, 1, v6, 0 ;＋Ｄ
;指数部調整（さいご）
cy
co 1, v7, 0, 0, 0
cybreak
coend
co 1, v7, 0, 0, 3
varsel v7 : var 2, 0, 1, 0
varsel v2 : var 3, 0, 10, 0
coend
co 1, v7, 0, 0, 4
varsel v7 : var 1, 0, 1, 0
varsel v3 : var 0, 1, v2, 0
varsel v3 : var 5, 0, 10, 0
varsel v2 : var 4, 0, 10, 0
co 1, v3, 0, 5, 1
varsel v2 : var 1, 0, 1, 0
coend
coend
cyend
send
|
