Return to 空気抵抗がある場合の斜方投射を調べる準備

Python で空気抵抗がある場合の斜方投射を調べる

無次元化された解

空気抵抗がある場合の無次元化された解は，

\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}

詳細は以下のページ

パッケージの 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]:
# savefig した svg が 600 x 450 になるように 72 で割った値にしてみる
plt.rcParams["figure.figsize"] = (8.3334, 6.25)
# font size の設定
plt.rcParams["font.size"] = 12
In [3]:
# 角度 th は度で 
def X(T, b, th):
    rth = np.radians(th)
    return (1-np.exp(-b*T))/b * np.cos(rth)

# T -> x に。あとで root_scalar で使うため。
def Y(x, b, th, H):
    rth = np.radians(th)
    return H + (1-np.exp(-b*x))/b * np.sin(rth) + (1-b*x-np.exp(-b*x))/b**2

滞空時間と水平到達距離

$T = 0$ で打ち出された粒子が地面 $Y=0$ に落ちるまでの滞空時間 $T_1$ は $Y(T, b, \theta, H) = 0$ となる $T$ を数値的に解いて求める。

Python で方程式の数値解を求めるために，scipy.optimize.root_scalar() を使う。

この時間での水平方向の到達距離 $L$ は，

$$L = X(T_1, b, \theta)$$

In [4]:
from scipy.optimize import root_scalar

def T1(b, th, H):
    # どの範囲で root_scalar するか。ここでは 0.1 ~ 3.0 と広めに 
    ans = root_scalar(Y, bracket = [0.1, 3], args=(b, th, H))
    return ans.root

def L(b, th, H):
    return X(T1(b, th, H), b, th)

$H = 0$ の場合の軌道

空気抵抗の係数 $b=0.5$ の場合に，$\theta = 45^{\circ}$ 前後の角度での斜方投射の軌道をグラフにしてみる。

In [5]:
# θ = 45° のグラフ
# 必要なパラメータの数値設定
th1 = 45
b1 = 0.5
H1 = 0

# 時間 T の範囲。滞空時間 T1 まで。
T_list = np.linspace(0, T1(b1, th1, H1), 100)

# 凡例の数値の桁数を揃える
# 上付下付添字は $ で挟んで LaTeX 形式で
key="$θ_0$ = %2d°" % th1

plt.plot(X(T_list, b1, th1), Y(T_list, b1, th1, H1), label = key)

# グラフのオプション設定例
# グリッドをドットで
plt.grid(linestyle='dotted')

# 座標軸のラベルとタイトル
plt.xlabel("X")
plt.ylabel("Y")
plt.title("b = %.1f, H = %.1f の場合の斜方投射" % (b1, H1))

# 横軸縦軸の表示範囲
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();

In [6]:
# 複数のグラフを一度に描く例
b1 = 0.5
H1 = 0

# 角度 th1 を変えて以下のように羅列すれば描けるが...

# 時間 T の範囲。滞空時間 T1 まで。
T_list = np.linspace(0, T1(b1, 45, H1), 100)
# 凡例の数値の桁数を揃える
# 上付下付添字は $ で挟んで LaTeX 形式で
key="$θ_0$ = %2d°" % 45
plt.plot(X(T_list, b1, 45), Y(T_list, b1, 45, H1), 
         label = key, linewidth = 0.8)

# 時間 T の範囲。滞空時間 T1 まで。
T_list = np.linspace(0, T1(b1, 44, H1), 100)
# 凡例の数値の桁数を揃える
# 上付下付添字は $ で挟んで LaTeX 形式で
key="$θ_0$ = %2d°" % 44
plt.plot(X(T_list, b1, 44), Y(T_list, b1, 44, H1), 
         label = key, linewidth = 0.8)

# オプションの設定部分は同じ
# グリッドをドットで
plt.grid(linestyle='dotted')

# 座標軸のラベルとタイトル
plt.xlabel("X")
plt.ylabel("Y")
plt.title("b = %.1f, H = %.1f の場合の斜方投射" % (b1, H1))

# 横軸縦軸の表示範囲
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();

In [7]:
# 複数のグラフを一度に描く場合は...
# 以下のように for ループで繰り返し処理させればいいかも。
thmax = 45
thmin = 38

