企业官网建设流程全解析

简介：这份压缩包为非线性振动周期解求解提供了一套完整的谐波平衡法 MATLAB 实现，面向机械、航空航天、土木等专业的工程学习者、研究生及科研人员，用于处理力和位移呈非线性关系的复杂振动问题。包内共 14 个 m 文件，全部为 MATLAB 脚本，分别承担非线性力项构造、外部激励输入、泰勒展开后的线性化方程组组装、谐波系数求解以及响应合成等任务，整体体积仅 6KB，便于快速加载与调试。已有 619 人学习使用，常被用作课堂教学演示和科研初期验证工具。借助该程序，用户只需修改系统参数，即可得到包含振幅与相位的近似周期解，进而观察多稳态、分岔等非线性特征；同时，源码的模块结构也便于读者深入理解谐波平衡法从建模到数值实现的完整流程。

1. 谐波平衡法：非线性振动周期解里最省力气的频域路线

拿到一个叫 NLvibration 的压缩包，很多人的第一反应是去翻里面的求解器和文件结构，但真正值钱的往往是包名背后那条技术线：非线性振动、周期解、谐波平衡法。扫频实验中常见的共振峰“跳跃”现象——升频和降频扫出来的幅值曲线不一样——用时域积分复现常常要跑几百个周期，一次扫描以分钟计；谐波平衡法把微分方程压成一组关于傅里叶系数的非线性代数方程，同一问题通常零点几秒就能扫出一条频响曲线。这套思路适合透平叶片、减振器、MEMS 谐振器和转子系统的工程师与研究生：手里有一组运动微分方程，想知道系统在不同激励频率下的周期响应长什么样，哪些解在实际物理过程中留得住。

2. 谐波平衡法的数学骨架：把 Duffing 方程改写成频域残差方程组

工程里大量非线性振动问题都能化简成 Duffing 方程或其变体，所以后面的代码和参数都围绕它来写。先把方程亮出来：

m x'' + c x' + k x + α x³ = f cos(ωt)

2.1 为什么选 Duffing 方程做基准

m 是等效质量，c 是阻尼，k 是线性刚度，α 是三次非线性项的系数，右边是简谐激励。αx³ 可以来自大挠度薄板的几何非线性、磁悬浮轴承的力-位移关系，也可以来自超弹性减振器的恢复力。加了这一项之后，响应不再和激励一一对应：当激励频率扫过固有频率附近时，等效刚度会随振幅变化，共振峰向高频或低频方向弯曲，一个激励频率可能对应三个幅值解，其中两个是实际系统里能观测到的，另一个夹在中间。扫频实验里出现“跳跃”，根子就在这里。

做谐波平衡法之前，先明确一个基本态度：这个方法不是把原微分方程丢掉，而是假设周期解存在，再把“解满足方程”这种强条件，换成“残差在所有傅里叶基函数上的投影为零”的弱条件。弱条件换来的是问题规模大幅度缩小——一个需要步进积分的二阶常微分方程，变成一组十几个未知数的非线性代数方程。

2.2 未知量和残差的频域展开：从微分方程到代数方程

把响应 x(t) 写成有限傅里叶级数：

x(t) = a0/2 + Σ_{k=1..H} [ a_k cos(kωt) + b_k sin(kωt) ]

H 是谐波截断阶数，也是这个方法里唯一需要人工设置的“网格”参数。它决定解里允许出现多少倍频成分。对硬弹簧 Duffing，基频激励下解的主要成分是 cos(ωt) 和 sin(ωt)，但三次方非线性会把能量耦合到 3 倍频、5 倍频，所以 H=1 只能得到共振峰的骨架，H=3 以上才能看到超谐波对峰形的修正。

把展开代进运动方程，把一切移到同一边，得到时域残差：

r(t) = m x'' + c x' + k x + α x³ - f cos(ωt)

理论上，如果 x(t) 是精确周期解，r(t) 恒为零。但有限截断下它不可能恒为零，于是强迫残差在周期上正交，即每个频率分量投影为零：

∫_0^T r(t) cos(kωt) dt = 0，∫_0^T r(t) sin(kωt) dt = 0，k = 0,1,...,H

T = 2π/ω。这样得到 2H+1 个方程和同样多的未知数。方程里每一项都是傅里叶系数的多项式，这就是“谐波平衡”这个名字的由来：按频率分量逐个摆平。

2.3 解析展开与数值变换：两条路线怎么选

