谐波平衡法求解非线性振动周期解的完整流程与代码实现

简介：非线性振动中周期解的求解常采用谐波平衡法，一份配套的MATLAB代码包可以帮助研究者与学习者快速验证算法、观察动态响应。该代码包面向机械、航空航天、土木工程等专业方向，适合非线性振动课程、结构动力学分析以及相关科研预研。谐波平衡法将周期解近似展开为基频整数倍谐波的叠加，通过非线性项展开与线性化求各阶幅值相位；包内代码围绕这一思路编写，涵盖主程序、非线性力定义、激励函数、响应计算、矩阵线性化等模块，用户可直接运行主程序获取周期解，也可通过调整参数或替换非线性项适配不同系统。压缩包共14个文件，全部为.m脚本与函数，整体仅6KB，体积小、结构清晰，便于逐段阅读算法流程；各模块耦合度低，支持独立修改与调用。该资源已有619人学习下载，可作为课堂演示、课题研究及二次开发的实用起点。

1. 谐波平衡法是求nonlinear vibration周期解最不依赖初值的一条路

很多人在算非线性振动（nonlinear vibration）的周期解时，第一反应是用Runge-Kutta把方程从零时刻积分到“看起来稳定”。实际上你会遇到两个麻烦：瞬态过程太长，步长取不好会导致高频分量被数值耗散吃掉；而在强非线性下，系统可能存在多个周期解，时域积分只能给出一张“吸引子快照”，你根本不知道在同一个激励频率下旁边还藏着另一个解。谐波平衡法的思路完全相反：它先把周期解假设成一组Fourier级数，再把原微分方程投影到各个谐波项上，把求ODE初值的问题变成求解一组代数方程。NLvibration这类小工具的核心就是这个过程，它不依赖时域积分初值，也不怕解跳到别的分支上，只要代数残差能收敛，周期解就摆在明处。对做转子动力学、MEMS谐振器、电力电子振荡器的人，这一招几乎绕不过去。本文直接讲清楚谐波平衡法的推导、参数设置、强非线性延拓和最终验证，每一步都给能跑的代码。

2. 从杜芬方程开始，用谐波平衡法把微分方程折成代数方程

2.1 周期解为什么能用Fourier截断而不丢主特征

一个周期为 T 的稳态解，如果满足Dirichlet条件，总可以写成Fourier级数：

x(t)=a0+∑_{k=1}^{∞} (ak cos(kωt)+bk sin(kωt))

谐波平衡法的核心假设是：对于工程里常见的阻尼系统，高频谐波在能量上占比极小，保留到 N 阶就够。这里的 N 在NLvibration里通常叫“谐波数（Harmonics）”或“截断阶数”。N 不能拍脑袋定，要看非线性项的强度。以杜芬方程为例：

x'' + 2ζx' + x + ε x^3 = f cos(ωt)

当 ε 较小时，响应主要由基波控制，N 取 1 或 2 就够；但当 ε 大到 1 以上，三次非线性会耦合出明显的 3 倍频成分，N 至少要取到 3，否则共振峰位置和跳跃点会完全算错。判断标准很直接：算完后看最后一个谐波系数的幅值，如果它和最大系数相比超过 1%，就把 N 加 1。

2.2 写出残差与谐波平衡条件

将周期解截断到 N 阶，待定系数有 2N+1 个（a0 以及各阶 ak, bk）。把假设解代入运动方程，方程左边会得到一个时间函数 R(t)，它不再恒等于零，而是包含一大堆高次谐波。谐波平衡条件就是让 R(t) 在 F0, cos(kωt), sin(kωt) 这些基函数上的投影为零：

∫_0^T R(t) dt = 0 ∫_0^T R(t) cos(kωt) dt = 0 (k=1..N) ∫_0^T R(t) sin(kωt) dt = 0 (k=1..N)

这一组积分方程等价于：把 R(t) 自身也做Fourier展开，让它的前 N 阶谐波系数全部置零。对于多项式非线性，x^3 的Fourier系数是原始系数的卷积，用手算很繁琐，但用符号计算可以自动展开。

2.2.1 用SymPy验证三次非线性的谐波耦合