for th in range(45, 38-1, -1):
    # 時間 T の範囲。滞空時間 T1 まで。
    T_list = np.linspace(0, T1(b1, th, H1), 100)

    # 凡例の数値の桁数を揃える
    # 上付下付添字は $ で挟んで LaTeX 形式で
    key="$θ_0$ = %2d°" % th

    plt.plot(X(T_list, b1, th), Y(T_list, b1, th, H1), 
             label = key, linewidth = 0.8)

# グラフのオプション設定例
# グリッドをドットで
plt.grid(linestyle='dotted')

# 座標軸のラベルとタイトル
plt.xlabel("X")
plt.ylabel("Y")
plt.title("b = %.1f, H = %.1f の場合の斜方投射" % (b1, H1))

# 横軸縦軸の表示範囲
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();

落下点付近を拡大表示します。

In [8]:
for th in range(45, 38-1, -1):
    # 時間 T の範囲。滞空時間 T1 まで。
    T_list = np.linspace(0, T1(b1, th, H1), 100)

    # 凡例の数値の桁数を揃える
    # 上付下付添字は $ で挟んで LaTeX 形式で
    key="$θ_0$ = %2d°" % th

    plt.plot(X(T_list, b1, th), Y(T_list, b1, th, H1), 
             label = key, linewidth = 0.8)

# グラフのオプション設定例
# グリッドをドットで
plt.grid(linestyle='dotted')

# 座標軸のラベルとタイトル
plt.xlabel("X")
plt.ylabel("Y")
plt.title("b = %.1f, H = %.1f の場合の斜方投射" % (b1, H1))

# 表示範囲を狭くして拡大表示する
plt.xlim(0.66, 0.68)
plt.ylim(0, 0.01)

# 判例を表示
plt.legend();

$b = 0.5$ の場合，最大水平到達距離となるのは $\theta$ が 39° から 40° くらいの時であることがわかります。これをもう少し定量的に求めてみます。

水平到達距離 L

定量的に調べるために，水平到達距離 L(b, th, H) をグラフにしてみる。

上で定義した L(b, th, H) では，引数 th の部分にリスト th_list を入れることができないので，少し工夫してみます。

In [9]:
th_list = np.linspace(38, 42, 100)

# L(b1, th_list, H1) とはいかないのでちょっと工夫
# th_list のそれぞれの要素 th = th_list[i] に対する
# Lの値 L(b1, th, H1) を要素にもつリスト Ll を作成する
Ll = []
for th in th_list:
    tmp = L(b1, th, H1)
    Ll.append(tmp)
    
# グリッドをドットで
plt.grid(linestyle='dotted')
# 座標軸のラベルとタイトル
plt.xlabel("打ち出し角度 $θ_0$ (°)")
plt.ylabel("規格化された水平到達距離 L")
plt.title("b = %.1f, H = %.1f の場合の打ち出し角度と水平到達距離" % (b1, H1))

plt.plot(th_list, Ll);

最大水平到達距離を数値微分で求める

$L(b, \theta, H)$ が最大となるのは $\displaystyle \frac{\partial L}{\partial \theta} = 0$ となる角度のときです。

$L(b, \theta, H)$ は解析的な関数として与えられていないので，ここでは以下のような数値微分で求めることにします。

$$\frac{\partial L}{\partial \theta} \simeq \frac{L(b, \theta+h, H)-L(b, \theta-h, H)}{2 h}$$

In [10]:
# 数値微分: h は小さいほどよい，というわけではない
# th -> x に
def dL(x, h):
    return (L(b1, x+h, H1)-L(b1, x-h, H1))/(2*h)
In [11]:
# 方程式の数値解
# dL = 0 となる角度を h を変えて求める
for h in [1.e-3, 1.e-4, 1.e-5, 1.e-6, 1.e-7]:
    print(root_scalar(dL, bracket = [38, 42], args=(h)).root)

39.487609707424156
39.48760970413319
39.4876097064069
39.48760985949289
39.487609219856495
In [12]:
# 以後は h = 1e-5 として計算することにしよう 
h = 1e-5
thmax = root_scalar(dL, bracket = [20, 45], args=(h)).root

print('b = %.1f, H = %.1f のとき' % (b1, H1), end = ' ')
print('L ば最大となる角度は %.6f°' % thmax)
print('    最大距離は %.6f' % L(b1, thmax, H1))