对残差做傅里叶投影，工程上有两条实现路线。第一条是纯解析展开：用三角函数恒等式把 αx³ 展开成各谐波系数相乘的形式。比如取 H=1，设 x = a1 cos(ωt)，则

x³ = a1³ cos³(ωt) = a1³ [ 3cos(ωt) + cos(3ωt) ] / 4

其中 3/4 a1³ 落在基频项上，相当于一个随幅值变化的等效刚度；cos(3ωt) 项则说明即使是单谐波假设，非线性也会向外辐射 3 倍频分量。解析展开在 H=1 时很漂亮，到 H=5 以上，交叉项多得可怕，维护成本指数上升。

第二条是数值路线：把 x(t) 和 r(t) 在一个周期内均匀采样 N 个点，用离散傅里叶变换或直接做数值积分投影，取出各频率分量的幅值作为残差向量。这个做法不关心非线性项的具体形式，只要能算出时域值就行。NLvibration 这类包里常见的谐波平衡求解器，多数走的也是数值路线。它的代价是失去了显式表达式带来的物理直觉，换来的是能用同一套代码处理干摩擦、间隙、迟滞这类根本写不出解析表达的强非线性。参数上只需要记住一个原则：采样点数 N 要明显大于 2H+1，工程上我习惯取 N = 4(2H+1) 或直接取 64、128 这类 2 的幂，避免高频分量折叠回基频区域造成混叠。

数值路线还有个隐藏优点：求傅里叶系数时可以用上一个频率点的解做初值，天然支持扫频延续。这一点到第 3 章写代码时会非常明显地体现出来。新手最容易犯的错误是一上来就开 H=7，然后得到一组残差全为零、但对扰动非常脆弱的虚假周期解。所以谐波平衡法的第一原则是：先用 H=1 把骨架解出来，确认主谐波的位置和跳跃区，再逐步加谐波。

3. 最小可运行实现：用 Python 给 NLvibration 搭一套谐波平衡求解核心

下面用 Python 实现一套最小但能直接跑的谐波平衡求解核心。无量纲 Duffing 方程取为：

x'' + 2ζx' + x + βx³ = f cos(ωt)

m=k=1，阻尼比 ζ、非线性系数 β、激励幅值 f 都是给定参数。物理时间 t 换成无量纲相位 φ=ωt，那么 d/dt = ω d/dφ，微分算子要在代码里乘上对应的 ω 幂次。

3.1 固定频率下的谐波平衡求解函数：先跑通一组系数

import numpy as np from scipy.optimize import fsolve def build_basis(H, N): phi = np.linspace(0.0, 2.0*np.pi, N, endpoint=False) cos_basis = [np.ones(N)] + [np.cos(k*phi) for k in range(1, H+1)] sin_basis = [np.zeros(N)] + [np.sin(k*phi) for k in range(1, H+1)] return phi, cos_basis, sin_basis def x_from_coeff(coeff, cos_basis, sin_basis, H): N = len(cos_basis[0]) x = coeff[0] * np.ones(N) for k in range(1, H+1): x = x + coeff[k] * cos_basis[k] + coeff[H+k] * sin_basis[k] return x def hb_residual(coeff, omega, zeta, beta, f, cos_basis, sin_basis, H): phi = np.linspace(0.0, 2.0*np.pi, len(cos_basis[0]), endpoint=False) x = x_from_coeff(coeff, cos_basis, sin_basis, H) xp = np.zeros_like(x) xpp = np.zeros_like(x) for k in range(1, H+1): xp += -k * coeff[k] * sin_basis[k] + k * coeff[H+k] * cos_basis[k] xpp += -k*k * coeff[k] * cos_basis[k] - k*k * coeff[H+k] * sin_basis[k] r = omega*omega*xpp + 2.0*zeta*omega*xp + x + beta*x**3 - f*np.cos(phi) res = np.zeros(2*H+1) res[0] = np.mean(r) for k in range(1, H+1): res[k] = 2.0*np.mean(r * np.cos(k*phi)) res[H+k] = 2.0*np.mean(r * np.sin(k*phi)) return res def solve_hb(omega, zeta=0.02, beta=1.0, f=1.0, H=5, N=64, guess=None): phi, cos_basis, sin_basis = build_basis(H, N) if guess is None: lin = f / np.sqrt((1.0 - omega**2)**2 + (2.0*zeta*omega)**2) guess = np.zeros(2*H+1) guess[1] = lin # 用线性共振幅值做骨架初值 sol, info, ier, msg = fsolve( lambda c: hb_residual(c, omega, zeta, beta, f, cos_basis, sin_basis, H), guess, full_output=True, xtol=1e-10) if ier != 1: raise RuntimeError(f"omega={omega} 未收敛: {msg}") return sol