import sympy as sp omega, t, eps = sp.symbols('omega t eps', positive=True) a1, b1, a3, b3 = sp.symbols('a1 b1 a3 b3', real=True) # 假设只保留基波和3次谐波，且系统无常数项 x = (a1*sp.cos(omega*t) + b1*sp.sin(omega*t) + a3*sp.cos(3*omega*t) + b3*sp.sin(3*omega*t)) x3 = sp.expand(x**3) # 提取cos(omega*t)的系数 proj_cos1 = sp.integrate(x3*sp.cos(omega*t), (t, 0, 2*sp.pi/omega)) / sp.pi print("cos(omega t) in x^3:", sp.simplify(proj_cos1))

这段代码做的事情是把三次非线性项在基波上的投影解出来。integrate是Fourier投影的解析实现，除以 pi 是因为完备基的内积归一化。输出里会看到a1**3和a1*a3之类的项，这说明基波系数和被截断的高阶系数相互耦合。实际手算很容易漏掉这种耦合，这也是谐波平衡法“折成代数方程”后必须用计算机处理的原因。

2.3 谐波平衡条件与虚功率平衡的等价性

另一种理解方式：把残差 R(t) 乘以每个基函数后在周期内积分，物理上就是在周期内对“虚位移”做功为零。所以谐波平衡法也叫“Galerkin法”或“Ritz平均法”。在NLvibration这类工具里，它底层就是一个非线性最小二乘问题：

minimize ||R_vector(X)||_2^2

其中 X 是所有谐波系数的向量，R_vector 是上述投影残差组成的向量。注意不要直接对这个平方和做梯度下降，因为强非线性下目标函数高度非凸，梯度法容易卡在局部极小。正确做法是直接对 R_vector(X)=0 用Newton迭代，这样收敛才是二阶的。

2.3.1 NLvibration中“谐波数”参数如何影响方程维度

表格里的维度是每个激励频率点上的规模。谐振频率扫描时，如果频率点有 200 个，N 取 3，就要解 200 次 7×7 的非线性方程组，这个计算量微秒级，完全不是瓶颈。真正的瓶颈是 Newton 迭代里每一步都要重新计算 Jacobian 对 N 的依赖：N 增加一阶，Jacobian 的填充量大约从 O(N^2) 涨到 O(N^3)。所以在NLvibration里如果看到“运算时间暴涨”，先查是不是把基波和常数项之外的谐波数开到了 10 以上，而不是查CPU。

推导到这里，谐波平衡法已经从“把一个微分方程变成一堆积分等式”落实到了“求解一组 F(X)=0”。下面用具体代码把这个求解过程跑通。

3. 用NLvibration/Python实现周期解的数值求解：最少代码跑通Duffing

3.1 离散谐波平衡的Newton迭代框架

实现谐波平衡法的关键是写出残差函数 F(X)。X 的定义要统一：把 cos(kωt) 系数和 sin(kωt) 系数按 k 从小到大排成一维向量，常数项放在最前面。运动方程是：

x'' + 2ζ x' + x + ε x^3 = f cos(ωt)

假设截断到 N 阶，x 用 X 组合得到。把 x''、x'、x 和非线性项 x^3 都投影到基函数上，得到 F(X)。用 scipy.optimize.fsolve 或手写Newton。注意 fsolve 默认用前向差分求 Jacobian，当 N 较大时数值 Jacobian 误差明显，建议用numeric_jacobian或者手写解析 Jacobian。为了演示可复现，我先给一份手写 Newton 的代码，它不依赖任何ODE积分器，只依赖 numpy。

