「高さ h からの斜方投射の問題を陰関数定理を使って解いてみる」のページで使っているので。
モジュールの import
In [1]:
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np
# グラフを SVG で Notebook にインライン表示
%config InlineBackend.figure_formats = ['svg']
plt.rcParams['mathtext.fontset'] = 'cm'
解の規格化
「高さ h からの斜方投射の問題を陰関数定理を使って解いてみる」にまとめているように，初期条件を $t = 0$ で
$$x = 0, \quad y = h, \quad v_x = v_0 \cos\theta, \quad v_y = v_0 \sin \theta$$
としたときの解は
\begin{eqnarray}
x(t, \theta) &=& v_0 \cos\theta\cdot t \\
y(t, \theta) &=& h + v_0 \sin\theta\cdot t -\frac{1}{2} g t^2
\end{eqnarray}
最大水平到達距離を与える $\theta_{\rm m}$ は
$$\theta_{\rm m} = \tan^{-1} \sqrt{\frac{v_0^2}{v_0^2 + 2 g h}}$$
そのときの滞空時間 $\tau_{\rm m}$ は
$$\tau_{\rm m} = \frac{v_0}{g \sin \theta_{\rm m}}$$
これらを系に特徴的な時間および長さで規格化すると，（要は $v_0 \rightarrow 1, \ g \rightarrow 1$ とすればよい）
\begin{eqnarray}
x(t, \theta) &\Rightarrow& \cos\theta\cdot t \\
y(t, \theta) &\Rightarrow& h + \sin\theta\cdot t -\frac{1}{2} t^2 \\
\theta_{\rm m} &\Rightarrow& \tan^{-1} \sqrt{\frac{1}{1 + 2 h}}\\
\tau_{\rm m} &\Rightarrow& \frac{1}{\sin \theta_{\rm m}}
\end{eqnarray}
関数の定義
$x(t, \theta), \ y(t, \theta, h), \ \theta_{\rm m}(h), \ \tau_{\rm m}(h)$ の定義。
In [2]:
def x(t, theta):
return np.cos(theta) * t
def y(t, theta, h):
return h + np.sin(theta) * t - t**2/2
def thetam(h):
return np.arctan(np.sqrt(1/(1 + 2*h)))
def taum(h):
return 1/np.sin(thetam(h))
グラフ例 1
In [3]:
fig = plt.figure(figsize=[6.4, 4.8])
ax = fig.add_subplot(aspect='equal')
#fig.tight_layout()
for i in range(0, 6):
h = 0.2*(5 - i)
thm = np.degrees(thetam(h))
Label = r'$h=%.1f,\ \theta_{\rm m}=%.1f^{\circ}$' % (h, thm)
t = np.linspace(0, taum(h))
plt.plot(x(t, thetam(h)), y(t, thetam(h), h), label=Label)
# 横軸縦軸の表示範囲
plt.xlim(0, 1.8)
plt.ylim(0, 1.2)
# グリッド（格子線）
plt.grid(ls=':')
# 軸ラベル，タイトル
plt.xlabel('$x$', fontsize=12)
plt.ylabel('$y$', fontsize=12)
plt.title('高さ $h$ からの斜方投射の最大到達距離', fontsize=12)
# x軸 y軸
plt.axhline(0, color='gray', ls = '--', lw=1)
plt.axvline(0, color='gray', ls = '--', lw=1);
plt.legend(fontsize=9);
グラフ例 2
In [4]:
fig = plt.figure(figsize=[6.4, 4.8])
ax = fig.add_subplot(aspect='equal')
fig.tight_layout()
# 軌道
h = 0.8
t = np.linspace(0, taum(h))
plt.plot(x(t, thetam(h)), y(t, thetam(h), h), c='blue')
# 高さ h
plt.plot([0, 0], [0, h], c='k')
ax.text(0.04, h/2, r"$h$",
fontsize="x-large", ha="center", va="center", c='k')
# 最大水平到達距離 ell_m
lm = x(taum(h),thetam(h))
plt.plot([0, lm], [0, 0], c='tab:red')
ax.text(lm/2-0.1, 0.09,
r"$\ell_{\rm m}=\frac{v_0^2}{g} \sqrt{1 + \frac{2gh}{v_0^2}}$",
fontsize="x-large", ha="center", va="center", c='tab:red')
# 初期位置
plt.scatter([0],[h], c='blue')
# 初速度ベクトル
v0x = 0.3*np.cos(thetam(h))
v0y = 0.3*np.sin(thetam(h))
plt.quiver([0], [h], [v0x], [v0y], color="blue", lw = 2, zorder=10,
angles='xy', scale_units='xy', scale=1
)
ax.text(v0x+0.04, h+v0y+0.04, r'$v_0$',
fontsize="x-large", ha="center", va="center", c='blue')
plt.plot([0,0.13],[h,h],lw=0.5,ls='--',c='blue')
# 角度 theta_i 部分
arci = patches.Arc(
xy=(0, h), width=0.2, height=0.2,
theta1=0, theta2=np.degrees(thetam(h)), ec="blue", lw=1)
ax.add_patch(arci)
ax.text(0.11, h+0.01,
r"$\theta_{\rm m} = \tan^{-1}\sqrt{\frac{v_0^2}{v_0^2 + 2 g h}}$",
fontsize="x-large", c='blue')
# 角度 theta_f 部分
arcf = patches.Arc(
xy=(lm, 0), width=0.2, height=0.2,
theta1=90+np.degrees(thetam(h)), theta2=180, ec="blue", lw=1)
ax.add_patch(arcf)
ax.text(1.23, 0.05, r"$\theta_{\rm f} = \frac{\pi}{2} - \theta_{\rm m}$",
fontsize="x-large", c='blue')
# 横軸縦軸の表示範囲
plt.xlim(-0.02, 1.8)
plt.ylim(-0.01, 1.2)
# タイトル
plt.title(r'高さ $h$ からの斜方投射：最大水平到達距離 $\ell_{\rm m}$ の場合の軌道')
# 目盛設定
ax.axis(False);
# x軸 y軸
plt.axhline(0, color='gray', ls = '--', lw=1, zorder=0)
plt.axvline(0, color='gray', ls = '--', lw=1, zorder=0);
グラフ例 3
In [5]:
fig = plt.figure(figsize=[6.4, 4.8])
ax = fig.add_subplot(aspect='equal')
fig.tight_layout()
# 軌道
h = 0.8
t = np.linspace(0, taum(h))
plt.plot(x(t, thetam(h)), y(t, thetam(h), h), c='blue',lw=2)
# 高さ h
#plt.plot([0, 0], [0, h], c='k')
ax.text(0, h/2, r"$h$",
fontsize="xx-large", ha="center", va="center", c='k')
plt.quiver([0], [0.45], [0], [h-0.45], color="k", zorder=10, width=0.003,
angles='xy', scale_units='xy', scale=1
)
plt.quiver([0], [0.35], [0], [-0.35], color="k", zorder=10, width=0.003,
angles='xy', scale_units='xy', scale=1
)
# 最大水平到達距離 ell_m
lm = x(taum(h),thetam(h))
#plt.plot([0, lm], [0, 0], c='tab:red')
ax.text(lm/2-0.1, 0,
r"$\ell_{\rm m}=\frac{v_0^2}{g} \sqrt{1 + \frac{2gh}{v_0^2}}$",
fontsize="xx-large", ha="center", va="center", c='tab:red')
plt.quiver([0.98], [0], [lm-0.98], [0], color="tab:red", zorder=10, width=0.003,
angles='xy', scale_units='xy', scale=1
)
plt.quiver([0.43], [0], [-0.43], [0], color="tab:red", zorder=10, width=0.003,
angles='xy', scale_units='xy', scale=1
)
# 初期位置
plt.scatter([0],[h], c='blue')
# 初速度ベクトル
v0x = 0.3*np.cos(thetam(h))
v0y = 0.3*np.sin(thetam(h))
plt.quiver([0], [h], [v0x], [v0y], color="blue", zorder=10,
angles='xy', scale_units='xy', scale=1
)
ax.text(v0x+0.04, h+v0y+0.04, r'$v_0$',
fontsize="x-large", ha="center", va="center", c='blue')
plt.plot([0,0.13],[h,h],lw=0.5,ls='--',c='blue')
# 角度 theta_i 部分
arci = patches.Arc(
xy=(0, h), width=0.2, height=0.2,
theta1=0, theta2=np.degrees(thetam(h)), ec="blue", lw=1)
ax.add_patch(arci)
ax.text(0.11, h+0.01,
r"$\theta_{\rm m} = \tan^{-1}\sqrt{\frac{v_0^2}{v_0^2 + 2 g h}}$",
fontsize="x-large", c='blue')
# 角度 theta_f 部分
arcf = patches.Arc(
xy=(lm, 0), width=0.2, height=0.2,
theta1=90+np.degrees(thetam(h)), theta2=180, ec="blue", lw=1)
ax.add_patch(arcf)
ax.text(1.23, 0.05, r"$\theta_{\rm f} = \frac{\pi}{2} - \theta_{\rm m}$",
fontsize="x-large", c='blue')
# 横軸縦軸の表示範囲
plt.xlim(-0.05, 1.8)
plt.ylim(-0.1, 1.05)
# タイトル
plt.title(r'高さ $h$ からの斜方投射：最大水平到達距離 $\ell_{\rm m}$ の場合の軌道')
# 目盛設定
ax.axis(False);
# x軸 y軸
plt.axhline(0, color='gray', ls = '--', c='lightgray', lw=0.5, zorder=0)
plt.axvline(0, color='gray', ls = '--', c='lightgray', lw=0.5, zorder=0);
