如何数值求解含参量驱动的受迫振子微分方程以捕捉反共振响应？

阿华AIGC实验室

2026-4-14

首先，咱们先理清你当前方法的问题所在，再一步步解决如何捕捉反共振的数值响应：

为什么你的PSD方法看不到反共振？

你现在的代码是固定了一个驱动频率Og=2，然后计算该系统的功率谱密度（PSD）。但反共振并不是单个驱动频率下的频谱特征——它是当你扫动驱动频率f时，系统在驱动频率处的稳态振幅出现的一个凹陷。打个比方，就像共振是某个频率下振幅突然变大，反共振则是某个频率下振幅意外变小，这只能通过「扫频实验」才能观察到：即对一系列不同的驱动频率f，分别计算系统的稳态响应振幅，再把振幅随f的变化画出来，才能看到那个凹陷。

正确的数值方法：频率扫频+稳态振幅提取

要捕捉反共振，你需要按以下步骤操作：

1. 确定驱动频率的扫描范围

首先，根据你的解析推导估算反共振频率的大致位置（比如靠近f0的某个比例值），然后设定一个覆盖该位置的f区间（比如从f0/2到2f0），取足够多的采样点（50-100个点就能得到平滑的曲线）。

2. 对每个驱动频率f，求解ODE至稳态

对每个f：

初始化位移和速度的初始值（比如都设为0，因为我们只关心最终的稳态响应，初始瞬态会被阻尼衰减掉）
积分ODE足够长的时间，确保初始瞬态完全消失（可以通过观察x(t)的曲线判断，或者根据阻尼系数Γ估算衰减时间）
截取积分结果的最后一段（稳态段）用于后续分析

3. 提取稳态响应在驱动频率f处的振幅

稳态下，x(t)的主导响应成分是频率为f的简谐振动（因为驱动项是sin(ft)）。你可以用两种可靠方法提取振幅：

曲线拟合法：用A*sin(ft + φ)拟合稳态段的x(t)数据，直接得到振幅A（这种方法精度高，能过滤掉参量驱动带来的小频率分量干扰）
傅里叶变换法：对稳态段做FFT，找到对应频率f的分量振幅（注意要保证FFT的频率分辨率足够高，避免频率偏移误差）

4. 绘制振幅-频率曲线

把每个f对应的振幅A画成曲线，就能清晰看到反共振对应的凹陷了。

改进后的代码示例

下面是基于你的需求修改的完整代码，实现了频率扫频和稳态振幅提取：

import numpy as np
from matplotlib import pyplot as plt
from scipy.optimize import curve_fit
from scipy.integrate import solve_ivp  # 用更高效的自适应步长求解器

# 系统参数
Gamma = 0.6       # 阻尼系数
f0 = 3            # 固有频率
Lambda = 0.5      # 参量驱动强度
F = 1             # 驱动强度

# 驱动频率扫描范围（覆盖可能的反共振点）
f_values = np.linspace(1, 4, 60)
steady_amplitudes = []

# 拟合用的正弦函数模板
def fit_sine(t, A, phi):
    return A * np.sin(2 * np.pi * f_current * t + phi)

# 对每个驱动频率进行计算
for f_current in f_values:
    # 积分总时间（足够长以到达稳态）
    tf = 100
    t_eval = np.linspace(0, tf, 100000)
    
    # 定义微分方程
    def ode_system(t, y):
        x, x_dot = y
        # 计算二阶导数x''
        x_double_dot = (-2 * Gamma * x_dot 
                       - (2 * np.pi * f0)**2 * (1 + Lambda * np.sin(2 * np.pi * 2 * f_current * t)) * x 
                       + F * np.sin(2 * np.pi * f_current * t))
        return [x_dot, x_double_dot]
    
    # 求解ODE
    sol = solve_ivp(ode_system, [0, tf], [0, 0], t_eval=t_eval)
    x_data = sol.y[0]
    t_data = sol.t
    
    # 截取稳态段（取最后20秒的数据）
    steady_idx = np.where(t_data >= tf - 20)[0][0]
    x_steady = x_data[steady_idx:]
    t_steady = t_data[steady_idx:]
    
    # 拟合提取振幅
    popt, _ = curve_fit(fit_sine, t_steady, x_steady)
    amplitude = np.abs(popt[0])
    steady_amplitudes.append(amplitude)

# 绘制结果
plt.figure(figsize=(10, 6))
plt.plot(f_values, steady_amplitudes, 'o-', color='#2ecc71', linewidth=2, markersize=4, label='数值稳态振幅')
plt.xlabel('驱动频率 f', fontsize=12)
plt.ylabel('稳态响应振幅', fontsize=12)
plt.title('含参量驱动的受迫振子响应曲线（反共振凹陷）', fontsize=14)
plt.legend(fontsize=10)
plt.grid(alpha=0.3)
plt.show()