ODEINT求解塔正弦激励运动方程：异常振荡问题排查求助

1. 振荡的物理原因：频率成分叠加 #

你看到的“异常振荡”本质是正常物理现象而非代码错误。根据你使用的台北101参数，塔的固有频率约为：

omega0 = np.sqrt(k1/m1)
f0 = omega0/(2*np.pi)
print(f"固有频率：{f0:.3f} Hz")  # 输出约0.147 Hz

你设置的驱动频率为50Hz，远高于固有频率，此时受迫振动的响应是**瞬态响应（固有频率的低频振动）与稳态响应（驱动频率的高频振动）**的叠加：

• 瞬态响应是系统因初始条件产生的自由振动，会随阻尼逐渐衰减
• 稳态响应是驱动频率下的受迫振动，因驱动频率远高于固有频率，振幅很小

两者叠加后就会呈现出“低频大振动上叠加高频小振荡”的视觉效果，这是受迫阻尼振子的典型响应。

2. 代码中的可优化点 #

虽然核心逻辑正确，但代码存在几个可以改进的地方：

(1) 不必要的超大时间点数 #

你设置了N=1000000，生成100万个时间点会大幅降低绘图效率。odeint是自适应步长求解器，内部不需要这么多点，插值到过多点反而可能引入微小误差。建议将N改为10000，既保证绘图平滑，又提升运行速度。

(2) 方程转换的冗余性 #

你的entry函数返回-omega²*Em*np.sin(omega*t)，随后在dX中又取负得到目标激励项，完全可以简化为直接在dX中计算，减少函数调用冗余：

def dX(X,t):
    x1, dx1 = X
    f = 50
    omega = 2*np.pi*f
    excitation = omega**2 * 1 * np.sin(omega*t)
    dv1_dt = - (k1/m1)*x1 - (f1/m1)*dx1 + excitation
    return [dx1, dv1_dt]

(3) 模块导入的冗余 #

你导入了from math import *，但代码中并未使用math模块的任何函数，反而可能与numpy的函数产生命名冲突（比如math.sin和np.sin），建议删除这行导入。

3. 验证稳态响应的方法 #

如果你想确认稳态响应的平滑性，可以延长仿真时间（比如duration=1000），让瞬态响应完全衰减，此时曲线会只剩下高频的稳态响应，呈现平滑的正弦曲线。

优化后的完整代码示例：

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

# 参数
m1 = 264000000.  # kg
f1 = 5000000.     # kg/s
k1 = 225000000.   # N/m
x0 = 0.           # 初始位移
dx0 = 0.          # 初始速度
duration = 1000   # 延长仿真时间
N = 10000         # 合理的时间点数
time = np.linspace(0, duration, N)

# 定义微分方程组
def dX(X, t):
    x1, dx1 = X
    f_drive = 50  # 驱动频率
    omega = 2 * np.pi * f_drive
    # 激励项对应原方程右侧
    excitation = omega**2 * 1 * np.sin(omega * t)
    dv1_dt = - (k1/m1)*x1 - (f1/m1)*dx1 + excitation
    return [dx1, dv1_dt]

# 求解
result = odeint(dX, [x0, dx0], time)

# 绘图
plt.plot(time, result[:, 0])
plt.xlabel('时间 (s)')
plt.ylabel('位移 (m)')
plt.title('塔在正弦激励下的位移响应')
plt.show()

总结 #

你看到的振荡是驱动频率远高于固有频率时的正常物理响应，代码核心逻辑正确，优化上述点后可以更高效地验证结果。

内容的提问来源于stack exchange，提问作者ASTARDANGER