企业官网建设流程全解析

简介：面向非线性振动研究与教学场景的MATLAB代码包，针对多稳态、混沌、分岔等非线性特征，围绕谐波平衡法求解周期解给出完整实现思路。资源共14个文件，全部为m脚本，压缩包仅6KB，包含主程序、力函数与响应求解、矩阵线性化等独立模块，并带多个测试函数，便于按步骤调用、调试和扩展。已有619人浏览学习，可用于机械、航空航天等领域的非线性振动分析与课程实验。代码从非线性方程定义、谐波展开、线性化到各阶振幅相位求解和叠加，将谐波平衡法计算流程模块化拆分；使用者可修改系统参数、调整谐波阶数，快速获得近似周期解并观察幅频关系与动态响应，支持对比不同初始状态下系统是否收敛至不同周期解，对科研验证、教学演示或理解复杂振动现象都很有帮助。

1. 谐波平衡法：非线性振动周期解分析的正确打开方式

做非线性振动稳态分析的人应该都有过这种经历：一个含立方刚度的减振器模型，用数值积分从零初值开始跑，瞬态过程拖了几百个激励周期还没有进稳态，算一条幅频曲线要在每个频率点重复这个过程，半天时间就没了。NLvibration 这个程序包正是用来改变这件事的——它用谐波平衡法（Harmonic Balance Method, HBM）把周期解问题从时域积分变成频域代数方程组，直接假设响应是有限阶傅里叶级数，然后在频域里让每个谐波的系数满足力平衡。求解结果是一组谐波幅值系数，展开即可得到整个稳态周期响应，还能进一步分析幅频曲线、跳跃现象和稳定性，适合转子动力学、减振器设计与微谐振器建模的工程师和研究生。

2. NLvibration 程序包解剖：从 Duffing 方程到最小可跑通代码

2.1 从 Duffing 方程看谐波平衡法的基本流程

非线性振动里最常用的基准模型是 Duffing 方程：

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

其中 m、c、k 分别是质量、阻尼和线性刚度，α 是立方非线性系数。这个方程能同时体现硬弹簧软弹簧效应、共振峰弯曲、多解与跳跃这些非线性现象，几乎所有谐波平衡法的入门代码都会先拿它做验证。

谐波平衡法的核心假设是：稳态周期解可以写成有限项傅里叶级数：

x(t) = a₀ + Σ_{k=1}^{H} [a_k cos(kωt) + b_k sin(kωt)]

把这一形式代回 Duffing 方程，x³ 仍然是一个周期函数，可以被展开成相同的谐波形式。接下来让 cos(kωt) 和 sin(kωt) 的系数在等式两侧分别相等，就得到 2H（或 2H+1，如果保留常数项）个代数方程。原来的二阶常微分方程，被转换成一个需要联立求解的非线性代数方程组。

NLvibration 这类谐波平衡程序包本质上只做三件事：把每个谐波系数排列成未知向量，把非线性项映射到频域，对非线性代数方程组做牛顿迭代。线性项部分甚至不需要迭代，因为它对每个谐波独立，可以用一个块对角矩阵直接表达。

2.2 NLvibration 的数据流与模块划分

常见做法是把程序拆成三个职责清晰的模块。第一个是系统描述模块，负责承载质量、阻尼、刚度、非线性系数和激励参数；第二个是谐波平衡核心模块，负责时域与频域交替变换（AFT）、残差组装和雅可比矩阵；第三个是求解器模块，负责牛顿迭代、扫频延续和稳定性判断。解压 NLvibration 后，源码文件通常也是按这三个职责组织的，后处理脚本单独放。

求解器输出的是每个频率点上的傅里叶系数，后处理再把这些系数换算成一阶幅值、相位和总谐波失真。稳定性标记由 Floquet 乘子计算得到，稳定段和不稳定段在幅频曲线上用不同符号区分。

2.3 最小可跑通示例：AFT 频时变换的代码

非线性项 x³ 在频域里没有简单的闭合形式，工程上最常见的做法是 AFT（Alternating Frequency-Time）方法：先由谐波系数在时域重建 x(t)，在时域计算非线性力，再用 FFT 把结果映射回频域。下面是完整可运行的单点求解脚本，直接复制就能看到一阶幅值的结果。