b = 0.5, H = 0.0 のとき L ば最大となる角度は 39.487610°
    最大距離は 0.679421
In [13]:
# あらためて使う数値をまとめておく。
H1 = 0
b1 = 0.5
# 到達距離が最大となる打ち出し角度
thmax = root_scalar(dL, bracket = [20, 45], args=(h)).root

# 時間 T の範囲。滞空時間 T1 まで。
T_list = np.linspace(0, T1(b1, thmax, H1), 100)

# 凡例の数値の桁数を揃える
# 上付下付添字は $ で挟んで LaTeX 形式で
key="b =%.1f, $θ_0$ =%.4f°, $L_{max}$ =%.4f" % (b1, thmax, L(b1, thmax, H1))

plt.plot(X(T_list, b1, thmax), Y(T_list, b1, thmax, H1), label = key)

# グラフのオプション設定例
# グリッドをドットで
plt.grid(linestyle='dotted')

# 座標軸のラベルとタイトル
plt.xlabel("X")
plt.ylabel("Y")
plt.title("H = %.1f の場合の斜方投射の最大水平到達距離" % H1)

# 横軸縦軸の表示範囲
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$$

ヒント：

b1 = 0.5 の場合は上記のように，

b1 = 0.5
thmax = ...
T_list = ...
key = ...
plt.plot(...)

としたわけだから，b1 の値がリスト bl の要素の場合は

for b1 in bl:
    thmax = ...
    T_list = ...
    key = ...
    plt.plot(...)

のようにしてやればよい。

解答例

In [14]:
bl = [0.5, 1.0, 1.5, 2.0]
In [15]:
H1 = 0

# b1 = 0.5 の場合の上記の例を
for b1 in bl:
    thmax = root_scalar(dL, bracket = [20, 45], args=(h)).root
    # 時間 T の範囲。滞空時間 T1 まで。
    T_list = np.linspace(0, T1(b1, thmax, H1), 100)

    # 凡例の数値の桁数を揃える
    # 上付下付添字は $ で挟んで LaTeX 形式で
    key="b =%.1f,  $θ_0$ =%.4f°,  $L_{max}$ =%.4f" % (b1, thmax, L(b1, thmax, H1))

    plt.plot(X(T_list, b1, thmax), Y(T_list, b1, thmax, H1), label = key)

# グリッドをドットで
plt.grid(linestyle='dotted')

# 座標軸のラベルとタイトル
plt.xlabel("X")
plt.ylabel("Y")
plt.title("空気抵抗がある場合の斜方投射の最大水平到達距離")

# 横軸縦軸の表示範囲
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();

参考までに，グラフの上下左右枠に主目盛・副目盛をつけて見栄えをサイエンティフィックにしたい場合は，以下のように。

In [16]:
# 先に主目盛・副目盛の設定
#グラフの上下左右に目盛線を付ける。
f = plt.figure()
ax = f.add_subplot(111)
ax.yaxis.set_ticks_position('both')
ax.xaxis.set_ticks_position('both')

# x の大目盛を 1 刻みに
plt.xticks(np.arange(0, 0.9, 0.1))
# y の大目盛を 0.5 刻みに
plt.yticks(np.arange(0., 0.5, 0.1))
# 小目盛も表示
plt.minorticks_on()
# tick の向きを内側に
plt.tick_params(which ="both",direction = 'in')

# あとは同じ
H1 = 0

for b1 in bl:
    thmax = root_scalar(dL, bracket = [20, 45], args=(h)).root
    # 時間 T の範囲。滞空時間 T1 まで。
    T_list = np.linspace(0, T1(b1, thmax, H1), 100)

    # 凡例の数値の桁数を揃える
    # 上付下付添字は $ で挟んで LaTeX 形式で
    key="b =%.1f,  $θ_0$ =%.4f°,  $L_{max}$ =%.4f" % (b1, thmax, L(b1, thmax, H1))

    plt.plot(X(T_list, b1, thmax), Y(T_list, b1, thmax, H1), label = key)

# グリッドをドットで
plt.grid(linestyle='dotted')

# 座標軸のラベルとタイトル
plt.xlabel("X")
plt.ylabel("Y")
plt.title("空気抵抗がある場合の斜方投射の最大水平到達距離")

