運動方程式 $\bm{F} = m\bm{a}$ はニュートン力学の中核ですが、この式を「知っている」だけでは物理の問題は解けません。大切なのは、具体的な状況に応じて運動方程式を立て、それを微分方程式として解くことです。

本記事では、質点の運動方程式を直交座標系と極座標系の両方で定式化し、放物運動・等速円運動・単振り子という3つの代表的な問題を通じて、方程式の立て方と解き方を丁寧に解説します。最後に Python を使った数値シミュレーションも行います。

本記事の内容

• 運動方程式を「解く」とは何か
• 直交座標系での運動方程式（$x$, $y$, $z$ 成分）
• 極座標系での運動方程式（$r$, $\theta$ 成分の導出）
• 具体例1: 放物運動（斜方投射の解析解）
• 具体例2: 等速円運動と向心力
• 具体例3: 単振り子の運動方程式（非線形と線形近似）
• Python による斜方投射と単振り子のシミュレーション

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

• ニュートンの運動の3法則をわかりやすく解説

運動方程式を「解く」とは何か

ニュートンの第2法則は、質量 $m$ の質点に合力 $\bm{F}$ が作用するとき、

$$ m\frac{d^2\bm{r}}{dt^2} = \bm{F} $$

と表されます。ここで $\bm{r}(t)$ は質点の位置ベクトルです。

この式は2階常微分方程式ですから、初期位置 $\bm{r}(0) = \bm{r}_0$ と初期速度 $\dot{\bm{r}}(0) = \bm{v}_0$ の2つの初期条件を与えれば、原理的に $\bm{r}(t)$ を一意に決定できます。

つまり「運動方程式を解く」とは、力 $\bm{F}$ と初期条件から質点の位置 $\bm{r}(t)$ を時刻 $t$ の関数として求めることです。解析的に閉じた形で解ける場合もあれば、数値的に解く必要がある場合もあります。

解法の一般的な手順は次のとおりです。

1. 座標系の選択: 問題の対称性に合った座標系を選ぶ
2. 力の洗い出し: 質点に作用するすべての力を特定する
3. 成分への分解: 運動方程式を各座標成分に分解する
4. 微分方程式を解く: 解析的または数値的に解く
5. 初期条件の適用: 積分定数を初期条件で決定する

座標系の選び方ひとつで、方程式が劇的に簡単になることがあります。直線運動には直交座標系、中心力の問題には極座標系が適しているのです。

直交座標系での運動方程式

定式化

3次元の直交座標系 $(x, y, z)$ では、位置ベクトルは

$$ \bm{r} = x\,\bm{e}_x + y\,\bm{e}_y + z\,\bm{e}_z $$

と表されます。ここで $\bm{e}_x$, $\bm{e}_y$, $\bm{e}_z$ は各軸方向の単位ベクトルです。直交座標系では単位ベクトルの方向が時間によらず一定なので、加速度は単純に各成分の2階時間微分で得られます。

$$ \bm{a} = \frac{d^2\bm{r}}{dt^2} = \ddot{x}\,\bm{e}_x + \ddot{y}\,\bm{e}_y + \ddot{z}\,\bm{e}_z $$

力 $\bm{F}$ も同様に成分分解すると、

$$ \bm{F} = F_x\,\bm{e}_x + F_y\,\bm{e}_y + F_z\,\bm{e}_z $$

運動方程式 $m\bm{a} = \bm{F}$ を成分ごとに書くと、

$$ m\ddot{x} = F_x $$

$$ m\ddot{y} = F_y $$

$$ m\ddot{z} = F_z $$

となります。直交座標系では、$x$, $y$, $z$ 方向の運動方程式が独立に（互いに結合せずに）得られるという利点があります。

極座標系での運動方程式

2次元極座標の定義

中心力の問題など、原点からの距離と角度で記述するのが自然な場合には、極座標 $(r, \theta)$ を使います。直交座標との関係は

$$ x = r\cos\theta, \quad y = r\sin\theta $$

です。

単位ベクトルの時間変化

極座標の肝は、単位ベクトル $\bm{e}_r$（動径方向）と $\bm{e}_\theta$（角度方向）が時間とともに回転する点です。

$$ \bm{e}_r = \cos\theta\,\bm{e}_x + \sin\theta\,\bm{e}_y $$

$$ \bm{e}_\theta = -\sin\theta\,\bm{e}_x + \cos\theta\,\bm{e}_y $$

これらを時間で微分すると、

