Return to 高さ h からの斜方投射の最大到達距離を求める準備

SymPy で高さ h からの斜方投射の最大到達距離を求める

地面から高さ $h$ の地点から空気抵抗なしの斜方投射を行うと，初速度の大きさを一定とした場合，水平方向の到達距離が最大となるのは打ち上げ角度（仰角）が何度のときかを Python の SymPy を使って（かつ SciPy や NumPy を使わずに）求める。

運動方程式の解

初期条件を (t = 0) で

$$x(0) = 0, \quad y(0) = h, \quad v_x(0)  = v_0 \cos\theta, \quad v_y(0) = v_0 \sin \theta$$

としたときの解は

\begin{eqnarray}
x(t) &=& v_0 \cos\theta\cdot t \\
y(t) &=& h + v_0 \sin\theta\cdot t – \frac{1}{2} g t^2 \\
v_x(t) &=& v_0 \cos\theta \\
v_y(t) &=& v_0 \sin\theta – g t
\end{eqnarray}

無次元化された解

この系に特徴的な時間 $\displaystyle \tau \equiv \frac{v_0}{g}$ および長さ $\displaystyle \ell \equiv v_0 \tau = \frac{v_0^2}{g}$ で解を無次元化すると，

\begin{eqnarray}
T &\equiv& \frac{t}{\tau} \\
H &\equiv& \frac{h}{\ell} \\
X &\equiv& \frac{x}{\ell} = \cos\theta\cdot T \\
Y&\equiv& \frac{y}{\ell} = H+ \sin\theta\cdot T – \frac{1}{2} T^2 \\
V_x &\equiv& \frac{d\bar{x}}{d\bar{t}} = \cos\theta \\
V_y &\equiv& \frac{d\bar{y}}{d\bar{t}} = \sin\theta – T
\end{eqnarray}

無次元化された滞空時間と水平到達距離

無次元化された滞空時間 $T_1$ は

$$T_1 = \sin\theta + \sqrt{\sin^2\theta + 2 H}$$

この時間での水平方向の無次元化された到達距離 $L$ は

$$L = X(T_1) = \cos\theta\cdot T_1$$

詳細は以下のページ

必要なモジュールの import

必要なモジュールを import します。

In [1]:
from sympy.abc import *
from sympy import * 
# SymPy Plotting Backends (SPB) でグラフを描く
from spb import *

# グラフを SVG で Notebook にインライン表示
%config InlineBackend.figure_formats = ['svg']

関数の定義

無次元化された解 $X(T), Y(T)$ および滞空時間 $T_1$ と水平到達距離 $L$ を関数として定義します。

In [2]:
# 角度 th は度で 
def X(T, th):
    theta = th * pi/180
    return cos(theta) * T

def Y(T, th, H):
    theta = th * pi/180
    return H + sin(theta)*T - T**2/2

def T1(th, H):
    theta = th * pi/180
    return sin(theta) + sqrt(sin(theta)**2 + 2*H)

def L(th, H):
    theta = th * pi/180
    return cos(theta) * T1(th, H)

$H = 0$ の場合の軌道

$\theta = 45^{\circ}$ 前後の角度での斜方投射の軌道をグラフにしてみる。

まずは，とりあえずのグラフ：

SPB で媒介変数表示のグラフを描くには，

plot_parametric(x(t), y(t), (t, tmin, tmax))

デフォルトで use_cm = True になっているので，カラーマップが不要なら False に。

In [3]:
# θ = 45° のグラフ
th = 45
H = 0

plot_parametric(X(T, th), Y(T, th, H), (T, 0, T1(th, H)));

以下のようなオプションを設定して，もう少し体裁を整えます。

  • $\theta$ の値を凡例に
  • カラーマップは不要 use_cm = False
  • 座標軸のラベルとグラフのタイトル
  • 横軸縦軸の表示範囲
  • 縦軸横軸のアスペクト比
  • 座標軸の目盛の間隔
In [4]:
# θ = 45° のグラフ
th = 45
H = 0

# 凡例
key = 'θ = %2d°' % th

p=plot_parametric(
    X(T, th), Y(T, th, H), (T, 0, T1(th, H)), key, 
    legend = True,
    use_cm = False, 
    xlabel = "X", ylabel = "Y",
    title = "H = 0 の場合の斜方投射",
    xlim = (0, 1.1), ylim = (0, 0.5), 
    aspect = "equal", show = False);

