6.3追跡曲線 - 驚異のアニヲタ社会復帰の予備

MGSでスネークがスティンガーミサイルを打って、リキッドやらソリダスやらメタルギアと闘うのだが、発射したミサイルがどうやって目標物にあたるか、その弾道と捕捉時間を求めてみよう。

$t=0$ においてミサイルは原点 $(0,\hspace{3}0)$ から、点 $(a,\hspace{3}b)$ にいる飛行機（カサッカ）に向かって飛ぶとする。MGS2でカサッカに乗っていたのはオタコンじゃね？とかそこ突っ込まない！
カサッカは定速度 $v_A$ で、x軸に平行に飛ぶと仮定する。
ミサイルは定速度 $v_m$ を持つとして、その座標を $(x_m,\hspace{3}y_m)$ とする。
ミサイルは点 $M$ で、カサッカに向かって最短距離で飛ぼうとする。つまり、点 $M$ でのミサイル軌道の接線方向に、カサッカが常にいるものと考えられる。この接戦の方程式は、
$\begin{align}y-y_m&=\frac{dy_m}{dx_m}(x-x_m)\\&=\frac{dy_m}{dt}\frac{dt}{dx_m}(x-x_m)\end{align}$
である。カサッカの座標 $A(x_a,\hspace{3}b)$ はこの直線上にあり、 $x_A=a+v_A t$ だから
$b-y_m=\frac{\dot{y}_m}{\dot{x}_m}(a+v_A t -x_m)$
と、速度の関係式
$v_m ^2=\dot x_m^2+\dot y_m^2$
を満たす。
ミサイルの軌道は、境界条件 $x_m(0)=y_m(0)=0$ を満たす、
$\dot x_m(b-y)=\dot y_m(a+v_A t-x_m)$ …(4)
$\dot x_m^2+\dot y_m^2=v_m ^2$ …(5)
を解けばよい。mは煩雑なので今後省略するとして、(4)は
$\frac{dx}{dy}(b-y)=a+v_At-x$
とかけるので、これを $t$ で微分すると
$\frac{d^2x}{dy^2}\frac{dy}{dt}(b-y)-\frac{dx}{dy}\frac{dy}{dt}=v_A-\frac{dx}{dt}$
となる。一方、(5)から
$\frac{dy}{dt}=\frac{v_m}{\sqrt{1+(\frac{dx}{dy})^2}}$
になるので、両者を組み合わせて、 $c=\frac{v_A}{v_m}$ とすると、
$\frac{d^2x}{dt^2}(b-y)=c\sqrt{1+(\frac{dx}{dy})^2}$
という、非線型第2階微分方程式が得られる。 $p=\frac{dx}{dy}$ という変換を用いると、 $\frac{dp}{dy}=\frac{d^x}{dy^2}$ であり、
$\frac{dp}{dy}(b-y)=c\sqrt{1+p^2}$ となる。これは変数分離で積分することが可能で、
$\int \frac{1}{\sqrt{1+p^2}}dp=c\int \frac{1}{b-y}dy$
$\log(p+\sqrt{1+p^2})=-c\log(b-y)+K$
となる。初期条件を考えると、 $t=0$ で
$y=0$
$p=\frac{dx}{dy}=\frac{a}{b}=d$
である。簡単のため、 $f=d+\sqrt{1+d^2}$ とおくと、
$K=\log(d+\sqrt{1+d^2})+c\log b=\log(fb^c)$
となるので、戻して
$\log(p+\sqrt{1+p^2})=-c\log(b-y)+\log(fb^c)$
$\frac{dx}{dy}=p=\frac{1}{2}\{ \frac{fb^c}{(b-y)^c}-\frac{(b-y)^c}{fb^c}\}$
$x=\frac{1}{2}\{ \frac{fb^c}{(c-1)(b-y)^{c-1}}+\frac{(b-y)^{c+1}}{(c+1)fb^c}\}+L$
と、積分までがんばればなる。 $L$ は、初期条件 $x=y=0$ より得られ、
$L=\frac{b\{ (f^2+1)c+f^2-1\}}{2f(1-c^2)}$
であるから、ミサイルの軌道は
　
$x=\frac{1}{2}\{ \frac{fb^c}{(c-1)(b-y)^{c-1}}+\frac{(b-y)^{c+1}}{(c+1)fb^c}\}+\frac{b\{ (f^2+1)c+f^2-1\}}{2f(1-c^2)}$
　
である。大変。
　
ミサイルがカサッカに当たるときの座標と時刻を考える。このとき、 $x=a+v_At$ 、 $y=b$ になるから、本中では謎の代数計算を行い、
$t=\frac{\sqrt{a^2+b^2} + a\frac{v_A}{v_m}}{v_m(1-(\frac{v_A}{v_m})^2)}$
となる。

a <- 5000
b <- 3000
d <- a / b
f <- d + sqrt(1 + d^2)
vA <- 1000
vm <- 2000
vAm <- vA / vm #cの代わり

#微分方程式は、距離が高度の関数になっていることに注意する。
missile <- function(y){
	L <- (b*((f^2 + 1)*vAm + f^2 - 1)) / (2*f*(1 - vAm^2) )
	leftbunshi <- (b -y)^(vAm +1)
	leftbunbo <- (vAm +1)*f*(b^vAm)
	
	rightbunshi <- f*(b^vAm)*((b - y)^(1 - vAm))
	rightbunbo <- 1 - vAm
	
	return(1/2 * ((leftbunshi / leftbunbo) - (rightbunshi / rightbunbo)) + L)
}

y.altitude <- seq(0, 3000, length=10000) #高度
x.dist <- missile(y.altitude) #距離
xlim <- range(x.dist, na.rm=TRUE)
ylim <- range(y.altitude)
plot(x.dist, y.altitude, type="l", xlab="distance", ylab="altitude", xlim=xlim, ylim=ylim)
points(a, b, cex=1, pch=16, col=2)
segments(x0=a, y0=b, x1=xlim[2], lty=2, col=2)

bakuhatsu <- (sqrt(a^2 + b^2) + a*vAm) / (vm*(1 - vAm^2)) #両者がぶつかる時刻