$$ \frac{d\bm{e}_r}{dt} = (-\sin\theta\,\bm{e}_x + \cos\theta\,\bm{e}_y)\dot{\theta} = \dot{\theta}\,\bm{e}_\theta $$

$$ \frac{d\bm{e}_\theta}{dt} = (-\cos\theta\,\bm{e}_x – \sin\theta\,\bm{e}_y)\dot{\theta} = -\dot{\theta}\,\bm{e}_r $$

速度と加速度の導出

位置ベクトルは $\bm{r} = r\,\bm{e}_r$ なので、速度は積の微分法則で求められます。

$$ \begin{align} \bm{v} = \dot{\bm{r}} &= \dot{r}\,\bm{e}_r + r\,\frac{d\bm{e}_r}{dt} \\ &= \dot{r}\,\bm{e}_r + r\dot{\theta}\,\bm{e}_\theta \end{align} $$

第1項 $\dot{r}\,\bm{e}_r$ は動径方向の速度、第2項 $r\dot{\theta}\,\bm{e}_\theta$ は接線方向の速度です。

加速度はもう一度微分します。

$$ \begin{align} \bm{a} = \dot{\bm{v}} &= \frac{d}{dt}\left(\dot{r}\,\bm{e}_r + r\dot{\theta}\,\bm{e}_\theta\right) \\ &= \ddot{r}\,\bm{e}_r + \dot{r}\,\frac{d\bm{e}_r}{dt} + \dot{r}\dot{\theta}\,\bm{e}_\theta + r\ddot{\theta}\,\bm{e}_\theta + r\dot{\theta}\,\frac{d\bm{e}_\theta}{dt} \\ &= \ddot{r}\,\bm{e}_r + \dot{r}\dot{\theta}\,\bm{e}_\theta + \dot{r}\dot{\theta}\,\bm{e}_\theta + r\ddot{\theta}\,\bm{e}_\theta + r\dot{\theta}(-\dot{\theta}\,\bm{e}_r) \\ &= (\ddot{r} – r\dot{\theta}^2)\,\bm{e}_r + (r\ddot{\theta} + 2\dot{r}\dot{\theta})\,\bm{e}_\theta \end{align} $$

極座標での運動方程式

力を動径方向と角度方向に分解して $\bm{F} = F_r\,\bm{e}_r + F_\theta\,\bm{e}_\theta$ と書くと、$m\bm{a} = \bm{F}$ の各成分は

$$ m(\ddot{r} – r\dot{\theta}^2) = F_r $$

$$ m(r\ddot{\theta} + 2\dot{r}\dot{\theta}) = F_\theta $$

となります。

$r$ 方向の方程式に現れる $-mr\dot{\theta}^2$ は遠心力に対応する項です。回転する座標系から見ると、質点は原点から遠ざかろうとする力を受けているように見えます。

$\theta$ 方向の方程式に現れる $2m\dot{r}\dot{\theta}$ はコリオリの力に対応する項です。これは $r$ が変化しながら回転する質点に現れる見かけの力です。

具体例1: 放物運動（斜方投射）

問題設定

質量 $m$ の質点を、地面（$y = 0$）から初速 $v_0$、仰角 $\alpha$ で投射します。空気抵抗は無視し、重力加速度 $g$ のみが鉛直下向きに作用するものとします。

運動方程式の立式

$x$ 軸を水平方向（投射方向）、$y$ 軸を鉛直上向きにとると、作用する力は重力のみなので

$$ F_x = 0, \quad F_y = -mg $$

運動方程式は各成分で

$$ m\ddot{x} = 0 $$

$$ m\ddot{y} = -mg $$

解の導出

$x$ 方向の方程式を解きます。$m$ で割って

$$ \ddot{x} = 0 $$

1回積分すると

$$ \dot{x}(t) = C_1 $$

初期条件 $\dot{x}(0) = v_0\cos\alpha$ より $C_1 = v_0\cos\alpha$ なので

$$ \dot{x}(t) = v_0\cos\alpha $$

もう1回積分して

$$ x(t) = v_0 t\cos\alpha + C_2 $$

初期条件 $x(0) = 0$ より $C_2 = 0$ なので

$$ x(t) = v_0 t\cos\alpha $$

次に $y$ 方向の方程式を解きます。

$$ \ddot{y} = -g $$

1回積分すると

$$ \dot{y}(t) = -gt + C_3 $$

初期条件 $\dot{y}(0) = v_0\sin\alpha$ より $C_3 = v_0\sin\alpha$ なので

$$ \dot{y}(t) = v_0\sin\alpha – gt $$

もう1回積分して

