放物運動は古典力学の最も基本的な問題の一つですが、空気抵抗を加えた途端に解析解が得られなくなり、数値計算が必要になります。実際の弾道計算、スポーツの球技解析、ロケットの軌道予測など、空気抵抗は無視できない要素です。
本記事の内容
- 真空中の放物運動（解析解）
- 速度に比例する抵抗（線形抵抗）
- 速度の2乗に比例する抵抗（2次抵抗）
- 最適発射角の変化
- Pythonでの弾道シミュレーション
前提知識
真空中の放物運動
初速 $v_0$、仰角 $\alpha$ のとき:
$$ x(t) = v_0\cos\alpha \cdot t, \quad y(t) = v_0\sin\alpha \cdot t – \frac{g}{2}t^2 $$
飛距離: $R = \frac{v_0^2\sin 2\alpha}{g}$、最大高さ: $H = \frac{v_0^2\sin^2\alpha}{2g}$
空気抵抗を考慮した運動方程式
2次抵抗（高速域）
$$ \bm{F}_{\text{drag}} = -\frac{1}{2}C_D \rho A |\bm{v}| \bm{v} $$
$C_D$: 抗力係数、$\rho$: 空気密度、$A$: 断面積
運動方程式:
$$ m\ddot{x} = -\frac{b}{m}v \cdot v_x, \quad m\ddot{y} = -mg – \frac{b}{m}v \cdot v_y $$
ここで $b = C_D\rho A/2$, $v = \sqrt{v_x^2 + v_y^2}$
Pythonでの実装
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
g = 9.8
m = 0.145 # 野球ボール [kg]
r = 0.0366 # 半径 [m]
A = np.pi * r**2
Cd = 0.3
rho = 1.225 # 空気密度 [kg/m^3]
b = 0.5 * Cd * rho * A
def projectile_drag(t, state):
x, y, vx, vy = state
v = np.sqrt(vx**2 + vy**2)
drag = b * v / m
return [vx, vy, -drag*vx, -g - drag*vy]
def hit_ground(t, state):
return state[1]
hit_ground.terminal = True
hit_ground.direction = -1
v0 = 40 # m/s
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# 異なる角度での比較
for alpha_deg in [15, 30, 45, 60]:
alpha = np.radians(alpha_deg)
# 空気抵抗あり
sol = solve_ivp(projectile_drag, [0, 20],
[0, 0, v0*np.cos(alpha), v0*np.sin(alpha)],
t_eval=np.linspace(0, 20, 1000), events=hit_ground, max_step=0.01)
axes[0].plot(sol.y[0], sol.y[1], '-', label=f'$\\alpha={alpha_deg}°$（抵抗あり）')
# 空気抵抗なし
t_f = 2*v0*np.sin(alpha)/g
t_nd = np.linspace(0, t_f, 200)
axes[0].plot(v0*np.cos(alpha)*t_nd, v0*np.sin(alpha)*t_nd-0.5*g*t_nd**2,
'--', alpha=0.4)
axes[0].set_title(f'弾道比較（実線:抵抗あり, 破線:なし, $v_0$={v0}m/s）')
axes[0].set_xlabel('x [m]')
axes[0].set_ylabel('y [m]')
axes[0].legend(fontsize=8)
axes[0].grid(True, alpha=0.3)
axes[0].set_ylim(0, None)
# 最適発射角の探索
angles = np.linspace(5, 80, 76)
ranges_drag = []
ranges_no_drag = []
for alpha_deg in angles:
alpha = np.radians(alpha_deg)
sol = solve_ivp(projectile_drag, [0, 20],
[0, 0, v0*np.cos(alpha), v0*np.sin(alpha)],
events=hit_ground, max_step=0.01)
ranges_drag.append(sol.y_events[0][0][0] if sol.y_events[0].size > 0 else 0)
ranges_no_drag.append(v0**2*np.sin(2*np.radians(alpha_deg))/g)
opt_drag = angles[np.argmax(ranges_drag)]
axes[1].plot(angles, ranges_no_drag, 'b--', label='抵抗なし')
axes[1].plot(angles, ranges_drag, 'r-', linewidth=2, label='抵抗あり')
axes[1].axvline(x=45, color='b', linestyle=':', alpha=0.5)
axes[1].axvline(x=opt_drag, color='r', linestyle=':', alpha=0.5)
axes[1].set_title(f'飛距離 vs 発射角\n最適角: 抵抗なし=45°, 抵抗あり≈{opt_drag:.0f}°')
axes[1].set_xlabel('発射角 [°]')
axes[1].set_ylabel('飛距離 [m]')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
# 速度の時間変化
alpha = np.radians(40)
sol = solve_ivp(projectile_drag, [0, 20],
[0, 0, v0*np.cos(alpha), v0*np.sin(alpha)],
t_eval=np.linspace(0, 10, 500), events=hit_ground, max_step=0.01)
speed = np.sqrt(sol.y[2]**2 + sol.y[3]**2)
axes[2].plot(sol.t, speed, 'b-', linewidth=2, label='速さ $|v|$')
axes[2].plot(sol.t, sol.y[2], 'r-', label='$v_x$')
axes[2].plot(sol.t, sol.y[3], 'g-', label='$v_y$')
axes[2].axhline(y=0, color='k', linewidth=0.5)
axes[2].set_title('速度の時間変化（$\\alpha=40°$）')
axes[2].set_xlabel('時間 [s]')
axes[2].set_ylabel('速度 [m/s]')
axes[2].legend()
axes[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('projectile_motion.png', dpi=150, bbox_inches='tight')
plt.show()
空気抵抗により最適発射角は45°から小さい方向にずれ、飛距離も大幅に減少することが確認できます。
まとめ
- 真空中: $R = v_0^2\sin(2\alpha)/g$、最適角 $= 45°$
- 空気抵抗（2次）があると最適角は $40°$ 前後に低下
- 空気抵抗付きでは解析解が得られず数値計算が必要
scipy.integrate.solve_ivpで弾道シミュレーション可能
次のステップとして、以下の記事も参考にしてください。