import numpy as np def duffing_hbm(omega=1.2, eps=0.5, zeta=0.05, f=0.3, N=3): """ 谐波平衡法求解杜芬方程 x''+2*zeta*x'+x+eps*x^3=f*cos(omega*t) 返回谐波系数向量 X, 以及残差范数, 最后一项系数幅值 """ n = 2 * N + 1 # 未知数个数: a0, a1..aN, b1..bN X = np.zeros(n) # 初始猜测: 线性解的基波余弦项 # 对线性化系统 x''+2*zeta*x'+x = f*cos(wt) # 稳态幅值 A = f / sqrt((1-w^2)^2 + (2*zeta*w)^2) # 相位偏移写在余弦项和正弦项的系数里 denom = (1 - omega**2)**2 + (2*zeta*omega)**2 A = f / np.sqrt(denom) # 当 omega 接近1时，相位接近-pi/2 phi = np.arctan2(-2*zeta*omega, 1-omega**2) X[1] = A * np.cos(phi) # a1 X[N+1] = A * np.sin(phi) # b1 (顺序: 索引0是a0, 后面接a们，再后面接b们) def projection_basis(k, order): # 返回 cos(k*w*t) 或 sin(k*w*t) 在一个周期上的采样向量 t = np.linspace(0, 2*np.pi/omega, 2048, endpoint=False) if order == 0: return np.cos(k*omega*t) else: return np.sin(k*omega*t) def residual(X): # 用连续时间采样计算残差 R(t)，再投影到基函数 t = np.linspace(0, 2*np.pi/omega, 2048, endpoint=False) # 重构 x(t) x = X[0] * np.ones_like(t) for k in range(1, N+1): x += X[k] * np.cos(k*omega*t) + X[N+k] * np.sin(k*omega*t) # 速度与加速度 xdot = np.zeros_like(t) xddot = np.zeros_like(t) for k in range(1, N+1): xdot += - X[k]*k*omega*np.sin(k*omega*t) + X[N+k]*k*omega*np.cos(k*omega*t) xddot += - X[k]*(k*omega)**2*np.cos(k*omega*t) - X[N+k]*(k*omega)**2*np.sin(k*omega*t) R = xddot + 2*zeta*xdot + x + eps*x**3 - f*np.cos(omega*t) # 投影 F = np.zeros(n) F[0] = np.mean(R) # 常数项投影 for k in range(1, N+1): F[k] = 2*np.mean(R*np.cos(k*omega*t)) F[N+k] = 2*np.mean(R*np.sin(k*omega*t)) return F # Newton 迭代 for it in range(50): F = residual(X) normF = np.linalg.norm(F, ord=np.inf) if normF < 1e-10: break # 数值 Jacobian (中心差分) J = np.zeros((n, n)) h = 1e-6 for j in range(n): Xp = X.copy(); Xp[j] += h Xm = X.copy(); Xm[j] -= h J[:, j] = (residual(Xp) - residual(Xm)) / (2*h) # 解线性方程 dX = np.linalg.solve(J, -F) X += dX # 阻尼Newton: 如果残差变大就折半 while np.linalg.norm(residual(X), ord=np.inf) > normF: dX *= 0.5 X = X - dX # 注意: 重新计算 if np.linalg.norm(dX) < 1e-14: break last_amp = np.hypot(X[N], X[2*N]) if N > 1 else 0.0 return X, np.linalg.norm(residual(X), ord=np.inf), last_amp

3.2 代码背后的三个关键参数

第一个关键参数是omega，激励频率。它决定了基函数的周期。谐波平衡法只在固定的 omega 下求解，相当于扫频时每个频率点都是独立求根。第二个是zeta，阻尼比。阻尼太大时高阶谐波会被压制，N 可以取小；阻尼接近零时，共振峰很尖锐，Newton 迭代容易从峰的一侧跳到另一侧，需要在下一章讲延拓。第三个是eps，非线性系数。eps 为 0 时方程退化为线性，系统只有单一解，谐波平衡法直接退化成频响函数；eps 增大后，共振峰向右侧弯曲（硬弹簧特性），同时出现多解区间，所以它才是整个代码里最需要关注的值。

代码中的2048个采样点是对残差做数值投影。这个点数的选择也有讲究：因为余弦和正弦函数正交性依赖周期，如果频率点取得不正好是周期端点，就会泄漏。这里用np.linspace(0, 2*pi/omega, 2048, endpoint=False)正好覆盖一个整数周期，可以避免泄漏。如果点数太少，比如 128 点，在 N 大于等于 5 时，高频投影会出现明显混叠，导致 Newton 迭代在达到机器精度前就停滞。点数也不用太多，2048 对双精度浮点和 10 阶以内的谐波已经绰绰有余。

3.3 从线性解起步为什么是可靠的第一个猜测