import numpy as np # Duffing 系统参数 M = 1.0 C = 0.02 K = 1.0 ALPHA = 0.5 F0 = 0.3 OMEGA = 1.2 # 激励频率 H = 5 # 最高谐波次数 N = 512 # 时域采样点数 # X = [a1, b1, a2, b2, ..., aH, bH] def time_signal(X, omega): t = np.linspace(0, 2 * np.pi, N, endpoint=False) / omega x = np.zeros(N) for k in range(1, len(X) // 2 + 1): a, b = X[2 * k - 2], X[2 * k - 1] x += a * np.cos(k * omega * t) + b * np.sin(k * omega * t) return t, x def nonlinear_freq(X, alpha): _, x = time_signal(X, OMEGA) Fnl = alpha * x ** 3 Fk = np.fft.fft(Fnl) / N Y = np.zeros_like(X) for k in range(1, len(X) // 2 + 1): Y[2 * k - 2] = 2 * Fk[k].real # 余弦分量 Y[2 * k - 1] = -2 * Fk[k].imag # 正弦分量 return Y def residual(X): R = np.zeros_like(X) for k in range(1, H + 1): w = k * OMEGA D = np.array([[K - M * w ** 2, C * w], [-C * w, K - M * w ** 2]]) R[2 * k - 2:2 * k] = D @ X[2 * k - 2:2 * k] R += nonlinear_freq(X, ALPHA) R[0] -= F0 # 外激励只作用在一阶余弦方程上 return R # 以线性系统解作为初值，避免收敛到零解 X0 = np.zeros(2 * H) X0[0] = F0 / (K - M * OMEGA ** 2) X = X0.copy() for it in range(30): R = residual(X) if np.linalg.norm(R) < 1e-11: break J = np.zeros((2 * H, 2 * H)) eps = 1e-7 for j in range(2 * H): Xp, Xm = X.copy(), X.copy() Xp[j] += eps Xm[j] -= eps J[:, j] = (residual(Xp) - residual(Xm)) / (2 * eps) X += np.linalg.solve(J, -R) print("一阶幅值:", round(np.hypot(X[0], X[1]), 8)) print("各阶幅值:", [round(np.hypot(X[2 * k], X[2 * k + 1]), 8) for k in range(H)])

逻辑说明：残差由三部分构成，线性动力学刚度矩阵贡献、非线性力贡献和外激励贡献。线性部分逐谐波组装，D 矩阵把第 k 阶余弦系数和正弦系数耦合在一起；非线性部分通过 AFT 得到；外激励 F cos(ωt) 只出现在一阶余弦方程里，所以只减在 R[0] 上。

参数说明：H=5 对应 10 个未知数，是弱非线性系统的保守选择。N=512 是 2 的幂，FFT 效率最高且采样密度足够，即使 x³ 产生最高 15 阶谐波也不会混叠。牛顿迭代里的 eps=1e-7 是固定绝对差分步长，对量级为 1 的变量合适；如果响应幅值很大或很小，这个值需要改成相对步长，具体做法放在第 3 章。

3. 周期解核心方程与截断阶数：把谐波平衡变成可收敛的迭代

3.1 谐波平衡方程组的完整组装

第 2 章的代码能跑通，但要真正调好参数，需要理解线性动力学刚度矩阵的来历。对第 k 阶谐波，设 x_k = a_k cos(kωt) + b_k sin(kωt)，代入 m x''、c x'、k x 后分别提取 cos 和 sin 的系数：

cos 方程：(k - m(kω)²) a_k + (c kω) b_k sin 方程：(-c kω) a_k + (k - m(kω)²) b_k

所以单个谐波对应的 2×2 矩阵是：

D_k = [[k - m(kω)², c kω], [-c kω, k - m(kω)²]]

注意副对角线的符号来自速度项的导数，c x' 在投影时会交叉耦合余弦和正弦系数，写错符号会导致共振频率偏移。整个线性部分是一个块对角矩阵，每个块就是 D_k，它不随迭代变化，可以在扫频前预先算好。R = D X + F_nl(X) - F_ext 就是完整的谐波平衡残差方程。

