無次元化された解を書くと
\begin{eqnarray}
X &=& \cos\theta\cdot T \\
Y &=& H+ \sin\theta\cdot T – \frac{1}{2} T^2 \\
\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
Python で必要なパッケージを import します。
In [1]:
# NumPy も使います
import numpy as np
# Matplotlib でグラフを描きます
import matplotlib.pyplot as plt
# グラフを SVG で Notebook にインライン表示させる設定
%config InlineBackend.figure_formats = ['svg']
In [2]:
# 角度 th は度で
def X(th, T):
rth = np.radians(th)
return np.cos(rth) * T
def Y(th, T, H):
rth = np.radians(th)
return H + np.sin(rth)*T - T**2/2
In [3]:
def T1(th, H):
rth = np.radians(th)
return np.sin(rth) + np.sqrt(np.sin(rth)**2 + 2*H)
def L(th, H):
rth = np.radians(th)
return np.cos(rth) * T1(th, H)
$H = 0$ の場合の軌道
$\theta = 45^{\circ}$ 前後の角度での射方投射の軌道をグラフにしてみる。
In [4]:
# θ = 45° のグラフ
th = 45
H = 0
# T=0 から滞空時間 T1 までを 100 等分
T_list = np.linspace(0, T1(th, H), 100+1)
key = 'θ0 = 45°'
plt.plot(X(th, T_list), Y(th, T_list, 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() を使う例
plt.xticks(np.arange(0, 1.1, 0.1))
# 要素数を指定して np.linspace() を使う例
plt.yticks(np.linspace(0, 0.5, 5+1))
# 45°が45°に見えるように
# 縦軸横軸の比を 'equal' に。
plt.gca().set_aspect('equal')
# 判例を表示
plt.legend();
In [5]:
# 複数のグラフを一度に描く場合は...
key = 'θ0 = 47°'
plt.plot(X(47, T_list), Y(47, T_list, H), label = key)
key = 'θ0 = 46°'
plt.plot(X(46, T_list), Y(46, T_list, H), label = key)
key = 'θ0 = 45°'
plt.plot(X(45, T_list), Y(45, T_list, H), label = key)
key = 'θ0 = 44°'
plt.plot(X(44, T_list), Y(44, T_list, H), label = key)
key = 'θ0 = 43°'
plt.plot(X(43, T_list), Y(43, T_list, 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))
# 45°が45°に見えるように
# 縦軸横軸の比を 'equal' に。
plt.gca().set_aspect('equal')
# 判例を表示
plt.legend();
In [6]:
# 複数のグラフを一度に描く場合，
# for ループで繰り返し処理したらいいかも
for th in range(47, 42, -1):
# 数値を判例に入れる例
key = 'θ0 = '+str(th)+'°'
plt.plot(X(th, T_list), Y(th, T_list, 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();
In [7]:
# 着地点付近を拡大表示して 45° が最大到達距離となることを確認
for th in range(47, 42, -1):
key = 'θ0 = '+str(th)+'°'
plt.plot(X(th, T_list), Y(th, T_list, H), label = key)
plt.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
plt.xlabel("X")
plt.ylabel("Y")
plt.title("H = 0 の場合の射方投射")
# 横軸縦軸の表示範囲を調整する
plt.xlim(0.996, 1.002)
plt.ylim(0, 0.001)
# 判例を表示
plt.legend();
$H = 0.1$ の場合の軌道
In [8]:
H = 0.1
for th in range(47, 37, -1):
T_list = np.linspace(0, T1(th, H), 100+1)
key = 'θ0 = '+str(th)+'°'
plt.plot(X(th, T_list), Y(th, T_list, H), label = key)
plt.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
plt.xlabel("X")
plt.ylabel("Y")
plt.title("H = 0.1 の場合の射方投射")
# 横軸縦軸の表示範囲を調整する
plt.xlim(0, 1.4)
plt.ylim(0, 0.7)
plt.xticks(np.arange(0, 1.4, 0.1))
plt.gca().set_aspect('equal')
# 判例を表示
plt.legend();
In [9]:
# 着地点付近を拡大表示
H = 0.1
for th in range(47, 37, -1):
T_list = np.linspace(0, T1(th, H), 100+1)
key = 'θ0 = '+str(th)+'°'
plt.plot(X(th, T_list), Y(th, T_list, H), label = key)
plt.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
plt.xlabel("X")
plt.ylabel("Y")
plt.title("H = 0.1 の場合の射方投射")
# 横軸縦軸の表示範囲を調整する
plt.xlim(1.08, 1.1)
plt.ylim(0, 0.001)
# 判例を表示
plt.legend();
水平到達距離が最大となるのは $\theta=45^{\circ}$ のときではない！ことがわかるだろう。
水平到達距離 L
定量的に調べるために，水平到達距離 L(th, H) をグラフにしてみる。
In [10]:
H = 0.1
th_list = np.linspace(38, 48, 100)
plt.plot(th_list, L(th_list, H))
plt.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
plt.xlabel("打ち出し角度 (°)")
plt.ylabel("規格化された水平到達距離")
plt.title("H = 0.1 の場合の打ち出し角度と水平到達距離")
# 横軸縦軸の表示範囲を調整する
plt.xlim(38, 48);
L が最大となる角度を数値微分で求める
数値微分して $\displaystyle \frac{dL(\theta, H)}{d\theta} = 0$ となる角度 $\theta$ を求めてみる。
In [11]:
# 数値微分として scipy.misc.derivative() があるが
# deprecated らしい。
# とりあえず，以下のように定義しておく。
# 数値微分: h は小さいほどよい，というわけではない
# scipy.optimize.root_scalar を使うときは
# x の関数として定義
def dL(x, H, h):
return (L(x+h, H) - L(x-h, H))/(2*h)
In [12]:
# 方程式の数値解
from scipy.optimize import root_scalar
H = 0.1
# h を変えてみる
for h 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, h)).root)
In [13]:
# h はこの値でいってみる
h = 1.e-5
thmax = root_scalar(dL, bracket = [38, 48], args=(H, h)).root
print('H = %3.1f のとき，L が最大となる角度は %9.6f °' % (H, thmax))
print(' 最大到達距離は %9.6f' % L(thmax, H))
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 [14]:
def Lu(u, H):
return (u+sqrt((1+2*H)*u**2 + 2*H))/(1 + u**2)
In [15]:
# Python で微分するために SymPy を import
from sympy import *
In [16]:
# sympy で微分等に使う変数は symbols() で定義
x, u, H = symbols('x u h')
# .subs() が大事。これをつけないとどうなるかな。
def dLu(x, H):
return diff(Lu(u, H), u).subs(u, x)
# 微分できたか確認
dLu(x, H)
Out[16]:
In [17]:
H = 0.1
# 微分がゼロとなる u は root_scalar() で数値的に求める
# dLu(x)=0 となる x を [0, 1] の範囲で求める
# x 以外の変数は args に
ans = root_scalar(dLu, bracket = [0, 1], args=(H)).root
# u = tan(theta) より theta = arctan(u)
# さらに theta を度になおす
thmax = np.degrees(np.arctan(ans))
print('H = %3.1f のとき，L が最大となる角度は %15.12f °' % (H, thmax))
print(' 最大到達距離は %15.12f' % L(thmax, H))
In [18]:
# 上の結果を使ってグラフを描く
# 時間 T の範囲。滞空時間 T1 まで。
T_list = np.linspace(0, T1(thmax, H), 100)
# 凡例の数値の桁数を揃える
sH = '%3.1f' % H
sthmax = '%6.3f' % thmax
sLmax = '%6.3f' % L(thmax, H)
key = 'H = ' + sH + ', θ = ' + sthmax +', Lmax = ' + sLmax
plt.plot(X(thmax, T_list), Y(thmax, T_list, 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$
In [19]:
Hl = [float(i/10) for i in range(2,11,2)]
Hl
Out[19]:
解答例
In [20]:
# 到達距離が最大となる角度を追加するリスト
ansl = []
thml = []
for H in Hl:
ans = root_scalar(dLu, bracket = [0, 1], args=(H)).root
thmax = np.degrees(np.arctan(ans))
# 解をリスト ansl, thml に追加する
ansl.append(ans)
thml.append(thmax)
In [21]:
# 上の線から描きたいので reversed() に
for i in reversed(range(len(Hl))):
H = Hl[i]
thmax = thml[i]
T_list = np.linspace(0, T1(thmax, H), 100)
# 凡例の数値の桁数を揃える
sH = '%3.1f' % H
sthmax = '%6.3f' % thmax
sLmax = '%6.3f' % L(thmax, H)
key = 'H = ' + sH + ', θ = ' + sthmax +', Lmax = ' + sLmax
plt.plot(X(thmax, T_list), Y(thmax, T_list, H), label = key)
# グラフのオプション設定例
# グリッドをドットで
plt.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
plt.xlabel("X")
plt.ylabel("Y")
plt.title("高さ H からの斜方投射の最大水平到達距離")
# 横軸縦軸の表示範囲，目盛，アスペクト
plt.xlim(0, 2.0)
plt.ylim(0, 1.2)
plt.xticks(np.arange(0, 2.1, 0.2))
plt.yticks(np.arange(0, 2.1, 0.2))
plt.gca().set_aspect('equal')
# 判例を表示
plt.legend();
問題 2
高さ $H$ を横軸に，そのときの水平到達距離が最大となる打ち出し角度 $\theta$ を縦軸にしたグラフを描け。
In [22]:
# リスト Hl を横軸，thml を縦軸に
plt.scatter(Hl, thml);
このままだとあまりにショボいので， なめらかな曲線のグラフにしてみよ，というのが問題。
横軸 $H$ の範囲に関して，Hl は（以前のリストと区別するために Hl2 として）以下のように細かくしてみよう。
In [23]:
# 0 ~ 1 区間を 50 等分する
Hl2 = np.linspace(0, 1, 50+1)
Hl2
Out[23]:
解答例 1
上記の例の Hl の部分を Hl2 に変え，thml にかわる thml2 を求めて
plt.plot(Hl2, thml2)
とすればよいだろう。
In [24]:
ansl2 = []
thml2 = []
for H in Hl2:
ans = root_scalar(dLu, bracket = [0, 1], args=(H)).root
thmax = np.degrees(np.arctan(ans))
# 解をリスト ansl2, thml2 に追加する
ansl2.append(ans)
thml2.append(thmax)
# リスト Hl2 を横軸，thml2 を縦軸に
plt.plot(Hl2, thml2)
# グラフのオプション設定例
# グリッドをドットで
plt.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
plt.xlabel("規格化された高さ H")
plt.ylabel("打出し角度 (°)")
plt.title("高さ H からの斜方投射の水平到達距離を最大にする打出し角度")
# 横軸縦軸の表示範囲，目盛
plt.xlim(0, 1)
plt.ylim(26, 45)
plt.xticks(np.linspace(0,1,11))
plt.yticks(np.linspace(26,46,11));
解答例 2
$H$ を与えると水平到達距離が最大となる打出し角度 $\theta$ を与える関数を作成してみる。
In [25]:
def thetamax(H):
ans = root_scalar(dLu, bracket = [0, 1], args=(H)).root
thmax = np.degrees(np.arctan(ans))
return thmax
In [26]:
# H = 0 のときの角度が 45°となることを確認
thetamax(0)
Out[26]:
しかし，これでは thetamax(Hl2) のように引数にリストを入れるとエラーになるので，引数はリスト決め打ちで以下のような関数も定義してみる。
独り言：
この手間は，ひとえに
matplotlib.pyplot.plot()が陽関数 $y = f(x)$ をそのままグラフにする機能がなく，必ずplt.plot([x0, x1, ...], [y0, y1, ...])のようにリストを入れる必要があることに起因していると思うのだが…
In [27]:
def thetamaxl(Hl):
# 引数がリストの場合
thmaxl = []
for H in Hl:
ans = root_scalar(dLu, bracket = [0, 1], args=(H)).root
thmax = np.degrees(np.arctan(ans))
thmaxl.append(thmax)
return thmaxl
In [28]:
thetamaxl([0,1])
Out[28]:
In [29]:
plt.plot(Hl2, thetamaxl(Hl2))
# グラフのオプション設定例
# グリッドをドットで
plt.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
plt.xlabel("規格化された高さ H")
plt.ylabel("打出し角度 (°)")
plt.title("高さ H からの斜方投射の水平到達距離を最大にする打出し角度")
# 横軸縦軸の表示範囲，目盛
plt.xlim(0, 1)
plt.ylim(26, 46)
plt.xticks(np.linspace(0,1,11))
plt.yticks(np.linspace(26,46,11));
plt.savefig('phtosha11.svg');
参考：水平到達距離が最大になる角度
実は別ページで書いているように，規格化された高さ $H$ が与えられたとき，水平到達距離が最大となる角度 $\theta_{\rm max}$ は，以下の式で与えられることがわかっている。
$$\theta_{\rm max} = \arctan \frac{1}{\sqrt{1+2H}}$$
In [30]:
def themax(H):
th = np.arctan(1/np.sqrt(1 + 2*H))
return np.degrees(th)
あえて2つの関数をグラフにしてみると，完全に重なっていることがわかる。
In [31]:
plt.plot(Hl2, thetamaxl(Hl2), linewidth=3)
plt.plot(Hl2, themax(Hl2), linewidth=1);
問題 3
高さ $H$ を横軸に，そのときの最大水平到達距離を縦軸にしたグラフを描け。
In [32]:
# Hl と ansl に対応する Lu(u, H) の値をリストにする
Lmaxl = []
for i in range(len(Hl)):
Lmax = Lu(ansl[i], Hl[i])
Lmaxl.append(Lmax)
plt.scatter(Hl, Lmaxl);
このままだとあまりにショボいので， なめらかな曲線のグラフにしてみよ，というのが問題。
解答例
In [33]:
# Hl2 と ansl2 に対応する Lu(u, H) の値をリストにする
Lmaxl2 = []
for i in range(len(Hl2)):
Lmax = Lu(ansl2[i], Hl2[i])
Lmaxl2.append(Lmax)
plt.plot(Hl2, Lmaxl2)
# グラフのオプション設定例
# グリッドをドットで
plt.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
plt.xlabel("規格化された高さ H")
plt.ylabel("規格化された最大水平到達距離")
plt.title("高さ H からの斜方投射の最大水平到達距離")
# 横軸縦軸の表示範囲，目盛
plt.xlim(0, 1)
plt.ylim(1, 2)
plt.xticks(np.linspace(0,1,11))
plt.yticks(np.linspace(1,2,11));
plt.savefig('phtosha13.svg');
参考：最大水平到達距離
実は別ページで書いているように，規格化された高さ $H$ が与えられたときの最大水平到達距離は，規格化された量 $L_{\rm m}$ で以下の式で与えられることがわかっている。
$$L_{\rm m} = \sqrt{1+2H}$$
In [34]:
def Lmax(H):
return np.sqrt(1+2*Hl2)
あえて2つの関数をグラフにしてみると，完全に重なっていることがわかる。
In [35]:
plt.plot(Hl2, Lmaxl2, linewidth=3)
plt.plot(Hl2, Lmax(Hl2), linewidth=1);
plt.savefig('phtosha14.svg');
問題 4
$h = 5\,\mbox{(m)}$ の高さから速さ $v_0 = 10\,\mbox{(m/s)}$ で斜方投射したときの最大水平到達距離は何 $\mbox{m}$ か。またそのときの打出し角度は何度か。
ヒント：
$h = 5\,\mbox{(m)}$ は規格化された高さ $H$ ではいくらか，また規格化された最大到達距離 $L_{\rm max}(H)$ を次元をもった量に直すといくらか，ということ。
\begin{eqnarray}
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 [36]:
h = 5
v0 = 10
g = 9.80665
# 規格化された高さ H
H = h * g/v0**2
# Lu が最大，つまり dLu がゼロとなる u を求める
ans = root_scalar(dLu, bracket = [0, 1], args=(H)).root
xmax = v0**2/g * Lu(ans, H)
print('高さ h = %2d m のときの最大水平到達距離は' % h, end=' ')
print('%9.6f m' % xmax)
thmax = np.degrees(np.arctan(ans))
print('そのときの打ち出し角度は %9.6f °' % thmax)