用MATLAB实现分数阶振动模型：粘弹性阻尼与短记忆法求解指南

搞机械振动的人，手里那把整数阶模型有时候真的不够用。你按达朗贝尔原理老老实实写出 (m\ddot{x}+kx=0)，算出的固有频率和实验对得上，可一旦材料换成橡胶、黏弹性阻尼器，或者你去看高分子复合梁的衰减曲线，理论解跟实测数据就明显分家——衰减速率按时间变化、频响曲线也拖个尾巴。把二阶导数替换成Caputo分数阶导数，让“阶次”从2连续变化到1.8甚至1.5，一套方程就能把这个“中间状态”的动力学行为描述出来。我这次就把经典的双耦合弹簧质量系统改造成分数阶形式，用MATLAB从零写求解器，做了几组参数实验，把时域响应、频谱、参数影响一条龙跑通了。项目不大，思路却很有代表性：你会看到分数阶导数怎么进入运动方程、GL短记忆法怎么离散、初始条件怎么修正，以及阶次 (\alpha) 对振动衰减和频率到底有什么影响。适合正在做分数阶建模、振动仿真、或准备用MATLAB写数值算法的朋友参考。

1. 分数阶振动模型：为什么要把“2”改成1.8

1.1 经典双耦合模型的局限在哪里

经典双耦合弹簧质量系统大家很熟：两个质量块 (m_1)、(m_2)，各自通过弹簧 (k_1)、(k_2) 接到固定端，中间再串一个耦合弹簧 (k_3)。忽略阻尼时，运动方程是：

[ \begin{cases} m_1\ddot{x}_1+(k_1+k_3)x_1-k_3x_2=0\[2mm] m_2\ddot{x}_2+(k_2+k_3)x_2-k_3x_1=0 \end{cases} ]

这个模型的好处是解析解完整，两个固有频率写得很漂亮；坏处是它只适合描述“理想弹簧振子”这种几乎没有能量耗散的场景。真实工程里的结构，材料往往带粘弹性，应力不仅取决于当前应变，还和过去一段时间内的应变历史有关。你把一块橡胶拉长再松开，它不会像理想弹簧那样干净利落地回到原位，而是带着明显的滞后和残余松弛。经典整数阶方程里，要描述这种效果必须额外加阻尼项 (c\dot{x})，但实际测量会发现，单纯的粘性阻尼也不够用——阻尼力往往和频率有关，低频段和高频段的衰减特性差别很大。这时候就需要一个更灵活的工具，让“导数阶次”本身成为可调参数。

1.2 Caputo导数补上了什么

Caputo分数阶导数的定义是：

[ {}^C_{t_0}D^\alpha_t x(t)=\frac{1}{\Gamma(n-\alpha)}\int_{t_0}^{t}\frac{x^{(n)}(\tau)}{(t-\tau)^{\alpha-n+1}},d\tau ]

其中 (n-1<\alpha\le n)。当 (\alpha=1) 时它退化为普通一阶导数，当 (\alpha=2) 时退化为普通二阶导数，所以分数阶模型天然兼容经典模型。真正关键的是积分核 ((t-\tau)^{\alpha-n+1})：它把过去每一时刻的变化量按权重累计到当前，(\alpha) 越小，久远历史的权重衰减越慢，“记忆”越强。

工程上选Caputo定义而不是Riemann-Liouville定义，有个非常实际的原因：Caputo导数的初值条件跟整数阶完全一致，只需要给 (x(0)) 和 (\dot{x}(0))。这意味着实验里测到的初始位移、初始速度可以直接代入方程，不用去算那些物理意义模糊的分数阶积分初值。对做实验的人而言，这是决定性的优势。

1.3 双耦合系统的分数阶运动方程

把经典方程里的二阶导数换成Caputo分数阶导数，就得到：

[ \begin{cases} m_1,{}^C!D^{\alpha_1}_t x_1+(k_1+k_3)x_1-k_3x_2=0\[2mm] m_2,{}^C!D^{\alpha_2}_t x_2+(k_2+k_3)x_2-k_3x_1=0 \end{cases} ]