非线性方程求根不像线性方程，Newton 迭代必须给一个靠得住的X0。线性解起步是最自然的：先把非线性项去掉，得到线性频响函数。在线性系统里，x(t) 的振幅和相位可以解析给出，把它作为 N 阶谐波解的初始猜测，在中等非线性强度下 Newton 一般三到五步就能收敛。如果 epsilon 较大，线性解作为初始猜测可能落在牛顿法的收敛域之外，这时 NLvibration 的做法是“增量加载”：先把 eps 设成 0.1 跑一遍，收敛后把结果作为 eps=0.2 的初值，逐步升到目标值。

3.3.1 收敛失败时先看残差曲线而不是先调初值

很多人在谐波平衡法不收敛时，第一反应是改初始猜测。实际上更有效的做法是画残差函数R(t)的时域曲线。如果残差在单个周期内呈现光滑波动，但投影后的 F 范数降不下去，这说明 N 截断不够，只增大谐波阶数即可。如果残差曲线呈锯齿状，则是采样点数不足或基函数内积泄漏。如果残差在某些时间段尤其大，且 N 增加后残差峰值没有下降，那问题一定出在非线性项投影符号上——比如把 x^3 的系数符号写反，或者漏了常数项。建议在调试时把residual(X)返回的 F 也返回时域残差 R，打印几个典型时刻的值。

4. 非线性振动分析中的强非线性问题：弧长延拓与多解追踪

4.1 为什么共振区附近牛顿法会跳变

接近共振峰时，谐波平衡方程组的解曲线在幅值-频率平面上呈现 S 形。S 形的上下两个分支是稳定解，中间分支是不稳定解。用固定频率点做逐点扫频时，Newton 迭代的结果会在某个频率点上突然从低幅值分支跳到高幅值分支，这不是程序bug，而是因为牛顿法是在找“最近的根”，而 S 形区域里同一个频率下存在三个根，初始猜测决定了它收敛到哪一个。要完整画出这条 S 形曲线，必须沿着解曲线本身推进，而不是沿着频率推进。

4.2 用伪弧长延拓扫频的落地方式

伪弧长延拓的思想是把频率 omega 也当成未知量，引入一个弧长参数 s，额外添加一个约束方程。NLvibration 的常见实现如下：

X = X_prev + ds * tangent_X omega = omega_prev + ds * tangent_omega

然后对扩展后的方程组加上球面约束：

||X - X_new||^2 + (omega - omega_new)^2 - ds^2 = 0

这里的ds是步长。步长不能固定，要加上自适应逻辑：本轮 Newton 迭代超过 8 次才收敛，就把 ds 减半；少于 3 次收敛，则下一轮扩大 1.5 倍。弧长延拓的另一个好处是能自然通过转向点（saddle-node），因为约束方程让迭代方向始终沿着解曲线走。

# 伪弧长延拓的核心步骤（伪代码，省略Jacobian组装） # 已知点 (X0, w0)，切向量 (dX0, dw0)，步长 ds，预测 X_pred = X0 + ds * dX0 w_pred = w0 + ds * dw0 # 校正：用Newton法求解增广残差 # 其中残差 F(X,w)=0 是谐波平衡残差，新增约束 # C(X,w)= (X-X0).T*(X-X0) + (w-w0)**2 - ds**2 = 0 # 每轮迭代求解 (2N+2) 维线性方程组 J_aug = np.block([ [J_hbm, dF_dw], [2*(X-X0), 2*(w-w0)] ]) # 然后解 J_aug @ delta = -residual_aug

J_hbm是谐波平衡残差对 X 的 Jacobian，dF_dw是残差对 omega 的偏导。dF_dw 的解析式来自运动方程里 omega 只出现在 cos(omega t) 和 sin(omega t) 的自变量中，以及激励项 f cos(omega t) 的频率位置。不要把 dF_dw 用差分近似，因为靠近转向点时差分误差会导致切向量方向反号。

4.2.1 弧长延拓结果如何判定多解区间

实际操作中，我会用向上扫频和向下扫频各跑一次，把两条幅值曲线画在同一张图上。两张图在共振峰附近围出的滞后环就是多解区间。这个区间边界正好对应 S 形曲线的两个转向点。如果只用单方向扫频，你永远不会意识到那个跳跃其实包含两个稳定的周期解和一个不稳定的周期解。

4.3 周期解的稳定性判断：Floquet理论还是简谐判据

