简介这份MATLAB谐波平衡法实现资源面向从事非线性电路、电力电子装置及动力系统频域分析的工程师、科研人员和学生针对手工推导谐波平衡方程组繁琐、难以快速验证的问题围绕将时域非线性方程经傅立叶级数展开、在频域离散为代数方程组并迭代求解的完整流程展开。压缩包整体67KB共58个文件主体为53个m脚本另含README说明、License、CITATION等配置目录按核心求解、工具函数、三维扩展、测试用例、参数设置等模块清晰划分。其中核心求解模块覆盖谐波平衡迭代、稳态响应与幅值曲线计算工具函数提供FFT/IFFT变换、谐波打包/拆分、雅可比矩阵与Floquet稳定性分析等支撑测试模块包含对平衡求解、激励赋值、频响函数等算例的验证脚本便于对照检查实现正确性。已有111人学习下载整套代码兼顾可读性与扩展性适合希望结合具体实现深入理解非线性频域分析细节并在此基础上开展二次开发或教学研究的MATLAB使用者。1. 谐波平衡法不玄学在 MATLAB 里它就是求一组傅里叶系数做非线性周期振动分析时最常见的尴尬是数值积分要跑几十上百个周期才能进入稳态然后你才发现步长不够结果全是假的。谐波平衡法不一样它直接假设稳态解是有限项傅里叶级数把微分方程扔进残差里再要求残差在每个谐波分量上为零。落到 MATLAB 中谐波平衡法的实现就是「搭好残差函数交给 fsolve 解一组非线性代数方程」几秒钟就能拿到频响点做多解和跳跃分析也方便。这套方法适合非线性振动、电路振荡、传动系统周期响应也适合新手在 MATLAB 里跑通第一个频域稳态解。下面不抄理论推导直接给可复现的代码和调参经验重点讲清楚标题里的 Implementation 到底写在哪个文件里。2. 从时域到频域谐波平衡法的数学结构与三条 MATLAB 落地方案2.1 残差到底是谁的残差一个周期解怎么变成一组代数方程拿一个最简单的单自由度非线性系统做例子[ m x c x k x f_{nl}(x, x) F \cos \omega t ]假设它存在稳态周期解周期是 (T2\pi/\omega)。谐波平衡法不把时间历程积分出来而是直接把解写成有限项傅里叶级数[ x(t) a_0 \sum_{k1}^{N} \left( a_k \cos k\omega t b_k \sin k\omega t \right) ]这里的 (N) 是截断阶数。把形式代入原方程后方程不再逐点等于零而是产生一个时域残差[ r(t) m x c x k x f_{nl}(x,x) - F \cos \omega t ]我们要做的是选一组谐波基 ({1, \cos\omega t, \sin\omega t, \cos2\omega t, \sin2\omega t, \dots})让残差 (r(t)) 在所有基上投影为零。这是 Galerkin 意义上的“平衡”。投影之后原来的微分方程变成 (2N1) 个非线性代数方程未知数就是 (a_0, a_1,\dots,a_N, b_1,\dots,b_N)。这是谐波平衡法最核心的一步把边值型周期问题转化成零根问题。从离散角度看如果我们在一个周期内取 (M) 个时刻 (t_j)把 (r(t_j)) 算出来再分别乘上 (\cos k\omega t) 和 (\sin k\omega t) 做数值平均本质上就是离散傅里叶投影。(M) 足够大时这个数值投影和解析积分投影几乎一样。这个「时域采样 投影取系数」的写法是 MATLAB 代码能把谐波平衡法写短的关键后面第 3 章的完整代码就是按它来的。如果激励频率 (\omega) 已知未知量就是 (2N1) 个如果 (\omega) 本身也要解出来比如做自由振动或幅频关系曲线就需要把它加进未知量同时补一个附加条件通常是固定某一谐波的相位或幅值否则方程少一个约束求解器会犯迷糊。这里有三件事必须在写函数之前定下来谐波截断数 (N)、时域采样点数 (M)、基频 (\omega) 是已知还是未知。(N) 决定你能捕捉到多少非线性成分(M) 决定数值投影有没有混叠(\omega) 是否未知决定你解的是单个响应点还是一整条频响曲线。很多人一上来把这三个数全拍脑袋定然后回头追着 fsolve 报错调试方向就错了。2.2 三条 MATLAB 落地方案怎么选fsolve、手写 Newton 还是 AFT我见过不少人对“MATLAB 中谐波平衡法的实现”有一种误解以为必须先手推出一大堆三角函数的积分表达式再写成脚本。其实常见做法有三条按落地难度从低到高排第一条残差函数 fsolve。把上面 (2N1) 个投影写成一个函数 (R(X))(X) 是傅里叶系数向量直接调fsolve(R, X0)。优点是代码短非线性项怎么奇怪都无所谓因为只需要在函数里算一遍时域值缺点是 fsolve 默认用数值差分 Jacobian自由度大时每次迭代要做 (2N1) 次残差计算速度会下来。第二条解析 Jacobian 自写 Newton 迭代。把 (R(X)) 对 (X) 的导数解析算出来写成 Newton 迭代[ X_{k1} X_k - J(X_k)^{-1} R(X_k) ]优点是迭代次数少在强非线性、接近分叉点时也能控制缺点是要推 Jacobian尤其遇到摩擦、间隙、分段非线性时表达式会非常长、非常容易错。一般用来做小规模高精度求解或者给 fsolve 提供Jacobian选项让求解器不再做数值差分。第三条AFT-HBM也就是混合时频谐波平衡。求残差时不在频域硬算非线性项而是把 (X) 映射回时域在采样点上完整算出非线性力再用 FFT 或内积取回频谱系数作为残差交给求解器。强非线性、迟滞、间隙这类问题基本都靠这条因为非线性项可以是一个函数求值你不用为每种非线性推一次积分公式。下表是我自己的选型依据方案什么时候用优点代价fsolve 数值差分第一次验证、N≤7代码短改方程快残差计算次数多慢解析 Jacobian Newton强非线性、需要高精度收敛快行为可控推导和调试成本高AFT 时频混合分段、迟滞、复杂非线性非线性项免推导要处理 FFT 系数对应关系我的习惯是先按第三种的简化版——时域采样加内积投影——把整个流程跑通因为非线性项只是函数求值替换成本最低如果系统低频、弱非线性、谐波数少再试直接投影两者的结果应该一致只有自由度很大、数值 Jacobian 病态时才值得花时间写解析 Jacobian。大部分非线性振动问题 (N) 在 3 到 7 之间未知量最多十几二十个fsolve 完全扛得住所以方案一不是教学玩具而是真的能用来扫频的方案。3. 跑通最小实现残差函数、傅里叶基与 fsolve 的完整代码3.1 用单自由度三次非线性系统做靶子先确定方程和参数为了不把篇幅浪费在物理背景上用最常见的单自由度三次非线性系统作为测试对象[ x c x k x \alpha x^3 F_0 \cos \omega t ]参数我习惯这样设(m1)(c0.05)(k1)(\alpha0.5)(F_01.0)激励频率 (\omega1.2)。这个组合的响应并不弱三次项会把激励频率的 3 倍成分带出来刚好可以用来验证谐波数 (N) 的作用。质量归一化之后方程就变成[ x -c x - k x - \alpha x^3 F_0 \cos \omega t ]选 (N3)未知量有 7 个(a_0,a_1,a_2,a_3,b_1,b_2,b_3)。如果系统是完美的反对称非线性(a_0) 和一些偶次项理论上接近零但开始时不要假设任何系数为零让求解器自己判断这样能看出初值偏移和收敛问题。时域采样点 (M) 要先于代码定好。因为 (x^3) 会把最高 (3N) 次谐波带回时域采样点至少要比 (3N) 高才不容易混叠。我惯用 (M4N1)这里就是 13 个点。实际算出来之后可以把 (M) 翻倍看结果是否不变这是一个便宜又必要的收敛性检查。3.2 构造残差函数直接投影法的 MATLAB 代码下面这段是谐波平衡法实现的核心它把傅里叶系数向量 (X) 映射成时域位移、速度和加速度组装出残差再投影回各个谐波基。代码可以直接复制保存为hbm_residual.m。function R hbm_residual(X, w, N, M, m, c, k, alpha, F0) % X [a0; a1; ...; aN; b1; ...; bN] % 返回 2*N1 个残差fsolve 会令它等于 0 % 1. 周期内均匀采样去掉末端重复点 t linspace(0, 2*pi/w, M1); t t(1:end-1); % 2. 构造傅里叶基矩阵 C0 ones(M, 1); Ccos zeros(M, N); Csin zeros(M, N); for h 1:N Ccos(:, h) cos(h * w * t); Csin(:, h) sin(h * w * t); end % 3. 从系数映射到时域位移、速度、加速度 x C0 * X(1) Ccos * X(2:N1) Csin * X(N2:2*N1); xd -Ccos * (w * (1:N) .* X(N2:2*N1)) ... Csin * (w * (1:N) .* X(2:N1)); xdd -Ccos * ((w * (1:N)).^2 .* X(2:N1)) ... - Csin * ((w * (1:N)).^2 .* X(N2:2*N1)); % 4. 时域残差把微分方程左边减右边 r m * xdd c * xd k * x alpha * x.^3 - F0 * cos(w * t); % 5. 离散投影对常数项、余弦、正弦分别做内积 dt (2*pi/w) / M; R(1) dt * sum(r); for h 1:N R(h1) dt * sum(r .* cos(h * w * t)); R(N1h) dt * sum(r .* sin(h * w * t)); end R R(:); end逻辑上要注意三点。第一X(2:N1)放余弦系数X(N2:2*N1)放正弦系数顺序必须一致否则速度和加速度表达式会错位。第二速度表达式里余弦项贡献负正弦正弦项贡献正余弦符号不要搞反。第三最后投影时用矩形积分因为 (M) 足够大时它和真实内积的误差低于我们关心的收敛阈值如果你追求更高精度可以换成 Simpson 积分但非线性项存在时意义不大不如先检查混叠。3.3 让 fsolve 输出解再用 ode45 做一致性验证残差函数写好后剩下的事就是给它一个初值。这个初值不能是全零否则对某些参数 fsolve 可能停在零解上。我一般把基波余弦系数设为 1也就是假设系统在线性区幅值量级为 1然后让求解器自己修正相位。% 主脚本 N 3; M 4*N 1; % 时域采样点经验值 w 1.2; m 1; c 0.05; k 1; alpha 0.5; F0 1.0; X0 zeros(2*N1, 1); X0(2) 1.0; % 猜 a1 1 opts optimoptions(fsolve, ... Display, iter, ... Algorithm, trust-region, ... FunctionTolerance, 1e-10, ... OptimalityTolerance, 1e-12); % 用函数句柄固定参数 fun (X) hbm_residual(X, w, N, M, m, c, k, alpha, F0); Xsol fsolve(fun, X0, opts); % 提取幅值 a0 Xsol(1); a Xsol(2:N1); b Xsol(N2:2*N1); amp1 sqrt(a(1)^2 b(1)^2); amp3 sqrt(a(3)^2 b(3)^2); fprintf(a0 %.4e, 基波幅值 %.4e, 3次谐波幅值 %.4e\n, ... a0, amp1, amp3);这里我通常把Algorithm设为trust-region。fsolve默认的trust-region-dogleg在残差较小的零根问题上表现不错但当你把谐波平衡残差写成这种非最小二乘形式时trust-region对边界和缩放更稳。具体用哪个要看报错没有绝对。验证不能只看 fsolve 说收敛就收工。常见做法是拿 ode45 做一次时域仿真跑足够多周期取后半段稳定波形做 FFT对比基波和三次谐波幅值。% 用 ODE45 做对照 odefun (t, y) [y(2); -c*y(2) - k*y(1) - alpha*y(1)^3 F0*cos(w*t)]; [tod, yod] ode45(odefun, [0, 400*2*pi/w], [0.1; 0]); % 取最后 20% 的时域响应忽略瞬态 nt round(0.8 * length(tod)); yt yod(nt:end, 1); tt tod(nt:end); % 用 FFT 提取基波和三次谐波幅值 fs 1 / (tt(2) - tt(1)); Y fft(yt); freqs (0:length(Y)-1) * fs / length(Y); Amp 2 * abs(Y(2:4)) / length(Y); % 近似取前三项这段验证代码的关键是初值别给得太离谱否则 ode45 也要跑很久才进稳态。如果你发现 fsolve 解出的幅值和 ode45 对不上优先怀疑 (M) 不够而不是怀疑 HBM 原理有问题。4. 从单点解到频响曲线谐波数、初值与频率延拓的配合4.1 谐波截断的代价为什么 N3 不够用(N3) 能算不代表任何工况都能用。三次非线性项 (x^3) 在频域会把基波的三倍频成分带出来如果系统本身有超谐波共振比如固有频率大概在三倍激励频率附近那么 (N3) 刚好能捕捉到一部分但如果系统里还有迟滞、间隙或摩擦非线性可能产生连续频谱成分截断误差会直接污染前面的低阶系数。我一般会根据系统非线性最高次数来估 (N)。对于多项式型非线性最高幂次是 (p)至少要取 (N \ge p) 才能覆盖最主要的超谐波。对于强非线性比如干摩擦、分段刚度谐波成分衰减很慢(N5) 起步、(N7) 不新鲜个别场景要到 11 或 15。(N) 越大fsolve 的未知量越多Jacobian 数值差分越贵所以不要无脑加大。一个实用的判断方法是“谐波幅值衰减比”。算出 (a_N,b_N) 后如果最后一个谐波的幅值比基波低三四个数量级通常认为截断足够如果最后一个谐波幅值还在基波的 1% 以上必须加 (N)。这个检查应该写进脚本代替肉眼判断。4.2 初值决定你拿到哪一支解双稳态和跳跃怎么处理三次硬化弹簧的频响曲线是弯的在某个频率区间内存在三个解上支、下支和中支通常不稳定。fsolve 是局部收敛算法初值往哪给它就往哪个分支跑。这不是 bug而是你必须要利用的行为。如果你只想要某一个频率点的响应初值可以从线性幅值公式给[ a_1 \approx \frac{F_0}{\sqrt{(k - \omega^2)^2 (c\omega)^2}} ]代入参数后得到 (a_1) 的近似值再微调相位。更稳的做法是设计参数延拓从远离共振的低频或高频端开始以很小的频率步长扫过去每一步用上一步的解作为初值。这样上支和下支都能连续追踪不会因为初值乱跳而丢失分支。遇到跳跃点时fsolve 可能残差下降一半就卡住说明当前频率处解分支发生了鞍结分岔需要改用弧长延拓或者干脆把 (\omega) 也当成未知量。4.3 一个不调包也能用的扫频延拓模板下面这段是轻量级频率延拓没有用专业 continuation 工具箱只靠“上一解做初值”的思路。N 5; % 稍微提高谐波阶数 M 4*N 1; m 1; c 0.05; k 1; alpha 0.5; F0 1.0; wList 0.4:0.02:2.2; % 先粗扫后面可在跳跃处加密 Xprev zeros(2*N1, 1); Xprev(2) F0 / (k - 0.4^2); % 低频线性近似 store zeros(length(wList), 2*N1); for i 1:length(wList) w wList(i); fun (X) hbm_residual(X, w, N, M, m, c, k, alpha, F0); opts optimoptions(fsolve, Display, off, ... FunctionTolerance, 1e-10, OptimalityTolerance, 1e-12); Xsol fsolve(fun, Xprev, opts); store(i, :) Xsol; Xprev Xsol; % 关键把上一频率解传给下一个 end amp1 sqrt(store(:,2).^2 store(:,N3).^2); plot(wList, amp1, o-);这个模板的参数有三处值得调。一是频率步长共振点附近必须加密到 0.005 以下不然解会从上支直接跳到下支二是初值更新方式如果步长略大导致不收敛可以用两三个先前解的线性外推做初值if i 3 Xguess store(i-1,:) ... (store(i-1,:) - store(i-2,:)) * ... (wList(i)-wList(i-1)) / (wList(i-1)-wList(i-2)); end三是每步收敛后检查最后一项谐波幅值。如果它不再比前一项小很多说明这个频率点需要更多谐波可能需要把 (N) 加一档重新扫。这个扫频过程能直接看到硬化弹簧的共振峰向右弯曲以及跳跃点附近数值解的跳变验证 HBM 的多解能力。5. 避坑指南MATLAB 谐波平衡法最容易翻车的四个环节5.1 收敛性翻车初值、对称性和算法设置现象fsolve 每步残差都在下降但迭代到上限还没收敛。原因多半是初值离目标解太远尤其是强非线性段从线性近似初值出发可能落进一个局部平台区。解决思路不是硬撑 fsolve 的最大迭代次数而是从弱激励开始延拓先把 (F_0) 设成 0.1算收敛后把解作为 (F_00.2) 的初值逐步升到目标载荷。这套载荷延拓和频率延拓配合起来几乎能解决所有初值敏感问题。现象对称系统求解时Jacobian 接近奇异fsolve 提示近似奇异矩阵。原因是对称系统下偶次谐波系数理论上是零如果你仍然把它们全部设为未知量残差函数对这些方向的导数接近零Jacobian 就病态。解决方法是利用对称性删掉偶次项或者增加一个相位固定条件比如令 (b_10)把未知量从 (2N1) 降到 (2N)消除刚体自由度。注意不要直接砍变量要连带改残差函数里的傅里叶基矩阵。现象扫频到跳跃点附近fsolve 突然不收敛但两侧都能算。原因是在鞍结分岔点附近解的 Jacobian 奇异频率方向成为解的切线方向固定 (\omega) 的方程在这个点不可逆。简单解决是缩小频率步长让初值更接近前一步解更彻底的做法是让 (\omega) 也变成未知量补一个弧长约束方程也就是做最简弧长延拓。对大多数单自由度问题缩小步长已经够用。5.2 数值实现翻车采样点数、时间缩放和残余混叠现象把 (N) 从 3 加到 5低频系数反而出现明显变化甚至出现很小的直流分量。原因不是方程变了而是时域采样点 (M) 没跟着加。非线性项计算会产生最高 (pN) 次谐波如果 (M) 小于这个最高频率的两倍混叠会把高频能量折回低频污染前几阶投影。解决方法是按非线性最高幂次决定采样点我一般用 (M \ge (pN \times 2 1))对三次非线性就是 (M \ge 6N1)保守起见取 (8N) 也没问题。加 (M) 后如果结果不变说明混叠已经压住了。现象(\omega) 较大时基函数矩阵病态fsolve 经常失败。原因是时间 (t) 直接参与了 ( \cos(h \omega t))高频让采样点分布在一个剧烈振荡的周期上投影矩阵的条件数随 (\omega) 增长。解决方式是做时间无量纲化令 (\tau \omega t)基函数变成 (\cos(h \tau)) 和 (\sin(h \tau))(\tau) 的采样区间固定为 ([0,2\pi])与 (\omega) 无关(\omega) 只出现在速度、加速度表达式的系数里。这样既改善条件数也让扫频时残差函数结构保持稳定。现象同一个频率用矩形积分和更高精度积分算出来的幅值差 1%反复确认不是收敛问题。原因是矩形积分对高频残差分量衰减慢内积精度不够。解决方法是把残差投影里的四边形换到integral或者直接把投影那一步改成fft取系数。如果只是想快速对比最简单的办法是把 (M) 翻倍看结果是否进入一个稳定值。数值积分精度不足不是算法问题是离散化层数不够优先加 (M)。6. 进阶用 FFT 把非线性项自动化让谐波平衡法跑通任意方程组如果每换一个非线性模型都要手改投影循环很快就烦了。到这一步我会把残差函数里的“内积投影”换成“FFT 取前 N 阶系数”。思想是时域残差 (r(t)) 本身已经包含了所有非线性信息直接用fft(r_t)/M得到离散傅里叶系数再取前 (N) 阶的实部和虚部构成残差向量。这样非线性项只要能在时域算出来就能塞进同一个模板。function R hbm_afft_residual(X, w, N, M, m, c, k, alpha, F0) % 与上一版相同的未知量 X但投影用 FFT tau linspace(0, 2*pi, M1); tau tau(1:end-1); t tau / w; % 实际时间也可以直接用 tau 构造基矩阵 Ccos zeros(M, N); Csin zeros(M, N); for h 1:N Ccos(:,h) cos(h * tau); Csin(:,h) sin(h * tau); end x X(1) Ccos * X(2:N1) Csin * X(N2:2*N1); xd -w * (Ccos * ((1:N).*X(N2:2*N1)) ... - Csin * ((1:N).*X(2:N1))); xdd -w^2 * (Ccos * ((1:N).^2.*X(2:N1)) ... Csin * ((1:N).^2.*X(N2:2*N1))); r m*xdd c*xd k*x alpha*x.^3 - F0*cos(tau); fhat fft(r) / M; % 离散傅里叶系数fhat(1) 是直流实部对应 cos负虚部对应 sin R [real(fhat(1)); real(fhat(2:N1)); -imag(fhat(2:N1))]; end要点有两个。第一tau被固定在 ([0,2\pi])所以基矩阵只在扫频前构造一次复用给所有频率适合批量扫频第二(X) 的速度和加速度表达式必须乘上 (\omega) 和 (\omega^2)因为时域导数是对 (t) 求的而基函数是对 (\tau) 写的。这样改完之后非线性项从alpha*x.^3换成任意函数fnl(x,xd)都不影响后面的残差组装。我最后想说一个自己踩过的坑早期做某跨平台系统的周期解分析时我为了省时间只取 (N3)(M) 也没配够结果算出来的高次谐波在倍频处出现一条假峰我还以为是系统真有那个共振点后来把采样点加到 (M8N1) 才复原真实曲线。从那以后我每换一组参数都强制跑一遍“(M) 翻倍判稳”的检查宁可选保守的采样点数也不省那个几毫秒的矩阵运算。谐波平衡法的实现并不难难的是你要知道结果什么时候可信。希望帮到你。本文还有配套的精品资源点击获取