平面内の運動
今回は月面で下図のような斜方投射実験を行なうことにします。
初速度
力学においては二次元運動(平面運動)を
を成分ごとに書き表すことです。左辺の
物体にはたらく力のベクトルは
と表せます (重力は
となります。両辺から
となります。(6) は加速度がないので等速直線運動を表しています。(7) は前回に扱った落下運動と全く同じ式です。ただし、原点から鉛直斜方に打ち上げているので初期条件が異なります。
初速度ベクトルは三角関数を使って
のように表せます。時刻
と表されます。
という式が得られます。(9) と (10) から
すなわち、
斜方投射シミュレーション
解析解は得られましたが、プログラミングの練習のために、今回も敢えてオイラー法で数値解を求めてみましょう。前回記事 で作成したファイルに続けてコードを記述してください。最初に初速度と角度を指定して時間と座標、速度の配列データを得る 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
秒速
# 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)
最高到達高度は
# In[8]
print("最高到達高度: {:.3f}[m]".format(max(y)))
print("水平到達距離: {:.3f}[m]".format(max(x)))
# 最高到達高度: 213.414[m]
# 水平到達距離: 486.328[m]
時速
水平到達距離を最大にする仰角
仰角を大きくとればボールの滞空時間は長くなり、逆に小さくすれば速度の水平成分が大きくなります。この兼ね合いによって水平到達距離は決まります。では、どのくらいの角度で投射すればボールを一番遠くへ飛ばせるでしょうか? 水平到達距離を最大にする仰角を数値計算で探してみましょう。まずは大体の様子を見るために、角度を大まかに刻んで軌跡をプロットしてみます。
# 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()
出力された図を見ると、おおよそ
# 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]
確かに仰角
[付録] 投射物体のアニメーション
おまけです。ゲーム作成用ライブラリ 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()
コメント
Pythonの情報を探していてこのサイトにたどり着きました。どの記事も分かりやすく大変ためになりました。ありがとうございます。
コメントありがとうございます。
そう言っていただけると私も嬉しいです。
これからも Python に関する色々な情報を提供していくつもりです。
下記は誤植と思われますので、ご確認ください。
(1), (4), (5), (6), (7) 式の分母で、dt → dt^2
直しました。
ありがとうございます。