这里的 (\alpha_1)、(\alpha_2) 可以不同，表示两个质量块连接的阻尼材料特性不一样。 (\alpha_1=\alpha_2=2) 时回到经典无阻尼模型； (\alpha_1=\alpha_2) 略小于2时，系统相当于自带“分数阶阻尼”——不需要额外加 (c\dot{x}) 项，能量耗散已经被包含在分数阶算子内部。这个性质的物理直觉是：分数阶算子在频域里同时调制幅值和相位，既改变等效刚度，又引入阻尼角，所以一个参数能同时影响频率和衰减。

我最初看到这个方程时也有点怀疑：就改一个阶次，真能描述粘弹性阻尼吗？但后续仿真结果说明，(\alpha) 从2降到1.6左右的响应，和实验里粘弹性材料自由衰减曲线的趋势非常像，这比硬凑一个阻尼系数要自然得多。

2. 数值求解与MATLAB实现

2.1 GL离散与短记忆原理

Caputo分数阶微分方程一般没有通用解析解，数值上最常用的路径是Grünwald-Letnikov（GL）离散。GL形式的分数阶导数可以写成：

[ D^\alpha x(t_k) \approx h^{-\alpha}\sum_{j=0}^{k} w_j^{(\alpha)} x(t_k-jh) ]

其中系数 (w_j^{(\alpha)}) 有一个非常漂亮的递推关系：

[ w_0=1,\qquad w_j=\left(1-\frac{\alpha+1}{j}\right)w_{j-1},\quad j=1,2,\dots ]

这个递推的好处是写代码时不用算Gamma函数，几行循环就能生成全部系数。还有个隐藏福利：当 (\alpha=2) 时，(w_0=1, w_1=-2, w_2=1)，之后全部为0，GL公式自动退化为中心二阶后向差分。这意味着同一套代码可以同时验证整数阶和分数阶，非常方便。

GL求和理论上需要从 (j=0) 加到 (j=k)，也就是要记住整个时间历程。好在系数 (w_j) 会随 (j) 增大逐渐趋近于0，尤其是阶次离整数不远的时候，衰减很快。于是就有了短记忆法：只保留最近 (L_m) 个点，忽略更早的历史贡献，计算量从 (O(N^2)) 降到 (O(NL_m))。

2.2 初始条件修正：最容易踩的坑

直接用GL公式去解Caputo方程，这里藏着一个大坑：GL离散天然逼近的是Riemann-Liouville导数的框架，而Riemann-Liouville导数与Caputo导数之间差了一个由初始条件引起的边界项。如果忽略这个差别，初值不为零时算出来的响应往往长期偏离正确解，而且怎么调步长都调不好。

修正方法不复杂。令：

[ P(t)=x(0)+\dot{x}(0)t ]

这是由初始位移和初始速度构成的线性多项式。把变量换成 (z(t)=x(t)-P(t))，可以证明：

[ {}^C!D^\alpha_t x(t)={}^{RL}!D^\alpha_t z(t) ]

也就是说，只要在GL求和时把每个历史点都换成 (x(t_j)-P(t_j))，就能用GL离散逼近Caputo导数。我习惯把这个修正理解为摄像里的“参考帧剔除”：先把系统状态搬到零点，算完相对运动，再把它搬回到真实位移上。不剔除，积分核会把初始状态的记忆反复叠加，结果自然飘掉。

2.3 隐式线性系统推导

把GL离散代入双耦合方程后，变量 (x_{1,k}) 和 (x_{2,k}) 会同时出现在当期项中。处理方法是把所有含当期项的部分集中到左边，整理成一个2×2线性方程组。

以质量块1为例，GL离散代入后得到：

[ h^{-\alpha_1}\left(x_{1,k}+\sum_{j\ge1}w_j^{(1)}z_{1,k-j}\right) +\frac{k_1+k_3}{m_1}x_{1,k}-\frac{k_3}{m_1}x_{2,k}=0 ]

两边乘 (h^{\alpha_1}/m_1) 并移项：