ax=p.ax
# ax を定義したら，あとは Matplotlib 流に
# SymPy だけでやる。NumPy に頼らない
# np.linspace() や np.arange() を使わない方針で
# x の目盛を 0.1 刻みに
ax.set_xticks([0.1*i for i in range(11)])
# y の目盛を 0.1 刻みに
ax.set_yticks([0.1*i for i in range(6)]);

$\theta$ の値を $43^{\circ}$ から $47^{\circ}$ まで $1^{\circ}$ 刻みで変えて，複数の曲線を一度に描く例。

なお，細かいことですが，グラフの最高点は $\theta$ の値が大きいときですので，先に $47^{\circ}$ から描きはじめて，$1^{\circ}$ ずつ減らしてみます。

（plot_parametric() の中で for などの繰り返し処理の仕方がわからないので… ）

In [5]:
def para(th, H):
    return X(T, th), Y(T, th, H), (T, 0, T1(th, H)), 'θ = %2d°' % th

H = 0
p=plot_parametric(
    para(47,H), para(46,H), para(45,H), para(44,H), para(43,H), 
    legend = True,
    use_cm = False, 
    xlabel = "X", ylabel = "Y",
    title = "H = %3.1f の場合の斜方投射" % H,
    xlim = (0, 1.1), ylim = (0, 0.5), 
    aspect = "equal", show = False);

ax=p.ax
# ax を定義したら，あとは Matplotlib 流に
# SymPy だけでやる。NumPy に頼らない
# np.linspace() や np.arange() を使わない方針で
# x の目盛を 0.1 刻みに
ax.set_xticks([0.1*i for i in range(11)])
# y の目盛を 0.1 刻みに
ax.set_yticks([0.1*i for i in range(6)]);

最大到達距離となる角度 $\theta$ を調べるために，着地点付近を拡大表示します。

In [6]:
H = 0
plot_parametric(
    para(47,H), para(46,H), para(45,H), para(44,H), para(43,H), 
    legend = True,
    use_cm = False, 
    xlabel = "X", ylabel = "Y",
    title = "H = %3.1f の場合の斜方投射" % H,
    # 着地点付近を拡大表示
    xlim = (0.99, 1.0), ylim = (0, 0.01));

上のグラフから，$H=0$ の場合には確かに $\theta=45^{\circ}$ の時に最大到達距離になっていることがわかります。

$H = 0.2$ の場合の軌道

今度は，$H=0.2$ としてやってみます。

In [7]:
H = 0.2
p=plot_parametric(
    para(47,H), para(46,H), para(45,H), para(44,H), para(43,H), 
    para(42,H), para(41,H), para(40,H), para(39,H), para(38,H), 
    legend = True,
    use_cm = False, 
    xlabel = "X", ylabel = "Y",
    title = "H = %3.1f の場合の斜方投射" % H,
    xlim = (0, 1.5), ylim = (0, 0.7), 
    aspect = "equal", show = False);

ax=p.ax
# ax を定義したら，あとは Matplotlib 流に
# SymPy だけでやる。NumPy に頼らない
# np.linspace() や np.arange() を使わない方針で
# x の目盛を 0.1 刻みに
ax.set_xticks([0.1*i for i in range(14)]);
# y の目盛を 0.1 刻みに
ax.set_yticks([0.1*i for i in range(7)]);

最大到達距離となる角度 $\theta$ を調べるために，着地点付近を拡大表示します。

In [8]:
H = 0.2
plot_parametric(
    para(47,H), para(46,H), para(45,H), para(44,H), para(43,H), 
    para(42,H), para(41,H), para(40,H), para(39,H), para(38,H), 
    legend = True,
    use_cm = False, 
    xlabel = "X", ylabel = "Y",
    title = "H = %3.1f の場合の斜方投射" % H,
    # 着地点付近を拡大表示
    xlim = (1.15, 1.19), ylim = (0, 0.02));

水平到達距離が最大となるのは $\theta=45^{\circ}$ のときではない！ことがわかるだろう。

水平到達距離 $L$

定量的に調べるために，水平到達距離 L(th, H) をグラフにしてみる。

In [9]:
H = 0.2
# th を初期化
var('th')
plot(L(th, H), (th, 38, 48), 
     xlabel = "打ち出し角度 (°)", 
     ylabel = "規格化された水平到達距離", 
     title = "H = %3.1f の場合の打ち出し角度と水平到達距離" % H);

上のグラフをみると，$H=0.2$ の場合に水平到達距離が最大となるのは $\theta=45^{\circ}$ のときではなく… どのくらい？

