地面から高さ $h$ の地点から空気抵抗なしの斜方投射を行うと，初速度の大きさを一定とした場合，水平方向の到達距離が最大となるのは打ち上げ角度（仰角）が何度のときかを Python の SciPy と NumPy を使って（かつ SymPy を使わずに）求める。

Python の Matplotlib によるグラフ作成には，plt.*** という関数のみを使った plt (pyplot) 流（pyplot インターフェースとも）と，ax.*** という関数のみを使った ax 流（オブジェクト指向インターフェースとも）の2つの方法があります。ここでは ax.*** のみを使ってグラフを描きます。

ax.*** のみを使ってグラフを作成する方法については，以下のページにまとめています。

Matplotlib でグラフ作成：ax 編

update: 2025.01.21

ライブラリの import

In [1]:

# 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, y (0) = h, v_{x} (0) = v_{0} \cos θ, v_{y} (0) = v_{0} \sin θ$

としたときの解は

$\begin{array}{rcl} x (t) & = & v_{0} \cos θ \cdot t \\ y (t) & = & h + v_{0} \sin θ \cdot t - \frac{1}{2} g t^{2} \end{array}$

無次元化された斜方投射の解

この系に特徴的な時間 $τ \equiv \frac{v_{0}}{g}$ および長さ $ℓ \equiv v_{0} τ = \frac{v_{0}^{2}}{g}$ で解を無次元化すると，

$\begin{array}{rcl} T & \equiv & \frac{t}{τ} \\ H & \equiv & \frac{h}{ℓ} \\ X & \equiv & \frac{x}{ℓ} = \cos θ \cdot T \\ Y & \equiv & \frac{y}{ℓ} = H + \sin θ \cdot T - \frac{1}{2} T^{2} \end{array}$

詳細は以下のページ

高さ h からの斜方投射の最大到達距離を求める準備

関数の定義

無次元化された解 $X (T, θ), Y (T, θ)$ および滞空時間 $T_{1}$ と水平到達距離 $L$ を関数として定義します。

In [2]:

# 角度 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$ の場合の軌道

$θ = 45^{\circ}$ 前後の角度での斜方投射の軌道をグラフにしてみる。

まずは，とりあえずのグラフ：

In [3]:

# 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));

以下のようなオプションを設定して，もう少し体裁を整えます。

$θ$ の値を凡例に
グリッド（格子線）をドットに
座標軸のラベルとグラフのタイトル
横軸縦軸の表示範囲
座標軸の目盛の間隔
縦軸横軸のアスペクト比

In [4]:

# 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();

$θ$ の値を $43^{\circ}$ から $47^{\circ}$ まで $1^{\circ}$ 刻みで変えて，複数の曲線を一度に描く例。

for で繰り返し処理させてみます。

なお，細かいことですが，グラフの最高点は $θ$ の値が大きいときですので，先に $47^{\circ}$ から描きはじめて， $1^{\circ}$ ずつ減らしてみます。

In [5]:

# 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();

最大到達距離となる角度 $θ$ を調べるために，着地点付近を拡大表示します。

In [6]:

# 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$ の場合には確かに $θ = 45^{\circ}$ の時に最大到達距離になっていることがわかります。

$H = 0.2$ の場合の軌道

今度は， $H = 0.2$ としてやってみます。

In [7]:

# 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();

最大到達距離となる角度 $θ$ を調べるために，着地点付近を拡大表示します。

In [8]:

# 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();

水平到達距離が最大となるのは $θ = 45^{\circ}$ のときではない！ことがわかるだろう。

水平到達距離 $L$

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

In [9]:

# 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$ の場合に水平到達距離が最大となるのは $θ = 45^{\circ}$ のときではなく… どのくらい？

L が最大となる角度を数値微分で求める

数値微分して $\frac{d L (θ, H)}{d θ} = 0$ となる角度 $θ$ を求めてみる。

$θ \to x$ として，

$\frac{d L (x, H)}{d x} ≃ \frac{L (x + δ) - L (x - δ)}{2 δ}, | δ | << 1$

In [10]:

# とりあえず，以下のように定義しておく。
# 数値微分: delta は小さいほどよい，というわけではない 
# scipy.optimize.root_scalar を使うときは
# x の関数として定義
def dL(x, H, delta=1e-5):
    return (L(x+delta, H) - L(x-delta, H))/(2*delta)

