简介：面向非线性振动周期解分析需求，这份MATLAB代码包聚焦谐波平衡法的工程实现，适合机械、航空航天等领域的工程师与研究生快速上手周期解计算方法。压缩包共14个文件，全部为.m脚本，整体仅6KB，涵盖主程序、系统方程定义、非线性力构造、谐波响应求解等模块，便于按流程调用与二次开发。已有619人学习下载。代码围绕非线性常微分方程，演示了将非线性项泰勒展开、按谐波频率保留对应分量、线性化代数方程组并求解各次谐波振幅与相位的完整流程；用户只需修改系统参数，即可获得近似周期解并分析其动态特性，也可结合具体工程问题扩展为多谐波叠加或更高阶逼近方案，对教学演示与课题验证均具有实用价值。

1. NLvibration.rar 到底是什么：谐波平衡法求非线性振动周期解的老牌工具箱

搞非线性振动的人，早晚会在工位屏幕上面对同一个场景：隔振系统里带了立方刚度，轴承间隙里带着干摩擦，时域积分跑了几百个周期幅值还在漂，而课题组老资料夹里静静躺着一个 NLvibration.rar。这个以谐波平衡法（HBM）为核心的非线性振动周期解工具箱，解决的是把 m x″ + c x′ + k x + f_nl(x, x′) = F cos ωt 这类方程变成一组关于傅里叶系数的代数方程，一次性求出稳态周期解，再把幅频曲线的“弯”完整扫出来。它适合正在做非线性动力学课程设计、转子动力学或隔振设计的人，也适合想快速验证自己写的 HBM 代码的从业者。

2. 从方程到周期解：谐波平衡法推导与 NLvibration 的求解框架

在解压 rar、一个个翻 .m 文件之前，先花二十分钟把谐波平衡法的原理在草稿纸上过一遍，后面 debug 会快得多。很多人在非线性振动课上学 HBM，记住的只有一句“设傅里叶级数，代进去，同谐波平衡”，但真正动笔时发现三个地方根本没人讲清：谐波截断到几阶、非线性力怎么在频域里表达、以及迭代不收敛时该怀疑谁。下面这套推导和实现，是我在复现多种 HBM 工具时确认过的最稳路径。

2.1 非线性振动为什么难算：频域视角下的周期解

线性系统为什么好办，因为它满足叠加原理：简谐激励下稳态响应只含同一个频率的分量，振幅和相位用复数代数就能算出来。非线性系统（nonlinear vibration）里一旦出现立方刚度、间隙、干摩擦这样的环节，响应会自产谐波，能量在高次谐波之间来回串，幅频曲线也不再是单值曲线，而是带有弯折和多解区的曲线。用时域积分去扫频，等于每次都从头忍受瞬态衰减，阻尼小的时候一个频率点要积分几百个周期，遇到强非线性还得把步长掐到很小，扫完一整条曲线往往耗掉一整晚。

谐波平衡法就是绕开这条路的：直接假定稳态周期解存在，并写成有限项傅里叶级数 x(t) = a₀ + Σ_{k=1}^{N}[a_k cos(kωt) + b_k sin(kωt)]，把这一串系数当作未知数。代入运动方程后，用三角函数在一个周期上的正交性对每个谐波投影，方程就变成了 2N+1 个代数方程。解这个代数方程组得到的是频域系数，时域波形随时可以重构出来。换句话说，HBM 把“初值问题”换成了“边值问题”，代价是牛顿迭代的初值、收敛性和矩阵病态开始主导成败。

普通 HBM 最大的难点在非线性力的投影：三次方在频域还能硬推，干摩擦、迟滞、间隙这类分段函数在频域基本写不出像样的闭式表达式。所以我在实现里一律走 AFT（Alternating Frequency-Time）路线，也就是先把频域系数逆变换回时域波形，在时域把非线性力算出来，再做一次 FFT，把力的各阶谐波系数抓出来。这个过程看起来多了一步，实际最省事、最不容易出错，也是 NLvibration 这类老工具包普遍采用的做法。