$L$ が最大値となる角度を解析的微分で求める

$u \equiv \tan\theta$ として水平到達距離 $L$ を $u$ で表すと

\begin{eqnarray}
L(u, H)
&=& \frac{1}{1+u^2} \left\{u + \sqrt{(1+2H) u^2 + 2H} \right\}
\end{eqnarray}

In [10]:
var('u H')
def Lu(u, H):
    return (u+sqrt((1+2*H)*u**2 + 2*H))/(1 + u**2)
Lu(u, H)
Out[10]:
$\displaystyle \frac{u + \sqrt{2 H + u^{2} \cdot \left(2 H + 1\right)}}{u^{2} + 1}$

$\displaystyle \frac{d Lu(u, H)}{du} = 0$ となる $u$ を solve() で求めます。

In [11]:
dL = diff(Lu(u, H), u)
sols = solve(dL, u)
sols
Out[11]:
[-sqrt(1/(2*H + 1)), sqrt(1/(2*H + 1))]

さすが SymPy! あっさりと答えを出してくれましたねぇ… などと安易に信じてはいけません！（私も最初はすっかり騙されました。）

最初の解が解になっていないことは，実際に dL に代入してみればわかります。

In [12]:
dLsubs0 = dL.subs(u, sols[0])
factor(dLsubs0)
Out[12]:
$\displaystyle – \frac{\left(2 H + 1\right) \left(- H \sqrt{2 H + 1} \sqrt{\frac{1}{2 H + 1}} – H + \frac{2 H}{2 H + 1} – 1 + \frac{1}{2 H + 1}\right)}{2 \left(H + 1\right)^{2}}$
In [13]:
simplify(dLsubs0*2*(1+H)**2/(1+2*H))
Out[13]:
$\displaystyle H \left(\sqrt{2 H + 1} \sqrt{\frac{1}{2 H + 1}} + 1\right)$

一方で，二つめのほうは代入すると，確かに解になっているようです。

In [14]:
dLsubs1 = dL.subs(u, sols[1])
factor(dLsubs1)
Out[14]:
$\displaystyle – \frac{\left(2 H + 1\right) \left(H \sqrt{2 H + 1} \sqrt{\frac{1}{2 H + 1}} – H + \frac{2 H}{2 H + 1} – 1 + \frac{1}{2 H + 1}\right)}{2 \left(H + 1\right)^{2}}$
In [15]:
simplify(dLsubs1*2*(1+H)**2/(1+2*H))
Out[15]:
$\displaystyle H \left(- \sqrt{2 H + 1} \sqrt{\frac{1}{2 H + 1}} + 1\right)$

（これ以上簡単化してくれないところも，SymPy の意固地なところです，全く。）

平方根を含む方程式の解を求めるときには，こっそり両辺を二乗したりするので，その際に本来の方程式の解ではない，うその解が紛れ込む可能性があります。

元の方程式 dL = 0 を確かに満たしている解かどうかを確認する必要があります。

題意より $u > 0$ ですから，$\displaystyle u = u_{\rm max} \equiv \frac{1}{\sqrt{1+ 2H}}$ のとき，$Lu(u, H)$ が極値（実際には最大値）を取ることがわかります。

In [16]:
def umax(H):
    return 1/sqrt(1+2*H)

このときの水平到達距離 $Lu(u_{\rm max}, H)$ は

In [17]:
var('H')
simplify(Lu(umax(H), H))
Out[17]:
$\displaystyle \sqrt{2 H + 1}$

つまり，高さ $H$ のときの最大水平到達距離は，

$$L_{\rm max} = \sqrt{1 + 2H}$$

となることがわかります。

また，後にグラフを描くために，

\begin{eqnarray}
u &=& \tan \theta \\
\therefore\ \ \theta_{\rm max} &=& \tan^{-1} u_{\rm max}
\end{eqnarray}

つまり，高さ $H$ に対して最大水平到達距離になる角度を返す関数を，以下のように定義しておきます。

In [18]:
def thmax(H):
    """ 高さ H からの斜方投射で最大到達距離になる角度を°で返す
    """
    return atan(umax(H))*180/float(pi)
In [19]:
thmax(0.2)
Out[19]:
$\displaystyle 40.2029658865698$
In [20]:
# 上の結果を使ってグラフを描く

H = 0.2

# 凡例の数値の桁数を揃える
key = 'H=%3.1f, θ=%6.3f, Lmax=%6.3f' % (H, thmax(H), L(thmax(H), H))

