地面から高さ $h$ の地点から空気抵抗なしの斜方投射を行うと，初速度の大きさを一定とした場合，水平方向の到達距離が最大となるのは打ち上げ角度（仰角）が何度のときかを Python の SciPy と NumPy を使って（かつ SymPy を使わずに）求める。
Python の Matplotlib によるグラフ作成には，plt.*** という関数のみを使った plt (pyplot) 流（pyplot インターフェースとも）と，ax.*** という関数のみを使った ax 流（オブジェクト指向インターフェースとも）の2つの方法があります。ここでは ax.*** のみを使ってグラフを描きます。
ax.*** のみを使ってグラフを作成する方法については，以下のページにまとめています。
update: 2025.01.21
ライブラリの 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']
運動方程式の解
初期条件を $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
\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
\end{eqnarray}
詳細は以下のページ
関数の定義
無次元化された解 $X(T, \theta), Y(T, \theta)$ および滞空時間 $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 X(T1(th, H), th)
$H = 0$ の場合の軌道
$\theta = 45^{\circ}$ 前後の角度での斜方投射の軌道をグラフにしてみる。
まずは，とりあえずのグラフ：
# ax を使う際の最初のおまじない
fig, ax = plt.subplots()
# θ = 45° のグラフ
th = 45
H = 0
# 滞空時間 T1
T_list = np.linspace(0, T1(th, H))
ax.plot(X(T_list, th), Y(T_list, th, H));
以下のようなオプションを設定して，もう少し体裁を整えます。
- $\theta$ の値を凡例に
- グリッド（格子線）をドットに
- 座標軸のラベルとグラフのタイトル
- 横軸縦軸の表示範囲
- 座標軸の目盛の間隔
- 縦軸横軸のアスペクト比
# ax を使う際の最初のおまじない
fig, ax = plt.subplots()
# θ = 45° のグラフ
th = 45
H = 0
# 滞空時間
T_list = np.linspace(0, T1(th, H))
# 凡例
key = 'θ = %2d°' % th
ax.plot(X(T_list, th), Y(T_list, th, H), label = key)
# グラフのオプション設定例
# グリッドをドットで
ax.grid(linestyle='dotted')
# 座標軸のラベル
ax.set_xlabel("X")
ax.set_ylabel("Y")
# グラフのタイトル
ax.set_title("H = 0 の場合の斜方投射")
# 横軸縦軸の表示範囲
ax.set_xlim(0, 1.1)
ax.set_ylim(0, 0.5)
# 座標軸の目盛の間隔を手動で設定する場合
# 間隔を指定して np.arange() を使う例
# 0 から 0.1 間隔で 1.1 未満（つまり 1.0）まで
ax.set_xticks(np.arange(0, 1.1, 0.1))
# 要素数を指定して np.linspace() を使う例
# (0, 0.5) 区間を 5 等分
ax.set_yticks(np.linspace(0, 0.5, 5+1))
# 45°が45°に見えるように
# 縦軸横軸のアスペクト比を 'equal' に。
ax.set_aspect('equal')
# 凡例を表示
ax.legend();
$\theta$ の値を $43^{\circ}$ から $47^{\circ}$ まで $1^{\circ}$ 刻みで変えて，複数の曲線を一度に描く例。
for で繰り返し処理させてみます。
なお，細かいことですが，グラフの最高点は $\theta$ の値が大きいときですので，先に $47^{\circ}$ から描きはじめて，$1^{\circ}$ ずつ減らしてみます。
# ax を使う際の最初のおまじない
fig, ax = plt.subplots()
H = 0
# 複数のグラフを一度に描く場合，
# for ループで繰り返し処理
for th in range(47, 42, -1):
# 滞空時間
T_list = np.linspace(0, T1(th, H))
# 凡例
key = 'θ = %2d°' % th
# plot でグラフを描く
ax.plot(X(T_list, th), Y(T_list, th, H), label = key)
ax.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_title("H = 0 の場合の斜方投射")
# 横軸縦軸の表示範囲と目盛
ax.set_xlim(0, 1.1)
ax.set_ylim(0, 0.5)
ax.set_xticks(np.arange(0, 1.1, 0.1))
ax.set_aspect('equal')
# 凡例を表示
ax.legend();
最大到達距離となる角度 $\theta$ を調べるために，着地点付近を拡大表示します。
# ax を使う際の最初のおまじない
fig, ax = plt.subplots()
H = 0
# 複数のグラフを一度に描く場合，
# for ループで繰り返し処理
for th in range(47, 42, -1):
# 滞空時間
T_list = np.linspace(0, T1(th, H))
# 凡例
key = 'θ = %2d°' % th
# plot でグラフを描く
ax.plot(X(T_list, th), Y(T_list, th, H), label = key)
ax.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_title("H = 0 の場合の斜方投射")
# 着地点付近を拡大表示
ax.set_xlim(0.99, 1.0)
ax.set_ylim(0, 0.01)
# 凡例を表示
ax.legend();
上のグラフから，$H=0$ の場合には確かに $\theta=45^{\circ}$ の時に最大到達距離になっていることがわかります。
$H = 0.2$ の場合の軌道
今度は，$H=0.2$ としてやってみます。
# ax を使う際の最初のおまじない
fig, ax = plt.subplots()
H = 0.2
for th in range(47, 37, -1):
T_list = np.linspace(0, T1(th, H))
key = 'θ = %2d°' % th
ax.plot(X(T_list, th), Y(T_list, th, H), label = key)
ax.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_title("H = 0.2 の場合の斜方投射")
# 横軸縦軸の表示範囲を調整する
ax.set_xlim(0, 1.5)
ax.set_ylim(0, 0.7)
ax.set_xticks(np.arange(0, 1.4, 0.1))
ax.set_aspect('equal')
# 凡例を表示
ax.legend();
最大到達距離となる角度 $\theta$ を調べるために，着地点付近を拡大表示します。
# ax を使う際の最初のおまじない
fig, ax = plt.subplots()
H = 0.2
for th in range(47, 37, -1):
T_list = np.linspace(0, T1(th, H))
key = 'θ = %2d°' % th
ax.plot(X(T_list, th), Y(T_list, th, H), label = key)
ax.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_title("H = 0.2 の場合の斜方投射")
# 着地点付近を拡大表示
ax.set_xlim(1.15, 1.19)
ax.set_ylim(0, 0.02)
# 凡例を表示
ax.legend();
水平到達距離が最大となるのは $\theta=45^{\circ}$ のときではない！ことがわかるだろう。
水平到達距離 $L$
定量的に調べるために，水平到達距離 L(th, H) をグラフにしてみる。
# ax を使う際の最初のおまじない
fig, ax = plt.subplots()
H = 0.2
th_list = np.linspace(38, 48)
ax.plot(th_list, L(th_list, H))
ax.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
ax.set_xlabel("打ち出し角度 (°)")
ax.set_ylabel("規格化された水平到達距離")
ax.set_title("H = 0.2 の場合の打ち出し角度と水平到達距離")
# 横軸縦軸の表示範囲を調整する
ax.set_xlim(38, 48)
ax.set_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=1e-5):
return (L(x+delta, H) - L(x-delta, H))/(2*delta)
$dL(x) = 0$ の解を $x$ について数値的に求めるには，scipy.optimize.root_scalar を使います。
H = 0.2
# delta がデフォルトの場合
print(root_scalar(dL, bracket = [38, 48], args=(H)).root)
# 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 thm(H):
""" 高さ H からの斜方投射で最大到達距離になる角度（°）
"""
ans = root_scalar(dL, bracket = [1, 89], args=(H)).root
return ans
定義した thm(H) を使って，H = 0.2 のときの最大到達距離となる軌道のグラフを描きます。
# 上の結果を使ってグラフを描く
# ax を使う際の最初のおまじない
fig, ax = plt.subplots()
H = 0.2
# 時間 T の範囲。滞空時間 T1 まで。
T_list = np.linspace(0, T1(thm(H), H), 100)
# 凡例の数値の桁数を揃える
key = 'H=%3.1f, θ=%6.3f°, Lmax=%6.3f' % (H, thm(H), L(thm(H), H))
ax.plot(X(T_list, thm(H)), Y(T_list, thm(H), H), label = key)
# グラフのオプション設定例
# グリッドをドットで
ax.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_title("高さ H からの斜方投射の最大水平到達距離")
# 横軸縦軸の表示範囲
ax.set_xlim(0, 1.2)
ax.set_ylim(0, 0.5)
ax.set_xticks(np.arange(0, 1.3, 0.1))
ax.set_aspect('equal')
# 凡例を表示
ax.legend();
○練習：$H$ を変えて最大水平到達距離を求める
高さ $H$ からの斜方投射の水平到達距離が最大となる打ち出し角度と，そのときの最大水平到達距離がわかるようなグラフを次の $H$ についてまとめて描け。
$H = 0.2, 0.4, 0.6, 0.8, 1.0$
Hs = np.linspace(0.2, 1, 5)
Hs
# ax を使う際の最初のおまじない
fig, ax = plt.subplots()
...
# グラフのオプション設定例
# グリッドをドットで
...
# 座標軸のラベルとタイトル
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_title("高さ H からの斜方投射の最大水平到達距離")
# 横軸縦軸の表示範囲
ax.set_xlim(0, 2.4)
ax.set_ylim(0, 1.5)
ax.set_aspect('equal')
# 凡例を表示
ax.legend();
○練習：高さ $H$ と打ち出し角度 $\theta_{\rm m}$ のグラフ
高さ $H$ を横軸に，そのときの水平到達距離が最大となる打ち出し角度 $\theta_{\rm m}$ を縦軸にしたグラフを描け。
# ax を使う際の最初のおまじない
fig, ax = plt.subplots()
# 横軸は H
# (0, 1) 区間
...
# 縦軸は thm(H)
...
...
# グラフのオプション設定例
# グリッドをドットで
ax.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
ax.set_title("水平到達距離が最大となる打ち出し角度");
ax.set_xlabel("規格化された高さ H")
ax.set_ylabel("打ち出し角度 (°)")
○練習：高さ $H$ と最大水平到達距離 $L_{\rm m}$ のグラフ
高さ $H$ を横軸に，そのときの最大水平到達距離 $L_{\rm m} = L(\theta_{\rm m}, H)$ を縦軸にしたグラフを描け。
# ax を使う際の最初のおまじない
fig, ax = plt.subplots()
...
...
...
# グラフのオプション設定例
# グリッドをドットで
ax.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
ax.set_title("最大水平到達距離");
ax.set_xlabel("規格化された高さ H")
ax.set_ylabel("規格化された高さ最大水平到達距離")
○練習：水平到達距離を求める
$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
...
...