$$ y(t) = v_0 t\sin\alpha – \frac{1}{2}gt^2 + C_4 $$

初期条件 $y(0) = 0$ より $C_4 = 0$ なので

$$ y(t) = v_0 t\sin\alpha – \frac{1}{2}gt^2 $$

軌跡の方程式

$x(t)$ の式から $t = \frac{x}{v_0\cos\alpha}$ を $y(t)$ に代入すると、時刻を消去できます。

$$ \begin{align} y &= v_0 \cdot \frac{x}{v_0\cos\alpha} \cdot \sin\alpha – \frac{1}{2}g\left(\frac{x}{v_0\cos\alpha}\right)^2 \\ &= x\tan\alpha – \frac{g}{2v_0^2\cos^2\alpha}x^2 \end{align} $$

これは $x$ の2次関数であり、放物線を表します。放物運動という名称はここに由来します。

到達距離と最高点

着地時刻は $y(t) = 0$ かつ $t > 0$ から求まります。

$$ v_0 t\sin\alpha – \frac{1}{2}gt^2 = 0 $$

$$ t\left(v_0\sin\alpha – \frac{1}{2}gt\right) = 0 $$

$t \neq 0$ の解は

$$ t_{\text{land}} = \frac{2v_0\sin\alpha}{g} $$

このとき、水平到達距離（飛距離）は

$$ R = x(t_{\text{land}}) = v_0 \cos\alpha \cdot \frac{2v_0\sin\alpha}{g} = \frac{v_0^2 \sin 2\alpha}{g} $$

$\sin 2\alpha$ は $\alpha = 45°$ のとき最大値 1 をとるので、飛距離は $\alpha = 45°$ で最大になります。

最高点の時刻は $\dot{y}(t) = 0$ から

$$ t_{\text{top}} = \frac{v_0\sin\alpha}{g} $$

最高点の高さは

$$ H = y(t_{\text{top}}) = \frac{v_0^2\sin^2\alpha}{2g} $$

具体例2: 等速円運動と向心力

問題設定

質点が半径 $R$ の円上を一定の角速度 $\omega$ で運動しているとします。

極座標での解析

等速円運動では $r = R$（一定）かつ $\dot{\theta} = \omega$（一定）なので、

$$ \dot{r} = 0, \quad \ddot{r} = 0, \quad \ddot{\theta} = 0 $$

極座標の運動方程式に代入すると、

$$ m(0 – R\omega^2) = F_r \quad \Longrightarrow \quad F_r = -mR\omega^2 $$

$$ m(R \cdot 0 + 2 \cdot 0 \cdot \omega) = F_\theta \quad \Longrightarrow \quad F_\theta = 0 $$

$F_r < 0$ は力が原点向き（動径方向内向き）であることを意味します。これが向心力です。

$$ |\bm{F}| = mR\omega^2 = m\frac{v^2}{R} $$

ここで $v = R\omega$ を用いました。

向心力がなくなると、質点は慣性の法則により接線方向に直進します。糸の付いたボールを円運動させているとき、糸を離すとボールが直線的に飛んでいくのはこのためです。

直交座標での確認

直交座標で等速円運動を表すと、

$$ x(t) = R\cos(\omega t), \quad y(t) = R\sin(\omega t) $$

加速度を計算すると

$$ \ddot{x} = -R\omega^2\cos(\omega t), \quad \ddot{y} = -R\omega^2\sin(\omega t) $$

加速度ベクトル $\bm{a} = (\ddot{x}, \ddot{y})$ は常に原点を向いており、その大きさは

$$ |\bm{a}| = R\omega^2 $$

と一定です。極座標で得た結果と一致しています。

具体例3: 単振り子の運動方程式

問題設定

長さ $l$ の質量の無視できる糸の先に質量 $m$ のおもりを吊るし、鉛直面内で振動させます。糸がたるまないとすると、おもりは半径 $l$ の円弧上を運動します。

運動方程式の導出

鉛直下向きからの角度を $\theta$ とすると、おもりの接線方向の加速度は $l\ddot{\theta}$ です。接線方向に作用する力は重力の接線成分 $-mg\sin\theta$（復元力なので負）だけなので、

$$ ml\ddot{\theta} = -mg\sin\theta $$

整理すると

$$ \ddot{\theta} + \frac{g}{l}\sin\theta = 0 $$

これが単振り子の運動方程式です。$\sin\theta$ が含まれているため、この方程式は非線形微分方程式です。一般には初等関数では解けず、楕円積分が必要になります。

微小振動の近似（線形化）