[ \left(1+\frac{h^{\alpha_1}(k_1+k_3)}{m_1}\right)x_{1,k} -\frac{h^{\alpha_1}k_3}{m_1}x_{2,k} = P_1(t_k) - \sum_{j\ge1}w_j^{(1)}z_{1,k-j} ]

这里右边 (P_1(t_k)) 就是初始条件修正中减掉再搬回的线性基准项。质量块2做同样操作后，得到完整的隐式格式：

[ \begin{bmatrix} 1+\frac{h^{\alpha_1}(k_1+k_3)}{m_1} & -\frac{h^{\alpha_1}k_3}{m_1}\[2mm] -\frac{h^{\alpha_2}k_3}{m_2} & 1+\frac{h^{\alpha_2}(k_2+k_3)}{m_2} \end{bmatrix} \begin{bmatrix} x_{1,k}\ x_{2,k} \end{bmatrix}

\begin{bmatrix} P_1(t_k)-\sum w_j^{(1)}z_{1,k-j}\ P_2(t_k)-\sum w_j^{(2)}z_{2,k-j} \end{bmatrix} ]

每一步都解这个2×2方程组。选隐式而不是显式，是因为显式格式在阶次接近整数、系统刚度偏大时容易振荡发散；隐式格式多写两行代码，但稳定性和步长宽容度好得多。

2.4 可直接运行的MATLAB代码

下面是我跑通的主程序，参数都放在文件头部，直接改数值就能复用。

%% 基于Caputo导数的双耦合弹簧质量系统自由振动仿真 % 初速度为零；初始位移激励；GL短记忆 + 初始条件修正 clear; clc; close all; %% 系统参数 m1 = 1.2; m2 = 0.8; k1 = 120; k2 = 80; k3 = 40; alpha1 = 1.8; alpha2 = 1.7; %% 数值参数 T = 6; % 时长 s h = 0.002; % 步长 s N = round(T/h); t = (0:N)*h; fs = 1/h; %% 初始条件 x10 = 0.05; x20 = -0.03; v10 = 0; v20 = 0; P1 = x10 + v10*t; P2 = x20 + v20*t; x1 = zeros(1,N+1); x1(1) = x10; x2 = zeros(1,N+1); x2(1) = x20; %% GL短记忆长度与系数 Lm = min(N, round(2/h)); % 2秒记忆 w1 = zeros(1,Lm+1); w1(1) = 1; w2 = zeros(1,Lm+1); w2(1) = 1; for j = 1:Lm w1(j+1) = (1 - (alpha1+1)/j) * w1(j); w2(j+1) = (1 - (alpha2+1)/j) * w2(j); end %% 隐式离散矩阵 A = [1 + h^alpha1*(k1+k3)/m1, -h^alpha1*k3/m1; -h^alpha2*k3/m2, 1 + h^alpha2*(k2+k3)/m2]; %% 主循环 for k = 2:N+1 s1 = 0; s2 = 0; for j = 1:min(k-1, Lm) idx = k - j; s1 = s1 + w1(j+1) * (x1(idx) - P1(idx)); s2 = s2 + w2(j+1) * (x2(idx) - P2(idx)); end rhs = A \ [P1(k) - s1; P2(k) - s2]; x1(k) = rhs(1); x2(k) = rhs(2); end %% 时域图 figure('Color','w','Position',[60 60 1000 650]); subplot(2,1,1); plot(t, x1*1000, 'b-', 'LineWidth',1.1); hold on; plot(t, x2*1000, 'r-', 'LineWidth',1.1); grid on; ylabel('位移 (mm)'); xlabel('t (s)'); legend('x_1','x_2','Location','best'); title('自由振动时程响应'); %% 频谱 Nf = 2^nextpow2(N+1); X1 = fft(x1 - mean(x1), Nf); X2 = fft(x2 - mean(x2), Nf); f = (0:Nf-1)*fs/Nf; subplot(2,1,2); plot(f, abs(X1), 'b-', 'LineWidth',1.1); hold on; plot(f, abs(X2), 'r-', 'LineWidth',1.1); grid on; xlim([0 10]); xlabel('频率 (Hz)'); ylabel('幅值'); legend('x_1','x_2'); title('FFT频谱');

