微分方程式モデルでPursuitCurve問題を解く:hawk-pigeon問題

微分方程式モデルでPursuitCurve問題:hawk-pigeon問題（鳩を追いかける鷹の追跡曲線）を解くメモ．

hawk-pigeonモデルとここでは呼んでいるが他の専門書では

商船を追いかける海賊船
飛んでいる飛行機を撃ち落とすための迎撃ミサイル

などの追跡曲線で紹介されている．

今回は以下をまとめる．

その微分方程式がどんな形になるか
どう解くか
capture time（捕捉時間）はどうなるか

参考：

環境

Windows10 ver2004（OSビルド 19041.330）
Docker for win v19.03.12
- jupyter/scipy-notebook

Pursuit Curveとは

Pursuit Curveとは日本語で追跡曲線のことである．ある追跡者(chaser)が環境の外部要因を受けながら対象(target)を追いかけるときの軌跡を求める問題である．

解析としては，初期位置や初期速度により追跡者が対象に到達可能化どうか，到達可能な場合そのときの捕捉時間(capture time)などを調べる．

対象物への追跡曲線の作り方・考え方としては，追跡曲線の接線が常に対象物と交わるにすることでその軌跡はいずれ接線と平行するようになり，十分な速度があれば対象物に到達できる．

問題設定

2次元平面$(x,y)$を考え，

鳩Q(target)が等速$v>0$でこの座標$(b=0,vt)$上を飛んでいるとする．
鷹P(chaser)は時間$t=0$で初期位置$(a>0, 0)$で鳩Qを見つけそれを等速$w>0$で追跡しようとする．

このときの鷹Pの追跡曲線を求める．

微分方程式を立てる

鷹の軌道$P(x_{h}, y_{h})$の接線の傾きは

$$ \frac{dy(x_{h})}{dx}= \frac{y - y_{h}}{x - x_{h}} $$

鳩Q$(b=0, vt)$と交わるので

$$ \frac{dy(x_h)}{dx} =\frac{y-y_h}{x-x_h}\\ = \frac{vt-y_h}{b-x_h}\\ = \frac{vt-y_h}{0-x_h}\\ = \frac{y_h-vt}{x_h}\\ $$

時間$t$の式に直す．

$\begin{eqnarray*} \frac{vt-y}{b-x} &=& \frac{dy}{dx}\\ vt-y &=& (b-x)\frac{dy}{dx}\\ vt &=& y + (b-x)\frac{dy}{dx}\\ \end{eqnarray*}$

$\begin{eqnarray} t = \frac{y}{v} - \frac{x-b}{v}\frac{dy}{dx}, (1) \end{eqnarray}$

合成速度$w$は

$$ w = \sqrt{\dot{x}^{2} + \dot{y}^{2}}\ = \sqrt{(\frac{dx}{dt})^{2} + (\frac{dy}{dt})^{2}} $$

時間tで積分する．

$$ wt = \int_{0}^{t} \sqrt{(\frac{dx}{dt'})^{2} + (\frac{dy}{dt'})^{2}}dt' $$

時間[0, t]でx軸で[a,a-x]まで移動したとすると

$$ wt = \int_{a}^{a-x} \sqrt{(dx')^{2} + (dy)^{2}}\\ = \int_{a}^{a-x} \sqrt{1 + (\frac{dy}{dx'})^{2}}dx' $$

$a-x=x'$として変数変換する.

