尧图网页设计 
简介面向非线性振动研究与教学场景的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₀ Σ_{k1}^{H} [a_k cos(kωt) b_k sin(kωt)]把这一形式代回 Duffing 方程x³ 仍然是一个周期函数可以被展开成相同的谐波形式。接下来让 cos(kωt) 和 sin(kωt) 的系数在等式两侧分别相等就得到 2H或 2H1如果保留常数项个代数方程。原来的二阶常微分方程被转换成一个需要联立求解的非线性代数方程组。NLvibration 这类谐波平衡程序包本质上只做三件事把每个谐波系数排列成未知向量把非线性项映射到频域对非线性代数方程组做牛顿迭代。线性项部分甚至不需要迭代因为它对每个谐波独立可以用一个块对角矩阵直接表达。2.2 NLvibration 的数据流与模块划分常见做法是把程序拆成三个职责清晰的模块。第一个是系统描述模块负责承载质量、阻尼、刚度、非线性系数和激励参数第二个是谐波平衡核心模块负责时域与频域交替变换AFT、残差组装和雅可比矩阵第三个是求解器模块负责牛顿迭代、扫频延续和稳定性判断。解压 NLvibration 后源码文件通常也是按这三个职责组织的后处理脚本单独放。模块职责关键函数输入输出系统描述set_system_parameters物理参数、激励幅值、扫频范围参数对象谐波平衡核心nonlinear_force_freq, residual谐波系数向量 X、频率 omega残差向量 R、雅可比矩阵 J求解器newton_solve, continuation初值 X0、频率序列各频率下的谐波系数后处理amplitude_phase, floquet谐波系数、系统参数幅值、相位、稳定性标记求解器输出的是每个频率点上的傅里叶系数后处理再把这些系数换算成一阶幅值、相位和总谐波失真。稳定性标记由 Floquet 乘子计算得到稳定段和不稳定段在幅频曲线上用不同符号区分。2.3 最小可跑通示例AFT 频时变换的代码非线性项 x³ 在频域里没有简单的闭合形式工程上最常见的做法是 AFTAlternating 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, endpointFalse) / 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] 上。参数说明H5 对应 10 个未知数是弱非线性系统的保守选择。N512 是 2 的幂FFT 效率最高且采样密度足够即使 x³ 产生最高 15 阶谐波也不会混叠。牛顿迭代里的 eps1e-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₀。这时未知量要扩成 2H1 维FFT 结果里的 DC 分量取 Fk[0].real不乘 2重建时域信号时额外加 a₀。很多程序包默认不含 a₀遇到不对称恢复力时结果会明显偏差这是第一个要检查的扩展点。3.2 谐波截断阶数与采样点数的关系H 的选取决定了周期解的逼近精度。弱非线性α x³ 的贡献小于线性项 10%时 H3 到 H5 足够强非线性下共振峰严重弯曲响应波形变成近似方波需要 H8 到 H15。工程上的判断方法是收敛后检查最高阶谐波幅值与主谐波幅值之比如果比值大于 1e-3说明截断过早需要加大 H 重新计算。N 的选取同样关键。AFT 方法里时域采样覆盖一个周期 T2π/ωFFT 的频率分辨率正好是 ω第 k 条谱线对应 kω。但 x³ 的非线性会把能量推到 3H 阶谐波如果 N 不够大这些高频分量会混叠回低频谱线污染结果。经验取值是 N 取 2 的高次幂且 N ≥ 8H比如 H5 时 N 至少 64但为了保险我一般直接用 N512。3.3 牛顿迭代参数初值、收敛判据与阻尼策略牛顿迭代是 HBM 求解器的主力。每步迭代中雅可比矩阵用中心差分近似代价是每次迭代要额外求 4H 次残差。对 10 个未知数来说不算大但对多自由度系统或高阶截断解析雅可比能省下大量计算时间。def numerical_jacobian(residual_fn, X, omega, eps_rel1e-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, tol1e-10, max_iter50): 带阻尼策略的牛顿迭代残差增大时自动减半步长 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 扫频参数步长、延续策略与双向扫描求幅频曲线时逐点独立求解每个频率太低效常见做法是延续法从低频端出发把上一个频率收敛到的解作为下一个频率的初值。因为相邻频率的系统参数变化很小这个初值通常落在牛顿迭代的收敛域内。扫频前先确定频率序列和初值参数推荐值说明omega 序列np.linspace(0.4, 2.6, 200)共振区附近需要更密可分段加密首个初值线性解 F0 / (k - mω²)低频端远离共振线性解足够好延续方式X_prev X下一个频率从当前解起步幅值提取a₁² b₁² 开根主谐波幅值用于画幅频曲线扫频方向有个关键坑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, tol1e-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, tol1e-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, endpointFalse) _, 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但输出的各阶幅值全部是零。原因x0 是 Duffing 方程的平衡点全零初值让牛顿法直接掉进这个平凡解。解决把初值改为线性系统解也就是 F0 / (k - mω²)或者从相邻频率的已收敛解延续过来。含有常数项的系统还需要给 a₀ 一个微小的非零初值比如 1e-6避免常数项通道也陷入零解。5.2 谐波截断阶数不足导致幅值偏差现象用小 H 计算时共振峰幅值比数值积分低超过 10%增加 H 后幅值明显变化。原因强非线性让高次谐波参与能量交换五阶甚至七阶谐波幅值不可忽略截断误差直接进入频域力平衡。解决先算一版 H3再算一版 H7比较主谐波幅值的变化量变化量小于 1% 才算收敛。如果响应包含次谐波共振比如响应周期是激励周期的 2 倍未知量必须按基频的一半展开单纯增大 H 是没用的。5.3 单向扫频导致跳跃路径丢失现象升频扫描在共振峰后幅值突然往下掉曲线中间缺一段降低频率步长也补不回来。原因多解区间的各分支被极限点隔开牛顿迭代的收敛域有限只能跟住起始分支无法跨过极限点。解决做双向扫频升频和降频的曲线合起来看如果还需要极限点之间的不稳定支要用伪弧长延拓不能再依赖普通延续法。5.4 FFT 采样点数不足导致谐波混叠现象N64、H8 时计算结果混乱增大 H 反而更乱。原因x³ 至少产生 3H 阶谐波64 个采样点不够容纳这些高频分量它们混叠回低频段污染前几阶谐波的系数。解决N 取 2 的高次幂且不小于 8H工程保险值直接取 N512。强非线性下先检查 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 稳定性判断加上才看清全貌。现在每次跑新系统都会先扫一版低分辨率全频率范围确认没有异常多解区后再决定哪一段需要加密或上伪弧长。希望帮到你。本文还有配套的精品资源点击获取 谐波平衡法剖析：从Duffing方程到幅频曲线与稳定性判断
简介面向非线性振动研究与教学场景的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₀ Σ_{k1}^{H} [a_k cos(kωt) b_k sin(kωt)]把这一形式代回 Duffing 方程x³ 仍然是一个周期函数可以被展开成相同的谐波形式。接下来让 cos(kωt) 和 sin(kωt) 的系数在等式两侧分别相等就得到 2H或 2H1如果保留常数项个代数方程。原来的二阶常微分方程被转换成一个需要联立求解的非线性代数方程组。NLvibration 这类谐波平衡程序包本质上只做三件事把每个谐波系数排列成未知向量把非线性项映射到频域对非线性代数方程组做牛顿迭代。线性项部分甚至不需要迭代因为它对每个谐波独立可以用一个块对角矩阵直接表达。2.2 NLvibration 的数据流与模块划分常见做法是把程序拆成三个职责清晰的模块。第一个是系统描述模块负责承载质量、阻尼、刚度、非线性系数和激励参数第二个是谐波平衡核心模块负责时域与频域交替变换AFT、残差组装和雅可比矩阵第三个是求解器模块负责牛顿迭代、扫频延续和稳定性判断。解压 NLvibration 后源码文件通常也是按这三个职责组织的后处理脚本单独放。模块职责关键函数输入输出系统描述set_system_parameters物理参数、激励幅值、扫频范围参数对象谐波平衡核心nonlinear_force_freq, residual谐波系数向量 X、频率 omega残差向量 R、雅可比矩阵 J求解器newton_solve, continuation初值 X0、频率序列各频率下的谐波系数后处理amplitude_phase, floquet谐波系数、系统参数幅值、相位、稳定性标记求解器输出的是每个频率点上的傅里叶系数后处理再把这些系数换算成一阶幅值、相位和总谐波失真。稳定性标记由 Floquet 乘子计算得到稳定段和不稳定段在幅频曲线上用不同符号区分。2.3 最小可跑通示例AFT 频时变换的代码非线性项 x³ 在频域里没有简单的闭合形式工程上最常见的做法是 AFTAlternating 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, endpointFalse) / 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] 上。参数说明H5 对应 10 个未知数是弱非线性系统的保守选择。N512 是 2 的幂FFT 效率最高且采样密度足够即使 x³ 产生最高 15 阶谐波也不会混叠。牛顿迭代里的 eps1e-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₀。这时未知量要扩成 2H1 维FFT 结果里的 DC 分量取 Fk[0].real不乘 2重建时域信号时额外加 a₀。很多程序包默认不含 a₀遇到不对称恢复力时结果会明显偏差这是第一个要检查的扩展点。3.2 谐波截断阶数与采样点数的关系H 的选取决定了周期解的逼近精度。弱非线性α x³ 的贡献小于线性项 10%时 H3 到 H5 足够强非线性下共振峰严重弯曲响应波形变成近似方波需要 H8 到 H15。工程上的判断方法是收敛后检查最高阶谐波幅值与主谐波幅值之比如果比值大于 1e-3说明截断过早需要加大 H 重新计算。N 的选取同样关键。AFT 方法里时域采样覆盖一个周期 T2π/ωFFT 的频率分辨率正好是 ω第 k 条谱线对应 kω。但 x³ 的非线性会把能量推到 3H 阶谐波如果 N 不够大这些高频分量会混叠回低频谱线污染结果。经验取值是 N 取 2 的高次幂且 N ≥ 8H比如 H5 时 N 至少 64但为了保险我一般直接用 N512。3.3 牛顿迭代参数初值、收敛判据与阻尼策略牛顿迭代是 HBM 求解器的主力。每步迭代中雅可比矩阵用中心差分近似代价是每次迭代要额外求 4H 次残差。对 10 个未知数来说不算大但对多自由度系统或高阶截断解析雅可比能省下大量计算时间。def numerical_jacobian(residual_fn, X, omega, eps_rel1e-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, tol1e-10, max_iter50): 带阻尼策略的牛顿迭代残差增大时自动减半步长 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 扫频参数步长、延续策略与双向扫描求幅频曲线时逐点独立求解每个频率太低效常见做法是延续法从低频端出发把上一个频率收敛到的解作为下一个频率的初值。因为相邻频率的系统参数变化很小这个初值通常落在牛顿迭代的收敛域内。扫频前先确定频率序列和初值参数推荐值说明omega 序列np.linspace(0.4, 2.6, 200)共振区附近需要更密可分段加密首个初值线性解 F0 / (k - mω²)低频端远离共振线性解足够好延续方式X_prev X下一个频率从当前解起步幅值提取a₁² b₁² 开根主谐波幅值用于画幅频曲线扫频方向有个关键坑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, tol1e-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, tol1e-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, endpointFalse) _, 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但输出的各阶幅值全部是零。原因x0 是 Duffing 方程的平衡点全零初值让牛顿法直接掉进这个平凡解。解决把初值改为线性系统解也就是 F0 / (k - mω²)或者从相邻频率的已收敛解延续过来。含有常数项的系统还需要给 a₀ 一个微小的非零初值比如 1e-6避免常数项通道也陷入零解。5.2 谐波截断阶数不足导致幅值偏差现象用小 H 计算时共振峰幅值比数值积分低超过 10%增加 H 后幅值明显变化。原因强非线性让高次谐波参与能量交换五阶甚至七阶谐波幅值不可忽略截断误差直接进入频域力平衡。解决先算一版 H3再算一版 H7比较主谐波幅值的变化量变化量小于 1% 才算收敛。如果响应包含次谐波共振比如响应周期是激励周期的 2 倍未知量必须按基频的一半展开单纯增大 H 是没用的。5.3 单向扫频导致跳跃路径丢失现象升频扫描在共振峰后幅值突然往下掉曲线中间缺一段降低频率步长也补不回来。原因多解区间的各分支被极限点隔开牛顿迭代的收敛域有限只能跟住起始分支无法跨过极限点。解决做双向扫频升频和降频的曲线合起来看如果还需要极限点之间的不稳定支要用伪弧长延拓不能再依赖普通延续法。5.4 FFT 采样点数不足导致谐波混叠现象N64、H8 时计算结果混乱增大 H 反而更乱。原因x³ 至少产生 3H 阶谐波64 个采样点不够容纳这些高频分量它们混叠回低频段污染前几阶谐波的系数。解决N 取 2 的高次幂且不小于 8H工程保险值直接取 N512。强非线性下先检查 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 稳定性判断加上才看清全貌。现在每次跑新系统都会先扫一版低分辨率全频率范围确认没有异常多解区后再决定哪一段需要加密或上伪弧长。希望帮到你。本文还有配套的精品资源点击获取