代码逻辑说明：build_basis 在无量纲相位 φ∈[0,2π) 上采样 N 个点，cos_basis 和 sin_basis 的第 k 个元素分别对应 cos(kφ) 与 sin(kφ)，这样后面可以用同一个索引访问。x_from_coeff 把傅里叶系数还原成时域采样点，注意 a0 项用的是系数原值而不是 a0/2，因为投影里均值项用的就是 a0。hb_residual 先做频域解析求导，再做时域非线性计算，最后投影回每个谐波，这种“频域求导、时域算非线性、再投影”的做法是数值谐波平衡法的标准结构。

参数说明：H=5 对多数三次非线性问题足够，N=64 对应 2H+1=11 个未知数，采样点远超最低要求，不会出现混叠；xtol=1e-10 让 fsolve 收敛到很高精度，但代价是强非线性时可能报错，实际遇到报错可以先降到 1e-8。初值用线性共振幅值 guess[1]=lin 是关键：如果猜全零，三次非线性项在零点附近导数很小，牛顿法容易原地踏步。

3.2 扫频初值延续：让曲线翻过共振峰

固定频率能跑通后，下一步就是扫频。最简单也最有效的扫频方式是“延续法”：把所有谐波解按频率从小到大串起来，每算完一个点，就用它当下一频率点的初值。这比每个频率都重新猜初值可靠得多。

def sweep_continuation(omegas, zeta=0.02, beta=1.0, f=1.0, H=5, N=64): sols = [] guess = None for omega in omegas: try: sol = solve_hb(omega, zeta, beta, f, H, N, guess) sols.append(sol) guess = sol # 沿扫频方向传递初值 except RuntimeError: sols.append(None) # 该点真没找到，先留空位 return sols omegas = np.linspace(0.5, 2.0, 301) sols = sweep_continuation(omegas) amp = [np.sqrt(s[1]**2 + s[H+1]**2) for s in sols if s is not None]

这段代码把频响曲线一次性画出来时，会发现在共振峰附近曲线出现一个竖直的断裂带。这不是 bug，而是数学模型本身有多解：反向扫频（从 2.0 到 0.5）得到的是另一条分支。把两条扫频曲线叠在一起，中间包起来的区域就是多解区。

延续初值解决了“初值怎么给”的问题，但代价是把前一个点的误差也带了过来。工程上我一般会做升频和降频两次扫描，若两条曲线在多解区之外的区域重合良好，就说明延续链没有累积出问题。注意 fsolve 失败时，不能简单地把 guess 置空重新用线性初值。更好的做法是保留上一次成功的 guess，并在日志里记录这次 omega，方便排查。

3.3 伪弧长延拓：扫频断了才需要上的重型工具

延续法到强非线性时会遇到真正的麻烦：共振峰弯曲过度，频响曲线在极值点处掉头，频率 ω 不再是曲线坐标的单值函数。此时升频中断，降频也断，中间那段 S 形分支完全拿不到。伪弧长延拓就是为这种时刻准备的。

它的思路是把 ω 也当成未知数，给整个解向量加一个额外的约束方程“下一步必须在切向量方向的超平面上”。预测步写成：

# y = [a0, a1..aH, b1..bH, omega] # 预测: y_pred = y_current + ds * v # v 是雅可比矩阵零空间里的切向量，和当前解曲线的走向一致 # 校正: 同时求解残差方程 r(y_next)=0 与弧长约束 # (y_next - y_current)^T * v - ds = 0

实现上需要先对雅可比矩阵做一次 QR 分解拿到 v，再让 fsolve 同时处理扩维后的方程组。完整代码比这里的扫频循环长得多，我第一次实现时花了整整一下午调试切向量符号。所以建议路径是这样：先跑通 3.1 和 3.2 的延续法，等确认模型本身没问题、确实有 S 形分支需求时，再上伪弧长延拓。伪弧长延拓的步长 ds 一般取频率网格间距的 0.5～1 倍，太大会让校正阶段收敛困难，太小会浪费大量计算。

4. 频响曲线和 Floquet 乘子：怎么判断求出来的周期解真的成立

谐波平衡法解出一组系数并画出漂亮的频响曲线，只完成了一半工作。剩下的一半是回答一个问题：这组代数解，实际系统里到底待得住吗？