求出了周期解不等于它物理上能出现。NLvibration 里一般在得到谐波系数后计算单值矩阵（Monodromy），通过 Floquet 特征乘子实部是否穿过 +1 来判断。对于单自由度系统有个更快的办法：把周期解代入变分方程：

δx'' + 2ζ δx' + (1 + 3ε x(t)^2) δx = 0

在周期解基础上做小扰动。如果 x(t) 的幅值在多个周期内衰减，则稳定；否则不稳定。在扫频延拓过程中，观察 Jacobian 矩阵行列式是否改变符号是判断转向点的常用技巧，但在转向点处 Jacobian 奇异，行列式过零不能直接当作稳定性翻转。更稳妥的是跟踪单值矩阵最大特征乘子的模长变化。

5. 验证周期解：把谐波平衡结果交给时域积分做交叉检查

5.1 用 RK4/odeint 对比一个周期内的漂移

谐波平衡法给出的是周期解系数，要验证它是否正确，最直接的方法是把它作为初始条件扔给时域积分器，积分若干周期，看轨迹是否还停留在原始解的附近。以 scipy 的solve_ivp为例：

from scipy.integrate import solve_ivp def duffing_rhs(t, y, omega, zeta, eps, f): # 状态向量 y = [x, v] return [y[1], f*np.cos(omega*t) - 2*zeta*y[1] - y[0] - eps*y[0]**3] # 假设 hbm_X 是上面谐波平衡法得到的系数向量，N=3 omega_val = 1.2 t_span = (0, 80*2*np.pi/omega_val) # 从谐波平衡解重构初始状态 t0 = 0 x0 = hbm_X[0] + sum(hbm_X[k]*np.cos(k*omega_val*t0) + hbm_X[N+k]*np.sin(k*omega_val*t0) for k in range(1, N+1)) v0 = sum(-hbm_X[k]*k*omega_val*np.sin(k*omega_val*t0) + hbm_X[N+k]*k*omega_val*np.cos(k*omega_val*t0) for k in range(1, N+1)) sol = solve_ivp(duffing_rhs, t_span, [x0, v0], args=(omega_val, zeta, eps, f), rtol=1e-10, atol=1e-10) # 比较最后一个周期与谐波平衡解的形态 t_span_end = np.linspace(sol.t[-200], sol.t[-1], 200)

这段代码的验证逻辑是：如果谐波平衡解是正确的周期解，把它当初始状态，积分 80 个周期后轨迹应该几乎不漂移。观察最后 200 个采样点的幅值变化：如果相对误差小于 1e-6，说明谐波截断充分；如果漂移明显，则意味着该周期解在动力学上不稳定，即使代数上满足谐波平衡，在实验中也不会出现。

5.2 验证时注意的三个坑

第一个坑是积分总时长要足够长，否则瞬态衰减没有完，你会把暂态漂移误判成解不稳定。我一般会做双保险：分别积分 20 个周期和 80 个周期，比较末尾一个周期的位移幅值。第二个坑是积分器容差设置太低。谐波平衡法本身可以到机器精度，但如果 RK45 的容差只放到 1e-6，可能把高阶谐波的误差放大。建议至少rtol=1e-10。第三个坑是初值不能用“在某个时刻的瞬时位移和速度”组合，而必须保证该初值严格位于重构的周期轨上。如果代码里补一个周期内的位移重构和 N+1 个等距点的采样对比，就能同时检查重构函数没有相位偏移。

5.3 一个省事的小技巧：残差能量百分比

在 NLvibration 输出结果里，除了谐波系数，还应该输出一个能量残差指标：

eta = sqrt(sum_{k=N+1}^{2N} (proj_k)^2) / sqrt(sum_{k=1}^{N} (proj_k)^2)

其中 proj_k 是把运动方程残差投影到第 k 阶基函数上的幅值。这个指标告诉你被截断掉的高阶谐波里还藏着多少残余能量。eta 小于 0.01 时，说明当前 N 已经足够；eta 在 0.01 到 0.05 之间时，结果还能用，但稳定性边界会有一点误差；eta 大于 0.05 时，增加 N 之前先检查采样点数和非线性项投影是否写对。把这个指标和时域交叉验证一起打到结果里，比只贴一条幅值曲线更能说服自己：谐波平衡法求出的周期解既代数可解，又物理可达。

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