地面から高さ $h$ の地点から空気抵抗なしの斜方投射を行うと，初速度の大きさを一定とした場合，水平方向の到達距離が最大となるのは打ち上げ角度（仰角）が何度のときかを Python の SciPy と NumPy を使って（かつ SymPy を使わずに）求める。
運動方程式の解
初期条件を (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 します。
# NumPy を使います
import numpy as np
# Matplotlib でグラフを描きます
import matplotlib.pyplot as plt
# 方程式の数値解
from scipy.optimize import root_scalar
# グラフを SVG で Notebook にインライン表示させる設定
%config InlineBackend.figure_formats = ['svg']
関数の定義
無次元化された解 $X(T), Y(T)$ および滞空時間 $T_1$ と水平到達距離 $L$ を関数として定義します。
# 角度 th は度で
def X(T, th):
theta = np.radians(th)
return np.cos(theta) * T
def Y(T, th, H):
theta = np.radians(th)
return H + np.sin(theta)*T - T**2/2
def T1(th, H):
theta = np.radians(th)
return np.sin(theta) + np.sqrt(np.sin(theta)**2 + 2*H)
def L(th, H):
theta = np.radians(th)
return np.cos(theta) * T1(th, H)
$H = 0$ の場合の軌道
$\theta = 45^{\circ}$ 前後の角度での斜方投射の軌道をグラフにしてみる。
まずは，とりあえずのグラフ：
# θ = 45° のグラフ
th = 45
H = 0
# 滞空時間 T1 を 100 等分
T_list = np.linspace(0, T1(th, H), 100+1)
plt.plot(X(T_list, th), Y(T_list, th, H));
以下のようなオプションを設定して，もう少し体裁を整えます。
- $\theta$ の値を凡例に
- グリッド（格子線）をドットに
- 座標軸のラベルとグラフのタイトル
- 横軸縦軸の表示範囲
- 座標軸の目盛の間隔
- 縦軸横軸のアスペクト比
# θ = 45° のグラフ
th = 45
H = 0
# 滞空時間 T1 を 100 等分
T_list = np.linspace(0, T1(th, H), 100+1)
# 凡例
key = 'θ = %2d°' % th
# plot でグラフを描く
plt.plot(X(T_list, th), Y(T_list, th, H), label = key)
# グラフのオプション設定例
# グリッドをドットで
plt.grid(linestyle='dotted')
# 座標軸のラベル
plt.xlabel("X")
plt.ylabel("Y")
# グラフのタイトル
plt.title("H = 0 の場合の斜方投射")
# 横軸縦軸の表示範囲
plt.xlim(0, 1.1)
plt.ylim(0, 0.5)
# 座標軸の目盛の間隔を手動で設定する場合
# 間隔を指定して np.arange() を使う例
# 0 から 0.1 間隔で 1.1 未満（つまり 1.0）まで
plt.xticks(np.arange(0, 1.1, 0.1))
# 要素数を指定して np.linspace() を使う例
# (0, 0.5) 区間を 5 等分
plt.yticks(np.linspace(0, 0.5, 5+1))
# 45°が45°に見えるように
# 縦軸横軸のアスペクト比を 'equal' に。
plt.gca().set_aspect('equal')
# 凡例を表示
plt.legend();
$\theta$ の値を $43^{\circ}$ から $47^{\circ}$ まで $1^{\circ}$ 刻みで変えて，複数の曲線を一度に描く例。
for で繰り返し処理させてみます。
なお，細かいことですが，グラフの最高点は $\theta$ の値が大きいときですので，先に $47^{\circ}$ から描きはじめて，$1^{\circ}$ ずつ減らしてみます。
H = 0
# 複数のグラフを一度に描く場合，
# for ループで繰り返し処理
for th in range(47, 42, -1):
# 滞空時間 T1 を 100 等分
T_list = np.linspace(0, T1(th, H), 100+1)
# 凡例
key = 'θ = %2d°' % th
# plot でグラフを描く
plt.plot(X(T_list, th), Y(T_list, th, H), label = key)
plt.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
plt.xlabel("X")
plt.ylabel("Y")
plt.title("H = 0 の場合の斜方投射")
# 横軸縦軸の表示範囲と目盛
plt.xlim(0, 1.1)
plt.ylim(0, 0.5)
plt.xticks(np.arange(0, 1.1, 0.1))
plt.gca().set_aspect('equal')
# 凡例を表示
plt.legend();
最大到達距離となる角度 $\theta$ を調べるために，着地点付近を拡大表示します。
H = 0
# 複数のグラフを一度に描く場合，
# for ループで繰り返し処理
for th in range(47, 42, -1):
# 滞空時間 T1 を 100 等分
T_list = np.linspace(0, T1(th, H), 100+1)
# 凡例
key = 'θ = %2d°' % th
# plot でグラフを描く
plt.plot(X(T_list, th), Y(T_list, th, H), label = key)
plt.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
plt.xlabel("X")
plt.ylabel("Y")
plt.title("H = 0 の場合の斜方投射")
# 着地点付近を拡大表示
plt.xlim(0.99, 1.0)
plt.ylim(0, 0.01)
# 凡例を表示
plt.legend();
plt.savefig("Nhtosha04.svg");
上のグラフから，$H=0$ の場合には確かに $\theta=45^{\circ}$ の時に最大到達距離になっていることがわかります。
$H = 0.2$ の場合の軌道
今度は，$H=0.2$ としてやってみます。
H = 0.2
for th in range(47, 37, -1):
T_list = np.linspace(0, T1(th, H), 100+1)
key = 'θ = %2d°' % th
plt.plot(X(T_list, th), Y(T_list, th, H), label = key)
plt.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
plt.xlabel("X")
plt.ylabel("Y")
plt.title("H = 0.2 の場合の斜方投射")
# 横軸縦軸の表示範囲を調整する
plt.xlim(0, 1.5)
plt.ylim(0, 0.7)
plt.xticks(np.arange(0, 1.4, 0.1))
plt.gca().set_aspect('equal')
# 凡例を表示
plt.legend();
最大到達距離となる角度 $\theta$ を調べるために，着地点付近を拡大表示します。
H = 0.2
for th in range(47, 37, -1):
T_list = np.linspace(0, T1(th, H), 100+1)
key = 'θ = %2d°' % th
plt.plot(X(T_list, th), Y(T_list, th, H), label = key)
plt.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
plt.xlabel("X")
plt.ylabel("Y")
plt.title("H = 0.2 の場合の斜方投射")
# 着地点付近を拡大表示
plt.xlim(1.15, 1.19)
plt.ylim(0, 0.02)
# 凡例を表示
plt.legend();
水平到達距離が最大となるのは $\theta=45^{\circ}$ のときではない！ことがわかるだろう。
水平到達距離 $L$
定量的に調べるために，水平到達距離 L(th, H) をグラフにしてみる。
H = 0.2
th_list = np.linspace(38, 48, 101)
plt.plot(th_list, L(th_list, H))
plt.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
plt.xlabel("打ち出し角度 (°)")
plt.ylabel("規格化された水平到達距離")
plt.title("H = 0.2 の場合の打ち出し角度と水平到達距離")
# 横軸縦軸の表示範囲を調整する
plt.xlim(38, 48)
plt.xticks(np.arange(38, 48, 1));
上のグラフをみると，$H=0.2$ の場合に水平到達距離が最大となるのは $\theta=45^{\circ}$ のときではなく… どのくらい？
L が最大となる角度を数値微分で求める
数値微分して $\displaystyle \frac{dL(\theta, H)}{d\theta} = 0$ となる角度 $\theta$ を求めてみる。
$\theta \rightarrow x$ として，
$$\frac{dL(x, H)}{dx} \simeq \frac{L(x + \delta) – L(x – \delta)}{2 \delta}, \quad |\delta| << 1$$
# とりあえず，以下のように定義しておく。
# 数値微分: delta は小さいほどよい，というわけではない
# scipy.optimize.root_scalar を使うときは
# x の関数として定義
def dL(x, H, delta):
return (L(x+delta, H) - L(x-delta, H))/(2*delta)
$dL(x) = 0$ の解を $x$ について数値的に求めるには，scipy.optimize.root_scalar を使います。
H = 0.2
# delta を変えてみる
for delta in [1.e-3, 1.e-4, 1.e-5, 1.e-6, 1.e-7]:
# x 以外の変数は args に
print(root_scalar(dL, bracket = [38, 48], args=(H, delta)).root)
上の結果から，delta = 1.e-5 にして以下のような関数を定義します。
def thmax(H):
""" 高さ H からの斜方投射で最大到達距離になる角度（°）
"""
delta = 1.e-5
ans = root_scalar(dL, bracket = [1, 89], args=(H, delta)).root
return ans
thmax(0.2)
定義した thmax(H) を使って，H = 0.2 のときの最大到達距離となる軌道のグラフを描きます。
# 上の結果を使ってグラフを描く
H = 0.2
# 時間 T の範囲。滞空時間 T1 まで。
T_list = np.linspace(0, T1(thmax(H), H), 100)
# 凡例の数値の桁数を揃える
key = 'H=%3.1f, θ=%6.3f, Lmax=%6.3f' % (H, thmax(H), L(thmax(H), H))
plt.plot(X(T_list, thmax(H)), Y(T_list, thmax(H), H), label = key)
# グラフのオプション設定例
# グリッドをドットで
plt.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
plt.xlabel("X")
plt.ylabel("Y")
plt.title("高さ H からの斜方投射の最大水平到達距離")
# 横軸縦軸の表示範囲
plt.xlim(0, 1.2)
plt.ylim(0, 0.5)
plt.xticks(np.arange(0, 1.3, 0.1))
plt.gca().set_aspect('equal')
# 凡例を表示
plt.legend();
問題 1
高さ $H$ からの斜方投射の水平到達距離が最大となる打ち出し角度と，そのときの最大水平到達距離がわかるようなグラフを次の $H$ についてまとめて描け。
$H = 0.2, 0.4, 0.6, 0.8, 1.0$
解答例
Hl = np.linspace(0.2, 1, 5)
Hl
for H in reversed(Hl):
# 時間 T の範囲。滞空時間 T1 まで。
T_list = np.linspace(0, T1(thmax(H), H), 100)
# 凡例の数値の桁数を揃える
key = 'H=%3.1f, θ=%6.3f, Lmax=%6.3f' % (H, thmax(H), L(thmax(H), H))
plt.plot(X(T_list, thmax(H)), Y(T_list, thmax(H), H), label = key)
# グラフのオプション設定例
# グリッドをドットで
plt.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
plt.xlabel("X")
plt.ylabel("Y")
plt.title("高さ H からの斜方投射の最大水平到達距離")
# 横軸縦軸の表示範囲
plt.xlim(0, 2.4)
plt.ylim(0, 1.5)
#plt.xticks(np.arange(0, 1.3, 0.1))
plt.gca().set_aspect('equal')
# 凡例を表示
plt.legend();
問題 2
高さ $H$ を横軸に，そのときの水平到達距離が最大となる打ち出し角度 $\theta$ を縦軸にしたグラフを描け。
解答例
# 横軸は H
# (0, 1) 区間を 100等分
Hlist = np.linspace(0, 1, 101)
# 縦軸は thmax(H)
thmaxlist = []
for H in Hlist:
thmaxlist += [thmax(H)]
plt.plot(Hlist, thmaxlist)
# グラフのオプション設定例
# グリッドをドットで
plt.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
plt.xlabel("高さ H")
plt.ylabel("θmax (°)")
plt.title("高さ H からの斜方投射の最大水平到達距離を与える角度");
問題 3
高さ $H$ を横軸に，そのときの最大水平到達距離を縦軸にしたグラフを描け。
解答例
# 横軸は H
# (0, 1) 区間を 100等分
Hlist = np.linspace(0, 1, 101)
# 縦軸は L(thmax(H), H)
thmaxlist = []
for H in Hlist:
thmaxlist += [thmax(H)]
plt.plot(Hlist, L(thmaxlist, Hlist))
# グラフのオプション設定例
# グリッドをドットで
plt.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
plt.xlabel("高さ H")
plt.ylabel("最大水平到達距離 L")
plt.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)$
解答例
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))