4.1 残差为零不代表解成立：Floquet 乘子是什么

残差为零只能说明这个周期解在数学上满足截断后的方程。真实系统里，噪声、数值误差、外部扰动随时会给解加一个小扰动。如果扰动在几个周期内被系统吞掉，这个周期解在实验里能复现；如果扰动被逐周期放大，它哪怕数学上存在，实操中也只会看到系统逃向别的响应形态。判断依据在非线性振动里叫 Floquet 乘子。

做法很直接：在周期解 x0(t) 上叠一个微小扰动量 δ(t)，把运动方程线性化。对 Duffing 方程，变分方程是：

δ'' + 2ζδ' + [1 + 3β x0²(t)] δ = 0

这是一个线性常微分方程，但括号里的刚度项随 x0² 周期性变化。取一个周期起点，给 δ 和 δ' 各设一个线性无关的初值，积分一个周期，就能得到一个 2×2 的单值矩阵 M。它的两个特征值就是 Floquet 乘子。每个乘子的模长把“扰动放大或缩小的年度结算单”写得很清楚：两个乘子模长都小于 1，扰动在周期内被压缩，解能被实验观测；只要有一个模长大于 1，扰动会被持续放大，这个解在物理过程中留不住。乘子模长越接近 1，系统越容易对参数变化表现出剧烈的响应切换，也常常对应频响曲线上的分岔点。

4.2 单值矩阵怎么算：两个初始扰动方向的周期积分

计算 M 不需要额外实现新的求解器，用 scipy.integrate.solve_ivp 就可以。下面的函数从前面的谐波解出发，直接返回两个 Floquet 乘子：

from scipy.integrate import solve_ivp def floquet_multipliers(coeff, omega, zeta, beta, H, N=512): phi = np.linspace(0.0, 2.0*np.pi, N, endpoint=True) cos_basis = [np.ones(N)] + [np.cos(k*phi) for k in range(1, H+1)] sin_basis = [np.zeros(N)] + [np.sin(k*phi) for k in range(1, H+1)] x_base = coeff[0]*np.ones(N) for k in range(1, H+1): x_base = x_base + coeff[k]*cos_basis[k] + coeff[H+k]*sin_basis[k] def rhs(tau, z): x0 = np.interp(tau, phi, x_base) kcoef = 1.0 + 3.0*beta*x0**2 z1, z1p, z2, z2p = z return [z1p, -2.0*zeta/omega*z1p - kcoef/omega**2*z1, z2p, -2.0*zeta/omega*z2p - kcoef/omega**2*z2] sol = solve_ivp(rhs, [0.0, 2.0*np.pi], [1.0, 0.0, 0.0, 1.0], t_eval=[2.0*np.pi], rtol=1e-9, atol=1e-12) z = sol.y[:, -1] M = np.array([[z[0], z[2]], [z[1], z[3]]]) return np.linalg.eigvals(M)

代码说明：初始条件 [1,0,0,1] 实际上是两个独立扰动方向的列向量，第一个方向给 δ=1、δ'=0，第二个方向给 δ=0、δ'=1，最后从 z 的四个分量里分别取出第一列和第二列的终点，拼成单值矩阵 M。注意微分方程里的阻尼项和刚度项都除以 omega²，因为相位 φ 和物理时间 t 差一个 ω 的缩放，这一处最容易抄错。

运行这个函数时，如果发现乘子模长恰好落在 1 附近，建议把 N 提高再算一次，并检查谐波解本身有没有收敛到足够的残差精度。Floquet 乘子对周期解的小误差非常挑剔，谐波解残差只有 1e-4 量级时，乘子误差可能大到影响判断。所以工程上真正的顺序是：先确认谐波解残差达标，再谈 Floquet 乘子。

4.3 三件套验收：残差、乘子、时域对照一起上

只靠一个指标定结论容易翻车。我给初做谐波平衡法的人定了三条验收线：第一，残差 RMS 至少要小于激励幅值的 1e-6；第二，Floquet 乘子模长要明确小于 1，而不是靠“差不多”。第三，挑共振峰两侧和跳跃区附近的点，把谐波重构的 x(t) 与时域四阶 Runge-Kutta 积分结果叠在一起看图，两个波形在 2～3 个周期里肉眼重合，才敢把这条频响曲线拿去做报告。

