谐波平衡法求解动力学方程。 线性非线性方程均可进行求解。弹簧振子突然开始蹦迪这事儿在非线性动力学里还真有可能发生。今天咱们不聊数值积分那套老办法直接上手谐波平衡法——这玩意儿处理周期响应就跟玩音游似的能直接抓住振动节奏的命门。先看个硬核例子Duffing方程。这货长得像标准弹簧振子但自带叛逆属性def duffing_equation(x, t, alpha, beta, delta, gamma, omega): dx1 x[1] dx2 -delta*x[1] - alpha*x[0] - beta*x[0]**3 gamma*np.cos(omega*t) return [dx1, dx2]注意那个三次方项betax*3这就是非线性的罪魁祸首。传统ODE解法在这会算到怀疑人生特别是想找稳态周期解的时候。谐波平衡法的骚操作在于直接假设解是傅里叶级数。比如我们猜解长这样N 3 # 取到三次谐波 t np.linspace(0, 2*np.pi, 100) cos_terms [np.cos(n*omega*t) for n in range(N1)] sin_terms [np.sin(n*omega*t) for n in range(1, N1)]接下来把假设的解代入原方程用最小二乘法逼着方程成立。实际操作时得构造一个超大的系数矩阵A np.zeros((2*N1, 2*N1)) for i in range(2*N1): for j in range(2*N1): # 处理各阶谐波的交叉乘积项 # 这里需要展开非线性项的谐波平衡关系 ...具体填矩阵元素时会发现非线性项引发的灾难——三次方项会产生交叉乘积项。这时候需要祭出三角恒等式大法把cos³ωt拆成3cosωt/4 cos3ωt/4这种形式。谐波平衡法求解动力学方程。 线性非线性方程均可进行求解。当系数矩阵填完后用最小二乘法求解coeffs np.linalg.lstsq(A, F, rcondNone)[0]解出来的系数直接对应各阶谐波的幅值。想看效果对比数值积分结果from scipy.integrate import odeint t np.linspace(0, 100, 5000) sol odeint(duffing_equation, [0, 0], t, args(1, 5, 0.2, 8, 1.2))把谐波平衡解和数值解画在同一张图上会发现前几个周期可能对不上——这很正常谐波平衡法抓的是稳态解。等瞬态响应衰减后两条曲线会完美重合。重点来了当遇到强非线性时比如beta100传统的谐波平衡可能需要手动调参。试试把N设为5然后在构造矩阵时处理五阶交叉项。这时候代码里会出现让人眼花缭乱的索引计算# 处理三次非线性项的三阶交叉 for m in range(-N, N1): for n in range(-N, N1): for p in range(-N, N1): if m n p k: # 计算对应项的系数贡献 ...这种时候就会明白为什么有些论文里的公式长得像外星文字——好在代码实现时可以用张量运算优化不然嵌套循环能跑出咖啡凉。最后奉劝各位虽然谐波平衡法能优雅地处理周期解但碰到混沌系统还是得乖乖回去做数值积分。毕竟这方法的核心假设就是存在周期性稳态解当系统本身不按套路出牌时再优雅的数学技巧也得跪。