谐波平衡法求解非线性振动周期解的完整流程与代码实现

简介非线性振动中周期解的求解常采用谐波平衡法一份配套的MATLAB代码包可以帮助研究者与学习者快速验证算法、观察动态响应。该代码包面向机械、航空航天、土木工程等专业方向适合非线性振动课程、结构动力学分析以及相关科研预研。谐波平衡法将周期解近似展开为基频整数倍谐波的叠加通过非线性项展开与线性化求各阶幅值相位包内代码围绕这一思路编写涵盖主程序、非线性力定义、激励函数、响应计算、矩阵线性化等模块用户可直接运行主程序获取周期解也可通过调整参数或替换非线性项适配不同系统。压缩包共14个文件全部为.m脚本与函数整体仅6KB体积小、结构清晰便于逐段阅读算法流程各模块耦合度低支持独立修改与调用。该资源已有619人学习下载可作为课堂演示、课题研究及二次开发的实用起点。1. 谐波平衡法是求nonlinear vibration周期解最不依赖初值的一条路很多人在算非线性振动nonlinear vibration的周期解时第一反应是用Runge-Kutta把方程从零时刻积分到“看起来稳定”。实际上你会遇到两个麻烦瞬态过程太长步长取不好会导致高频分量被数值耗散吃掉而在强非线性下系统可能存在多个周期解时域积分只能给出一张“吸引子快照”你根本不知道在同一个激励频率下旁边还藏着另一个解。谐波平衡法的思路完全相反它先把周期解假设成一组Fourier级数再把原微分方程投影到各个谐波项上把求ODE初值的问题变成求解一组代数方程。NLvibration这类小工具的核心就是这个过程它不依赖时域积分初值也不怕解跳到别的分支上只要代数残差能收敛周期解就摆在明处。对做转子动力学、MEMS谐振器、电力电子振荡器的人这一招几乎绕不过去。本文直接讲清楚谐波平衡法的推导、参数设置、强非线性延拓和最终验证每一步都给能跑的代码。2. 从杜芬方程开始用谐波平衡法把微分方程折成代数方程2.1 周期解为什么能用Fourier截断而不丢主特征一个周期为 T 的稳态解如果满足Dirichlet条件总可以写成Fourier级数x(t)a0∑_{k1}^{∞} (ak cos(kωt)bk sin(kωt))谐波平衡法的核心假设是对于工程里常见的阻尼系统高频谐波在能量上占比极小保留到 N 阶就够。这里的 N 在NLvibration里通常叫“谐波数Harmonics”或“截断阶数”。N 不能拍脑袋定要看非线性项的强度。以杜芬方程为例x 2ζx x ε x^3 f cos(ωt)当 ε 较小时响应主要由基波控制N 取 1 或 2 就够但当 ε 大到 1 以上三次非线性会耦合出明显的 3 倍频成分N 至少要取到 3否则共振峰位置和跳跃点会完全算错。判断标准很直接算完后看最后一个谐波系数的幅值如果它和最大系数相比超过 1%就把 N 加 1。2.2 写出残差与谐波平衡条件将周期解截断到 N 阶待定系数有 2N1 个a0 以及各阶 ak, bk。把假设解代入运动方程方程左边会得到一个时间函数 R(t)它不再恒等于零而是包含一大堆高次谐波。谐波平衡条件就是让 R(t) 在 F0, cos(kωt), sin(kωt) 这些基函数上的投影为零∫_0^T R(t) dt 0 ∫_0^T R(t) cos(kωt) dt 0 (k1..N) ∫_0^T R(t) sin(kωt) dt 0 (k1..N)这一组积分方程等价于把 R(t) 自身也做Fourier展开让它的前 N 阶谐波系数全部置零。对于多项式非线性x^3 的Fourier系数是原始系数的卷积用手算很繁琐但用符号计算可以自动展开。2.2.1 用SymPy验证三次非线性的谐波耦合import sympy as sp omega, t, eps sp.symbols(omega t eps, positiveTrue) a1, b1, a3, b3 sp.symbols(a1 b1 a3 b3, realTrue) # 假设只保留基波和3次谐波且系统无常数项 x (a1*sp.cos(omega*t) b1*sp.sin(omega*t) a3*sp.cos(3*omega*t) b3*sp.sin(3*omega*t)) x3 sp.expand(x**3) # 提取cos(omega*t)的系数 proj_cos1 sp.integrate(x3*sp.cos(omega*t), (t, 0, 2*sp.pi/omega)) / sp.pi print(cos(omega t) in x^3:, sp.simplify(proj_cos1))这段代码做的事情是把三次非线性项在基波上的投影解出来。integrate是Fourier投影的解析实现除以 pi 是因为完备基的内积归一化。输出里会看到a1**3和a1*a3之类的项这说明基波系数和被截断的高阶系数相互耦合。实际手算很容易漏掉这种耦合这也是谐波平衡法“折成代数方程”后必须用计算机处理的原因。2.3 谐波平衡条件与虚功率平衡的等价性另一种理解方式把残差 R(t) 乘以每个基函数后在周期内积分物理上就是在周期内对“虚位移”做功为零。所以谐波平衡法也叫“Galerkin法”或“Ritz平均法”。在NLvibration这类工具里它底层就是一个非线性最小二乘问题minimize ||R_vector(X)||_2^2其中 X 是所有谐波系数的向量R_vector 是上述投影残差组成的向量。注意不要直接对这个平方和做梯度下降因为强非线性下目标函数高度非凸梯度法容易卡在局部极小。正确做法是直接对 R_vector(X)0 用Newton迭代这样收敛才是二阶的。2.3.1 NLvibration中“谐波数”参数如何影响方程维度谐波数 N未知量个数含常数项残差方程个数需要求解的线性系统规模1333×32555×53777×7N2N12N1(2N1)×(2N1)表格里的维度是每个激励频率点上的规模。谐振频率扫描时如果频率点有 200 个N 取 3就要解 200 次 7×7 的非线性方程组这个计算量微秒级完全不是瓶颈。真正的瓶颈是 Newton 迭代里每一步都要重新计算 Jacobian 对 N 的依赖N 增加一阶Jacobian 的填充量大约从 O(N^2) 涨到 O(N^3)。所以在NLvibration里如果看到“运算时间暴涨”先查是不是把基波和常数项之外的谐波数开到了 10 以上而不是查CPU。推导到这里谐波平衡法已经从“把一个微分方程变成一堆积分等式”落实到了“求解一组 F(X)0”。下面用具体代码把这个求解过程跑通。3. 用NLvibration/Python实现周期解的数值求解最少代码跑通Duffing3.1 离散谐波平衡的Newton迭代框架实现谐波平衡法的关键是写出残差函数 F(X)。X 的定义要统一把 cos(kωt) 系数和 sin(kωt) 系数按 k 从小到大排成一维向量常数项放在最前面。运动方程是x 2ζ x x ε x^3 f cos(ωt)假设截断到 N 阶x 用 X 组合得到。把 x、x、x 和非线性项 x^3 都投影到基函数上得到 F(X)。用 scipy.optimize.fsolve 或手写Newton。注意 fsolve 默认用前向差分求 Jacobian当 N 较大时数值 Jacobian 误差明显建议用numeric_jacobian或者手写解析 Jacobian。为了演示可复现我先给一份手写 Newton 的代码它不依赖任何ODE积分器只依赖 numpy。import numpy as np def duffing_hbm(omega1.2, eps0.5, zeta0.05, f0.3, N3): 谐波平衡法求解杜芬方程 x2*zeta*xxeps*x^3f*cos(omega*t) 返回谐波系数向量 X, 以及残差范数, 最后一项系数幅值 n 2 * N 1 # 未知数个数: a0, a1..aN, b1..bN X np.zeros(n) # 初始猜测: 线性解的基波余弦项 # 对线性化系统 x2*zeta*xx f*cos(wt) # 稳态幅值 A f / sqrt((1-w^2)^2 (2*zeta*w)^2) # 相位偏移写在余弦项和正弦项的系数里 denom (1 - omega**2)**2 (2*zeta*omega)**2 A f / np.sqrt(denom) # 当 omega 接近1时相位接近-pi/2 phi np.arctan2(-2*zeta*omega, 1-omega**2) X[1] A * np.cos(phi) # a1 X[N1] A * np.sin(phi) # b1 (顺序: 索引0是a0, 后面接a们再后面接b们) def projection_basis(k, order): # 返回 cos(k*w*t) 或 sin(k*w*t) 在一个周期上的采样向量 t np.linspace(0, 2*np.pi/omega, 2048, endpointFalse) if order 0: return np.cos(k*omega*t) else: return np.sin(k*omega*t) def residual(X): # 用连续时间采样计算残差 R(t)再投影到基函数 t np.linspace(0, 2*np.pi/omega, 2048, endpointFalse) # 重构 x(t) x X[0] * np.ones_like(t) for k in range(1, N1): x X[k] * np.cos(k*omega*t) X[Nk] * np.sin(k*omega*t) # 速度与加速度 xdot np.zeros_like(t) xddot np.zeros_like(t) for k in range(1, N1): xdot - X[k]*k*omega*np.sin(k*omega*t) X[Nk]*k*omega*np.cos(k*omega*t) xddot - X[k]*(k*omega)**2*np.cos(k*omega*t) - X[Nk]*(k*omega)**2*np.sin(k*omega*t) R xddot 2*zeta*xdot x eps*x**3 - f*np.cos(omega*t) # 投影 F np.zeros(n) F[0] np.mean(R) # 常数项投影 for k in range(1, N1): F[k] 2*np.mean(R*np.cos(k*omega*t)) F[Nk] 2*np.mean(R*np.sin(k*omega*t)) return F # Newton 迭代 for it in range(50): F residual(X) normF np.linalg.norm(F, ordnp.inf) if normF 1e-10: break # 数值 Jacobian (中心差分) J np.zeros((n, n)) h 1e-6 for j in range(n): Xp X.copy(); Xp[j] h Xm X.copy(); Xm[j] - h J[:, j] (residual(Xp) - residual(Xm)) / (2*h) # 解线性方程 dX np.linalg.solve(J, -F) X dX # 阻尼Newton: 如果残差变大就折半 while np.linalg.norm(residual(X), ordnp.inf) normF: dX * 0.5 X X - dX # 注意: 重新计算 if np.linalg.norm(dX) 1e-14: break last_amp np.hypot(X[N], X[2*N]) if N 1 else 0.0 return X, np.linalg.norm(residual(X), ordnp.inf), last_amp3.2 代码背后的三个关键参数第一个关键参数是omega激励频率。它决定了基函数的周期。谐波平衡法只在固定的 omega 下求解相当于扫频时每个频率点都是独立求根。第二个是zeta阻尼比。阻尼太大时高阶谐波会被压制N 可以取小阻尼接近零时共振峰很尖锐Newton 迭代容易从峰的一侧跳到另一侧需要在下一章讲延拓。第三个是eps非线性系数。eps 为 0 时方程退化为线性系统只有单一解谐波平衡法直接退化成频响函数eps 增大后共振峰向右侧弯曲硬弹簧特性同时出现多解区间所以它才是整个代码里最需要关注的值。代码中的2048个采样点是对残差做数值投影。这个点数的选择也有讲究因为余弦和正弦函数正交性依赖周期如果频率点取得不正好是周期端点就会泄漏。这里用np.linspace(0, 2*pi/omega, 2048, endpointFalse)正好覆盖一个整数周期可以避免泄漏。如果点数太少比如 128 点在 N 大于等于 5 时高频投影会出现明显混叠导致 Newton 迭代在达到机器精度前就停滞。点数也不用太多2048 对双精度浮点和 10 阶以内的谐波已经绰绰有余。3.3 从线性解起步为什么是可靠的第一个猜测非线性方程求根不像线性方程Newton 迭代必须给一个靠得住的X0。线性解起步是最自然的先把非线性项去掉得到线性频响函数。在线性系统里x(t) 的振幅和相位可以解析给出把它作为 N 阶谐波解的初始猜测在中等非线性强度下 Newton 一般三到五步就能收敛。如果 epsilon 较大线性解作为初始猜测可能落在牛顿法的收敛域之外这时 NLvibration 的做法是“增量加载”先把 eps 设成 0.1 跑一遍收敛后把结果作为 eps0.2 的初值逐步升到目标值。3.3.1 收敛失败时先看残差曲线而不是先调初值很多人在谐波平衡法不收敛时第一反应是改初始猜测。实际上更有效的做法是画残差函数R(t)的时域曲线。如果残差在单个周期内呈现光滑波动但投影后的 F 范数降不下去这说明 N 截断不够只增大谐波阶数即可。如果残差曲线呈锯齿状则是采样点数不足或基函数内积泄漏。如果残差在某些时间段尤其大且 N 增加后残差峰值没有下降那问题一定出在非线性项投影符号上——比如把 x^3 的系数符号写反或者漏了常数项。建议在调试时把residual(X)返回的 F 也返回时域残差 R打印几个典型时刻的值。4. 非线性振动分析中的强非线性问题弧长延拓与多解追踪4.1 为什么共振区附近牛顿法会跳变接近共振峰时谐波平衡方程组的解曲线在幅值-频率平面上呈现 S 形。S 形的上下两个分支是稳定解中间分支是不稳定解。用固定频率点做逐点扫频时Newton 迭代的结果会在某个频率点上突然从低幅值分支跳到高幅值分支这不是程序bug而是因为牛顿法是在找“最近的根”而 S 形区域里同一个频率下存在三个根初始猜测决定了它收敛到哪一个。要完整画出这条 S 形曲线必须沿着解曲线本身推进而不是沿着频率推进。4.2 用伪弧长延拓扫频的落地方式伪弧长延拓的思想是把频率 omega 也当成未知量引入一个弧长参数 s额外添加一个约束方程。NLvibration 的常见实现如下X X_prev ds * tangent_X omega omega_prev ds * tangent_omega然后对扩展后的方程组加上球面约束||X - X_new||^2 (omega - omega_new)^2 - ds^2 0这里的ds是步长。步长不能固定要加上自适应逻辑本轮 Newton 迭代超过 8 次才收敛就把 ds 减半少于 3 次收敛则下一轮扩大 1.5 倍。弧长延拓的另一个好处是能自然通过转向点saddle-node因为约束方程让迭代方向始终沿着解曲线走。# 伪弧长延拓的核心步骤伪代码省略Jacobian组装 # 已知点 (X0, w0)切向量 (dX0, dw0)步长 ds预测 X_pred X0 ds * dX0 w_pred w0 ds * dw0 # 校正用Newton法求解增广残差 # 其中残差 F(X,w)0 是谐波平衡残差新增约束 # C(X,w) (X-X0).T*(X-X0) (w-w0)**2 - ds**2 0 # 每轮迭代求解 (2N2) 维线性方程组 J_aug np.block([ [J_hbm, dF_dw], [2*(X-X0), 2*(w-w0)] ]) # 然后解 J_aug delta -residual_augJ_hbm是谐波平衡残差对 X 的 JacobiandF_dw是残差对 omega 的偏导。dF_dw 的解析式来自运动方程里 omega 只出现在 cos(omega t) 和 sin(omega t) 的自变量中以及激励项 f cos(omega t) 的频率位置。不要把 dF_dw 用差分近似因为靠近转向点时差分误差会导致切向量方向反号。4.2.1 弧长延拓结果如何判定多解区间延拓方向omega 变化幅值变化判定结果从低频向高频持续增加幅值先升后跳降存在跳跃从高频向低频持续减小幅值先升后跳升存在跳跃两个方向扫出的幅值曲线不重合频率区间重叠幅值不同多解区间确认实际操作中我会用向上扫频和向下扫频各跑一次把两条幅值曲线画在同一张图上。两张图在共振峰附近围出的滞后环就是多解区间。这个区间边界正好对应 S 形曲线的两个转向点。如果只用单方向扫频你永远不会意识到那个跳跃其实包含两个稳定的周期解和一个不稳定的周期解。4.3 周期解的稳定性判断Floquet理论还是简谐判据求出了周期解不等于它物理上能出现。NLvibration 里一般在得到谐波系数后计算单值矩阵Monodromy通过 Floquet 特征乘子实部是否穿过 1 来判断。对于单自由度系统有个更快的办法把周期解代入变分方程δx 2ζ δx (1 3ε x(t)^2) δx 0在周期解基础上做小扰动。如果 x(t) 的幅值在多个周期内衰减则稳定否则不稳定。在扫频延拓过程中观察 Jacobian 矩阵行列式是否改变符号是判断转向点的常用技巧但在转向点处 Jacobian 奇异行列式过零不能直接当作稳定性翻转。更稳妥的是跟踪单值矩阵最大特征乘子的模长变化。5. 验证周期解把谐波平衡结果交给时域积分做交叉检查5.1 用 RK4/odeint 对比一个周期内的漂移谐波平衡法给出的是周期解系数要验证它是否正确最直接的方法是把它作为初始条件扔给时域积分器积分若干周期看轨迹是否还停留在原始解的附近。以 scipy 的solve_ivp为例from scipy.integrate import solve_ivp def duffing_rhs(t, y, omega, zeta, eps, f): # 状态向量 y [x, v] return [y[1], f*np.cos(omega*t) - 2*zeta*y[1] - y[0] - eps*y[0]**3] # 假设 hbm_X 是上面谐波平衡法得到的系数向量N3 omega_val 1.2 t_span (0, 80*2*np.pi/omega_val) # 从谐波平衡解重构初始状态 t0 0 x0 hbm_X[0] sum(hbm_X[k]*np.cos(k*omega_val*t0) hbm_X[Nk]*np.sin(k*omega_val*t0) for k in range(1, N1)) v0 sum(-hbm_X[k]*k*omega_val*np.sin(k*omega_val*t0) hbm_X[Nk]*k*omega_val*np.cos(k*omega_val*t0) for k in range(1, N1)) sol solve_ivp(duffing_rhs, t_span, [x0, v0], args(omega_val, zeta, eps, f), rtol1e-10, atol1e-10) # 比较最后一个周期与谐波平衡解的形态 t_span_end np.linspace(sol.t[-200], sol.t[-1], 200)这段代码的验证逻辑是如果谐波平衡解是正确的周期解把它当初始状态积分 80 个周期后轨迹应该几乎不漂移。观察最后 200 个采样点的幅值变化如果相对误差小于 1e-6说明谐波截断充分如果漂移明显则意味着该周期解在动力学上不稳定即使代数上满足谐波平衡在实验中也不会出现。5.2 验证时注意的三个坑第一个坑是积分总时长要足够长否则瞬态衰减没有完你会把暂态漂移误判成解不稳定。我一般会做双保险分别积分 20 个周期和 80 个周期比较末尾一个周期的位移幅值。第二个坑是积分器容差设置太低。谐波平衡法本身可以到机器精度但如果 RK45 的容差只放到 1e-6可能把高阶谐波的误差放大。建议至少rtol1e-10。第三个坑是初值不能用“在某个时刻的瞬时位移和速度”组合而必须保证该初值严格位于重构的周期轨上。如果代码里补一个周期内的位移重构和 N1 个等距点的采样对比就能同时检查重构函数没有相位偏移。5.3 一个省事的小技巧残差能量百分比在 NLvibration 输出结果里除了谐波系数还应该输出一个能量残差指标eta sqrt(sum_{kN1}^{2N} (proj_k)^2) / sqrt(sum_{k1}^{N} (proj_k)^2)其中 proj_k 是把运动方程残差投影到第 k 阶基函数上的幅值。这个指标告诉你被截断掉的高阶谐波里还藏着多少残余能量。eta 小于 0.01 时说明当前 N 已经足够eta 在 0.01 到 0.05 之间时结果还能用但稳定性边界会有一点误差eta 大于 0.05 时增加 N 之前先检查采样点数和非线性项投影是否写对。把这个指标和时域交叉验证一起打到结果里比只贴一条幅值曲线更能说服自己谐波平衡法求出的周期解既代数可解又物理可达。本文还有配套的精品资源点击获取