振幅が十分に小さいとき（$|\theta| \ll 1$）、テイラー展開

$$ \sin\theta = \theta – \frac{\theta^3}{6} + \cdots \approx \theta $$

を適用すると、運動方程式は

$$ \ddot{\theta} + \frac{g}{l}\theta = 0 $$

と線形化されます。$\omega_0 = \sqrt{g/l}$ と置くと

$$ \ddot{\theta} + \omega_0^2\,\theta = 0 $$

これは単振動の方程式です。一般解は

$$ \theta(t) = A\cos(\omega_0 t + \phi) $$

ここで $A$ は振幅、$\phi$ は初期位相であり、いずれも初期条件によって決まります。

周期は

$$ T = \frac{2\pi}{\omega_0} = 2\pi\sqrt{\frac{l}{g}} $$

となり、質量 $m$ に依存しないという有名な結果が得られます。振り子の等時性と呼ばれるこの性質は、ガリレオ・ガリレイが発見したとされています。

ただし、この結果は微小振動の近似のもとでのみ成り立ちます。振幅が大きくなると周期は長くなり、厳密な解は楕円積分を用いて表されます。

Python による数値シミュレーション

斜方投射の軌跡

まず、斜方投射の解析解を可視化します。複数の投射角度で軌跡を描き、$\alpha = 45°$ で飛距離が最大になることを確認しましょう。

import numpy as np
import matplotlib.pyplot as plt

# パラメータ
v0 = 20.0   # 初速 [m/s]
g = 9.8     # 重力加速度 [m/s^2]

# 複数の投射角度
angles_deg = [15, 30, 45, 60, 75]

plt.figure(figsize=(10, 6))

for alpha_deg in angles_deg:
    alpha = np.radians(alpha_deg)

    # 着地時刻
    t_land = 2 * v0 * np.sin(alpha) / g
    t = np.linspace(0, t_land, 300)

    # 解析解
    x = v0 * t * np.cos(alpha)
    y = v0 * t * np.sin(alpha) - 0.5 * g * t**2

    # 飛距離
    R = v0**2 * np.sin(2 * alpha) / g

    plt.plot(x, y, linewidth=2, label=f'$\\alpha = {alpha_deg}°$ (R = {R:.1f} m)')

plt.xlabel('x [m]', fontsize=12)
plt.ylabel('y [m]', fontsize=12)
plt.title('Projectile Motion for Various Launch Angles', fontsize=14)
plt.legend(fontsize=10)
plt.grid(True, alpha=0.3)
plt.xlim(0, None)
plt.ylim(0, None)
plt.gca().set_aspect('equal')
plt.tight_layout()
plt.savefig('projectile_motion.png', dpi=150, bbox_inches='tight')
plt.show()

単振り子の数値解（非線形 vs 線形近似）

次に、単振り子の運動方程式を scipy.integrate.solve_ivp で数値的に解きます。非線形の厳密な方程式と線形近似の両方を計算し、振幅が大きい場合に近似がどの程度ずれるかを観察します。

2階微分方程式を1階の連立方程式に変換します。$\omega = \dot{\theta}$ と置くと、

$$ \frac{d\theta}{dt} = \omega, \quad \frac{d\omega}{dt} = -\frac{g}{l}\sin\theta $$

線形近似の場合は $\sin\theta$ を $\theta$ に置き換えます。

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# パラメータ
g = 9.8    # 重力加速度 [m/s^2]
l = 1.0    # 振り子の長さ [m]
omega0 = np.sqrt(g / l)  # 固有角振動数

# 非線形の運動方程式（厳密）
def pendulum_nonlinear(t, state):
    theta, omega = state
    dtheta_dt = omega
    domega_dt = -(g / l) * np.sin(theta)
    return [dtheta_dt, domega_dt]

# 線形近似の運動方程式
def pendulum_linear(t, state):
    theta, omega = state
    dtheta_dt = omega
    domega_dt = -(g / l) * theta
    return [dtheta_dt, domega_dt]

# シミュレーション時間
t_span = (0, 10)
t_eval = np.linspace(0, 10, 1000)

# 2つの初期角度で比較
initial_angles_deg = [10, 60]

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