2.2 求解流程：残差函数、雅可比矩阵与牛顿迭代

这类工具箱的骨架，我一般按五步走：先把 m、c、k、非线性刚度和外激励幅值写进一个参数对象；选谐波阶数 N 和频率扫描区间；在起始频率给一组初值，通常从零解或同频率线性解起步；构造残差函数，用有限差分或解析式生成雅可比矩阵，做牛顿迭代；最后把上一个频率的解作为下一个频率的初值，继续扫。后两步是整个流程的心脏，也是每个人会写出不同风格的地方。

残差函数里有一个容易踩的细节：外激励 F cos ωt 只投影到一次余弦项上，如果你把方程左右除以质量后忘掉 F₀ 的系数，幅频曲线的峰值会整体错位；另外零次谐波 a₀ 对应静态偏移，对称 Duffing 系统理论上 a₀=0，但迭代中不能因为理论值是零就直接把它从未知数里删掉，否则数值上少一个自由度，多个方程之间互相牵制，收敛反而变差。下面是我常用的残差实现，未知量按 [a₀; a₁; b₁; …; a_N; b_N] 排布：

% residual_hbm.m % 构造谐波平衡法残差 R(u, w)，u 为傅里叶系数向量 % u = [a0; a1; b1; a2; b2; ...; aN; bN]，共 2N+1 个未知量 function R = residual_hbm(u, w, p, N) a0 = u(1); a = u(2:2:end); % cos(kwt) 系数 b = u(3:2:end); % sin(kwt) 系数 % 在一个周期内等间距采样，点数取 2 的幂 nT = 2^nextpow2(2*N + 1); t = (0:nT-1)' * (2*pi/(nT*w)); % 时域位移重构 x = a0 + cos(w*t*(1:N))*a + sin(w*t*(1:N))*b; fnl = p.k3 * x.^3; % 非线性力（这里为立方刚度） % 力在频域的 0 次与 1..N 次投影 fnl0 = mean(fnl); fnl_cos = zeros(N,1); fnl_sin = zeros(N,1); for j = 1:N fnl_cos(j) = 2*mean(fnl .* cos(j*w*t)); fnl_sin(j) = 2*mean(fnl .* sin(j*w*t)); end % 组装残差 R = zeros(2*N+1, 1); R(1) = p.k*a0 + fnl0; % 零次：静态力平衡 for j = 1:N % 线性部分：M x'' + C x' + K x 在频域的投影 Lc = (p.k - p.m*j^2*w^2)*a(j) + p.c*j*w*b(j); Ls = (p.k - p.m*j^2*w^2)*b(j) - p.c*j*w*a(j); R(2*j) = Lc + fnl_cos(j) - p.F0*(j==1); % 外力 F0*cos(wt) R(2*j+1) = Ls + fnl_sin(j); end end

这段代码先把频域系数还原到 nT 个时间采样点，nT 取不小于 2N+1 的 2 的幂，是为了保证后面用均值投影时采样定理不被击穿；然后用 mean 做零次谐波积分，用 2mean(fnl .cos(jwt)) 做余弦投影，系数 2 来自正交归一化。线性部分的频域形式是直接对 x″ 和 x′ 做微分得到的：j²ω² 项带负号，阻尼项 c*jω 从余弦投影和正弦投影两侧交错出现。外力只在 R(2) 里减去 F₀，对应 1 次余弦，其它谐波没有外力分量，这个位置写错，扫出来的共振峰幅值就全错了。

雅可比矩阵我用前向差分而不是解析式，原因很实际：当非线性力换一种形式（比如从 x³ 换成间隙函数）时，解析式要全部重推，而差分公式一行不用改。差分步长 h 取 1e-6，这是针对系数量级约在 10⁻³ 到 1 之间的典型解。如果直接套在数值量级已经到 10⁶ 的解上，h 要相应放大到 1e-3 量级，否则差分完全被舍入误差吃掉。