如果系统包含平方非线性项 β x² 或静态预载，响应里会出现非零常数项 a₀。这时未知量要扩成 2H+1 维，FFT 结果里的 DC 分量取 Fk[0].real，不乘 2，重建时域信号时额外加 a₀。很多程序包默认不含 a₀，遇到不对称恢复力时结果会明显偏差，这是第一个要检查的扩展点。

3.2 谐波截断阶数与采样点数的关系

H 的选取决定了周期解的逼近精度。弱非线性（α x³ 的贡献小于线性项 10%）时 H=3 到 H=5 足够；强非线性下共振峰严重弯曲，响应波形变成近似方波，需要 H=8 到 H=15。工程上的判断方法是：收敛后检查最高阶谐波幅值与主谐波幅值之比，如果比值大于 1e-3，说明截断过早，需要加大 H 重新计算。

N 的选取同样关键。AFT 方法里时域采样覆盖一个周期 T=2π/ω，FFT 的频率分辨率正好是 ω，第 k 条谱线对应 kω。但 x³ 的非线性会把能量推到 3H 阶谐波，如果 N 不够大，这些高频分量会混叠回低频谱线，污染结果。经验取值是 N 取 2 的高次幂且 N ≥ 8H，比如 H=5 时 N 至少 64，但为了保险我一般直接用 N=512。

3.3 牛顿迭代参数：初值、收敛判据与阻尼策略

牛顿迭代是 HBM 求解器的主力。每步迭代中，雅可比矩阵用中心差分近似，代价是每次迭代要额外求 4H 次残差。对 10 个未知数来说不算大，但对多自由度系统或高阶截断，解析雅可比能省下大量计算时间。

def numerical_jacobian(residual_fn, X, omega, eps_rel=1e-7): """中心差分数值雅可比，步长随变量量级自适应""" n = len(X) J = np.zeros((n, n)) for j in range(n): eps = eps_rel * max(1.0, abs(X[j])) Xp, Xm = X.copy(), X.copy() Xp[j] += eps Xm[j] -= eps J[:, j] = (residual_fn(Xp, omega) - residual_fn(Xm, omega)) / (2 * eps) return J def newton_solve(X0, omega, residual_fn, tol=1e-10, max_iter=50): """带阻尼策略的牛顿迭代，残差增大时自动减半步长""" X = X0.copy() for it in range(max_iter): R = residual_fn(X, omega) rn = np.linalg.norm(R) if rn < tol: return X, it, True J = numerical_jacobian(residual_fn, X, omega) dx = np.linalg.solve(J, -R) X_next = X + dx # 阻尼：如果新残差不降反升，说明步长过大，折半重试 while np.linalg.norm(residual_fn(X_next, omega)) > rn and np.linalg.norm(dx) > 1e-12: dx *= 0.5 X_next = X + dx X = X_next return X, max_iter, False

参数说明：eps_rel 是相对差分步长，取 1e-6 到 1e-7 之间比较稳妥。太小会让有限差分受浮点舍入误差主导，太大则雅可比偏离真实导数。收敛判据用的是双条件：残差范数小于 tol，或者步长 dx 已经小到不再改变解，两者满足其一即可退出。阻尼策略解决的是初值偏差较大时牛顿迭代发散的常见问题，虽然会多算几次残差，但比直接抛异常要实用得多。

提示：如果迭代次数超过 max_iter 但残差还在下降，不要把 max_iter 无限加大，先检查初值和差分步长，多数情况是初值落入了错误分支而不是收敛慢。

4. 幅频曲线与稳定性验证：扫频参数、Floquet 乘子和结果核验

4.1 扫频参数：步长、延续策略与双向扫描

求幅频曲线时，逐点独立求解每个频率太低效，常见做法是延续法：从低频端出发，把上一个频率收敛到的解作为下一个频率的初值。因为相邻频率的系统参数变化很小，这个初值通常落在牛顿迭代的收敛域内。扫频前先确定频率序列和初值：

扫频方向有个关键坑：Duffing 系统在硬弹簧条件下幅频曲线向右弯曲，共振峰附近存在多解区间。从低频往高频扫，会得到上面一支；从高频往低频扫，会得到下面一支。只做一个方向，中间会有一段"跳变"，看起来像计算结果缺失。正确做法是双向扫频，两条曲线合在一起，多解区间自然露出。