p=plot_parametric(
    X(T, thmax(H)), Y(T, thmax(H), H), (T, 0, T1(thmax(H), H)), key,
    legend = True,
    use_cm = False, 
    xlabel = "X", ylabel = "Y",
    title = "H = %3.1f の場合の斜方投射" % H,
    xlim = (0, 1.2), ylim = (0, 0.5), 
    aspect = "equal", show = False);

ax=p.ax
# ax を定義したら，あとは Matplotlib 流に
# SymPy だけでやる。NumPy に頼らない
# np.linspace() や np.arange() を使わない方針で
# x の目盛を 0.1 刻みに
ax.set_xticks([0.1*i for i in range(13)]);
# y の目盛を 0.1 刻みに
ax.set_yticks([0.1*i for i in range(6)]);

問題 1

高さ $H$ からの斜方投射の水平到達距離が最大となる打ち出し角度と，そのときの最大水平到達距離がわかるようなグラフを次の $H$ についてまとめて描け。

$H = 0.2, 0.4, 0.6, 0.8, 1.0$

解答例

In [21]:
def para(H):
    th = thmax(H)
    key = 'H=%3.1f, θ=%6.3f, Lmax=%6.3f' % (H, thmax(H), L(thmax(H), H))
    return X(T, th), Y(T, th, H), (T, 0, T1(th, H)), key

p=plot_parametric(
    para(1.0), para(0.8), para(0.6), para(0.4), para(0.2), 
    legend = True,
    use_cm = False, 
    xlabel = "X", ylabel = "Y",
    title = "高さ H からの斜方投射の最大水平到達距離",
    xlim = (0, 2.4), ylim = (0, 1.5), 
    aspect = "equal", show = False);

ax=p.ax
# ax を定義したら，あとは Matplotlib 流に
# SymPy だけでやる。NumPy に頼らない
# np.linspace() や np.arange() を使わない方針で
# x の目盛を 0.1 刻みに
ax.set_xticks([0.2*i for i in range(11)])
# y の目盛を 0.1 刻みに
ax.set_yticks([0.2*i for i in range(8)]);

問題 2

高さ $H$ を横軸に，そのときの水平到達距離が最大となる打ち出し角度 $\theta$ を縦軸にしたグラフを描け。

解答例

In [22]:
# H を初期化しておく
var('H')
plot(thmax(H), (H, 0, 1), 
     xlabel = "高さ H", 
     ylabel = "θmax (°)", 
     title = "高さ H からの斜方投射の最大水平到達距離を与える角度");

問題 3

高さ $H$ を横軸に，そのときの最大水平到達距離を縦軸にしたグラフを描け。

解答例

In [23]:
# 既に問題2でやっているが，念のため H を初期化
var('H')
plot(L(thmax(H), H), (H, 0, 1), 
     xlabel = "高さ H", 
     ylabel = "最大水平到達距離 L", 
     title = "高さ H からの斜方投射の最大水平到達距離");

問題 4

$h = 5\,\mbox{(m)}$ の高さから速さ $v_0 = 10\,\mbox{(m/s)}$ で斜方投射したときの最大水平到達距離は何 $\mbox{m}$ か。またそのときの打出し角度は何度か。

ヒント：

$h = 5\,\mbox{(m)}$ は規格化された高さ $H$ ではいくらか，また規格化された最大到達距離 $L_{\rm max}(H)$ を次元をもった量に直すといくらか，ということ。

\begin{eqnarray}
\ell &=& \frac{v_0^2}{g} \\
H &=& \frac{h}{\ell} = \frac{g h}{v_0^2} \\
x &=& \ell X = \frac{v_0^2}{g} X \\
x_{\rm max} &=& \frac{v_0^2}{g} L_{\rm max}(H)
\end{eqnarray}

重力加速度 $g = 9.80665 \,(\mbox{m}/\mbox{s}^2)$

解答例

In [24]:
h = 5
v0 = 10
g = 9.80665

# 系に特徴的な長さ
ell = v0**2/g

# 規格化された高さ H
H = h/ell

# ell をかけて次元をもった量に変換
xmax = ell * L(thmax(H), H)

print('高さ h = %2d m のときの最大水平到達距離は' % h, end=' ')
print('%9.6f m' % xmax)

print('そのときの打ち出し角度は %9.6f °' % thmax(H))

高さ h =  5 m のときの最大水平到達距離は 14.351088 m
そのときの打ち出し角度は 35.395688 °