% jacobian_hbm.m % 用前向差分近似雅可比矩阵，避免手推复杂非线性力的导数 function J = jacobian_hbm(u, w, p, N) R0 = residual_hbm(u, w, p, N); m = length(u); J = zeros(m, m); h = 1e-6; for i = 1:m up = u; up(i) = up(i) + h; J(:, i) = (residual_hbm(up, w, p, N) - R0) / h; end end

牛顿迭代的退出条件我同时看残差范数和修正量范数，两个条件多半先满足的是后者。注意这里我用的是 inf 范数，在低频段 a₀ 接近零、高频段高阶谐波接近零时，inf 范数比二范数更早、更明确地暴露某个分量没收敛。

% newton_hbm.m % 牛顿迭代解单个频率下的周期解系数 function [u, ok] = newton_hbm(u0, w, p, N) u = u0; ok = false; for iter = 1:50 R = residual_hbm(u, w, p, N); if norm(R, inf) < 1e-9 ok = true; break; end J = jacobian_hbm(u, w, p, N); du = -J \ R; u = u + du; if norm(du, inf) < 1e-10 ok = true; break; end end end

扫频时每步把上一频率的解直接搬过来当初值，这叫简单继续法，便宜但容易在弯折区域翻车。如果是第一次跑通，用这个版本先验证整体逻辑，再决定要不要上伪弧长。下面的 warning 分支在未收敛时保留上一步解，是为了不中断整个 for 循环，扫完后你可以在图上看到断点和台阶，再针对性加密频率步长。

% hbm_sweep.m % 对频率列表做正向扫描，简单继续法提供初值 function [ws, U] = hbm_sweep(p, N, w_list, u_guess) ws = zeros(size(w_list)); U = zeros(2*N+1, length(w_list)); u0 = u_guess; for i = 1:length(w_list) w = w_list(i); [u, ok] = newton_hbm(u0, w, p, N); if ~ok warning('w=%.3f 迭代未收敛，继续沿用上一步解', w); u = u0; end U(:, i) = u; ws(i) = w; u0 = u; end end

3. 跑通最小算例：NLvibration 的调用方式与周期解重构

3.1 文件结构梳理与主函数入口

把 NLvibration.rar 解压后，第一件事不是读 README，而是用 dir 列一遍文件。这类老工具箱的命名通常很随性，hbm、hb、nltool、solve 混着用，靠眼睛找不如靠规则找：主入口脚本多数以 demo、example、test 或 main 开头，参数集中写在开头几行，后面是长长的计算和画图。我的习惯是先跑通 demo，再去翻它调用的求解器函数，这样能先看到预期输出，回头读源码时心里有底，不至于一进来就被十几二十个函数绕晕。

这类老工具箱里还经常出现 global m c k F0 的写法，第一眼看能跑，换个参数就要去改全局变量，复用性很差。我一般会把参数全部收进一个结构体 p，所有函数只传 p 不传散点参数，这样后续做参数扫描时只需要在脚本里改 p，不需要碰求解器内部逻辑，也方便以后把代码改写成函数式模块。

3.2 最小算例：Duffing 振子幅频曲线

用立方刚度的 Duffing 振子做最小算例，跑上面那套函数，参数如下：质量 1 kg，阻尼 0.02 N·s/m，线性刚度 1 N/m，此时线性固有频率 ωₙ=1 rad/s；立方刚度 k₃=0.3 属于中等硬化特性，共振峰向高频方向弯曲；激励幅值 0.1 N，让响应幅值保持在 0.1 到 0.5 这个量级，既能体现非线性，又不会让迭代因为幅值太大而频繁碰壁。