$$ dx' = -dx\\ wt = \int_{0}^{x'} \sqrt{1 + (\frac{dy}{dx})^{2}}(-dx)\\ = - \int_{0}^{x'} \sqrt{1 + (\frac{dy}{dx})^{2}}dx $$

時間tの式に直すと

$\begin{eqnarray} t= -\frac{1}{w} \int_0^x \sqrt{1 + (\frac{dy}{dx'})^2}dx', (2) \end{eqnarray}$

(1)式と(2)式で等式を作ると，

$$ -\frac{1}{w} \int_{0}^{x} \sqrt{1 + (\frac{dy}{dx'})^{2}}dx' = \frac{y}{v} - \frac{x-b}{v}\frac{dy}{dx} $$

非線形微分方程式になる．

置換して解ける形に直す． $p(x) = \frac{dy}{dx}$として置換する．

$$ -\frac{1}{w} \int_{0}^{x} \sqrt{1 + p(x')^{2}}dx' = \frac{y}{v} - \frac{x-b}{v}p(x)\ $$

両辺を$x$で微分する．

$$ -\frac{1}{w} (\sqrt{1 + p(x)^{2}} - \sqrt{1 + p(a)^{2}}) = \frac{1}{v}\frac{dy}{dx} - \frac{x-b}{v}\frac{dp}{dx} - \frac{1}{v}p(x) $$

$dp/dx$の式に直す．

$\begin{eqnarray*} \frac{x-b}{v}\frac{dp}{dx} &=& \frac{1}{w} (\sqrt{1 + p(x)^2})\\ (x-b)\frac{dp}{dx} &=& \frac{v}{w} (\sqrt{1 + p(x)^2})\\ (x-b)\frac{dp}{dx} &=& u (\sqrt{1 + p(x)^2}), u=\frac{v}{w}\\ \end{eqnarray*}$

微分方程式を解く

変数分離可能な形になったので解くことができる．

変数分離に直す．

$$ \frac{dp}{\sqrt{1 + p(x)^{2}}} = \frac{u}{x-b} dx\ $$

$x>b$として積分する．

$$ \ln{[p + \sqrt{1 + p(x)^{2}}]} + C = u\ln{(x-b)}\ $$

初期条件$x=a, y=0, p(a)=dy/dx = 0, b=0 $として定数を求める．

$\begin{eqnarray*} \ln{[p + \sqrt{1 + p(x)^2}]} + C &=& u\ln{(x-b)}\\\\ \ln{[0 + \sqrt{1 + 0^2}]} + C &=& u\ln{(x-b)}\\\\ C &=& u\ln{(a-b)}, (a>b) \end{eqnarray*}$

定数を代入して

$\begin{eqnarray*} \ln{[p + \sqrt{1 + p(x)^2}]} + u\ln{(a-b)} &=& u\ln{(x-b)}\\\\ \ln{[p + \sqrt{1 + p(x)^2}]} - \ln{(x-b)^{u}}+\ln{(a-b)^u} &=& 0 \end{eqnarray*}$

変形する

$\begin{eqnarray*} 0 = \ln{[p + \sqrt{1 + p(x)^2}]} - \ln{(x-b)^u}+\ln{(a-b)^u}\\ 0 = \ln{[p + \sqrt{1 + p(x)^2}]} + \ln{\left(\frac{a-b}{x-b}\right)^u} \end{eqnarray*}$

$\ln{1}=0$より

$$ [p + \sqrt{1 + p(x)^{2}}] \cdot \left(\frac{a-b}{x-b}\right)^{u} = 1 $$

左辺をpの式に直す．

$$ [p + \sqrt{1 + p(x)^{2}}] = \left(\frac{a-b}{x-b}\right)^{-u}\ $$

$q = \left(\frac{a-b}{x-b}\right)^{-u}$とおいて

$\begin{eqnarray*} p + \sqrt{1+p^2} &=& q \\ \sqrt{1+p^2} &=& q - p \\ 1+p^2 &=& (q-p)^2\\ 1+p^2 &=& q^2 - 2pq + p^2\\ 2pq &=& q^2 -1 \\ \end{eqnarray*}$

$$ p = \frac{1}{2}\left( q - \frac{1}{q} \right) $$

pとqの変数をもとに戻すと

$$ \frac{dy}{dx} = \frac{1}{2}\left( \left(\frac{a-b}{x-b}\right)^{-u} - \left(\frac{a-b}{x-b}\right)^{u} \right) $$

積分すると

$\begin{eqnarray*} y(x) &=& \frac{1}{2}\left( \int \left(\frac{a-b}{x-b}\right)^{-u}dx - \int\left(\frac{a-b}{x-b}\right)^{u} dx \right) +C\\ &=& \frac{1}{2}\left( \left( \frac{(x-b)^{u+1}}{(u+1)(a-b)^{u}} \right) - \left(\frac{(a-b)^{u}(x-b)^{-u+1}}{(-u+1)}\right) \right) +C \end{eqnarray*}$

初期条件$x=a, y(a)=0$より，定数を求める．

$\begin{eqnarray*} 0 &=& \frac{1}{2}\left( \left( \frac{(a-b)^{u+1}}{(u+1)(a-b)^{u}} \right) - \left(\frac{(a-b)^{u}(a-b)^{-u+1}}{(-u+1)}\right) \right)+C\\ C &=& \frac{1}{2}\left( \left( \frac{(a-b)}{(-u+1)} \right) - \left(\frac{(a-b)}{(u+1)}\right) \right)\\ &=& \frac{1}{2}\left( \frac{(a-b)(1+u) - (a-b)(1-u)}{1-u^2} \right)\\ &=& \frac{1}{2}\left( \frac{2u(a-b)}{1-u^2} \right)\\ &=& \left( \frac{u(a-b)}{1-u^2} \right)\\ \end{eqnarray*}$

まとめると，

$\begin{eqnarray} y(x) &=& \frac{1}{2}\left( \left( \frac{(x-b)^{u+1}}{(u+1)(a-b)^{u}} \right) - \left(\frac{(a-b)^{u}(x-b)^{-u+1}}{(-u+1)}\right) \right) +C\\ C&=& \left( \frac{u(a-b)}{1-u^2} \right)\\ \end{eqnarray}$

鳩Qの初期x座標が$b=0$とすると

$\begin{eqnarray} y(x) &=& \frac{1}{2}\left( \left( \frac{x^{u+1}}{(u+1)a^{u}} \right) - \left(\frac{a^{u}x^{-u+1}}{(-u+1)}\right) \right) +C\\ C&=& \left( \frac{ua}{1-u^2} \right)\\ \end{eqnarray}$

可視化と分析

鳩と鷹の速度$v, w$を変えて導出した式を使って可視化する．

$\begin{eqnarray} y(x) &=& \frac{1}{2}\left( \left( \frac{x^{u+1}}{(u+1)a^{u}} \right) - \left(\frac{a^{u}x^{-u+1}・・}{(-u+1)}\right) \right) +C\\ C&=& \left( \frac{ua}{1-u^2} \right)\\ \end{eqnarray}$

import numpy as np
import matplotlib.pyplot as plt

v = 1
w = 1.5
u = v/w
a = 1
C = ( u * a / (1-u**2) )
# print(C)

x = np.linspace(0, a, 1000)
y = .5 * ((x**(u+1)/((u+1)*a**(u)) ) - ((x**(-u+1)*(a)**u)/((-u+1)))) + C
cap_t = a/(w*(1-u**2))
print(f"cap_t = {cap_t}")

fig = plt.figure(figsize=(5,5))
plt.plot(x, y, label=f"$v/w={u:.3}$")

v = 1
w = 0.999
u = v/w
C = ( u * a / (1-u**2) )
y = .5 * ((x**(u+1)/((u+1)*a**(u)) ) - ((x**(-u+1)*(a)**u)/((-u+1)))) + C
plt.plot(x, y, label=f"$v/w={u:.3}$")
cap_t = a/(w*(1-u**2))
print(f"cap_t = {cap_t}")

v = 1.5
w = 1
u = v/w
C = ( u * a / (1-u**2) )
y = .5 * ((x**(u+1)/((u+1)*a**(u)) ) - ((x**(-u+1)*(a)**u)/((-u+1)))) + C
plt.plot(x, y, label=f"$v/w={u:.3}$")
cap_t = a/(w*(1-u**2))
print(f"cap_t = {cap_t}")

plt.xlim(-.1, 2)
plt.ylim(-.1, 2)
plt.grid(True)
plt.legend()
plt.show()

速度比3パターンを比較する．

$u < w$: 鳩の横軸$b=0$に近づくのが早く鳩に到達する．
$u=w$: yの式から計算不可能，つまり鳩に追いつかない．
$u>w$: 計算は可能だが鳩の横軸$b=0$に近づくのが遅く，横軸にたどり着いたとしても垂直直線y軸から速度差により追いつけない．

capture time（捕捉時間）を求める

capture timeを求める．

$x=b, y=vt$

$\begin{eqnarray*} y(x=b) &=& 0 +C\\\\ vt &=& 0 + C\\\\ vt &=& C\\\\ t &=& \frac{C}{v}\\\\ &=& \frac{1}{v}\frac{u(a-b)}{1-u^2}\\\\ &=& \frac{(a-b)}{w(1-u^2)} \end{eqnarray*}$

$b=0$とすると

$\begin{eqnarray} t &=& \frac{(a-b)}{w(1-u^2)} \end{eqnarray}$

$1-u^{2} > 0 \to 1 > u^{2} \to w^{2} > v^{2} \to w > v$のとき，計算可能．つまり鳩に追いつく．
$1-u^{2} \leq 0 \to 1 \leq u^{2} \to w^{2} \leq v^{2} \to w \leq v$のときcapture timeはマイナスになり計算不可，つまり鷲は鳩に追いつけない．

はしくれエンジニアもどきのメモ

情報系技術・哲学・デザインなどの勉強メモ・備忘録です。

微分方程式モデルでPursuitCurve問題を解く:hawk-pigeon問題

微分方程式モデルでPursuitCurve問題を解く:hawk-pigeon問題

環境

Pursuit Curveとは

問題設定

微分方程式を立てる

微分方程式を解く

可視化と分析

capture time（捕捉時間）を求める