現実の振動系には必ず減衰（エネルギー散逸）が存在し、外部から強制力が加わることもあります。減衰振動と強制振動は、機械振動の設計、地震応答解析、電気回路のフィルタ設計など、工学のあらゆる分野で重要です。
特に共振現象は、橋の崩壊（タコマナローズ橋事件）から楽器の音色まで、劇的な効果をもたらします。本記事では、減衰振動の3つのケースと強制振動の定常解を、微分方程式の解法とともに丁寧に導出します。
本記事の内容
- 減衰振動の運動方程式と3つのケース
- 強制振動と定常解の導出
- 共振現象の数学的記述
- Pythonでの可視化
前提知識
減衰振動
運動方程式
速度に比例する減衰力 $F_d = -c\dot{x}$ を加えると：
$$ m\ddot{x} + c\dot{x} + kx = 0 $$
$\gamma = c/(2m)$（減衰係数）、$\omega_0 = \sqrt{k/m}$（固有角振動数）とおくと：
$$ \begin{equation} \ddot{x} + 2\gamma\dot{x} + \omega_0^2 x = 0 \end{equation} $$
特性方程式
$x = e^{\lambda t}$ を仮定すると：
$$ \lambda^2 + 2\gamma\lambda + \omega_0^2 = 0 \implies \lambda = -\gamma \pm \sqrt{\gamma^2 – \omega_0^2} $$
$\gamma$ と $\omega_0$ の大小関係で3つのケースに分かれます。
ケース1: 不足減衰（$\gamma < \omega_0$）
減衰振動数 $\omega_d = \sqrt{\omega_0^2 – \gamma^2}$ とおくと：
$$ \begin{equation} x(t) = Ae^{-\gamma t}\cos(\omega_d t + \phi) \end{equation} $$
振動しながら指数関数的に減衰します。最もよく見られるケースです。
ケース2: 臨界減衰（$\gamma = \omega_0$）
$$ \begin{equation} x(t) = (C_1 + C_2 t)e^{-\gamma t} \end{equation} $$
振動せずに最も速く平衡位置に戻ります。ドアクローザーやショックアブソーバーの理想的な設計です。
ケース3: 過減衰（$\gamma > \omega_0$）
$$ \begin{equation} x(t) = C_1 e^{\lambda_1 t} + C_2 e^{\lambda_2 t} \end{equation} $$
$\lambda_1, \lambda_2$ はともに負の実数。振動せずにゆっくりと平衡に戻ります。
減衰比（ダンピング比）
$$ \zeta = \frac{\gamma}{\omega_0} = \frac{c}{2\sqrt{mk}} $$
$\zeta < 1$: 不足減衰、$\zeta = 1$: 臨界減衰、$\zeta > 1$: 過減衰。
強制振動
運動方程式
外力 $F(t) = F_0\cos\Omega t$ が加わる場合：
$$ \begin{equation} \ddot{x} + 2\gamma\dot{x} + \omega_0^2 x = \frac{F_0}{m}\cos\Omega t \end{equation} $$
定常解の導出
十分時間が経った後の定常解（過渡解は減衰して消える）を求めます。
$x_p(t) = X\cos(\Omega t – \delta)$ と仮定し、代入すると：
$$ \begin{equation} X = \frac{F_0/m}{\sqrt{(\omega_0^2 – \Omega^2)^2 + (2\gamma\Omega)^2}} \end{equation} $$
位相遅れ：
$$ \begin{equation} \tan\delta = \frac{2\gamma\Omega}{\omega_0^2 – \Omega^2} \end{equation} $$
共振
振幅 $X$ が最大になる周波数は、$\frac{dX}{d\Omega} = 0$ から：
$$ \Omega_{\text{res}} = \sqrt{\omega_0^2 – 2\gamma^2} \quad (\gamma < \omega_0/\sqrt{2} \text{ のとき}) $$
減衰が小さいとき $\Omega_{\text{res}} \approx \omega_0$ で、このとき振幅は $X_{\max} \approx F_0/(2m\gamma\omega_0)$ と非常に大きくなります。
Pythonでの実装
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
omega_0 = 2 * np.pi # 固有角振動数
t = np.linspace(0, 10, 2000)
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# 左上: 3種の減衰振動
zeta_values = [0.1, 0.5, 1.0, 2.0]
labels = ['不足減衰 (ζ=0.1)', '不足減衰 (ζ=0.5)', '臨界減衰 (ζ=1.0)', '過減衰 (ζ=2.0)']
for zeta, label in zip(zeta_values, labels):
gamma = zeta * omega_0
def eom(state, t):
x, v = state
return [v, -2*gamma*v - omega_0**2 * x]
sol = odeint(eom, [1.0, 0.0], t)
axes[0,0].plot(t, sol[:, 0], linewidth=2, label=label)
axes[0,0].axhline(y=0, color='k', linewidth=0.5)
axes[0,0].set_xlabel('時間 [s]', fontsize=12)
axes[0,0].set_ylabel('変位 x', fontsize=12)
axes[0,0].set_title('減衰振動の3ケース', fontsize=14)
axes[0,0].legend(fontsize=9)
axes[0,0].grid(True, alpha=0.3)
# 右上: 共振曲線（振幅-周波数応答）
Omega = np.linspace(0.01, 3 * omega_0, 500)
F0_m = 1.0
for zeta in [0.05, 0.1, 0.3, 0.5, 1.0]:
gamma = zeta * omega_0
X = F0_m / np.sqrt((omega_0**2 - Omega**2)**2 + (2*gamma*Omega)**2)
axes[0,1].plot(Omega / omega_0, X, linewidth=2, label=f'ζ={zeta}')
axes[0,1].set_xlabel(r'$\Omega / \omega_0$', fontsize=12)
axes[0,1].set_ylabel('振幅 X', fontsize=12)
axes[0,1].set_title('共振曲線（周波数応答）', fontsize=14)
axes[0,1].legend(fontsize=10)
axes[0,1].grid(True, alpha=0.3)
axes[0,1].set_xlim(0, 3)
# 左下: 位相-周波数応答
for zeta in [0.05, 0.1, 0.3, 0.5, 1.0]:
gamma = zeta * omega_0
delta = np.arctan2(2*gamma*Omega, omega_0**2 - Omega**2)
axes[1,0].plot(Omega / omega_0, np.degrees(delta), linewidth=2, label=f'ζ={zeta}')
axes[1,0].set_xlabel(r'$\Omega / \omega_0$', fontsize=12)
axes[1,0].set_ylabel('位相遅れ δ [°]', fontsize=12)
axes[1,0].set_title('位相-周波数応答', fontsize=14)
axes[1,0].legend(fontsize=10)
axes[1,0].grid(True, alpha=0.3)
axes[1,0].set_xlim(0, 3)
# 右下: 強制振動の時間応答（過渡+定常）
zeta = 0.1
gamma = zeta * omega_0
Omega_force = omega_0 * 0.9
F0 = 1.0
def forced_eom(state, t):
x, v = state
return [v, -2*gamma*v - omega_0**2*x + F0*np.cos(Omega_force*t)]
sol = odeint(forced_eom, [0.0, 0.0], t)
axes[1,1].plot(t, sol[:, 0], 'b-', linewidth=1.5, label='過渡+定常応答')
# 定常解のみ
X_ss = F0 / np.sqrt((omega_0**2 - Omega_force**2)**2 + (2*gamma*Omega_force)**2)
delta_ss = np.arctan2(2*gamma*Omega_force, omega_0**2 - Omega_force**2)
x_steady = X_ss * np.cos(Omega_force * t - delta_ss)
axes[1,1].plot(t, x_steady, 'r--', linewidth=1.5, label='定常解のみ', alpha=0.7)
axes[1,1].set_xlabel('時間 [s]', fontsize=12)
axes[1,1].set_ylabel('変位 x', fontsize=12)
axes[1,1].set_title(f'強制振動 (ζ={zeta}, Ω/ω₀=0.9)', fontsize=14)
axes[1,1].legend(fontsize=10)
axes[1,1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
まとめ
- 減衰振動: $\ddot{x} + 2\gamma\dot{x} + \omega_0^2 x = 0$
- 3ケース: 不足減衰（$\zeta<1$, 振動減衰）、臨界減衰（$\zeta=1$, 最速収束）、過減衰（$\zeta>1$, 緩やかな収束）
- 強制振動の定常解: $X = F_0/m / \sqrt{(\omega_0^2-\Omega^2)^2 + (2\gamma\Omega)^2}$
- 共振: $\Omega \approx \omega_0$ で振幅が最大化。減衰が小さいほど鋭い共振
次のステップとして、以下の記事も参考にしてください。
- 万有引力とケプラーの法則
- 2次遅れ系の応答 — 制御工学での類似概念