简介面向非线性振动周期解分析需求这份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₀ Σ_{k1}^{N}[a_k cos(kωt) b_k sin(kωt)]把这一串系数当作未知数。代入运动方程后用三角函数在一个周期上的正交性对每个谐波投影方程就变成了 2N1 个代数方程。解这个代数方程组得到的是频域系数时域波形随时可以重构出来。换句话说HBM 把“初值问题”换成了“边值问题”代价是牛顿迭代的初值、收敛性和矩阵病态开始主导成败。普通 HBM 最大的难点在非线性力的投影三次方在频域还能硬推干摩擦、迟滞、间隙这类分段函数在频域基本写不出像样的闭式表达式。所以我在实现里一律走 AFTAlternating 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]共 2N1 个未知量 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*N1, 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*(j1); % 外力 F0*cos(wt) R(2*j1) Ls fnl_sin(j); end end这段代码先把频域系数还原到 nT 个时间采样点nT 取不小于 2N1 的 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*N1, 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 end3. 跑通最小算例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*N1, 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 给一个 2N1 长的零向量对多数弱非线性系统够用但强非线性系统上零初值几乎必挂这一点在第 5 章会专门展开。3.3 周期解重构从傅里叶系数回到时域波形算完幅频曲线别急着交差把某个频率的系数代回时域看一眼波形很多问题在这一步就现形了。比如取 w1.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 谐波数怎么选N1、N3 还是 N17N 决定内存和迭代代价也决定解的截断误差。对对称 Duffing 这类只有奇次谐波的问题N5 实际只用到了 1、3、5 三个频点N 取偶数时偶次谐波项虽然会参与迭代但解出来基本等于零白白浪费算力。如果你的系统带有明显非对称性例如单侧限位、含预载荷的间隙偶次谐波会真的长大那时就别再只取奇数。系统类型建议 N说明弱非线性、小振幅13基波为主3 次谐波只做校正中等立方非线性59保证 3、5 次谐波误差小于 1%干摩擦/间隙/迟滞1321波形含拐点高阶谐波衰减慢强非线性 冲击31 以上配合时域分段解析才可靠纯 HBM 有风险给一个自检办法跑完 N 阶解之后把 N 改成 N2比较两次结果的基波幅值 A₁相对变化小于 1% 就说明截断够了。另外看最高阶谐波系数在解向量里的占比如果 a_N、b_N 相对 a₁ 还在 1% 以上说明高频能量没有被压住波形在拐点附近必然失真。 注意不要把谐波阶数和 FFT 采样点数混为一谈N 是代数方程里的截断阶数nT 是时域采样点数两者独立但 nT 必须不小于 2N1否则投影误差会反过来污染牛顿迭代。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 ≥ 2N1采样数不够导致的投影误差足够破坏牛顿迭代这条在 N≥15 时特别容易踩。5. 谐波平衡法避坑NLvibration 调试中的 5 个典型问题与排查这一章是血泪经验区。下面五条我都亲眼见过而且每条都曾让整条幅频曲线看起来像某个同学喝多了画出来的。每条按「现象 → 原因 → 解决」直接说。5.1 现象一共振峰附近残差不降反涨现象扫频到峰值一侧牛顿迭代连续 20 步残差不减反增最终退出标记为未收敛。原因雅可比矩阵在峰附近条件数升高前向差分计算出的梯度方向和真实下降方向已经偏了迭代陷入振荡而不是收敛。解决先用阻尼牛顿做两步试探修正量长度超过当前解的量级时强制减半再把步长减半重新扫这段区间。别一上来就加大 N那个方向在峰附近帮不上忙只会让矩阵更大、迭代更慢。5.2 现象二低阶谐波解和时域积分对不上现象N1 时幅频曲线光滑但取几个频率点和四阶 Runge-Kutta 的稳态结果一对比差 5% 以上。原因N1 忽略了三次方带来的倍频能量实际响应的 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)把原方程在解附近线性化得到带周期系数的变分方程把一个单位扰动矩阵从 t0 积分到 tT2π/ω得到的单值矩阵的特征值就是 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(i1) - 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% 就回头查谐波阶数和采样点数。这套双验证跟了我的谐波平衡法好几年唯一一次把解算错是忘了检查时域重构里的采样间隔后来把这条检查写成固定注释挂在重构函数头里就再没犯过。希望帮到你。本文还有配套的精品资源点击获取
