2. 阻尼单摆的衰减振荡与能量耗散分析

模型与参数设定

在无阻尼单摆方程基础上加入与角速度成正比的线性阻尼项，得到本节使用的二阶非线性常微分方程：

θ'' + c·θ' + (g/L)·sin θ = 0

其中 θ 为摆角（rad），c 为阻尼系数（1/s），g = 9.81 m/s²，L 为摆长。取 L = 1.0 m，则固有角频率 ω₀ = √(g/L) ≈ 3.132 rad/s。初始条件统一取 θ(0) = 0.6 rad、θ'(0) = 0，即从静止释放。

本节对比三组阻尼系数：c = 0.05（弱阻尼）、c = 0.5（中等阻尼）、c = 2.0（过阻尼）。判别依据是阻尼比 ζ = c/(2ω₀)：ζ < 1 为欠阻尼，ζ = 1 为临界阻尼，ζ > 1 为过阻尼。三组参数分别对应 ζ ≈ 0.008、0.080、0.319，均处于欠阻尼区间，但衰减速度差异显著，便于观察能量耗散速率的变化。

数值积分实现

将二阶方程化为一阶方程组，用四阶龙格–库塔法（RK4）积分。状态向量 y = [θ, ω]，其中 ω = θ'：

import numpy as np

g, L = 9.81, 1.0

def deriv(y, c):
    theta, omega = y
    return np.array([omega, -c*omega - (g/L)*np.sin(theta)])

def rk4(y0, c, dt, n_steps):
    y = np.array(y0, dtype=float)
    traj = np.empty((n_steps+1, 2))
    traj[0] = y
    for k in range(n_steps):
        k1 = deriv(y, c)
        k2 = deriv(y + 0.5*dt*k1, c)
        k3 = deriv(y + 0.5*dt*k2, c)
        k4 = deriv(y + dt*k3, c)
        y = y + dt/6.0*(k1 + 2*k2 + 2*k3 + k4)
        traj[k+1] = y
    return traj

dt, T = 0.005, 40.0
n_steps = int(T/dt)
for c in (0.05, 0.5, 2.0):
    traj = rk4([0.6, 0.0], c, dt, n_steps)
    # 后续用于能量计算与相图绘制

时间步长 dt = 0.005 s 远小于最快振荡周期（约 2 s），保证积分精度。总时长 40 s 足以让弱阻尼情形也出现明显衰减。

能量计算与耗散率

单摆的机械能由动能与势能组成（以最低点为势能零点）：

E = 0.5*m*L**2*omega**2 + m*g*L*(1 - np.cos(theta))

取 m = 1 kg。无阻尼时 E 守恒；加入阻尼后，能量耗散率满足 dE/dt = -c·m·L²·ω²，恒为负值，说明系统持续失能。数值上可用相邻时刻能量差除以 dt 估计瞬时耗散率，并与解析式对照验证积分正确性。

衰减振荡与相轨迹收敛

下图给出三组阻尼系数下的关键行为对比：左侧为摆角时间序列的包络衰减，右侧为相平面轨迹向原点（稳定平衡点 θ = 0, ω = 0）收敛的过程。阻尼越大，包络衰减越快，相轨迹螺旋收缩越迅速。

结果判读与验证

从仿真输出中提取三项指标进行判读：

表中能量占比为按上述参数积分后的预期量级，实际数值随步长与积分时长略有差异，应以程序输出为准。验证耗散系统收敛性的直接判据是：相轨迹在足够长时间后是否单调趋近原点，且能量序列是否单调不增（允许数值噪声引起的微小波动）。若能量出现持续上升，应检查阻尼项符号是否写反或步长是否过大。