for idx, theta0_deg in enumerate(initial_angles_deg):
    theta0 = np.radians(theta0_deg)
    y0 = [theta0, 0.0]  # 初期角度 theta0, 初期角速度 0

    # 非線形（厳密）
    sol_nonlinear = solve_ivp(
        pendulum_nonlinear, t_span, y0,
        method='RK45', t_eval=t_eval, rtol=1e-10
    )

    # 線形近似
    sol_linear = solve_ivp(
        pendulum_linear, t_span, y0,
        method='RK45', t_eval=t_eval, rtol=1e-10
    )

    # 解析解（線形近似）
    theta_analytical = theta0 * np.cos(omega0 * t_eval)

    ax = axes[idx]
    ax.plot(sol_nonlinear.t, np.degrees(sol_nonlinear.y[0]),
            'b-', linewidth=2, label='Nonlinear (exact)')
    ax.plot(sol_linear.t, np.degrees(sol_linear.y[0]),
            'r--', linewidth=1.5, label='Linear approx.')
    ax.plot(t_eval, np.degrees(theta_analytical),
            'g:', linewidth=1.5, label='Analytical (linear)')
    ax.set_xlabel('Time [s]', fontsize=12)
    ax.set_ylabel('$\\theta$ [deg]', fontsize=12)
    ax.set_title(f'Simple Pendulum ($\\theta_0 = {theta0_deg}°$)', fontsize=13)
    ax.legend(fontsize=10)
    ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('simple_pendulum.png', dpi=150, bbox_inches='tight')
plt.show()

# 周期の比較を表示
print("=== 周期の比較 ===")
T_linear = 2 * np.pi / omega0
print(f"線形近似の周期: T = 2π√(l/g) = {T_linear:.4f} s")

for theta0_deg in [10, 30, 60, 90]:
    theta0 = np.radians(theta0_deg)
    sol = solve_ivp(
        pendulum_nonlinear, (0, 20), [theta0, 0.0],
        method='RK45', t_eval=np.linspace(0, 20, 10000), rtol=1e-12
    )
    # ゼロクロス（正→負）を検出して周期を推定
    theta_vals = sol.y[0]
    crossings = []
    for i in range(1, len(theta_vals)):
        if theta_vals[i-1] > 0 and theta_vals[i] <= 0:
            # 線形補間でゼロクロスの時刻を推定
            t_cross = sol.t[i-1] + (sol.t[i] - sol.t[i-1]) * theta_vals[i-1] / (theta_vals[i-1] - theta_vals[i])
            crossings.append(t_cross)
    if len(crossings) >= 2:
        T_numerical = 2 * (crossings[1] - crossings[0])
        error = (T_numerical - T_linear) / T_linear * 100
        print(f"θ₀ = {theta0_deg:3d}° → T = {T_numerical:.4f} s (線形近似からの誤差: {error:+.2f}%)")

上記のコードを実行すると、次のことが確認できます。

• 初期角度 10 度の場合: 非線形と線形近似の結果はほぼ完全に一致します。微小角度の近似が十分に有効であることがわかります
• 初期角度 60 度の場合: 時間が経つにつれて非線形と線形近似の位相がずれていきます。非線形の場合の周期が長くなるため、振動の山と谷の位置が徐々にずれるのです
• 周期の比較: 初期角度が大きくなるにつれて、実際の周期が線形近似の周期 $T = 2\pi\sqrt{l/g}$ からずれていく様子が数値で確認できます

まとめ

本記事では、質点の運動方程式の立て方と解法について解説しました。

• 運動方程式を解くとは: 力と初期条件から、位置 $\bm{r}(t)$ を時刻の関数として求めることです。座標系の選択が問題を簡単にする鍵となります
• 直交座標系: 運動方程式が $x$, $y$, $z$ 成分に独立に分解でき、直線的な運動の解析に適しています
• 極座標系: 加速度に遠心力項 $-r\dot{\theta}^2$ とコリオリ項 $2\dot{r}\dot{\theta}$ が現れます。中心力の問題や円運動の解析に自然な座標系です
• 放物運動: 水平と鉛直が独立に解け、軌跡は放物線になります。飛距離は $\alpha = 45°$ で最大です
• 等速円運動: 向心力 $F = mR\omega^2 = mv^2/R$ が常に中心向きに必要です
• 単振り子: 非線形方程式 $\ddot{\theta} + (g/l)\sin\theta = 0$ を微小角度で線形化すると単振動になり、周期 $T = 2\pi\sqrt{l/g}$ が得られます

次のステップとして、以下のトピックもぜひ学んでみてください。

• 単振動とバネ-質点系の運動方程式は、運動方程式の最も基本的な応用です
• エネルギー保存則を導入すると、運動方程式を解かずに運動の定性的な性質を知ることができます
• 解析力学（ラグランジュ力学）では、座標系に依存しない一般的な定式化が可能になります