# 升频扫描 omega_range = np.linspace(0.4, 2.6, 200) amp_up = np.zeros(len(omega_range)) X_prev = np.zeros(2 * H) X_prev[0] = F0 / (K - M * omega_range[0] ** 2) for i, w in enumerate(omega_range): X, _, ok = newton_solve(X_prev, w, residual, tol=1e-10) amp_up[i] = np.hypot(X[0], X[1]) X_prev = X # 降频扫描：从高频端反向走，初值同样用线性解 omega_rev = omega_range[::-1] amp_down = np.zeros(len(omega_rev)) X_prev[0] = F0 / (K - M * omega_rev[0] ** 2) for i, w in enumerate(omega_rev): X, _, ok = newton_solve(X_prev, w, residual, tol=1e-10) amp_down[i] = np.hypot(X[0], X[1]) X_prev = X

参数说明：升频和降频两个方向用各自的线性解起步，是为了绕过多解区间对初值的吸引。如果两个方向扫出的曲线在某一频率段不重合，那一段就是多解区，跳跃点就在边界上。要继续深挖多解分支内部的结构，就需要第 6 章的伪弧长延拓。

4.2 稳定性判断：Floquet 乘子的计算

HBM 给出的是周期解，但这个解在物理上不一定能实现。判断周期解是否稳定，需要对解施加一个小扰动，观察扰动在一个周期后的发展。把 Duffing 方程在周期解附近线性化，得到变分方程：

δx'' + (c/m) δx' + (k + 3α x(t)²)/m δx = 0

其中 x(t) 是 HBM 解重建的时域位移。状态转移矩阵 Φ(T) 的特征值就是 Floquet 乘子，乘子模最大值大于 1 时周期解不稳定。

def floquet_multipliers(X, omega): """计算 Floquet 乘子，返回最大模""" T = 2 * np.pi / omega t = np.linspace(0, T, 1000, endpoint=False) _, x = time_signal(X, omega) dt = t[1] - t[0] Phi = np.eye(2) for i in range(len(t) - 1): k_t = K + 3 * ALPHA * x[i] ** 2 A = np.array([[0, 1], [-k_t / M, -C / M]]) Phi = Phi @ (np.eye(2) + A * dt) # 一阶近似，步长足够小 mu = np.linalg.eigvals(Phi) return np.max(np.abs(mu))

逻辑说明：变分矩阵 A 在一个周期内随时间变化，所以用欧拉逐步逼近积分状态转移矩阵。1000 个采样点的步长约为 T/1000，对 Duffing 这类系统精度足够。乘子最大模大于 1.0 时标记为不稳定，实际计算时建议用 1.0001 作为阈值，避免数值误差引起误判。

4.3 与数值积分对比的验证流程

HBM 结果是近似解，必须和时域数值积分结果做交叉验证。验证时一个实用技巧是直接用 HBM 解重建的位移和速度作为数值积分的初始条件，这样瞬态很短，通常几十个周期就能进入稳态，省掉从零初值开始的漫长等待。

from scipy.integrate import solve_ivp def duffing_state(t, y): x, v = y return [v, (F0 * np.cos(omega * t) - C * v - K * x - ALPHA * x ** 3) / M] # 从 HBM 解重建初始位移和速度 T = 2 * np.pi / omega t_cycle, x_hbm = time_signal(X, omega) x0_val = x_hbm[0] v0_val = (x_hbm[1] - x_hbm[-1]) / (2 * T / N) # 中心差分近似初速度

数值积分跑 100 个周期，丢弃前 80 个周期的瞬态，对最后 20 个周期做 FFT，提取一阶幅值与 HBM 对比。一阶幅值相对误差小于 1e-3 视为通过。如果误差偏大，优先检查 H 和 N，然后检查激励频率是否落在极限点附近，极限点附近对初值极其敏感，数值积分和 HBM 都可能各自收敛到不同分支。

5. 谐波平衡法避坑指南：5 个高频问题和排查流程

5.1 牛顿迭代收敛到平凡零解

