简介面向分数阶系统与混沌动力学研究者的Matlab工具包专注非相称分数阶自治连续时间系统的李雅普诺夫指数计算。代码基于Caputo导数建模提供主函数与辅助函数覆盖系统模型定义、分数阶微分方程数值求解到李雅普诺夫指数提取的完整流程无需额外工具箱即可直接运行。压缩包体积仅3KB共包含两个m文件结构精简、层次清晰便于科研人员阅读、复用与二次开发。该代码包已有五十三人学习适用于自动化控制、信号处理、电路设计等涉及分数阶动态特性的场景。借助此代码研究者可以快速评估非相称分数阶系统是否具有正的李雅普诺夫指数进而判定系统的混沌行为并为稳定性分析、混沌控制及系统优化设计提供量化依据是相关领域理论验证与工程落地的实用工具。1. 非相称FO系统的LE计算为什么常规代码会失效做分数阶混沌的LELyapunov指数时Matlab里能找到的现成脚本九成只支持相称系统所有状态变量的导数阶次共用一个 q。可换成非相称分数阶non-commensurate FO系统后同样的代码跑出来的LE谱要么整体偏移要么连正负号都不可靠。问题出在三个地方——记忆系数只有一套、变分方程各分量仍然使用同一个阶次、步长幂次 h^q 被写死成标量。相称代码里这些假设通常是隐式的不报错所以结果错了也找不到原因。这篇文章给出一个可以直接运行的Matlab方案用Grunwald-LetnikovGL离散化同时求解状态方程和变分方程通过QR分解累积LE谱并给出验证、参数设置和三个常见坑。适合正在写分数阶混沌论文的研究生以及想把整数阶混沌系统迁移到分数阶场景的工程技术人员。2. 非相称分数阶系统LE从Caputo定义到变分方程2.1 相称与非相称到底差在哪分数阶系统的常见建模形式是Caputo型导数方程组D^{q_i} x_i(t) f_i(x_1, x_2, ..., x_n), i 1, 2, ..., n当所有 q_i 都等于同一个 q 时系统是相称的commensurate只要存在 i≠j 使得 q_i≠q_j系统就是非相称的non-commensurate。表面上看只是把阶次向量从标量换成数组但离散化时差异是结构性的。以GL离散化为例相称系统只需要维护一条记忆系数序列 w_j(q)非相称系统则要为每个状态变量各维护一条因为每个分量的记忆核衰减速度不同。更麻烦的是变分方程。建立Lyapunov指数需要同时求解扰动传播方程其第 i 个分量的导数阶次是 q_i 而不是统一的 q这导致Jacobian矩阵乘上扰动向量之后必须按行乘以不同的步长幂次 h^{q_i}。这两处改动在数学上都不难但大部分现成代码没有把它们实现出来。下表对比了相称与非相称系统在数值实现上的关键差异处理对象相称系统 q_i q非相称系统 q_i 不全相等GL记忆系数1套w(q)n套w(q_i)i1..n状态更新步长因子h^q 统一乘到右端第i分量乘 h^{q_i}变分方程离散化整个向量用 h^qJacobian×扰动后按行用 h^{q_i}稳定性与收敛阶与经典分数阶ODE一致各方向收敛不同步需缩短步长“把 q 从标量改成向量”不是对相称代码做简单替换而是要把上述四个维度全部改成逐分量处理否则LE谱必然失真。2.2 变分方程与QR分解LE谱怎么从轨迹里算出来设原始系统 D^{q_i} x_i f_i(x)。对初始扰动 δx 做线性化得到变分方程D^{q_i} δx_i Σ_{j1}^n J_{ij}(x(t)) δx_j其中 J(x) 是系统右端 f 的Jacobian矩阵。注意这个方程组里不同 δx_i 分量的导数阶次仍然各不相同这是非相称LE计算的数学难点。联立原始系统和变分方程之后对扰动矩阵 V每列是一个扰动方向做数值积分。为防止所有扰动方向在最大李雅普诺夫指数方向上坍缩需要周期性执行Gram-Schmidt正交化——工程上直接用QR分解实现。经典做法是每走一小步把当前的 V 做一次QR分解V Q R把 R 的对角元素绝对值取对数并累加每一步把 V 重置为 Q继续积分。最终得到的是λ_i ≈ (1/T) Σ log | R_ii |这个流程和整数阶混沌系统的Wolf算法非常接近区别在于每一步积分V时用到的记忆系数和步长幂次都按非相称阶次逐行处理。2.3 为什么不能直接套用整数阶LE代码整数阶LE代码和分数阶之间隔着两个不可忽略的因素记忆效应和步长幂次。整数阶系统没有历史记忆项当前时刻的导数只依赖当前状态分数阶系统则依赖从 0 到当前时刻的全部轨迹而且阶次越低记忆衰减越慢。这意味着变分方程的离散化不能写成简单的 δx_{k1} δx_k h·J·δx_k必须还原成带卷积和的GL格式。另一个容易被忽略的问题是非相称系统各分量的 h^{q_i} 不同步长稍大一点不同分量间的误差就会失去平衡。许多从整数阶迁移过来的代码直接在残差项里漏掉历史卷积或者只对状态方程补了卷积、变分方程仍用整数阶格式得到的结果看起来像那么回事拿相称退化系统一验就露馅。因此一个专门针对非相称系统的LE计算函数是必要的下面给出完整实现。3. Matlab实现基于GL离散化的非相称FO LE代码3.1 GL记忆系数初始化每个阶次一条系数链GL离散化的核心是把 Caputo 分数阶导数近似为加权历史差值之和D^q f(t) ≈ h^{-q} Σ_{j0}^{∞} (-1)^j C(q, j) f(t - jh)定义系数 w_j (-1)^j C(q, j)代入组合数定义可以得到递推关系w_0 1, w_j w_{j-1} · (j - 1 - q) / j在Matlab里按列存储不同阶次的系数链。实际代码中不会真正累加到无穷而是设置一个记忆长度 M只取最近 M 个历史时刻。下面这段代码初始化 n 套GL系数%% 参数区 q [0.90; 0.92; 0.94]; % 非相称阶次向量 n numel(q); % 系统维数 h 0.005; % 步长 M round(5 / h); % 记忆长度取最近5秒历史 %% GL系数矩阵每列对应一个q_i G zeros(M, n); for i 1:n G(1, i) 1; % w_0 1 for j 2:M G(j, i) G(j-1, i) * (j - 2 - q(i)) / (j - 1); end end % G(j,i) w_{j-1}(q_i)调用时G(1)是w0G(2)是w1依此类推这段代码把阶次向量 q 的每个分量分别生成一条记忆系数链。递推式里 (j - 2 - q(i)) 来自组合数递推 w_j / w_{j-1} (j-1-q)/j把数组索引偏移量也考虑进去。注意 q(i) 越接近 1系数衰减越快q(i) 越小历史记忆的影响持续越久需要足够大的 M 才能保证精度。3.2 状态方程和变分方程联立求解非相称系统的GL一步格式从下面的近似式出发D^q x(t) ≈ h^{-q} Σ_{j0}^{∞} w_j x(t - jh)经过整理得到显式更新公式。状态方程按行更新%% 主循环骨架假设X是n×N的状态历史矩阵 % X(:,1)是初始状态X(:,k)是当前时刻状态 % 状态更新第i个分量 for i 1:n conv 0; % 历史卷积累加项 L min(k, M-1); % 当前可用的记忆长度 for j 1:L conv conv G(j1, i) * X(i, k1-j); end X(i, k1) h^q(i) * fk(i) - conv; end这段代码的含义是当前状态的第 i 个分量等于 h^{q_i} 乘以系统右端项再减去过去 M 个时刻该分量乘上对应GL系数的加权和。注意乘在右端的是 h^{q_i} 而不是同一个 h 的固定幂这是非相称系统特有的处理方式。变分方程的离散化思路完全相同。把变分矩阵 Vn×n展平成长度 n^2 的向量存入历史矩阵 V_hist第 i 个分量的更新同样使用第 i 条GL系数链%% 变分方程更新 % vk是当前n×n变分矩阵Jk是当前JacobianV_hist是(n*n)×N历史 Fv Jk * vk; % 整数阶右端尚未乘h^q_i for i 1:n for m 1:n row (m - 1) * n i; % 展平后的行号 conv 0; L min(k, M-1); for j 1:L conv conv G(j1, i) * V_hist(row, k1-j); end V_hist(row, k1) h^q(i) * Fv(i, m) - conv; end end外层循环的 i 对应变分分量的导数阶次内层循环的 m 对应第 m 个扰动方向。每一次都要用各自的 q(i) 做卷积不能像相称系统那样用一个 G 向量统一处理。3.3 完整函数与调用示例把上述逻辑组合在一起加上QR分解和LE累积就是下面这个可运行函数。以非相称分数阶Chen系统为例function LEs LE_noncommensurate_fo() % 非相称分数阶Chen系统的Lyapunov指数谱计算 % 系统 % D^q1 x a(y - x) % D^q2 y (c-a)x - xz c y % D^q3 z xy - b z % 阶次 q 互不相同 - non-commensurate a 35; b 3; c 28; q [0.90; 0.92; 0.94]; n numel(q); T 150; h 0.005; N round(T/h); M round(5 / h); % 记忆长度 G zeros(M, n); for i 1:n G(1, i) 1; for j 2:M G(j, i) G(j-1, i) * (j - 2 - q(i)) / (j - 1); end end X zeros(n, N); % 状态历史 Vh zeros(n*n, N); % 变分矩阵展平历史 X(:,1) [1; 1; 1]; Vh(:,1) reshape(eye(n), n*n, 1); sum_logR zeros(n, 1); TT T * round(1/h) / (N-1); % 实际累计时间修正 for k 1:N-1 xk X(:,k); vk reshape(Vh(:,k), n, n); % Jacobian和系统右端 Jk [-a, a, 0; c-a-xk(3), c, -xk(1); xk(2), xk(1), -b]; fk [a*(xk(2)-xk(1)); (c-a)*xk(1) - xk(1)*xk(3) c*xk(2); xk(1)*xk(2) - b*xk(3)]; % 状态更新 for i 1:n conv 0; L min(k, M-1); for j 1:L conv conv G(j1, i) * X(i, k1-j); end X(i, k1) h^q(i) * fk(i) - conv; end % 变分方程更新 Fv Jk * vk; for i 1:n for m 1:n row (m-1)*n i; conv 0; L min(k, M-1); for j 1:L conv conv G(j1, i) * Vh(row, k1-j); end Vh(row, k1) h^q(i) * Fv(i,m) - conv; end end % QR分解累积 P reshape(Vh(:,k1), n, n); [Q, R] qr(P); s sign(diag(R)); Q Q * diag(s); R diag(s) * R; sum_logR sum_logR log(abs(diag(R))); Vh(:,k1) Q(:); % 用正交基替换扰动方向 end LEs sum_logR / T; fprintf(非相称FO Chen系统LE: (%.4f, %.4f, %.4f)\n, LEs); end调用方式就是直接运行LEs LE_noncommensurate_fo()。主循环里的三个关键参数——记忆长度 M、步长 h、总时长 T——直接影响LE谱的精度。M 控制历史信息保留量h 控制离散化误差T 控制统计平均的收敛程度。改变系统时只需要替换 fk、Jk 和初始条件注意状态维数 n 与 q 的长度一致即可。4. 参数设置、验证与三处必踩的坑4.1 阶次向量q的取值范围分数阶系统的稳定性与 q 的取值有强依赖关系。在0到1之间的阶次是最常见的场景此时系统有记忆衰减特性q_i 大于1时将出现超扩散行为数值格式的收敛性会随之变化。对非相称系统来说不仅要看每个 q_i 是否落在合理区间还要看各阶次间的差是否过大。下表给出常用的设置区间和建议场景q_i 常用区间说明分数阶混沌分析0.85 ~ 0.99混沌现象最丰富LE收敛也较快分数阶超混沌0.90 ~ 0.98需要至少两个正的LE阶次差不宜超过0.1分数阶控制0.5 ~ 0.9兼顾稳定裕度与响应速度退化验证所有 q_i 0.995应与整数阶系统LE近似一致需要特别提醒的是如果 q_i 之间的差超过0.2系统各分量的时间尺度会明显分离可能导致某个方向上的积分严重失真。此时必须缩小步长必要时对不同分量采用不同的积分步数但这样做会让代码复杂度上升一般先检查是否真的需要这么大的阶次差。4.2 步长h与记忆长度M的权衡GL显式格式的误差随步长增大而快速增长。对非相称系统最小阶次 q_min 决定收敛难度因为 h^{q_min} 衰减最慢记忆效应最强。经验做法是让 h 满足h 0.01 / max(1, max(abs(eig(Jk))))也就是根据Jacobian的最大特征值模来限定步长。对Chen系统这类混沌系统h 0.005 是比较稳妥的起点。M 的选择同样关键M round(5/h) 意味着保留最近5秒的历史对 q 0.85 的系统已经足够若 q 低于0.8记忆衰减变慢应把记忆时长提高到10秒以上。另一个容易忽视的是时间长度 T。LE谱需要足够长的统计平均时间才能收敛一般要求 T ≥ 100且至少包含数万个步长。判断是否收敛的方法是把 T 从100增加到200观察LE变化量是否小于1%。4.3 三个典型错误相称代码改参数时必出事把相称代码改成非相称时下面三个错误几乎必然会犯。第一个错误是记忆系数漏更新。相称代码里只有一条 G 向量改成非相称时只更新了状态方程变分方程仍然使用原来那条 G。这样算出的变分矩阵在量级上完全不对LE谱只会在极小步长下碰巧接近真实值。排查方式很简单打印 G(:,1) 和 G(:,2)两列不相等的状态下任何共用G的代码都有嫌疑。第二个错误是QR分解频率过高。每步都做QR分解虽然可行但会让GL记忆的历史值全部被正交化旋转。严格来说GL记忆项中的历史变分矩阵应该保持原始坐标每步强制替换成正交基会引入额外误差。工程上的折中是每 510 步做一次QR分解同时把变分矩阵的历史也做同样的基变换。更极端一点的做法是用短记忆让被旋转的历史只影响最近几步误差可忽略。上面示例代码采用每步QR是因为步长已经足够小如果发现LE收敛值随时间震荡优先改成每10步QR一次。第三个错误是初始瞬态没有排除。无论分数阶还是整数阶系统都需要一段瞬态过渡才能达到吸引子。GL方法自带记忆效应初始值的影响会持续很久。稳妥做法是从第50秒之后才开始累积 Q*R 的对数分解结果在此之前只积分不累积。把主循环改成if k * h 50 sum_logR sum_logR log(abs(diag(R))); end这样得到的LE谱不受初始条件影响重复计算的一致性也更好。4.4 用一个退化系统验证代码正确性拿到任何一套LE代码第一步不是直接算目标系统而是做退化验证。把 q 全部设为相同的0.995非相称系统退化为相称系统再把 h 缩小到0.001分数阶系统应当逼近整数阶Chen系统的动力学。此时参考LE值可以用MATLAB自带的整数阶同步计算或经典文献值对比通常为 (1.0, 0, -34) 这个量级。如果退化后LE谱与参考值偏差超过0.1说明GL系数递推或变分方程更新代码有误应当先修好再回到非相称场景。% 退化验证示例 q_test [0.995; 0.995; 0.995]; % 运行后应得到接近整数阶Chen系统的LE三要素5. 用LE谱判定非相称FO超混沌一个可复现的验证流程5.1 构造非相称阶超混沌系统并计算LE验证代码正确性之后就可以扩展到非相称超混沌系统。一个常见做法是把四维超混沌Lorenz系统改造成非相称阶次结构例如阶次向量取 q [0.95, 0.97, 0.99, 0.92]四个分量各自拥有不同的记忆核。计算LE谱时把本博客第三节的函数改一下 fk 和 Jk 的维度即可核心的GL系数初始化、变分方程逐行卷积和QR累积逻辑完全不需要动。下面给出Jacobian的一般计算思路% 四维非相称系统的Jacobian示例以超混沌Lorenz为例 % 系统右端 f 对 x 求偏导得到4×4矩阵 Jk zeros(4); Jk(1,1) -a; Jk(1,2) a; Jk(2,1) c - x(3); Jk(2,3) -x(1); Jk(3,1) x(2); Jk(3,2) x(1); Jk(3,3) -b; % 第四行和其余耦合项按实际系统填写算出的LE谱中正的指数个数就是判定混沌的重要依据一个正LE意味着混沌两个正LE意味着超混沌。若计算结果显示 λ1 0、λ2 0、λ3 ≈ 0、λ4 0则说明该系统在非相称阶次设置下确实处于超混沌状态。5.2 与分岔图互相交叉验证LE谱单算一次还不够可靠特别是非相称系统的记忆效应会让数值结果对初值敏感。建议把LE谱和分岔图结合使用固定除某个阶次之外的所有参数让该阶次从0.85扫描到0.99每步重新计算LE同时记录状态变量的极值绘制分岔图。当分岔图从周期窗口进入混沌区域时LE谱应当同步出现由负转正的跳变当分岔图中吸引子失去混沌特征时最大LE应当回落到0附近。这种交叉验证能让论文中的混沌判据更有说服力也能反过来暴露代码在特定阶次附近的数值不稳定问题。Q*R累积和分岔扫描的耗时问题也不能忽略。GL方法每步都要做卷积求和扫描30个阶次点时建议先用短记忆 M round(2/h) 做初筛找到混沌区间后再用 M round(5/h) 精算。这样既控制总运行时长又不牺牲最终结果的可靠性。对非相称系统这种“每个分量带着自己的记忆核”的结构能做到稳定重现、交叉验证一致LE代码才算真正落地可用。本文还有配套的精品资源点击获取