时域对照的具体做法是：用谐波解在周期末的状态 [x(0), x'(0)] 作为初值，用 solve_ivp 积分 30 个周期，然后取最后 2 个周期的波形和 x(t) 比较。这样做的好处是，即便谐波解本身乘子模长略大于 1，只要差异不大，短时间内波形不会漂移太远，依然能确认“这个解在局部是轨道意义上成立的”。如果 RK 积分很快就跑离谐波解，别急着怪 RK——先检查 Floquet 乘子，大概率乘子已经越过 1，谐波解的物理性本身有问题。

5. 谐波平衡法实战避坑：初值、谐波数、容差与模态选择的四个翻车现场

代码能在自己机器上跑通后，真正的工程问题才开始。这些年我见过最多、自己也踩过最多的是下面四个，按“现象、原因、解决”写清楚。

5.1 扫频曲线在共振峰附近断裂：初值链断了

现象：升频扫描从 0.5 开始一切顺利，到 0.92 附近突然 fsolve 报错，曲线断开；但降频扫描能顺利穿过同一区域。原因：这是经典的多解区入口。在跳跃频率附近，共振峰主分支弯回来，升频方向已经没有“连得上的解”，每个频率点用上一个解做初值也会掉进不合法的邻域。解决：先别上伪弧长延拓，升频和降频各扫一次，把两条分支叠在一起，确认断裂位置是否存在重叠解区；如果确实需要中间那条 S 形分支，才动用伪弧长延拓。

5.2 fsolve 明明给了初值却说迭代不前进：非线性项频域表达算错了

现象：初值已经在正确解附近，残差在 1e-2 量级，继续迭代却卡住，报 “The iteration is not making good progress”。原因：用解析谐波平衡时漏了非线性项的交叉频率组合。例如把 x³ 只展开成 3 倍频分量，忘了它还有 cos(ωt) 分量的贡献，导致残差里基频方程少了一项，雅可比矩阵结构不对。解决：放弃手写非线性频域表达，改用第 3 章的“时域算非线性、频域做投影”方案。这也是很多 NLvibration 包宁可少写几个解析函数也要保留时域投影的原因。

5.3 谐波数从 3 加到 7 反而出现锯齿波形：混叠在作怪

现象：H=3 时频响曲线光滑，H=7 时解出波形带锯齿，残差表面更小，但和时域积分对不上。原因：N 没跟着 H 增大。如果 N=32、H=7，未知数有 15 个，采样点只有 32，7 倍频分量几乎没有有效样本支撑，高频能量折叠回低频谱线形成混叠，方程反而被“假谐波”污染。解决：把 N 和 H 绑定，N 至少 4(2H+1)，工程上直接取 N=128、H 设到 5～7。改完再扫一遍，锯齿消失。

5.4 HBM 解和时域积分对不上：先查这四处

现象：谐波平衡法算出的共振峰位置比 RK 时域积分高了 5%，相位也对不上。原因通常不是方法本身，而是四个工程细节：一是 ω 的定义没对齐，代码里用的是圆频率而激励参数给的是 Hz，差 2π；二是无量纲化参数搞混，质量归一化时把刚度系数写错；三是时域积分初瞬态没耗尽，拿还没进入周期形态的波形和周期解比；四是响应里存在次谐波（比如 1/2 次谐波），而 H 只取了整数倍频。解决顺序：先打印 omega 的单位，再检查无量纲方程系数，然后把 RK 积分周期加长到 50～100 个周期，最后如果前三点都没问题，把 H 的范围放宽到包含半次谐波试试。注意最后一条意味着谐波平衡法的基频不再等于激励频率，而是取二者的最小公倍数周期。

6. 谐波平衡法的参数配方：一组让我少走弯路的默认设置

如果从零开始搭一套谐波平衡代码，我建议一组保守默认值。表格里每一行都经历过翻车教训，不是拍脑袋来的：

用这套配置处理硬弹簧 Duffing 时，升频和降频两条扫频曲线只会在跳跃区分离，其他区域重合得很好。这个现象本身就是最好的自检。我自己第一次做磁悬浮轴承转子时，图省事用了 N=32、H=3，结果边界附近怎么都对不上，花了两个星期才发现是采样点数太少导致高次谐波折叠，从那以后就一直保留这个配方，宁可慢一点，也不让混叠混进来。

最后说一句习惯：每个周期解我都会顺手存一份 H=1 的骨架解，再做高次谐波的修正。两套解差得多，说明非线性效应强，需要加倍检查；两套解几乎一样，说明这个问题本质上接近线性，可以更放心。这套习惯救过我好几次，希望帮到你。

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