现象：迭代正常结束，残差范数小于 1e-12，但输出的各阶幅值全部是零。原因：x=0 是 Duffing 方程的平衡点，全零初值让牛顿法直接掉进这个平凡解。解决：把初值改为线性系统解，也就是 F0 / (k - mω²)，或者从相邻频率的已收敛解延续过来。含有常数项的系统还需要给 a₀ 一个微小的非零初值，比如 1e-6，避免常数项通道也陷入零解。

5.2 谐波截断阶数不足导致幅值偏差

现象：用小 H 计算时共振峰幅值比数值积分低超过 10%，增加 H 后幅值明显变化。原因：强非线性让高次谐波参与能量交换，五阶甚至七阶谐波幅值不可忽略，截断误差直接进入频域力平衡。解决：先算一版 H=3，再算一版 H=7，比较主谐波幅值的变化量；变化量小于 1% 才算收敛。如果响应包含次谐波共振，比如响应周期是激励周期的 2 倍，未知量必须按基频的一半展开，单纯增大 H 是没用的。

5.3 单向扫频导致跳跃路径丢失

现象：升频扫描在共振峰后幅值突然往下掉，曲线中间缺一段；降低频率步长也补不回来。原因：多解区间的各分支被极限点隔开，牛顿迭代的收敛域有限，只能跟住起始分支，无法跨过极限点。解决：做双向扫频，升频和降频的曲线合起来看；如果还需要极限点之间的不稳定支，要用伪弧长延拓，不能再依赖普通延续法。

5.4 FFT 采样点数不足导致谐波混叠

现象：N=64、H=8 时计算结果混乱，增大 H 反而更乱。原因：x³ 至少产生 3H 阶谐波，64 个采样点不够容纳这些高频分量，它们混叠回低频段，污染前几阶谐波的系数。解决：N 取 2 的高次幂且不小于 8H，工程保险值直接取 N=512。强非线性下先检查 N 够不够，再讨论 H 够不够，顺序不要反。

5.5 数值雅可比差分步长选错导致迭代停滞

现象：牛顿迭代残差降到 1e-6 左右就卡住，加大迭代次数也没有改善。原因：固定差分步长选得太小，比如 1e-12，残差的浮点舍入噪声在雅可比里占了主导。解决：改用相对步长 eps_rel × max(1.0, |X_j|)，取 1e-6 到 1e-7；如果问题规模大，直接推导解析雅可比，对 Duffing 这类立方非线性，解析式很容易写。

6. 从周期解到极限点追踪：伪弧长延拓与多谐波进阶

6.1 伪弧长延拓：穿过极限点的可靠路径

普通延续法把上一个频率的解作为初值，但极限点处雅可比矩阵奇异，牛顿迭代无法逾越。伪弧长延拓把频率 ω 也当作未知量，在 (X, ω) 扩展空间里沿着解曲线的弧长方向预测和校正。预测步用前两步解的差作为切向，校正步在扩展空间中求解加约束的方程组：

# 预测：沿前两步方向外推 dx_dir = (X_prev - X_prev2) / ds_prev dw_dir = (omega_prev - omega_prev2) / ds_prev X_guess = X_prev + ds * dx_dir omega_guess = omega_prev + ds * dw_dir # 校正：扩展残差方程，把切向约束加入未知量 # [R(X, omega); (X - X_guess) @ dx_dir + (omega - omega_guess) * dw_dir] = 0

扩展方程的好处是即使原雅可比奇异，扩展雅可比通常仍然是满秩的，可以平滑通过极限点。做幅频曲线多解区间分析时，伪弧长延拓是标准工具，普通延续只适合远离多解区的简单扫频。

6.2 谐波自适应策略与激励频率扩展

外部激励包含多个频率分量时，基频应取各激励频率的最大公约数，谐波索引按这个基频的整数倍排列。收敛后检查最高阶谐波幅值与最大幅值之比，超过 1e-3 就自动增加 H 并重算，形成简单的自适应循环。我早年只靠单向扫频，碰到跳跃现象时一度以为是程序写错了，后来把双向扫频和 Floquet 稳定性判断加上才看清全貌。现在每次跑新系统，都会先扫一版低分辨率全频率范围，确认没有异常多解区后，再决定哪一段需要加密或上伪弧长。希望帮到你。

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