$d L (x) = 0$ の解を $x$ について数値的に求めるには，scipy.optimize.root_scalar を使います。

In [11]:

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)

40.20296588835251
40.2029658845539
40.20296588569512
40.20296588835251
40.20296588448801
40.20296551268753

上の結果から，delta = 1.e-5 にして以下のような関数を定義します。

In [12]:

def thm(H):
    """ 高さ H からの斜方投射で最大到達距離になる角度（°）
    """
    ans = root_scalar(dL, bracket = [1, 89], args=(H)).root    
    return ans

定義した thm(H) を使って，H = 0.2 のときの最大到達距離となる軌道のグラフを描きます。

In [13]:

# 上の結果を使ってグラフを描く
# 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$

In [14]:

Hs = np.linspace(0.2, 1, 5)
Hs

Out[14]:

array([0.2, 0.4, 0.6, 0.8, 1. ])

In [15]:

# 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$ と打ち出し角度 $θ_{m}$ のグラフ

高さ $H$ を横軸に，そのときの水平到達距離が最大となる打ち出し角度 $θ_{m}$ を縦軸にしたグラフを描け。

In [16]:

# 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_{m}$ のグラフ

高さ $H$ を横軸に，そのときの最大水平到達距離 $L_{m} = L (θ_{m}, H)$ を縦軸にしたグラフを描け。

In [17]:

# ax を使う際の最初のおまじない
fig, ax = plt.subplots()

...
...
...


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

# 座標軸のラベルとタイトル
ax.set_title("最大水平到達距離");
ax.set_xlabel("規格化された高さ H")
ax.set_ylabel("規格化された高さ最大水平到達距離")

○練習：水平到達距離を求める

$h = 5 (m)$ の高さから速さ $v_{0} = 10 (m/s)$ で斜方投射したときの最大水平到達距離は何 $m$ か。またそのときの打出し角度は何度か。

ヒント：

$h = 5 (m)$ は規格化された高さ $H$ ではいくらか，また規格化された最大到達距離 $L_{m a x} (H)$ を次元をもった量に直すといくらか，ということ。

$\begin{array}{rcl} ℓ & = & \frac{v_{0}^{2}}{g} \\ H & = & \frac{h}{ℓ} = \frac{g h}{v_{0}^{2}} \\ x & = & ℓ X = \frac{v_{0}^{2}}{g} X \\ x_{m a x} & = & \frac{v_{0}^{2}}{g} L_{m a x} (H) \end{array}$

重力加速度 $g = 9.80665 (m / s^{2})$

In [18]:

h = 5
v0 = 10
g = 9.80665

...
...

高さ h =  5 m のときの最大水平到達距離は **.****** m
そのときの打ち出し角度は **.******°

SciPy で高さ h からの斜方投射の最大到達距離を求める

ライブラリの import

運動方程式の解

無次元化された斜方投射の解

関数の定義

$H = 0$ の場合の軌道

$H = 0.2$ の場合の軌道

水平到達距離 $L$

L が最大となる角度を数値微分で求める

○練習： $H$ を変えて最大水平到達距離を求める

○練習：高さ $H$ と打ち出し角度 $θ_{m}$ のグラフ

○練習：高さ $H$ と最大水平到達距離 $L_{m}$ のグラフ

○練習：水平到達距離を求める

このページの目次

このセクションの構成

SciPy で高さ h からの斜方投射の最大到達距離を求める

ライブラリの import

運動方程式の解

無次元化された斜方投射の解

関数の定義

H=0 の場合の軌道

H=0.2 の場合の軌道

水平到達距離 L

L が最大となる角度を数値微分で求める

○練習：H を変えて最大水平到達距離を求める

○練習：高さ H と打ち出し角度 θm のグラフ

○練習：高さ H と最大水平到達距離 Lm のグラフ

○練習：水平到達距離を求める

このページの目次

このセクションの構成

$H = 0$ の場合の軌道

$H = 0.2$ の場合の軌道

水平到達距離 $L$

○練習： $H$ を変えて最大水平到達距離を求める

○練習：高さ $H$ と打ち出し角度 $θ_{m}$ のグラフ

○練習：高さ $H$ と最大水平到達距離 $L_{m}$ のグラフ