% demo_duffing.m % Duffing 振子幅频响应：HBM 最小算例 p = struct('m', 1.0, 'c', 0.02, 'k', 1.0, 'k3', 0.3, 'F0', 0.1); N = 5; % 谐波阶数 w_list = (0.5:0.01:1.8)'; % 扫频范围，单位 rad/s u_guess = zeros(2*N+1, 1); % 从零解起步 [ws, U] = hbm_sweep(p, N, w_list, u_guess); A1 = sqrt(U(2,:).^2 + U(3,:).^2); % 基波幅值 plot(ws, A1, 'b-', 'LineWidth', 1.5); xlabel('omega (rad/s)'); ylabel('A_1');

这段脚本里 w_list 用 0.01 rad/s 的步长从 0.5 扫到 1.8，覆盖了共振峰前后足够宽的范围。注意起始频率取 0.5 而不是 0，是因为零频率附近对应的是静态解，HBM 的频域假设在 ω→0 时会有数值退化；u_guess 给一个 2N+1 长的零向量，对多数弱非线性系统够用，但强非线性系统上零初值几乎必挂，这一点在第 5 章会专门展开。

3.3 周期解重构：从傅里叶系数回到时域波形

算完幅频曲线，别急着交差，把某个频率的系数代回时域看一眼波形，很多问题在这一步就现形了。比如取 w=1.05 时，如果画出的是一个多峰畸变波形，说明谐波数取低了；如果波形光滑正弦但幅值和时域积分对不上，多半是投影系数或外力项出了问题。

% 周期解重构：把傅里叶系数变成时域波形 w_index = find(abs(ws - 1.05) < 1e-6, 1); w = ws(w_index); u = U(:, w_index); nt = 256; t = (0:nt-1)' * (2*pi/(nt*w)); x = u(1) + cos(w*t*(1:N))*u(2:2:end) + sin(w*t*(1:N))*u(3:2:end); plot(t, x, 'b-');

这里重构采样点数取了 256 个，远超残差函数里的最少采样点数，画图足够光滑；t 覆盖 0 到 2π/w，正好一个周期。把重构波形和四阶 Runge-Kutta 积分结果的末端一个周期叠在一起画，是最直接的交叉验证。很多人只对数值、不对波形，结果幅值对了但相位差半个周期，这类错误在后续做稳定性判断时会被 Floquet 乘子放大暴露出来。

4. 谐波平衡法三大参数：谐波阶数、频率步长与收敛判据怎么定

HBM 在实际使用里的手感，基本由三个参数决定：谐波阶数 N、扫频步长和牛顿迭代容差。这三个参数之间互相牵制：N 给得越大，雅可比矩阵越大，单次牛顿迭代越贵，能容忍的扫频步长也越窄；步长加密之后，继续法的初值好，迭代次数反而少。下面把三条线的调参经验一次说清，也都是 NLvibration 这类包最常被改的参数位置。

4.1 谐波数怎么选：N=1、N=3 还是 N=17

N 决定内存和迭代代价，也决定解的截断误差。对对称 Duffing 这类只有奇次谐波的问题，N=5 实际只用到了 1、3、5 三个频点；N 取偶数时偶次谐波项虽然会参与迭代，但解出来基本等于零，白白浪费算力。如果你的系统带有明显非对称性，例如单侧限位、含预载荷的间隙，偶次谐波会真的长大，那时就别再只取奇数。

给一个自检办法：跑完 N 阶解之后，把 N 改成 N+2，比较两次结果的基波幅值 A₁，相对变化小于 1% 就说明截断够了。另外看最高阶谐波系数在解向量里的占比，如果 a_N、b_N 相对 a₁ 还在 1% 以上，说明高频能量没有被压住，波形在拐点附近必然失真。> 注意：不要把谐波阶数和 FFT 采样点数混为一谈，N 是代数方程里的截断阶数，nT 是时域采样点数，两者独立，但 nT 必须不小于 2N+1，否则投影误差会反过来污染牛顿迭代。

4.2 频率步长与初值继续法：如何把“弯”扫出来

