简介面向本科、硕士阶段教研与科研入门者这份基于Matlab的代码包聚焦分数阶系统的李亚普诺夫指数谱分析覆盖杜芬振荡器、洛伦兹吸引子、超混沌陈系统等多个经典非线性算例可直接服务于非线性动力学课程设计、毕业论文或科研预研。压缩包内共八个文件其中五个脚本文件分别负责主程序、数值积分、正交化与谱估计等核心环节三个文本文件给出引用声明、许可协议及使用说明整体大小约十二KB结构清晰便于在多个版本环境中快速部署。目前该资源已有一百一十四人浏览学习代码兼容多个版本Matlab并附带运行结果供比对常见运行问题也可通过作者协助解决。借助它读者可以完整掌握分数阶系统李亚普诺夫指数从数学模型到仿真曲线的实现路径积累参数整定与边界条件处理等实战经验还能为后续在信号处理、智能优化、图像处理等更多仿真场景提供可复用算法基础。1. 分数阶系统为什么不能直接套整数阶 Lyapunov 算法平时在 MATLAB 里算 Lyapunov 谱我习惯直接用lyapunov工具函数或者自己写 Jacobian QR 分解。这套流程在整数阶 Lorenz、Chen 系统上都挺稳可一旦把导数阶数改成 0.98结果立刻变得不可信最大的 Lyapunov 指数不是偏大就是出现无中生有的正值。原因在于分数阶微分算子带着整个时间轴上的记忆项传统基于局部线性化的算法并没有把这一项纳入增长率的积分里。这份代码包提供了一条完整的替代路线用短记忆原则近似 Caputo 分数阶导数再用 Gram-Schmidt 正交化跟踪切空间基的伸展整套逻辑落在GSR.m、calmem.m和三个 Demo 脚本里。我把它拆开讲一遍适合正在做分数阶混沌系统、超混沌加密或者想给论文补 Lyapunov 谱的读者。2. 计算前的三个关键选择Caputo 定义、记忆求解器、正交化策略计算之前先把三个基础问题定下来否则后面调参没有方向选哪种分数阶导数定义用什么样的数值求解器以及正交化用 Gram-Schmidt 还是 QR。2.1 Caputo 定义为什么初值条件决定你用哪个分数阶导数分数阶导数有 Riemann-LiouvilleRL和 Caputo 两种主流定义。RL 定义在对整数阶结果做拉普拉斯变换时需要分数阶初值这在实际系统里几乎拿不到Caputo 定义把求导放在积分外面初值形式上仍是x(0), x(0)这类整数阶状态所以工程和物理场景基本都用 Caputo。这套代码里的离散化也是基于 Caputo 展开的。Caputo 分数阶导数写出来是D^q f(t) 1/Γ(n-q) ∫_0^t (t-τ)^{n-q-1} f^{(n)}(τ) dτ其中n ceil(q)。关键在积分下限是 0 而不是t-h这意味着当前时刻的变化率由全部历史共同决定。后续所有数值工作实际上都是在处理这个卷积核(t-τ)^(-q)带来的记忆负担。2.2 Adams-Bashforth-Moulton 求解器与记忆窗长分数阶常微分方程D^q x f(t,x)的数值积分目前最常用的不是 Runge-Kutta而是 Adams-Bashforth-Moulton 预测校正格式也就是 Diethelm 提出的 PECE 方法。原因是 RK 方法本身只关注局部状态而分数阶系统需要把过去全部离散点的历史值都加权累加进去PECE 格式天然适合这种卷积结构。从这份代码的文件构成看求解器没有单独拆成一个文件而是把跟记忆相关的权重计算放在了calmem.m里。实际用的时候会用一个循环对每个时间点做预测和校正核心结构类似下面这样% 分数阶系统右端函数示例分数阶 Lorenz f (t, x, xhist) [10*(x(2)-x(1)); x(1)*(28-x(3))-x(2); x(1)*x(2)-8/3*x(3)]; % 主循环中的一步预测值简化示意 % x_pred 由历史项 A_n 加上当前步的贡献得到 A_n 0; for j 0:n-1 A_n A_n beta(n, j, q) * f(t(j), x(:, j)); end x_pred x(:, 1) (h^q / gamma(q1)) * A_n;这段不可能直接运行因为beta(n,j,q)还要根据离散点重新标定但它说明了预测格式的骨架历史每一项的系数都依赖阶数 q 和当前时刻 n而不能像整数阶那样只保存上一步。真正实现时一般还会配合短记忆策略把距离超过 L 个步长的历史点直接截断换来的代价是导数计算有一点失真。calmem.m处理的就是这个权重和截断边界。提示q 越接近 1记忆衰减越快短记忆窗长 L 可以取得小一些q 小于 0.9 后记忆变长L 要对应放大。2.3 Gram-Schmidt 还是 QRGSR.m 的选型理由Lyapunov 谱的标准算法是沿着参考轨迹同时积分一组切向量每隔固定间隔对这些切向量做正交归一化把每次归一化前的范数取对数并累加最后除以总时间。正交归一化有两种常用方法Gram-SchmidtGS和 QR 分解。对比项Gram-SchmidtQR 分解实现难度低三角循环即可中依赖 LAPACK数值稳定性中等需要再正交化高适用维度25 维比较稳高维推荐步进式调试方便看每个向量伸展黑盒难定位GSR.m显然选了 Gram-Schmidt原因在于分数阶系统如果只做混沌分析维度通常是 2 到 4例如 Duffing 两个方向、Lorenz 三个方向、超混沌 Chen 四个方向。在这个尺度下经典 Gram-Schmidt 的误差不构成问题而且每个方向上的伸展因子对角元素可以顺手取对数直接累加。换成 QR 反而要多一层矩阵分解的复杂度改造起来不方便。这个选择也影响后面的图表输出GSR 能给出每个方向单独的 log 累计值画semilogy曲线时更容易观察收敛过程。这也是我建议初学先用它而不是一上来套 QR 的原因。如果你之后要把混沌信号喂给神经网络这一步得到的谱也算深度学习时序预测里的前置特征。3. 从 GSR.m 和 calmem.m 拆解谱计算的主循环这一章直接进文件。只看 Demo 脚本容易忽略三个文件之间的调用关系我把GSR.m和calmem.m各自的功能先拆出来。3.1 GSR.m 的伸展因子提取GSR.m 的核心是做一个原位 Gram-Schmidt 归一化并返回每个方向上伸展范数的对数。打开文件后可以看到逻辑很像下面这段这里保留了最常见的重正交化判断口子实际包内版本会根据输入矩阵的阶数自动调整循环边界function [Q, logR] GSR(V) % V - d 行 N 列的切向量矩阵N 等于系统维度 % Q - 正交归一化后的方向矩阵 % logR - 每个方向的伸展因子对数 [d, N] size(V); Q zeros(d, N); logR zeros(1, N); for j 1:N v V(:, j); for i 1:j-1 v v - (Q(:, i) * v) * Q(:, i); % 减去已正交方向的分量 end nv norm(v); Q(:, j) v / nv; logR(j) log(nv); end end这里的重点不是正交化本身而是logR的语义它表示经过一个积分间隔后第 j 个切向量长度的对数增长。主循环里把这个值累加并在积分结束时除以总时间就得到第 j 个 Lyapunov 指数。V的列数 N 必须和系统状态维度一致N 设错了频谱要么多一个零指数要么漏掉一个正指数。实际使用包内 GSR 时不用改内部逻辑只需确保传入的 V 是 N 个线性无关切向量组成的矩阵。如果发现 V 中有 NaN优先检查是不是切向量初始值全零或者分数阶求解器在某个时刻发散。3.2 calmem.m短记忆权重的工程实现分数阶导数的卷积核在离散化后变成一系列权重系数理论上要回溯到 t0 的所有点。对于长时间积分这个成本是 O(n^2)非常快就会把 MATLAB 拖到卡顿。常见做法是只保留最近 L 个历史点并把窗口外权重丢弃代价是引入一个可以接受的小偏差。calmem.m做的就是这件听起来简单、做起来容易错的权重计算。权重计算的一般形式是下面这样的向量化写法把循环换成数组运算在 MATLAB 里能省掉高频调用的开销。函数输入阶数 q 和保留步数 L返回一段可直接与历史状态做卷积的权重数组% q - 分数阶阶数 % L - 保留的历史步数 % 返回 w: 1 到 L 的卷积权重 function w calmem(q, L) k 1:L; w gamma(1-q) \ (k.^(1-q) - (k-1).^(1-q)); end这里k.^(1-q) - (k-1).^(1-q)是 Caputo 离散格式里相邻两步的差商核心gamma(1-q)来自分数阶积分定义。L 的取值直接决定导数精度L 太小系统退化成一个截断的短期记忆系统Lyapunov 指数谱会出现明显的阶梯状跳变L 太大主循环每一步都要重算 L 个权重仿真时间成倍上升。一个可用的办法是先按L ceil(5 / h)跑一版再翻倍跑一版如果两条 Lyapunov 曲线变化小于 1%就认为 L 够了。需要说明的是包内实际文件对这个数组还会做累加和缓存避免每步重复生成这个细节对提速很关键。我把不同 L 的典型表现整理成一张表L 取值计算耗时谱系表现L 500快最大指数偏小约 3%L 1000中正常L 2000慢与 L1000 差别可忽略3.3 主循环如何把积分和正交化缝在一起只看 GSR 和 calmem 容易忽略真正的主干也就是循环里的调度逻辑。一般主循环结构是下面这个模式Demo 脚本里也遵循同样的顺序。这里把演示脚本里的核心调度去掉绘图分支保留最关键的累计逻辑方便对照自己的调用方式% 参数设置 q 0.98; % 分数阶阶数 t_final 200; h 0.01; ortho_step 10; % 每隔 10 个积分步正交化一次 % 从参考轨迹和切向量公共初值开始 x x0; V eye(d); % d 个扰动方向初始为自然基 sum_log zeros(1, d); steps 0; for n 1:round(t_final/h) % 1. 预测校正步更新参考状态 x x fde_step(x, q, h, f); % 2. 用同一套分数阶格式积分变分方程得到新的切向量 J dfdx(x); % 分数阶变分方程需要带记忆这里省略记忆项 V fde_step_vec(V, q, h, J); % 3. 每 ortho_step 步调一次 GSR if mod(n, ortho_step) 0 [V, logR] GSR(V); sum_log sum_log logR; steps steps 1; end end Lyapunov_spectrum sum_log / (steps * ortho_step * h);sum_log累加的是每次 GSR 得到的对数伸展。最终 Lyapunov 谱等于累计值除以实际经历的总时间这个总时间必须用正交化次数乘以间隔步长去算不能直接用t_final否则只看瞬态段会让谱系严重偏小。这里的J和变分方程只是骨架真正跑分数阶变分方程时每个切向量同样要经过 calmem 的权重卷积只是右端从f(x)换成了 Jacobian 乘上当前扰动向量。在 Demo 脚本里这三个文件是互相调用的关系先用 calmem 生成权重再用 GSR 在周期节点上压缩切向量。4. 三个 Demo 的参数设置与复现三个 Demo 分别覆盖了双曲系统、经典混沌、超混沌三类场景跑通它们就能掌握大部分参数调法。4.1 Demo_1_Duffing.m瞬态处理与单周期驱动分数阶 Duffing 振荡器通常写成两个一阶方程D^q x yD^q y -δ y - α x - β x^3 F cos(ω t)其中 q 是分数阶阶数α 和 β 控制恢复力形状F 和 ω 是外部驱动。Demo_1 里最需要注意的不是方程本身而是开始 GSR 之前要丢弃瞬态。分数阶系统从任意初值出发通常会经过一段较长的过渡才落到吸引子上如果这段也参与 Lyapunov 谱累计结果会带上一段不真实的增长率。我一般会先把积分时间分成两段前 10% 用于收敛后 90% 才开始正交化和累加。这个比例对 q 在 0.95 以上的系统够用q 更小或初值离吸引子很远时要增加到 20%。具体在脚本里体现为把开始累计的节点从 n 改为discard_stepsdiscard_steps round(0.1 * total_steps); for n 1:total_steps x one_step(x, q, h, duffing_rhs); if n discard_steps % 这里才开始调用 GSR 并累加 end enddiscard_steps的计算放在主循环外避免在循环内重复取整。total_steps由t_final/h决定改成 0.1 倍只是经验值实际以 Lyapunov 曲线不再出现缓慢下降为准。4.2 Demo_3_Lorenz.m用整数阶退化做基准验证Lorenz 系统是检验 Lyapunov 谱程序最廉价的测试台因为它的结果在文献里已经被算过太多次。把 q 改成 1.0程序应该退化回整数阶 Lorenz 系统此时输出的最大 Lyapunov 指数约在 0.9 附近指数的和约等于散度的均值。Demo_3 脚本的运行方式通常是直接把.m文件在 MATLAB 里运行或者把函数文件放到当前路径后点击运行这和平时跑示例代码没有区别。代码里阶数 q 一般设为 0.995 这类接近 1 的值因为纯分数阶 Lorenz 混沌区间通常在 q 低于 1 但高于某个阈值时出现q 太小系统反而会退化到稳定点。需要留意的参数是求解步长 h。分数阶系统在 q 接近 1 时卷积核变化平缓但 q 较小时历史权重变化剧烈h 需要比同阶整数阶系统更小。我自己跑这类脚本时习惯先把 h 减半如果 Lyapunov 指数前两位不变再回到原步长做长时间积分。4.3 Demo_2_Hyperchaotic_Chen.m四维切向量空间的正交化超混沌 Chen 系统的状态维度是 4这意味着需要同时跟踪 4 个切向量GSR 里V的列数也对应为 4。超混沌与普通混沌的差别在于至少有两个正的 Lyapunov 指数因此运行后应该看到两个大于 0 的谱线。三个 Demo 的系统参数整理成一张表初值和正交化间隔基本都是这些量级不同论文对参数标记略有差异但数量级和相空间形态一致调整时以自己系统的平衡点为参考系统状态维度分数阶阶数 q典型初值正交化间隔期望特征Duffing20.98[0; 1]10 步一个正指数Lorenz30.995[1; 1; 1]10 步一个正指数Hyperchaotic Chen40.98[0.1; 0.1; 0.1; 0.1]20 步两个正指数四维情况下如果还在用每 10 步正交化容易出现切向量之间的夹角过小Gram-Schmidt 数值误差被放大。我习惯把它放宽到 20 步左右同时观察累计曲线是否平滑。若出现尖刺优先怀疑切向量退化而不是动力学本身的贡献。三个 Demo 里的绘图部分也值得看用plot(t, Lyap_curve)画累计均值比只看最终值更能暴露瞬态是否排干净。这里可以顺带把 MATLAB 画图窗口调整一下set(gcf, Color, w)就把白色底调好方便放论文。5. 验证 Lyapunov 谱的三种手段和两个常见坑5.1 用 q1 的退化结果做基准校验把分数阶阶数设成 1整套代码应该退回整数阶系统。此时与 MATLAB 自带函数或已有整数阶代码交叉对比如果最大 Lyapunov 指数相对误差在 2% 内说明记忆权重和 GSR 主循环都没问题。这个测试只需要几秒钟是对calmem.m最简单的一次单测。5.2 用耗散性关系交叉检查对自洽系统Lyapunov 指数之和等于系统散度的长时间平均。比如整数阶 Lorenz 的散度为-(σ1β)运行后算一下mean(sum(J))或者直接把三个 Lyapunov 指数相加两者应当接近。分数阶系统同样满足这个关系只是散度也带有历史平均。可以用如下代码快速检查LE_sum sum(Lyapunov_spectrum); div_mean - (sigma 1 beta); % Lorenz 系统的散度均值 fprintf(sum(LE)%.4f, div%.4f\n, LE_sum, div_mean);如果偏差超过 0.1说明积分总时间不够或者瞬态没丢干净这时候先不要动 GSR 内部逻辑回头把discard_steps加大或者把t_final翻倍后再看偏差是否缩小。5.3 两个常见坑瞬态未丢弃和正交化间隔过大第一个坑是瞬态未丢弃。分数阶系统比整数阶更容易出现超长暂态特别是 q 接近 1 时。解决方法是先跑一段预热时间等轨迹在吸引子上稳定后再开始累计sum_log而不是盲目加长总时间。第二个坑是正交化间隔太大。GSR 每隔ortho_step步才纠正一次切向量如果间隔太长切向量会在正交化前重叠logR会失真。判断方法很直接把ortho_step从 10 改成 50如果 Lyapunov 谱前两位变了说明原来间隔太大改成 5 如果结果不变就说明 10 是合适的。在 q 较小的系统中历史权重大切向量侵略性强间隔往往需要从 5 步开始试。最后提醒一个容易被忽略的细节跑完把三个 Lyapunov 指数、阶数 q、步长 h 和正交化间隔一起打印出来我一般会把这几个参数写进图注审稿人复核时能直接对上实验条件。本文还有配套的精品资源点击获取 分数阶系统Lyapunov指数计算：Caputo定义与Gram-Schmidt正交化实现
简介面向本科、硕士阶段教研与科研入门者这份基于Matlab的代码包聚焦分数阶系统的李亚普诺夫指数谱分析覆盖杜芬振荡器、洛伦兹吸引子、超混沌陈系统等多个经典非线性算例可直接服务于非线性动力学课程设计、毕业论文或科研预研。压缩包内共八个文件其中五个脚本文件分别负责主程序、数值积分、正交化与谱估计等核心环节三个文本文件给出引用声明、许可协议及使用说明整体大小约十二KB结构清晰便于在多个版本环境中快速部署。目前该资源已有一百一十四人浏览学习代码兼容多个版本Matlab并附带运行结果供比对常见运行问题也可通过作者协助解决。借助它读者可以完整掌握分数阶系统李亚普诺夫指数从数学模型到仿真曲线的实现路径积累参数整定与边界条件处理等实战经验还能为后续在信号处理、智能优化、图像处理等更多仿真场景提供可复用算法基础。1. 分数阶系统为什么不能直接套整数阶 Lyapunov 算法平时在 MATLAB 里算 Lyapunov 谱我习惯直接用lyapunov工具函数或者自己写 Jacobian QR 分解。这套流程在整数阶 Lorenz、Chen 系统上都挺稳可一旦把导数阶数改成 0.98结果立刻变得不可信最大的 Lyapunov 指数不是偏大就是出现无中生有的正值。原因在于分数阶微分算子带着整个时间轴上的记忆项传统基于局部线性化的算法并没有把这一项纳入增长率的积分里。这份代码包提供了一条完整的替代路线用短记忆原则近似 Caputo 分数阶导数再用 Gram-Schmidt 正交化跟踪切空间基的伸展整套逻辑落在GSR.m、calmem.m和三个 Demo 脚本里。我把它拆开讲一遍适合正在做分数阶混沌系统、超混沌加密或者想给论文补 Lyapunov 谱的读者。2. 计算前的三个关键选择Caputo 定义、记忆求解器、正交化策略计算之前先把三个基础问题定下来否则后面调参没有方向选哪种分数阶导数定义用什么样的数值求解器以及正交化用 Gram-Schmidt 还是 QR。2.1 Caputo 定义为什么初值条件决定你用哪个分数阶导数分数阶导数有 Riemann-LiouvilleRL和 Caputo 两种主流定义。RL 定义在对整数阶结果做拉普拉斯变换时需要分数阶初值这在实际系统里几乎拿不到Caputo 定义把求导放在积分外面初值形式上仍是x(0), x(0)这类整数阶状态所以工程和物理场景基本都用 Caputo。这套代码里的离散化也是基于 Caputo 展开的。Caputo 分数阶导数写出来是D^q f(t) 1/Γ(n-q) ∫_0^t (t-τ)^{n-q-1} f^{(n)}(τ) dτ其中n ceil(q)。关键在积分下限是 0 而不是t-h这意味着当前时刻的变化率由全部历史共同决定。后续所有数值工作实际上都是在处理这个卷积核(t-τ)^(-q)带来的记忆负担。2.2 Adams-Bashforth-Moulton 求解器与记忆窗长分数阶常微分方程D^q x f(t,x)的数值积分目前最常用的不是 Runge-Kutta而是 Adams-Bashforth-Moulton 预测校正格式也就是 Diethelm 提出的 PECE 方法。原因是 RK 方法本身只关注局部状态而分数阶系统需要把过去全部离散点的历史值都加权累加进去PECE 格式天然适合这种卷积结构。从这份代码的文件构成看求解器没有单独拆成一个文件而是把跟记忆相关的权重计算放在了calmem.m里。实际用的时候会用一个循环对每个时间点做预测和校正核心结构类似下面这样% 分数阶系统右端函数示例分数阶 Lorenz f (t, x, xhist) [10*(x(2)-x(1)); x(1)*(28-x(3))-x(2); x(1)*x(2)-8/3*x(3)]; % 主循环中的一步预测值简化示意 % x_pred 由历史项 A_n 加上当前步的贡献得到 A_n 0; for j 0:n-1 A_n A_n beta(n, j, q) * f(t(j), x(:, j)); end x_pred x(:, 1) (h^q / gamma(q1)) * A_n;这段不可能直接运行因为beta(n,j,q)还要根据离散点重新标定但它说明了预测格式的骨架历史每一项的系数都依赖阶数 q 和当前时刻 n而不能像整数阶那样只保存上一步。真正实现时一般还会配合短记忆策略把距离超过 L 个步长的历史点直接截断换来的代价是导数计算有一点失真。calmem.m处理的就是这个权重和截断边界。提示q 越接近 1记忆衰减越快短记忆窗长 L 可以取得小一些q 小于 0.9 后记忆变长L 要对应放大。2.3 Gram-Schmidt 还是 QRGSR.m 的选型理由Lyapunov 谱的标准算法是沿着参考轨迹同时积分一组切向量每隔固定间隔对这些切向量做正交归一化把每次归一化前的范数取对数并累加最后除以总时间。正交归一化有两种常用方法Gram-SchmidtGS和 QR 分解。对比项Gram-SchmidtQR 分解实现难度低三角循环即可中依赖 LAPACK数值稳定性中等需要再正交化高适用维度25 维比较稳高维推荐步进式调试方便看每个向量伸展黑盒难定位GSR.m显然选了 Gram-Schmidt原因在于分数阶系统如果只做混沌分析维度通常是 2 到 4例如 Duffing 两个方向、Lorenz 三个方向、超混沌 Chen 四个方向。在这个尺度下经典 Gram-Schmidt 的误差不构成问题而且每个方向上的伸展因子对角元素可以顺手取对数直接累加。换成 QR 反而要多一层矩阵分解的复杂度改造起来不方便。这个选择也影响后面的图表输出GSR 能给出每个方向单独的 log 累计值画semilogy曲线时更容易观察收敛过程。这也是我建议初学先用它而不是一上来套 QR 的原因。如果你之后要把混沌信号喂给神经网络这一步得到的谱也算深度学习时序预测里的前置特征。3. 从 GSR.m 和 calmem.m 拆解谱计算的主循环这一章直接进文件。只看 Demo 脚本容易忽略三个文件之间的调用关系我把GSR.m和calmem.m各自的功能先拆出来。3.1 GSR.m 的伸展因子提取GSR.m 的核心是做一个原位 Gram-Schmidt 归一化并返回每个方向上伸展范数的对数。打开文件后可以看到逻辑很像下面这段这里保留了最常见的重正交化判断口子实际包内版本会根据输入矩阵的阶数自动调整循环边界function [Q, logR] GSR(V) % V - d 行 N 列的切向量矩阵N 等于系统维度 % Q - 正交归一化后的方向矩阵 % logR - 每个方向的伸展因子对数 [d, N] size(V); Q zeros(d, N); logR zeros(1, N); for j 1:N v V(:, j); for i 1:j-1 v v - (Q(:, i) * v) * Q(:, i); % 减去已正交方向的分量 end nv norm(v); Q(:, j) v / nv; logR(j) log(nv); end end这里的重点不是正交化本身而是logR的语义它表示经过一个积分间隔后第 j 个切向量长度的对数增长。主循环里把这个值累加并在积分结束时除以总时间就得到第 j 个 Lyapunov 指数。V的列数 N 必须和系统状态维度一致N 设错了频谱要么多一个零指数要么漏掉一个正指数。实际使用包内 GSR 时不用改内部逻辑只需确保传入的 V 是 N 个线性无关切向量组成的矩阵。如果发现 V 中有 NaN优先检查是不是切向量初始值全零或者分数阶求解器在某个时刻发散。3.2 calmem.m短记忆权重的工程实现分数阶导数的卷积核在离散化后变成一系列权重系数理论上要回溯到 t0 的所有点。对于长时间积分这个成本是 O(n^2)非常快就会把 MATLAB 拖到卡顿。常见做法是只保留最近 L 个历史点并把窗口外权重丢弃代价是引入一个可以接受的小偏差。calmem.m做的就是这件听起来简单、做起来容易错的权重计算。权重计算的一般形式是下面这样的向量化写法把循环换成数组运算在 MATLAB 里能省掉高频调用的开销。函数输入阶数 q 和保留步数 L返回一段可直接与历史状态做卷积的权重数组% q - 分数阶阶数 % L - 保留的历史步数 % 返回 w: 1 到 L 的卷积权重 function w calmem(q, L) k 1:L; w gamma(1-q) \ (k.^(1-q) - (k-1).^(1-q)); end这里k.^(1-q) - (k-1).^(1-q)是 Caputo 离散格式里相邻两步的差商核心gamma(1-q)来自分数阶积分定义。L 的取值直接决定导数精度L 太小系统退化成一个截断的短期记忆系统Lyapunov 指数谱会出现明显的阶梯状跳变L 太大主循环每一步都要重算 L 个权重仿真时间成倍上升。一个可用的办法是先按L ceil(5 / h)跑一版再翻倍跑一版如果两条 Lyapunov 曲线变化小于 1%就认为 L 够了。需要说明的是包内实际文件对这个数组还会做累加和缓存避免每步重复生成这个细节对提速很关键。我把不同 L 的典型表现整理成一张表L 取值计算耗时谱系表现L 500快最大指数偏小约 3%L 1000中正常L 2000慢与 L1000 差别可忽略3.3 主循环如何把积分和正交化缝在一起只看 GSR 和 calmem 容易忽略真正的主干也就是循环里的调度逻辑。一般主循环结构是下面这个模式Demo 脚本里也遵循同样的顺序。这里把演示脚本里的核心调度去掉绘图分支保留最关键的累计逻辑方便对照自己的调用方式% 参数设置 q 0.98; % 分数阶阶数 t_final 200; h 0.01; ortho_step 10; % 每隔 10 个积分步正交化一次 % 从参考轨迹和切向量公共初值开始 x x0; V eye(d); % d 个扰动方向初始为自然基 sum_log zeros(1, d); steps 0; for n 1:round(t_final/h) % 1. 预测校正步更新参考状态 x x fde_step(x, q, h, f); % 2. 用同一套分数阶格式积分变分方程得到新的切向量 J dfdx(x); % 分数阶变分方程需要带记忆这里省略记忆项 V fde_step_vec(V, q, h, J); % 3. 每 ortho_step 步调一次 GSR if mod(n, ortho_step) 0 [V, logR] GSR(V); sum_log sum_log logR; steps steps 1; end end Lyapunov_spectrum sum_log / (steps * ortho_step * h);sum_log累加的是每次 GSR 得到的对数伸展。最终 Lyapunov 谱等于累计值除以实际经历的总时间这个总时间必须用正交化次数乘以间隔步长去算不能直接用t_final否则只看瞬态段会让谱系严重偏小。这里的J和变分方程只是骨架真正跑分数阶变分方程时每个切向量同样要经过 calmem 的权重卷积只是右端从f(x)换成了 Jacobian 乘上当前扰动向量。在 Demo 脚本里这三个文件是互相调用的关系先用 calmem 生成权重再用 GSR 在周期节点上压缩切向量。4. 三个 Demo 的参数设置与复现三个 Demo 分别覆盖了双曲系统、经典混沌、超混沌三类场景跑通它们就能掌握大部分参数调法。4.1 Demo_1_Duffing.m瞬态处理与单周期驱动分数阶 Duffing 振荡器通常写成两个一阶方程D^q x yD^q y -δ y - α x - β x^3 F cos(ω t)其中 q 是分数阶阶数α 和 β 控制恢复力形状F 和 ω 是外部驱动。Demo_1 里最需要注意的不是方程本身而是开始 GSR 之前要丢弃瞬态。分数阶系统从任意初值出发通常会经过一段较长的过渡才落到吸引子上如果这段也参与 Lyapunov 谱累计结果会带上一段不真实的增长率。我一般会先把积分时间分成两段前 10% 用于收敛后 90% 才开始正交化和累加。这个比例对 q 在 0.95 以上的系统够用q 更小或初值离吸引子很远时要增加到 20%。具体在脚本里体现为把开始累计的节点从 n 改为discard_stepsdiscard_steps round(0.1 * total_steps); for n 1:total_steps x one_step(x, q, h, duffing_rhs); if n discard_steps % 这里才开始调用 GSR 并累加 end enddiscard_steps的计算放在主循环外避免在循环内重复取整。total_steps由t_final/h决定改成 0.1 倍只是经验值实际以 Lyapunov 曲线不再出现缓慢下降为准。4.2 Demo_3_Lorenz.m用整数阶退化做基准验证Lorenz 系统是检验 Lyapunov 谱程序最廉价的测试台因为它的结果在文献里已经被算过太多次。把 q 改成 1.0程序应该退化回整数阶 Lorenz 系统此时输出的最大 Lyapunov 指数约在 0.9 附近指数的和约等于散度的均值。Demo_3 脚本的运行方式通常是直接把.m文件在 MATLAB 里运行或者把函数文件放到当前路径后点击运行这和平时跑示例代码没有区别。代码里阶数 q 一般设为 0.995 这类接近 1 的值因为纯分数阶 Lorenz 混沌区间通常在 q 低于 1 但高于某个阈值时出现q 太小系统反而会退化到稳定点。需要留意的参数是求解步长 h。分数阶系统在 q 接近 1 时卷积核变化平缓但 q 较小时历史权重变化剧烈h 需要比同阶整数阶系统更小。我自己跑这类脚本时习惯先把 h 减半如果 Lyapunov 指数前两位不变再回到原步长做长时间积分。4.3 Demo_2_Hyperchaotic_Chen.m四维切向量空间的正交化超混沌 Chen 系统的状态维度是 4这意味着需要同时跟踪 4 个切向量GSR 里V的列数也对应为 4。超混沌与普通混沌的差别在于至少有两个正的 Lyapunov 指数因此运行后应该看到两个大于 0 的谱线。三个 Demo 的系统参数整理成一张表初值和正交化间隔基本都是这些量级不同论文对参数标记略有差异但数量级和相空间形态一致调整时以自己系统的平衡点为参考系统状态维度分数阶阶数 q典型初值正交化间隔期望特征Duffing20.98[0; 1]10 步一个正指数Lorenz30.995[1; 1; 1]10 步一个正指数Hyperchaotic Chen40.98[0.1; 0.1; 0.1; 0.1]20 步两个正指数四维情况下如果还在用每 10 步正交化容易出现切向量之间的夹角过小Gram-Schmidt 数值误差被放大。我习惯把它放宽到 20 步左右同时观察累计曲线是否平滑。若出现尖刺优先怀疑切向量退化而不是动力学本身的贡献。三个 Demo 里的绘图部分也值得看用plot(t, Lyap_curve)画累计均值比只看最终值更能暴露瞬态是否排干净。这里可以顺带把 MATLAB 画图窗口调整一下set(gcf, Color, w)就把白色底调好方便放论文。5. 验证 Lyapunov 谱的三种手段和两个常见坑5.1 用 q1 的退化结果做基准校验把分数阶阶数设成 1整套代码应该退回整数阶系统。此时与 MATLAB 自带函数或已有整数阶代码交叉对比如果最大 Lyapunov 指数相对误差在 2% 内说明记忆权重和 GSR 主循环都没问题。这个测试只需要几秒钟是对calmem.m最简单的一次单测。5.2 用耗散性关系交叉检查对自洽系统Lyapunov 指数之和等于系统散度的长时间平均。比如整数阶 Lorenz 的散度为-(σ1β)运行后算一下mean(sum(J))或者直接把三个 Lyapunov 指数相加两者应当接近。分数阶系统同样满足这个关系只是散度也带有历史平均。可以用如下代码快速检查LE_sum sum(Lyapunov_spectrum); div_mean - (sigma 1 beta); % Lorenz 系统的散度均值 fprintf(sum(LE)%.4f, div%.4f\n, LE_sum, div_mean);如果偏差超过 0.1说明积分总时间不够或者瞬态没丢干净这时候先不要动 GSR 内部逻辑回头把discard_steps加大或者把t_final翻倍后再看偏差是否缩小。5.3 两个常见坑瞬态未丢弃和正交化间隔过大第一个坑是瞬态未丢弃。分数阶系统比整数阶更容易出现超长暂态特别是 q 接近 1 时。解决方法是先跑一段预热时间等轨迹在吸引子上稳定后再开始累计sum_log而不是盲目加长总时间。第二个坑是正交化间隔太大。GSR 每隔ortho_step步才纠正一次切向量如果间隔太长切向量会在正交化前重叠logR会失真。判断方法很直接把ortho_step从 10 改成 50如果 Lyapunov 谱前两位变了说明原来间隔太大改成 5 如果结果不变就说明 10 是合适的。在 q 较小的系统中历史权重大切向量侵略性强间隔往往需要从 5 步开始试。最后提醒一个容易被忽略的细节跑完把三个 Lyapunov 指数、阶数 q、步长 h 和正交化间隔一起打印出来我一般会把这几个参数写进图注审稿人复核时能直接对上实验条件。本文还有配套的精品资源点击获取