月面における斜方投射

平面内の運動

今回は月面で下図のような斜方投射実験を行なうことにします。
 
初速度 $v_0$、角度 $\theta$ で打ち上げられたボールの軌道を求めることが目的です。
力学においては二次元運動（平面運動）を $x$ 方向と $y$ 方向の一次元運動に分解して考えることができるので、2 次元になったからといって特に何か新しいことを学ぶ必要はなく、前回までに使っていたアルゴリズムを、ほとんどそのまま適用できます。唯一行なうべき仕事は ベクトル 形式の運動方程式
 \[m\frac{d^2 \boldsymbol{r}}{dt^2}=\boldsymbol{F}\tag{1}\]
を成分ごとに書き表すことです。左辺の $\boldsymbol{r}$ は投射された物体の位置ベクトルです。
 \[\boldsymbol{r}=\begin{bmatrix}x\\y\end{bmatrix}\tag{2}\]
物体にはたらく力のベクトルは
 \[\boldsymbol{F}=\begin{bmatrix}0\\-mg\end{bmatrix}\tag{3}\]
と表せます (重力は $y$ 軸方向に下向きにはたらきます)。すなわち、運動方程式 (1) を成分ごとに書き下すと
 \[\begin{align*}&m\frac{d^2 x}{dt^2}=0\tag{4}\\[6pt]&m\frac{d^2 y}{dt^2}=-mg\tag{5}\end{align*}\]
となります。両辺から $m$ を取り除くと
 \[\begin{align*}&\frac{d^2 x}{dt^2}=0\tag{6}\\[6pt]&\frac{d^2 y}{dt^2}=-g\tag{7}\end{align*}\]
となります。(6) は加速度がないので等速直線運動を表しています。(7) は前回に扱った落下運動と全く同じ式です。ただし、原点から鉛直斜方に打ち上げているので初期条件が異なります。
 
初速度ベクトルは三角関数を使って
 \[\boldsymbol{v}_0=\begin{bmatrix}|\boldsymbol{v}_0|\cos\theta\\|\boldsymbol{v}_0|\sin\theta\end{bmatrix}\tag{8}\]
のように表せます。時刻 $t$ における物体の $x$ 座標は
 \[x=(|\boldsymbol{v}_0|\cos\theta)t\tag{9}\]
と表されます。$y$ 軸方向には初速度 $|\boldsymbol{v}_0|\sin\theta$ で打ち上げられて、徐々に失速しながらある高度まで達した後、今度は落下に転じます。実際、(7) を積分して初期条件を当てはめると、
 \[y=(|\boldsymbol{v}_0|\sin\theta)t-\frac{1}{2}gt^2\tag{10}\]
という式が得られます。(9) と (10) から $t$ を消去すると、ボールの軌道の式が得られます。
 \[y=(\tan\theta)x-\frac{g}{2|\boldsymbol{v}_0|^2\cos^2\theta}x^2\tag{10}\]
すなわち、$y$ は $x$ の二次関数であり、軌跡は上に凸の放物線となります。

斜方投射シミュレーション

解析解は得られましたが、プログラミングの練習のために、今回も敢えてオイラー法で数値解を求めてみましょう。前回記事 で作成したファイルに続けてコードを記述してください。最初に初速度と角度を指定して時間と座標、速度の配列データを得る projectile_motion() 関数を定義します。この関数の内部で前回作成した euler_method_2ord() を呼び出します。

# In[6]

# 投射関数
# t0:初期時刻
# v0:初速度の絶対値
# r0:投射地点の座標
# theta:仰角
# g:重力加速度
def projectile_motion(t0, v0, r0, tmax, theta, g=9.8, n=12):

    # 初速度のx成分とy成分
    vx0 = v0 * np.cos(theta)
    vy0 = v0 * np.sin(theta)

    # 投射地点の座標
    x0 = r0[0]
    y0 = r0[1]

    # x方向の運動方程式をオイラー法で解く
    t, vx, x = euler_method_2ord(t0, vx0, x0, tmax,
                                 lambda t,y,v:0, n=n)

    # y方向の運動方程式をオイラー法で解く
    t, vy, y = euler_method_2ord(t0, vy0, y0, tmax,
                                 lambda t,y,v:-g, n=n)

    return t, x, y, vx, vy

秒速 $30$ m (時速 $108$ km)、仰角 $60^{\circ}$ で打ち出されたボールの軌道を Matplotlib で描いてみます。projectile_motion() の引数 g は地球上を想定してデフォルトで 9.8 となっているので、1.62 を指定しておきます。

# In[7]

# 初期時刻
t0 = 0

# 初速度の絶対値
v0 = 30

# 投射地点の座標
r0 = (0, 0)

# 時刻の最大値
tmax = 150

# 投射角度(60°)
theta = np.pi / 3

# スライスでt,x,yのデータを取得
t, x, y = projectile_motion(t0, v0, r0, tmax, theta, g=1.62, n=8)[0:3]

# y>0のデータを抽出
t, x, y = t[y >= 0], x[y >= 0], y[y >= 0]