幅频曲线弯折的根本原因是鞍结分岔：在分岔点附近，雅可比矩阵接近奇异，解对频率极敏感，步长稍大就会从一条分支跳到另一条，看起来就像扫出一条断开的曲线。常见做法是把简单继续法升级成带线性外推的版本，下面这段代码就是扫频循环里换掉初值那几行：

% 扫频循环里，用线性外推代替“上一解直接当初值” if i > 2 slope = (U(:, i-1) - U(:, i-2)) / (ws(i-1) - ws(i-2)); u0 = U(:, i-1) + slope * (w_list(i) - ws(i-1)); else u0 = U(:, i-1); end

外推斜率用上一个频率点的位移增量除频率增量，得到解随频率的局部斜率。线性外推在分岔点附近容易过冲，外推初值解不出来时，要退回简单继续法并用更小步长。步长的基准不是绝对数值，而是共振峰宽度；一般取峰值半功率带宽的 1/20 到 1/50 作为初始步长。峰附近自动减半三次以上还失败，就该怀疑不是步长问题，而是分支本身已经转向，需要换初值方向或切到伪弧长继续法。

4.3 收敛判据与雅可比病态的三条处理经验

第一条，容差按无量纲化后的残差范数取，我一般设 1e-9。如果方程里质量和刚度数量级差太大，先做坐标归一化，否则容差会被某个大数分量带着跑，看似收敛实际残差各项差了几十个数量级。第二条，有限差分雅可比在接近分岔点时几乎必然病态，condition number 会飙到 1e12 以上，此时强制用 J\R 得到的是垃圾修正量；我的处理是先检查 cond(J)，超过 1e10 就切换到阻尼牛顿，把迭代步长乘一个 0.3 到 0.5 的衰减系数，宁慢勿崩。第三条，当某一频率的迭代发散了，不要急着增大 N，先看时域采样点数 nT 是否满足 nT ≥ 2N+1，采样数不够导致的投影误差足够破坏牛顿迭代，这条在 N≥15 时特别容易踩。

5. 谐波平衡法避坑：NLvibration 调试中的 5 个典型问题与排查

这一章是血泪经验区。下面五条我都亲眼见过，而且每条都曾让整条幅频曲线看起来像某个同学喝多了画出来的。每条按「现象 → 原因 → 解决」直接说。

5.1 现象一：共振峰附近残差不降反涨

现象：扫频到峰值一侧，牛顿迭代连续 20 步残差不减反增，最终退出标记为未收敛。原因：雅可比矩阵在峰附近条件数升高，前向差分计算出的梯度方向和真实下降方向已经偏了，迭代陷入振荡而不是收敛。解决：先用阻尼牛顿做两步试探，修正量长度超过当前解的量级时强制减半；再把步长减半重新扫这段区间。别一上来就加大 N，那个方向在峰附近帮不上忙，只会让矩阵更大、迭代更慢。

5.2 现象二：低阶谐波解和时域积分对不上

现象：N=1 时幅频曲线光滑，但取几个频率点和四阶 Runge-Kutta 的稳态结果一对比，差 5% 以上。原因：N=1 忽略了三次方带来的倍频能量，实际响应的 3 次谐波虽然幅值小，但和基波耦合后改变了共振峰的位置和峰值。解决：把 N 加到 5 再对比，用上一章的尾系数占比做判断；另外注意时域积分的稳态判据要取后几个周期的最大幅值，很多对不上的案例其实是对比了瞬态，这在 vibration 工程验证里是非常经典的对拍误区。

5.3 现象三：牛顿迭代出现 NaN，波形直接炸掉

现象：某次迭代后解向量出现 NaN，程序没有直接报错，继续往外传，后面频率的解也跟着全挂。原因：迭代中间某个试探解幅值过大，x³ 在双精度下溢出；或者时间向量构造时用了 linspace(0, T, nT)，末尾点重复导致投影基准差了一个采样间隔。解决：给牛顿迭代加一条 line search，修正量超过振幅参考值 0.5 倍就整体缩短；时间向量改用 (0:nT-1)*T/nT 这种不含重复端点的写法。NaN 一旦出现，及时用 isfinite 检查残差，立刻中止本轮迭代，比让它传播下去好救得多。

