等加速度系における、衝突までのステップ数計算

自機Sと、敵弾Bの衝突までの時間を求めたい。どちらも、座標と一定の半径(R_S, R_B)と一定速度or一定加速度を持っているものとする。

S・Bが等加速度運動・等速運動・静止のどれかの運動をしていて、少なくとも一方が等加速度運動をしている場合、相対座標は時間tの2次式で書かれる。ユークリッド距離空間では、相対座標の2乗和<=(衝突半径の和)^2なので、4次不等式を解く事になる。

つまり、
$\sum_{k\in K} (a_k\frac{t(t+1)}{2}+v_k t+x_k)^2 \leq (R_S+R_B)^2$
かつt>=0を満たす最小解を求める。(精度は整数レベル)
精度が整数レベル、というのは、1ステップ(自機や敵弾の座標に速度を足す)を時間で1としているため。整数解で十分なのだ。
ここでKは2次元だと{x,y}, 3次元だと{x,y,z}といった感じ。t(t+1)/2というのは、速度に加速度を足した後に座標に速度を足す操作をt回繰り返したと言う意味でこうなっている。

さて定式化したのに申し訳ないのだが、上で色々書いた変数とはお別れである。つまるところt=0での(左辺)-(右辺)の正負を調べた後、正であれば4次方程式のt>=0の最小解を求め、そうでなければ、t=0が解となる。よって4次方程式を高速に解ければよい。

最初、ベアストウ法を使ってみていた。これは、任意の多次方程式の虚数解も含んだすべての解を数値計算する(近似方法はニュートン法)優れものであるが、実際に動かしてみると、特殊な条件下では収束までにえらくステップ数がかかることがわかった。

下左側は、方程式 $0.01x^4-0.7x^3-65.2x^2+2703x+187440.9847=0$ をWikipediaに書いてある通りに実装して解いた場合の、4次の項での各ステップの(p,q)の座標である。p,qは、2次方程式でいうx^2+px+q=0の部分に相当していて、これを収束させればよい。
下右側は、下左側で見える直線上での動きをみたもの。収束点にたまたま近づいたから収束したという感じにもみえる。

結論から書くと、まず初期値の設定が問題で、 $ax^4+bx^3+cx^2+dx+e=0$ といった方程式では、通常は上位の3項を2次方程式とみなして初期値をはじき出す。ここで、eの絶対値がa,b,cに比べてかなり大きいと、初期値からかなり離れたところに解ができてしまい、収束しにくい。
また、逐次法なので、1ステップごとに進む幅は実装者任せである。細かく取りすぎると近づくのが遅くなり、大きく取りすぎると、近くなる点を通り過ぎて逆に遠くなってしまう場合が出る。

ということで、これはいかんといろいろ調べていたら、Strumの定理というものを見つけた。これは、任意の多次方程式の実数開区間(a,b)内の解の個数を求める。もちろん実数解しかわからず、重解も1個と数えてしまうが、今回の目的には特に影響はない。
これを、適当なTを決めて(0,T)で個数を求めて、1個以上であれば、(0,T/2)で個数を求めて、なければ(T/2,T)のほうにあると判断して(T/2,3T/4)を探し、あればさらに(0,T/4)で・・というふうに2分探索すればいいのではないかと考えた。実装してみると結構よさげに動いているようだ。

4次程度ではあんまり効果がないかもしれないが、(0,T)でM個解があったとき、それらの解は一様分布していると仮定して、2分探索よりももっと小さい範囲を次に指定して、効率のいい探索ができないかなーと考えたりしている。