# 斜方投射された物体の軌道を散布図にプロット
fig, ax = plt.subplots(figsize=(8, 4))
ax.set_title("Projectile motion on the moon", fontsize=16, pad=10)
ax.grid()
ax.set_xlabel("x [m]", fontsize=14, labelpad=8)
ax.set_ylabel("y [m]", fontsize=14, labelpad=8)
ax.scatter(x, y, color="blue", s=2)

最高到達高度は $213$ m、水平到達距離は $486$ m です：

# In[8]

print("最高到達高度: {:.3f}[m]".format(max(y)))
print("水平到達距離: {:.3f}[m]".format(max(x)))
# 最高到達高度: 213.414[m]
# 水平到達距離: 486.328[m]

時速 $108$ km はプロ野球選手が投げる球よりずっと遅い速度ですが、それでも月面ではこれだけの距離を飛びます。

水平到達距離を最大にする仰角

仰角を大きくとればボールの滞空時間は長くなり、逆に小さくすれば速度の水平成分が大きくなります。この兼ね合いによって水平到達距離は決まります。では、どのくらいの角度で投射すればボールを一番遠くへ飛ばせるでしょうか？　水平到達距離を最大にする仰角を数値計算で探してみましょう。まずは大体の様子を見るために、角度を大まかに刻んで軌跡をプロットしてみます。

# In[9]

# FigureとAxesの設定
fig, ax = plt.subplots(figsize=(8, 4))
ax.set_title("Projectile motion on the moon", fontsize=16, pad=10)
ax.grid()
ax.set_xlabel("x [m]", fontsize=14, labelpad=8)
ax.set_ylabel("y [m]", fontsize=14, labelpad=8)

# 0～π/2を9等分
theta = np.linspace(0, np.pi/2, 9)

# 仰角を少しずつ変えながら軌跡をプロット
for k in theta:
    t, x, y = projectile_motion(t0, v0, r0, tmax, k, g=1.62, n=8)[0:3]
    t, x, y = t[y >= 0], x[y >= 0], y[y >= 0]
    ax.plot(x, y, label="{}°".format(np.rad2deg(k)))

# 凡例を表示
ax.legend()

出力された図を見ると、おおよそ $45^{\circ}$ で水平到達距離が最大になっているようです。次に $45^{\circ}$ 付近を細かく刻んで水平到達距離が最大となる仰角を探索してみます。精度を上げるために、オイラー法の時間分割数 n を大きくとります。

# In[10]

x_max = 0
theta_max = 0

# 仰角π/4付近のデータ
theta = np.linspace(np.pi/4 - 0.165, np.pi/4 + 0.165, 2**5+1)

# 水平到達距離が最大となる仰角を探索
for k in theta:
    x, y = projectile_motion(t0, v0, r0, tmax, k, g=1.62, n=16)[1:3]
    x = x[y >= 0]
    x = np.max(x)
    if x > x_max:
        x_max = x
        theta_max = np.rad2deg(k)

print("水平到達距離を最大にする仰角: {:.2f}[deg]".format(theta_max))
# 水平到達距離を最大にする仰角: 45.00[deg]

確かに仰角 $45^{\circ}$ で水平到達距離が最大になるようです。正確には変数 theta の刻み幅が不確かさとなります。

[付録] 投射物体のアニメーション

おまけです。ゲーム作成用ライブラリ Pygame Zero を使って、斜方投射のアニメーション・プログラムを作成してみました。興味のある方は、Pygame Zero をインストールして、IDLE などで次のコードを記述したpyファイルを実行してみてください。100 ピクセル (表示されるウィンドウの縦幅の 1/5) が 1m に相当します。月面上のゆったりとしたボールの動きを観察できます。

# Projectile

import pgzrun
import math

# スクリーンのサイズを設定
WIDTH = 900
HEIGHT = 500

# 月面の重力加速度[m/s^2]
g = 1.62

# 描画スケール(1m当たりのピクセル)
scale = 100

# 投射物体クラス
class Projectile:

    def __init__(self, mass, v0, angle,
                 radius=0.1, color='RED', dt=1/60):
        self.mass = mass
        self.angle = math.radians(angle)
        self.radius = radius
        self.dt = dt
        self.color = color
        self.v0 = v0
        self.vx = v0 * math.cos(self.angle)
        self.vy = v0 * math.sin(self.angle)
        self.x = 0
        self.y = 0

    # 微小時間における投射体の位置変更メソッド
    def change_pos(self):
        self.x += self.vx * self.dt
        self.y += self.vy * self.dt
        self.vy += - g * self.dt
        self.position = (self.x, self.y)

    # 投射体の描画メソッド
    def motion(self):
        scale_x = 10 + scale * self.position[0]
        scale_y = 400 - scale * self.position[1]

        scale_position = (scale_x, scale_y)
        scale_radius = int(scale * self.radius)
        screen.draw.filled_circle(scale_position,
                                  scale_radius, self.color)

        # x座標が一定値を超えたらリセットして再描画
        if scale_x > 1200:
            self.x = 0
            self.y = 0
            self.vy = self.v0 * math.sin(self.angle)

# 投射物体オブジェクトを生成
# 質量1kg,初速度10m/s,仰角60°
ball = Projectile(1, 4, 60)

# スクリーンとオブジェクトの描画
def draw():
    screen.clear()
    ball.motion()

# スクリーン更新
def update():
    ball.change_pos()

pgzrun.go()