必要なモジュールの import
必要なモジュールを import します。
# NumPy の関数を使います
import numpy as np
# 方程式の数値解
from scipy.optimize import root_scalar
# Matplotlib でグラフを描きます
import matplotlib.pyplot as plt
# グラフを SVG で Notebook にインライン表示させる設定
%config InlineBackend.figure_formats = ['svg']
空気抵抗がある場合の斜方投射の解
空気抵抗がある場合の無次元化された斜方投射の解は，
\begin{eqnarray}
X &=& \frac{1-e^{-b T}}{b} \, \cos\theta \\
Y &=& H + \frac{1-e^{-b T}}{b} \, \sin\theta + \frac{1 – b T -e^{-b T}}{b^2}
\end{eqnarray}
# 角度 th は度で
def X(T, th, b):
theta = np.radians(th)
return (1-np.exp(-b*T))/b * np.cos(theta)
# T -> x に。あとで root_scalar で使うため。
def Y(x, th, H, b):
theta = np.radians(th)
return H + (1-np.exp(-b*x))/b * np.sin(theta) + (1-b*x-np.exp(-b*x))/b**2
滞空時間 $T_1$
$T = 0$ で打ち出された粒子が地面 $Y=0$ に落ちるまでの滞空時間 $T_1$ は以下を満たす。
\begin{eqnarray}
Y(T_1, \theta, H, b) &=& H + \frac{1-e^{-b T_1}}{b} \, \sin\theta + \frac{1 – b T_1 -e^{-b T_1}}{b^2} = 0
\end{eqnarray}
この式は $T_1$ について超越方程式であり，解析的に解を求めることはできない。
そこで Python で方程式の数値解を求めるために，scipy.optimize.root_scalar() を使う。
水平到達距離 $L$
$Y(T_1, \theta, H, b) =0$ の数値解 $T_1(\theta, H, b)$ が求まったら，この滞空時間の間の水平方向の到達距離 $L$ は，
$$L = X(T_1(\theta, H, b), \theta, b)$$
def T1(th, H, b):
# どの範囲で root_scalar するか。ここでは 0.1 ~ 3.0 と広めに
ans = root_scalar(Y, bracket = [0.1, 3], args=(th, H, b))
return ans.root
def L(th, H, b):
return X(T1(th, H, b), th, b)
$H = 0$ の場合の軌道
空気抵抗の係数 $b=0.5$ の場合に，$\theta = 45^{\circ}$ 前後の角度での斜方投射の軌道をグラフにしてみる。
とりあえずは $\theta = 45^{\circ}$ の場合を1本，グラフにしてみる。
th = 45
H = 0
b =0.5
# 時間 T の範囲。滞空時間 T1 まで。
T_list = np.linspace(0, T1(th, H, b), 100)
plt.plot(X(T_list, th, b), Y(T_list, th, H, b));
オプションを設定してそれらしいグラフに。
# θ = 45° のグラフ
# 必要なパラメータの数値設定
th = 45
H = 0
b = 0.5
# 時間 T の範囲。滞空時間 T1 まで。
T_list = np.linspace(0, T1(th, H, b), 100)
# 凡例の数値の桁数を揃える
key="θ = %2d°" % th
plt.plot(X(T_list, th, b), Y(T_list, th, H, b), label = key)
# グラフのオプション設定例
# グリッドをドットで
plt.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
plt.xlabel("X")
plt.ylabel("Y")
plt.title("b = %.1f, H = %.1f の場合の斜方投射" % (b, H))
# 横軸縦軸の表示範囲
plt.xlim(0, 0.8)
plt.ylim(0, 0.4)
plt.xticks(np.arange(0, 0.9, 0.1))
plt.yticks(np.arange(0, 0.5, 0.1))
# 縦横の比，アスペクト比
plt.gca().set_aspect('equal')
# 判例を表示
plt.legend();
$\theta$ の値を $38^{\circ}$ から $45^{\circ}$ まで $1^{\circ}$ 刻みで変えて，複数の曲線を一度に描く例。
for で繰り返し処理させてみます。
なお，細かいことですが，グラフの最高点は $\theta$ の値が大きいときですので，先に $45^{\circ}$ から描きはじめて，$1^{\circ}$ ずつ減らしてみます。
# for で 1 ずつ減らしていく例 1
for i in range(5, 0, -1):
print(i, end=' ')
# for で 1 ずつ減らしていく例 2
for i in reversed(range(6)):
print(i, end=' ')
H = 0
b = 0.5
# 複数のグラフを一度に描く場合は...
# 以下のように for ループで繰り返し処理させればよい。
for th in range(45, 37, -1):
# 時間 T の範囲。滞空時間 T1 まで。
T_list = np.linspace(0, T1(th, H, b), 100)
# 凡例の数値の桁数を揃える
key="θ = %2d°" % th
plt.plot(X(T_list, th, b), Y(T_list, th, H, b),
label = key, linewidth = 0.8)
# グラフのオプション設定例
# グリッドをドットで
plt.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
plt.xlabel("X")
plt.ylabel("Y")
plt.title("b = %.1f, H = %.1f の場合の斜方投射" % (b, H))
# 横軸縦軸の表示範囲
plt.xlim(0, 0.8)
plt.ylim(0, 0.4)
plt.xticks(np.arange(0, 0.9, 0.1))
plt.yticks(np.arange(0, 0.5, 0.1))
plt.gca().set_aspect('equal')
# 判例を表示
plt.legend();
落下点付近を拡大表示します。
H = 0
b = 0.5
for th in range(45, 37, -1):
# 時間 T の範囲。滞空時間 T1 まで。
T_list = np.linspace(0, T1(th, H, b), 100)
# 凡例の数値の桁数を揃える
key="θ = %2d°" % th
plt.plot(X(T_list, th, b), Y(T_list, th, H, b),
label = key, linewidth = 0.8)
# グラフのオプション設定例
# グリッドをドットで
plt.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
plt.xlabel("X")
plt.ylabel("Y")
plt.title("b = %.1f, H = %.1f の場合の斜方投射" % (b, H))
# 落下点付近を拡大表示
plt.xlim(0.66, 0.68)
plt.ylim(0, 0.01)
# 判例を表示
plt.legend();
$b = 0.5$ の場合，最大水平到達距離となるのは $\theta$ が 39° から 40° くらいの時であることがわかります。これをもう少し定量的に求めてみます。
水平到達距離 $L$ のグラフ
定量的に調べるために，水平到達距離 L(th, H, b) をグラフにしてみる。
上で定義した L(th, H, b) では，引数 th の部分にリスト th_list を入れることができないので，少し工夫してみます。
H = 0
b = 0.5
th_list = np.linspace(38, 42, 100)
L_list = []
for th in th_list:
L_list += [L(th, H, b)]
plt.plot(th_list, L_list);
# グリッドをドットで
plt.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
plt.xlabel("投射角度 θ (°)")
plt.ylabel("規格化された水平到達距離 L")
plt.title("b = %.1f, H = %.1f の場合の投射角度と水平到達距離" % (b, H));
最大水平到達距離となる角度 $\theta_{\rm max}$ を数値微分で求める
$L(\theta, H, b)$ が最大となるのは $\displaystyle \frac{d L}{d \theta} = 0$ となる角度のときです。
$L(\theta, H, b)$ は解析的な関数として与えられていないので，ここでは以下のような数値微分で求めることにします。
$$\frac{d L}{d \theta} \simeq \frac{L(\theta+\delta, H, b)-L(\theta-\delta, H, b)}{2 \delta}$$
$\displaystyle \frac{d L}{d \theta} = 0$ となる角度 $\theta$ を方程式の解を数値的に解く root_scalar() を使って求めます。
# 数値微分: delta は小さいほどよい，というわけではない
# th -> x に
def dL(x, H, b, delta):
return (L(x+delta, H, b)-L(x-delta, H, b))/(2*delta)
b = 0.5
H = 0
# 方程式の数値解
# dL = 0 となる角度を delta を変えて求める
for delta in [1e-3, 1e-4, 1e-5, 1e-6, 1e-7, 1e-8]:
print(root_scalar(dL, bracket = [38, 42],args=(H,b,delta)).root)
上の結果から，delta = 1e-4 にして以下のような関数を定義します。
def thmax(H, b):
""" 高さ H からの斜方投射で最大到達距離になる角度（°）
"""
delta = 1e-4
ans = root_scalar(dL, bracket=[10,45], args=(H,b,delta)).root
return ans
$H=0, b = 0.5$ の場合
最大水平到達距離と投射角度
H = 0
b = 0.5
print('b = %.1f, H = %.1f のとき' % (b, H))
print('L が最大となる角度 θmax は %.6f°' % thmax(H, b))
print('そのときの最大水平到達距離 Lmax は %.6f' % L(thmax(H, b), H, b))
軌道のグラフ
# あらためて使う数値をまとめておく。
H = 0
b = 0.5
thm = thmax(H, b)
# 時間 T の範囲。滞空時間 T1 まで。
T_list = np.linspace(0, T1(thm, H, b), 100)
# 凡例の数値の桁数を揃える
key="b = %.1f, θ = %.4f°, Lmax = %.4f" % (b, thm, L(thm, H, b))
plt.plot(X(T_list, thm, b), Y(T_list, thm, H, b), label = key)
# グラフのオプション設定例
# グリッドをドットで
plt.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
plt.xlabel("X")
plt.ylabel("Y")
plt.title("H = %.0f の場合の斜方投射の最大水平到達距離" % H)
# 横軸縦軸の表示範囲
plt.xlim(0, 0.8)
plt.ylim(0, 0.4)
plt.xticks(np.arange(0, 0.9, 0.1))
plt.yticks(np.arange(0, 0.5, 0.1))
plt.gca().set_aspect('equal')
# 判例を表示
plt.legend();
問題 1
空気抵抗がある場合の高さ $H=0$ からの斜方投射の水平到達距離が最大となる投射角度と，そのときの最大水平到達距離がわかるようなグラフを次の規格化された抵抗係数 $b$ の場合についてまとめて描け。
$$ b = 0.5, \ 1.0, \ 1.5, \ 2.0$$
解答例
b_list = [0.5, 1.0, 1.5, 2.0]
H = 0
for b in b_list:
# 到達距離が最大となる投射角度
thm = thmax(H, b)
# 時間 T の範囲。滞空時間 T1 まで。
T_list = np.linspace(0, T1(thm, H, b), 100)
# 凡例の数値の桁数を揃える
key="b = %.1f, θmax = %.4f°, Lmax = %.4f" % (b, thm, L(thm, H, b))
plt.plot(X(T_list, thm, b), Y(T_list, thm, H, b), label = key)
# グラフのオプション設定例
# グリッドをドットで
plt.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
plt.xlabel("X")
plt.ylabel("Y")
plt.title("H = %.1f の場合の斜方投射の最大水平到達距離" % H)
# 横軸縦軸の表示範囲
plt.xlim(0, 0.8)
plt.ylim(0, 0.4)
plt.xticks(np.arange(0, 0.9, 0.1))
plt.yticks(np.arange(0, 0.5, 0.1))
plt.gca().set_aspect('equal')
# 判例を表示
plt.legend();
問題 2
問題 1 で，$H = 0$ の場合，空気抵抗がないとき，つまり $b=0$ のときには $\theta=45^{\circ}$ のときに最大水平到達距離になるが，このままでは b=0 の場合のグラフを追加しようとするとエラーになるはずである。
その原因をつきとめ，改善し，$b = 0, 0.5, 1.0, 1.5$ の場合のグラフをまとめて描け。
問題 3
空気抵抗 $b = 0.5$ の場合の高さ $H$ からの斜方投射の水平到達距離が最大となる投射角度と，そのときの最大水平到達距離がわかるようなグラフを次の $H$ の場合についてまとめて描け。
$$ H = 0.0, 0.2, 0.4, 0.6, 0.8, 1.0$$
解答例
# NumPy に頼らない方法
H_list = [i/10 for i in range(0, 11, 2)]
H_list
# NumPy に頼る方法 1
H_list = np.linspace(0, 1, 6)
H_list
# NumPy に頼る方法 2
H_list = np.arange(0, 1.1, 0.2)
H_list
b = 0.5
for H in reversed(H_list):
# 到達距離が最大となる投射角度
thm = thmax(H, b)
# 時間 T の範囲。滞空時間 T1 まで。
T_list = np.linspace(0, T1(thm, H, b), 100)
# 凡例の数値の桁数を揃える
# 上付下付添字は $ で挟んで LaTeX 形式で
key="H =%.1f, θ =%.4f°, Lmax =%.4f" % (H, thm, L(thm, H, b))
plt.plot(X(T_list, thm, b), Y(T_list, thm, H, b), label = key)
# グラフのオプション設定例
# グリッドをドットで
plt.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
plt.xlabel("X")
plt.ylabel("Y")
plt.title("b = %.1f の場合の斜方投射の最大水平到達距離" % b)
# 横軸縦軸の表示範囲
plt.xlim(0, 1.8)
plt.ylim(0, 1.1)
plt.xticks(np.arange(0, 1.8, 0.2))
plt.yticks(np.arange(0, 1.1, 0.2))
plt.gca().set_aspect('equal')
# 判例を表示
plt.legend(prop={"size": "small"});