代码里两个细节建议关注下。第一，(P_1)、(P_2) 是按初始速度构造的线性基准，如果以后想改成有初速度的工况，只需要改v10、v20，修正逻辑不用动。第二，(A) 矩阵在定常线性系统里是常数矩阵，提前算一次就够了；若后面做参数扫描，每改一组参数重算一次即可。

2.5 验证：让α退回2，和ode45对齐

写数值代码最怕“看着合理但错得离谱”，所以第一件事就是验证。把上面程序里的alpha1=2; alpha2=2;，再用ode45解经典整数阶方程：

f = @(t, y) [y(3); y(4); -((k1+k3)*y(1) - k3*y(2))/m1; -(k3*y(2) - k3*y(1) + k2*y(2))/m2]; [t45, y45] = ode45(f, [0 T], [x10 x20 0 0]);

两个结果对比，曲线几乎重合。原因前面说过，(\alpha=2) 时GL系数自动变成 (1,-2,1)，离散格式退化到标准二阶后向差分，所以框架本身不会偏。验证通过后，再把 (\alpha) 改回1.8，看到的就是分数阶模型带来的真实差异。

3. 参数实验设计与结果分析

3.1 实验方案：扫阶次、扫耦合刚度

模型跑通以后，我用同一套代码做了三组数值实验。第一组是基准对比：整数阶 (\alpha=2) 对分数阶 (\alpha_1=1.8,\alpha_2=1.7)。第二组扫描阶次：(\alpha) 从1.4到1.9，步长0.1，两个质量块取相同阶次。第三组扫描耦合刚度：(k_3) 从20扫到60，步长10，阶次固定在1.8。

我记录的指标有三个：前两秒FFT谱峰对应的主频；自由响应包络衰减到初始峰值的50%所需时间；两个质量块之间能量交换的拍频周期。包络提取用的是MATLAB自带的hilbert取幅值，简单稳定。

3.2 阶次α对振动响应的显著影响

先看阶次扫描结果。以本组参数为例，整数阶模型两个主频大约在1.59Hz和2.15Hz，且几乎不衰减。随着 (\alpha) 减小，现象非常清晰：

规律很直白：阶次越低，等效阻尼越强，主频越低，衰减越快。背后的原理是分数阶算子在频域里相当于一个幅值和相位同时变化的“柔性算子”，阶次下降时高频成分被压制得更狠，系统的等效刚度下降，所以频率往下走；而记忆核让能量被持续耗散，所以衰减变快。这个趋势和我用粘弹性材料测过的自由衰减曲线从形态上是对得上的。

3.3 耦合刚度k3改变模态分离程度

再看 (k_3) 的影响：

(k_3) 越大，两个主频离得越远，拍频周期越短。这和整数阶耦合振动系统的结论一致，但分数阶版本多了一个变化：阶次越低，同等 (k_3) 下拍频周期越长，能量交换越慢。可以这样理解——分数阶阻尼消耗了一部分可交换的能量，两个质量块之间的“能量跷跷板”摆动幅度变小、速度变慢。

3.4 频谱与能量交换特征

FFT谱在整数阶情况下是两根又尖又高的谱线；改成分数阶后，谱峰变矮、变宽，位置略左移。谱峰变宽对应半功率带宽变大，这是阻尼增加的直接证据。时域里两质量块振幅包络呈周期性交替，一个块振幅大时另一个块振幅小；分数阶系统中这种交替依然存在，但每一轮交换都会损失一部分能量，所以包络整体呈衰减趋势。

4. 常见问题与调试技巧

4.1 发散或输出NaN

如果一运行结果就飞到天上，最常见原因是步长 (h) 太大。分数阶方程的离散核比较敏感，尤其阶次接近2时，后向差分的一阶精度会放大局部误差。我的建议：先把 (h) 设成0.001，稳定后再逐步放大看结果是否一致。第二个原因是短记忆长度 (L_m) 太短。阶次偏低比如1.4时，(w_j) 衰减较慢，2秒记忆可能不够，需要加大到5秒甚至更长。第三个原因很隐蔽：初始条件修正没做或者修正后的 (P_1,P_2) 没更新到循环里，导致常数基准项和实际初始状态不匹配。

