stormworksのLuaで二分法する話

始めに

ストームワークスでluaをやってるとたまに出てくる二分法のお話。
比較的プログラミング初級者向けの内容です。
そしてこの記事はこちらのアドベントカレンダーの記事です。

プログラムにおける二分法とは

二分法というのは、プログラムで数値解析を行う求根アルゴリズムのひとつです。
ひとつ例を挙げると、平方根があります。

luaなら math.sqrt(x) で求まりますが、それを使わずに √x を求めてください。

というときに使われる手法の1つが二分法です。
一番簡単かつ安定している求根アルゴリズムです。
あと文字数が少ないというメリットがあります。（ストワでは大事）

二分法の簡単な考え方

たとえば √2 なら二乗が2になる数値を二分しながら探します。
最初の範囲は1から2、中間は1.5です。
二乗が2になる数値（1.41421356）は1から1.5の間にあるので、次の範囲を1から1.5にします。
範囲1から1.5の中間は1.25になり、次の範囲は1.25から1.5になります。
その次は1.3625…と範囲の指定と中間の数値を繰り返して求めていけば行けばそのうち1.41421356に近づいて行きます。
これが二分法の原理です。
実装するうえで、 f(x)=0 の方が都合かいいのでそれで実装します。

二分法のプログラム

--a,b範囲
--maxTime最大繰返し回数
--eps許容誤差
--f指定する関数

function bisection(a, b, eps, f,maxTime)
	local i, s = 0, 0
	while math.abs(a - b) > eps and i<maxTime do
		i = i + 1
		s = (a + b) / 2
		if (f(a) * f(s) < 0) then --f(a)*f(s)<0のとき、その間にf(x)=0となるxがある
			b = s
		else
			a = s
		end
	end
    if i>=maxTime then s=0 end --error
	return s
end

これが二分法の関数です。
実際簡単に実装可能で、for文で実装する方法もあります。

ただし使用する場合条件があります。

• 範囲内に確実に0となる f(x) が ひとつだけ 存在すること
• f(x) は指定範囲内で 単調増加または単調減少である こと

根が二つ以上存在する場合は範囲を分割する必要があります。
許容誤差はepsです、ストワなら0.01くらいで充分だと思われます。
また間違った関数と範囲を指定すると無限ループになってしまうので最大繰り返し数(maxTime)を設けています。
luaは引数に関数を渡せるので、二分法を呼ぶときに f(x) を指定します

ストワで活用

いわゆる弾道計算器に使います。
弾道計算に関してはふるにん氏のこちらの動画がわかりやすいです。

以下のようなコードを使うことで距離500mにLACを射撃したときの時間 ballistitTime が求まります。

ballisticData={}
targetData={}
gunData = {
    --{v0:Muzzle Velocity, k:Drag, lifetime }
	{ 800,  1.5,  5}, -- MG
	{ 1000, 1.2,  5}, -- LAC
	{ 1000, 0.6,  5}, -- RAC
	{ 900,  0.3,  10} -- HAC }
}

function gunf(t)
	local E=1-math.exp(-ballisticData.k*t)
    
   return targetData.w^2+(targetData.h+ballisticData.g*t/ballisticData.k-ballisticData.g/ballisticData.k^2*E)^2-(ballisticData.v0/ballisticData.k)^2*E^2
end

function bisection(a, b, eps, f,maxTime)
	local i, s = 0, 0
	while math.abs(a - b) > eps and i<maxTime do
		i = i + 1
		s = (a + b) / 2
		if (f(a) * f(s) < 0) then
			b = s
		else
			a = s
		end
	end
    if i>=maxTime then s=0 end --error
	return s
end

function onTick()
    gunIndex=1--LAC
    ballisticData={g=30,v0=gunData[gunIndex][1],k=gunData[gunIndex][2],lifeTime=gunData[gunIndex][3]}
    --g:gravitational acceleration　30m/s^2
    targetData={w=500,h=0}--w:目標までの水平距離,h：目標までの垂直距離
    ballisticTime=bisection(0,ballisticData.lifeTime,0.01,gunf,50)

    debug.log(ballisticTime)
end

同様のコードを使用してこちらのマイコンは製作しています

弾道計算機は命中までの時間が一番重要で、必ず繰り返し処理が必要となる部分です。
その式が gunf(x) 関数です。
これにより命中時間が求まり、命中時間が求まると仰角や方角も計算出来るようになります。

以上、ストワで二分法を使う話でした