5.4 现象四：正向扫频漏掉了多解区间

现象：从低到高扫的频率响应曲线，在理论预测应该有三解的频率区间只有一支，翻过来从高到低扫变成另一支，中间那支怎么都见不到。原因：简单继续法始终沿一个方向推进，在鞍结分岔点行列式变号处，解分支发生了转向，频率不再是好的“推进坐标”。解决：把频率改成弧长参数做伪弧长继续；不想上复杂算法的话，至少做一次反向扫频，把两条分支拼起来，中间的失稳分支在后续 Floquet 分析里自然会被标出来。拼曲线时注意频率坐标要对齐，这是谐波平衡法实际工程里重复出现的问题。

5.5 现象五：干摩擦或迟滞系统结果一直振荡

现象：残差整体收敛，但每两个频率之间幅值跳来跳去，曲线毛刺严重。原因：干摩擦力的不连续造成时域波形有拐点，Fourier 截断到 N 后产生吉布斯现象，投影得到的力谐波和真实力的谐波差一个固定的过冲。解决：把 N 提高到 15 阶以上，时间采样点数跟着放大到 64 个以上；如果还抖，用光滑化的 tanh 模型近似干摩擦，代替严格的 sign 函数，代价是摩擦力的拐点被抹圆，等效摩擦功略微偏低，但曲线稳定性和迭代成功率都明显改善。

6. 验证周期解稳不稳：Floquet 乘子计算与我的固定做法

HBM 解出来的每条分支都能画出来，但工程上必须知道它是不是物理上真实存在的稳定周期解。标准做法是在周期解 x̄(t) 上叠加一个无穷小扰动 δ(t)，把原方程在解附近线性化，得到带周期系数的变分方程；把一个单位扰动矩阵从 t=0 积分到 t=T=2π/ω，得到的单值矩阵的特征值就是 Floquet 乘子。乘子模都小于 1，解渐近稳定；有乘子模大于 1，解不稳定；出现模等于 1 的位置，正好对应幅频弯折处的分岔点。

% floquet_stability.m % 快速估算给定周期解的 Floquet 乘子 function mu = floquet_stability(u, w, p, N) T = 2*pi/w; nt = 1024; tgrid = linspace(0, T, nt); Phi = eye(2); for i = 1:nt-1 a = u(2:2:end); b = u(3:2:end); x_t = u(1) + cos(w*tgrid(i)*(1:N))*a + sin(w*tgrid(i)*(1:N))*b; A = [0, 1; -(p.k + 3*p.k3*x_t^2)/p.m, -p.c/p.m]; dtau = tgrid(i+1) - tgrid(i); Phi = (eye(2) + A*dtau) * Phi; end mu = eig(Phi); end

这里的时间推进用一阶 Euler 做示意，真正要精确就换成 ode45 对矩阵列向量积分，但一阶 Euler 配合 1024 个采样点足够判断乘子模的大致位置。然后把每个频率的最大乘子模画在幅频曲线下方，小于 1 的区间标绿色，大于 1 的标红色，这张“幅频 + 稳定性”组合图是目前振动行业最常用的 HBM 后处理方式。对多自由度系统，原理完全一样，只是 Φ 从 2×2 涨到 2n×2n。我现在固定做两件事：第一，每一条扫出来的分支都跑一次 Floquet，把模大于 1 的段标出来，这些段在时域积分里永远看不到，最容易骗人；第二，在分岔点附近把步长加密到原来的 1/4，用正向扫频和反向扫频各来一遍，两条曲线分岔点位置差超过 1% 就回头查谐波阶数和采样点数。这套双验证跟了我的谐波平衡法好几年，唯一一次把解算错，是忘了检查时域重构里的采样间隔，后来把这条检查写成固定注释挂在重构函数头里，就再没犯过。希望帮到你。

本文还有配套的精品资源，点击获取