4.2 和ode45的结果对不齐

绝大多数情况是 (h) 不够小造成的。GL后向差分截断误差是一阶的，而 (ode45) 是四阶精度，两者直接对比时步长必须压得很低。另一种可能：有人试图把分数阶方程改写成“整数阶方程加阻尼”，然后让ode45去解这个替代系统。这个做法原则上是不成立的，因为分数阶算子的频响特性不能简单用几个整数阶项拼出来。如果只是想做验证，我建议找一个第三方分数阶求解器交叉验证，比如Garrappa公开的fde12，拿同一组参数对比趋势，只要变化形态一致就说明实现没问题。

4.3 计算慢怎么办

短记忆长度 (L_m) 直接决定耗时。我测试时2秒记忆已经能得到稳定结果，没必要保留全历史。如果要做参数扫描，先在 (h=0.005) 粗步长下扫一遍找趋势，锁定感兴趣区间后再缩小步长精算。另一种提速思路是把内层求和写成卷积形式，MATLAB的conv对一维数组非常擅长，复杂度远低于for循环。我试过把核心循环改成卷积后，(N=6000) 的情况下单次仿真时间从十几秒降到一秒以内。

4.4 从实测数据粗估阶次α

做实验时，阶次 (\alpha) 不是直接量出来的，要靠识别。我常用的土办法有两个。第一个是时域包络法：对自由衰减响应做Hilbert包络，在双对数坐标里拟合包络的衰减斜率。阶次越低，包络衰减斜率绝对值越大。第二个是频域法：从实验频响函数里读半功率带宽，算出等效阻尼比，再用不同 (\alpha) 的仿真结果做一张阻尼比-阶次对照表，反查阶次。这个方法不需要高深优化算法，误差大概在±0.05以内，工程预研阶段完全够用。

5. 工程扩展与我的实操心得

5.1 扩展到多自由度和非线性系统

双耦合系统只是个起点。把2×2矩阵推广到 (N\times N)，就变成多自由度分数阶振动系统，每步要解的线性方程组变大了一些，但GL求和和初始条件修正的框架不用变。如果弹簧是非线性的，比如Duffing型恢复力，那么 (A) 矩阵在每一步都会变化，需要在每个时间步内迭代求解。代码结构上仍然是“算历史项、组装矩阵、解方程”三步，只是组装部分从常数矩阵变成随状态更新的矩阵。

5.2 参数识别和分数阶控制方向

一旦阶次 (\alpha) 能从实验数据里识别出来，这个模型就能用于实际系统：可以用它做振动响应预测，也可以进一步设计分数阶PID控制器。分数阶PID比整数阶PID多两个自由度，对阻尼和刚度同时变化的系统适应性更强，这也是很多课题组把它用在做粘弹性隔振器控制上的原因。做控制仿真时，上面这套MATLAB求解器可以直接当作被控对象模型。

5.3 做真实实验时的几点建议

如果想拿实验数据验证这个模型，选材上建议优先用橡胶、聚氨酯这类粘弹性材料，它们比金属更容易展现出明显的分数阶特征。激振方式推荐锤击法或低频扫频，尽量避免初始速度过大，保证初速度条件可以被忽略。采样率至少500Hz以上，因为识别阶次需要用足够长的衰减段，采样率太低会把衰减尾巴的细节滤掉。后处理时先做低通滤波再提取包络，否则高频噪声会把代数值衰减趋势淹没。

我最初用GL法直接把Caputo方程当RL方程解，结果初始位移不为零时响应一直不对劲，后来才发现是初始条件修正没做。分数阶数值计算里，这种“看起来能用、一深究就翻车”的坑很多。你要是自己动手做，建议先用 (\alpha=2) 把代码验证一遍，再慢慢把阶次往下调；这样出现异常时，至少能分清是数值框架的问题，还是分数阶模型本身的动力学行为。整套流程跑通之后你会发现，分数阶模型其实不玄乎，它就是给经典力学方程加了一份“记忆档案”，让仿真结果更贴近真实材料的脾性。