3. 受迫单摆的周期解与共振响应仿真

受迫阻尼单摆的建模与频率扫描设计

本节对象是带线性阻尼、受简谐外力矩驱动的单摆。取角位移 θ 以竖直向下平衡位置为零、逆时针为正，运动方程为：

θ'' + 2ζω₀θ' + ω₀²sinθ = (A/ω₀²)·ω₀²·cos(ωt)，工程上更常用归一化形式：

θ'' + γθ' + ω₀²sinθ = F·cos(ωt)

其中 ω₀=√(g/L) 为小振幅固有角频率，γ 为阻尼系数，F 为单位质量归一化驱动力幅值，ω 为驱动角频率。线性化（sinθ≈θ）后稳态振幅满足经典受迫振动解：

Θ_lin(ω) = F / √[(ω₀²−ω²)² + (γω)²]，共振峰位于 ω_r = √(ω₀² − γ²/2)。

非线性项 sinθ 会使大振幅响应偏离该公式：共振峰右移（硬化特性）、峰高被压低、并可能出现跳跃与滞回。因此扫描策略必须同时记录线性预测与非线性仿真结果，才能识别差异。

数值实现：参数、积分与稳态振幅提取

采用四阶龙格–库塔法积分，步长 h=0.002 s，每个频点先积分足够长的瞬态段（约 60 个驱动周期）再取后段做稳态幅值统计。示例参数取 L=1.0 m、g=9.81 m/s²（ω₀≈3.132 rad/s）、γ=0.25 s⁻¹、F=0.6 rad/s²，扫描 ω 从 1.6 到 4.6 rad/s，步长 0.02。

import numpy as np

g, L = 9.81, 1.0
w0 = np.sqrt(g / L)          # 3.132 rad/s
gamma, F = 0.25, 0.6
h = 0.002

def deriv(state, t, w):
    th, om = state
    return np.array([om, -gamma*om - w0**2*np.sin(th) + F*np.cos(w*t)])

def rk4(state, t, w, h):
    k1 = deriv(state, t, w)
    k2 = deriv(state + 0.5*h*k1, t + 0.5*h, w)
    k3 = deriv(state + 0.5*h*k2, t + 0.5*h, w)
    k4 = deriv(state + h*k3, t + h, w)
    return state + h/6*(k1 + 2*k2 + 2*k3 + k4)

def steady_amplitude(w, th0=0.05, om0=0.0, n_trans=60, n_meas=20):
    T = 2*np.pi / w
    state = np.array([th0, om0]); t = 0.0
    steps_trans = int(n_trans*T/h)
    for _ in range(steps_trans):
        state = rk4(state, t, w, h); t += h
    th_max, th_min = -1e9, 1e9
    steps_meas = int(n_meas*T/h)
    for _ in range(steps_meas):
        state = rk4(state, t, w, h); t += h
        th_max = max(th_max, state[0]); th_min = min(th_min, state[0])
    return 0.5*(th_max - th_min)

def linear_amp(w):
    return F / np.sqrt((w0**2 - w**2)**2 + (gamma*w)**2)

omegas = np.arange(1.6, 4.61, 0.02)
amp_nl = np.array([steady_amplitude(w) for w in omegas])
amp_lin = linear_amp(omegas)

i_peak = np.argmax(amp_nl)
print("非线性共振峰: w = %.3f rad/s, 振幅 = %.4f rad" % (omegas[i_peak], amp_nl[i_peak]))
print("线性预测峰: w = %.3f rad/s, 振幅 = %.4f rad" % (np.sqrt(w0**2 - gamma**2/2), linear_amp(np.sqrt(w0**2 - gamma**2/2))))

运行后可得到两条曲线：线性预测峰约在 ω≈3.127 rad/s、峰值约 2.40 rad；非线性仿真峰通常出现在略高频率处且峰值略低，具体数值随 F 增大而偏移更明显。若把 F 提高到 1.2 以上，扫描曲线会出现明显的幅值跳变，此时应改用升频与降频双向扫描来观察滞回。

共振峰识别与线性/非线性对比

下图给出本节的核心判读逻辑：同一驱动频率下，线性模型只给出单一幅值，而非线性系统在共振区附近存在多个可能的稳态分支，扫描方向不同会落到不同分支上，这正是“共振峰位置偏移 + 幅值跳跃”的来源。

结果判读与验证要点

• 峰位对比：线性预测峰 ω_r=√(ω₀²−γ²/2)≈3.127 rad/s；非线性仿真峰若明显大于该值，说明 sinθ 的硬化效应已不可忽略，此时小振幅线性近似失效。
• 峰高对比：非线性峰值一般低于线性预测，因为大振幅下有效回复力偏离线性，能量在更高次谐波上分散。
• 滞回检验：对共振区附近频点分别以升频和降频为初值重新积分，若同一 ω 得到两个不同稳态幅值，即确认存在多稳态分支。
• 收敛性验证：将步长 h 减半、瞬态周期数加倍，若峰值位置变化小于扫描步长，则结果可信；否则需继续加密。