# 横軸縦軸の表示範囲
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

空気抵抗 $b = 0.5$ の場合の高さ $H$ からの斜方投射の水平到達距離が最大となる打ち出し角度と，そのときの最大水平到達距離がわかるようなグラフを次の $H$ の場合についてまとめて描け。

$$ H = 0.0, 0.2, 0.4, 0.6, 0.8, 1.0$$

ヒント：

H1 = 0.0 の場合は上記のように，

H1 = 0.0
thmax = ...
T_list = ...
key = ...
plt.plot(...)

としたわけだから，H1 の値がリスト Hl の要素の場合は

for H1 in Hl:
    thmax = ...
    T_list = ...
    key = ...
    plt.plot(...)

のようにしてやればいいだろう。また，以下では for H1 in reversed(Hl): としていますが，reversed() にしている理由はわかりますか？

解答例

In [17]:
Hl = [float(i/10) for i in range(0, 11, 2)]
Hl
Out[17]:
[0.0, 0.2, 0.4, 0.6, 0.8, 1.0]
In [18]:
b1 = 0.5

for H1 in reversed(Hl):
    thmax = root_scalar(dL, bracket = [10, 45], args=(h)).root
    # 時間 T の範囲。滞空時間 T1 まで。
    T_list = np.linspace(0, T1(b1, thmax, H1), 100)

    # 凡例の数値の桁数を揃える
    # 上付下付添字は $ で挟んで LaTeX 形式で
    key="H =%.1f, $θ_0$ =%.4f°, $L_{max}$ =%.4f" % (H1, thmax, L(b1, thmax, H1))

    plt.plot(X(T_list, b1, thmax), Y(T_list, b1, thmax, H1), label = key)

# グリッドをドットで
plt.grid(linestyle='dotted')

# 座標軸のラベルとタイトル
plt.xlabel("X")
plt.ylabel("Y")
plt.title("空気抵抗 b =%.1f の場合の斜方投射の最大水平到達距離" % b1)

# 横軸縦軸の表示範囲
plt.xlim(0, 2)
plt.ylim(0, 1.2)
plt.xticks(np.arange(0, 2.1, 0.2))
plt.yticks(np.arange(0, 1.3, 0.2))
plt.gca().set_aspect('equal')

# 判例を表示
plt.legend();

参考までに，グラフの上下左右枠に主目盛・副目盛をつけて見栄えをサイエンティフィックにしたい場合は，以下のように。

In [19]:
# 先に主目盛・副目盛の設定
#グラフの上下左右に目盛線を付ける。
f = plt.figure()
ax = f.add_subplot(111)
ax.yaxis.set_ticks_position('both')
ax.xaxis.set_ticks_position('both')

# x の大目盛を 1 刻みに
plt.xticks(np.arange(0, 0.9, 0.1))
# y の大目盛を 0.5 刻みに
plt.yticks(np.arange(0., 0.5, 0.1))
# 小目盛も表示
plt.minorticks_on()
# tick の向きを内側に
plt.tick_params(which ="both",direction = 'in')

# 以下は同じ
b1 = 0.5

for H1 in reversed(Hl):
    thmax = root_scalar(dL, bracket = [10, 45], args=(h)).root
    # 時間 T の範囲。滞空時間 T1 まで。
    T_list = np.linspace(0, T1(b1, thmax, H1), 100)

    # 凡例の数値の桁数を揃える
    # 上付下付添字は $ で挟んで LaTeX 形式で
    key="H =%.1f, $θ_0$ =%.4f°, $L_{max}$ =%.4f" % (H1, thmax, L(b1, thmax, H1))

    plt.plot(X(T_list, b1, thmax), Y(T_list, b1, thmax, H1), label = key)

# グリッドをドットで
plt.grid(linestyle='dotted')

# 座標軸のラベルとタイトル
plt.xlabel("X")
plt.ylabel("Y")
plt.title("空気抵抗 b =%.1f の場合の斜方投射の最大水平到達距離" % b1)

# 横軸縦軸の表示範囲
plt.xlim(0, 2)
plt.ylim(0, 1.2)
plt.xticks(np.arange(0, 2.1, 0.2))
plt.yticks(np.arange(0, 1.3, 0.2))
plt.gca().set_aspect('equal')

# 判例を表示
plt.legend();