企业官网建设流程全解析

把分数阶系统（FOS）的Lyapunov指数谱算出来，这事看起来小众，实际上坑特别多。先说结论：如果直接把整数阶系统的Lyapunov指数计算代码拿过来，套上分数阶Lorenz、4D分数阶Chen或者分数阶Duffing，大概率会得到一组“看起来像那么回事但根本不能用”的数字。原因在分数阶导数的记忆效应，它让切空间演化不再是一步一Jacobian那么简单。这篇内容我会把分数阶混沌系统的Lyapunov指数谱计算思路、三个典型系统的模型形式、以及一套在Matlab里可以直接跑的代码框架完整拆开讲。写给正在做FOS混沌分析、系统辨识、控制器稳定性验证的工程师和研究生，也适合刚接触分数阶非线性动力学但不想只看理论公式的人。

1. 内容整体设计与思路拆解

1.1 分数阶系统的Lyapunov指数谱为什么难算

先统一口径。这里说的分数阶系统，一般指用Caputo分数阶导数或Riemann-Liouville分数阶导数描述的非线性常微分方程组：

[ {}^C_0D_t^q x(t)=f(x(t),t), \quad 0<q\le 1 ]

从物理直觉上，分数阶系统比整数阶系统多了一个“历史记忆通道”。当前时刻的状态变化率不仅由当前状态决定，还受到整个时间历史的影响。这个特点让分数阶混沌系统在实际建模中很受欢迎，因为很多真实物理过程、黏弹性材料、生物组织、电化学系统都表现出类似的记忆和遗传特性。

但计算Lyapunov指数谱的时候，这个记忆特性恰恰成了大麻烦。经典Wolf算法和QR分解法都建立在一个前提上：系统是有限维常微分方程，切空间演化可以用Jacobian矩阵的瞬时值来推进。分数阶系统虽然形式上也是有限维状态，但本质上它是一个包含全局记忆的积分微分系统，状态演化方程里的卷积核不能随便丢掉。把整数阶Wolf算法直接搬过来，相当于强行把记忆效应切掉，得到的结果自然不靠谱。

我见过不止一个同学把整数阶Lorenz系统的谱计算代码改两个参数就套到分数阶Lorenz上，结果正Lyapunov指数飘到了三点几，比系统维度还离谱。这类问题不是代码bug，是算法前提不对。

1.2 为什么选择“分数阶预测-校正积分+修正Wolf算法”这个组合

目前文献里做FOS的Lyapunov指数谱，主流思路大体分三类：

第一类：用预测-校正法或分数阶Adams法对原系统做数值积分，然后把连续轨迹按照固定步长离散化，在离散时间序列上套用Wolf算法或Sano-Sawada算法。这种思路的好处是代码改动小，坏处是容易丢失分数阶系统的切空间演化信息，实际效果依赖步长和记忆长度。

第二类：把分数阶系统扩展成高阶整数阶系统，用频域展开来近似分数阶算子，再套整数阶Lyapunov指数算法。这类方法实现复杂，而且频域近似本身会引入相位误差，优点不明显。

第三类：在求解原系统的同时，把“扰动传播”也放在同一个分数阶求解器里推进，每隔固定时间做一次重正交化和QR分解。这类方法理论上更接近分数阶系统的变分结构，实践中也比较稳定。

我采用的策略是第三类的工程化版本：先用分数阶预测-校正法（PECE框架）对原系统和扰动基向量做同步推进，每个归一化间隔用QR分解提取谱值。这个方案既避免了直接套整数阶算法的结构错误，又不会让代码复杂到无法维护。

下面先快速过一遍必要的数学背景，直接给出能写进代码的表达形式。

2. 分数阶模型与Lyapunov指数谱的数学基础

2.1 Caputo导数与Grunwald-Letnikov离散

对于(0<q<1)，Caputo分数阶导数定义为：

[ {}^C_0D_t^q x(t)=\frac{1}{\Gamma(1-q)}\int_0^t \frac{x'(\tau)}{(t-\tau)^q}d\tau ]

这个定义的好处是初值条件可以直接用整数阶导数的初值(x(0))表示，工程应用里非常友好。进行数值计算时，不直接处理积分核，而是用Grunwald-Letnikov的离散形式来近似：

[ D^q x(t_n) \approx h^{-q}\sum_{j=0}^{n} \omega_j^{(q)} x(t_{n-j}) ]

其中权重系数(\omega_j^{(q)})满足递推关系：

[ \omega_0^{(q)}=1,\qquad \omega_j^{(q)}=\left(1-\frac{q+1}{j}\right)\omega_{j-1}^{(q)} ]

这个递推关系是实际写代码最常用的，因为它避免直接组合数计算，稳定且快。所有需要“记忆历史”的项，都被这些权重系数承载了。

2.2 Lyapunov指数谱的物理含义

Lyapunov指数谱描述的是相空间中相邻轨道的指数分离速率。设系统状态维度为(n)，在轨道附近取一套正交扰动基向量，经过足够长的时间，各方向的平均指数增长率就是Lyapunov指数：

[ \lambda_i=\lim_{t\to\infty}\frac{1}{t}\ln\frac{| \delta x_i(t)|}{|\delta x_i(0)|} ]

混沌系统的特征是至少存在一个正的Lyapunov指数；超混沌系统至少有两个正的Lyapunov指数。判断一个分数阶系统是周期、拟周期还是混沌，最稳妥的办法不是只靠相图，而是看完整谱的符号分布。

2.3 分数阶切空间推演时需要注意的点

对整数阶系统，扰动向量(\delta x(t))可以通过Jacobian矩阵线性化直接演化。分数阶系统同样可以在形式上写变分方程：

[ D^q \delta x = J(x(t)),\delta x,\qquad J(x)=\frac{\partial f}{\partial x} ]

由于Caputo导数对线性叠加是成立的，这个变分方程在形式上与原系统保持一致阶数。真正难的是：原系统的状态(x(t))本身带有记忆，而Jacobian矩阵随时间变化，扰动基向量的演化因此不仅依赖于当前时刻的Jacobian，还隐含地依赖过去所有时刻的轨道信息。也就是说，求解扰动基向量时，必须原系统和变分方程同步推进，不能先独立算完轨道，再用瞬时Jacobian去推扰动。

这个细节直接决定谱的精度，也是后面代码架构里的核心约束。

3. 三个典型系统的模型与参数选择

3.1 分数阶Lorenz系统

Lorenz系统是混沌理论里最常用的基准系统，分数阶形式为：

[ \begin{aligned} D^{q_1}x &= \sigma(y-x)\ D^{q_2}y &= x(\rho-z)-y\ D^{q_3}z &= xy-\beta z \end{aligned} ]

经典参数取(\sigma=10,\ \rho=28,\ \beta=8/3)，这个参数下整数阶Lorenz系统处于混沌状态。分数阶版本通常把三个阶数取相同值(q\in(0.9,1])，当(q)略小于1时依然能观察到混沌。

实际验证时，我习惯从(q=0.995)开始试，逐个往下扫到(q=0.95)左右。这样做的好处是便于排查程序问题：如果(q)接近1时算出来的谱与整数阶Lorenz标准值对不上，说明求解器或者谱估计算法有bug，不用等到分数阶记忆效应显著了才去调试。

3.2 4D分数阶Chen超混沌系统

Chen系统相比Lorenz系统多一个非对称项，分数阶形式为：

[ \begin{aligned} D^{q_1}x &= a(y-x)\ D^{q_2}y &= (c-a)x-xz+cy\ D^{q_3}z &= xy-bz \end{aligned} ]

为了得到超混沌行为，需要在Chen系统的基础上扩展一个状态变量(w)，构成四维系统。这里采用一种常用的反馈扩展形式：

[ \begin{aligned} D^{q_1}x &= a(y-x)\ D^{q_2}y &= dx-xz+cy-w\ D^{q_3}z &= xy-bz\ D^{q_4}w &= yz+rw \end{aligned} ]

典型参数为(a=35,\ b=3,\ c=12,\ d=7)，反馈参数(r)作为超混沌的可调因子。(r)的取值会显著影响谱的符号分布：偏小时系统处于混沌但只有单一正指数，偏大时系统可能出现两个正Lyapunov指数，进入超混沌状态。

初值一般取(x(0)=[1,1,1,1]^T)附近即可，混沌系统对初值不敏感于最终谱估计，但瞬态过程中的初值会影响需要烧掉的时长。

3.3 分数阶Duffing振荡器

Duffing振荡器是典型的非自治系统，也是机械振动分析里出现频率最高的混沌模型。整数阶标准形式为：

[ \ddot{x}+\delta\dot{x}+x^3-x=\gamma\cos(\omega t) ]

写成分数阶状态空间形式需要引入(y=\dot{x})，此时把二阶导数的记忆效应分配到两个一阶分数阶方程里，更常见的工程写法是把原状态方程直接推广为：

[ \begin{aligned} D^{q_1}x &= y\ D^{q_2}y &= -\delta y + x - x^3 + \gamma\cos(\omega t) \end{aligned} ]

由于时间(t)以显式方式进入方程，这是一个非自治分数阶系统。计算Lyapunov指数谱时，对非自治系统可以先把(t)也扩展成一个新状态(D^q s=1)，也可以直接保留显式时间项，二者结果一致，但写代码时需要注意求解器每一步传入的当前时间必须一致。

典型的混沌参数取(\delta=0.05,\ \gamma=0.3,\ \omega=1.2)。这个参数组合下，整数阶Duffing系统已经表现出双阱势能间的跳跃混沌行为，分数阶略调低后会看到失稳和周期窗口交替出现。

4. Matlab代码框架与核心实现

4.1 主函数：分数阶Lyapunov指数谱估计

下面给出的是我自己常用的计算框架，统一适用于上述三类系统。代码核心思路：使用分数阶预测-校正法推进状态，同时利用切空间变分方程推进扰动基向量，每隔固定时间间隔做一次QR分解。

function [lambda, LEhist] = fos_lyapunov_spectrum(fun, q, x0, opts) % fos_lyapunov_spectrum 计算分数阶系统的Lyapunov指数谱 % 输入fun为系统右端函数，格式 @(t, x) % q为分数阶阶数，0<q<=1，可传标量或与状态同维向量 % opts为求解参数，包含积分步长h、总时长T、瞬态时间tran和归一化间隔teval % % 输出lambda为最终Lyapunov指数谱 % LEhist为每步归一化后记录的谱收敛历史 n = numel(x0); x = x0(:); % 默认参数 if ~isfield(opts, 'h'), opts.h = 0.005; end if ~isfield(opts, 'T'), opts.T = 600; end if ~isfield(opts, 'tran'), opts.tran = 50; end if ~isfield(opts, 'teval'), opts.teval = 1.0; end h = opts.h; T = opts.T; tran = opts.tran; % 初始单位正交基向量 Q0 = eye(n); % 记录当前状态和切基向量 xcur = x; Qcur = Q0; % 累计对数增长 cumlog = zeros(n, 1); LEhist = []; tcur = 0; lastsave = 0; while tcur < T if tcur < tran % 瞬态阶段：只推进状态，不估计谱 xcur = fos_pce_step(fun, q, xcur, tcur, h); else % 谱估计阶段：状态与切基向量同步推进 tau_now = tcur - tran; % 先推进中心轨道一个步长 xnext = fos_pce_step(fun, q, xcur, tcur, h); % 对每个切向量分别推进变分方程 Qnew = zeros(n, n); for k = 1:n Qnew(:, k) = fos_pce_variation_step(fun, q, xcur, Qcur(:, k), tcur, h); end % QR分解实现Gram-Schmidt正交化 [Qu, Ru] = qr(Qnew, 0); % 处理符号，保证谱收敛稳定 sg = sign(diag(Ru)); sg(sg == 0) = 1; Qu = Qu .* sg'; Ru = Ru .* sg'; rd = diag(Ru); rd(rd < 1e-30 & rd > -1e-30) = 1e-30; cumlog = cumlog + log(abs(rd)); % 更新状态和切基向量 xcur = xnext; Qcur = Qu; if tcur - lastsave >= opts.teval lambda_est = cumlog / (tcur - tran); LEhist(:, end+1) = lambda_est; %#ok<AGROW> lastsave = tcur; end end tcur = tcur + h; end % 最终谱 lambda = cumlog / (T - tran); end

上面的主函数里我留了两个子函数没有展开：fos_pce_step和fos_pce_variation_step。它们分别负责原系统和变分方程的分数阶预测-校正推进。

4.2 子函数：分数阶预测-校正单步推进

单步推进实现的是分数阶PECE方法。假设所有阶数相同为(q)，使用历史状态数据来构造预测值。

function xnext = fos_pce_step(fun, q, xcur, tcur, h) % 单步预测-校正推进 % 为了保持代码可读性，这里采用简化记忆实现： % 将分数阶导数等效为历史权重积分，通过梯形型校正修正 % 取当前点作为初值预测 tpred = tcur + h; fcur = fun(tcur, xcur); % 预测值采用一阶欧拉型近似，结合阶数修正 xpred = xcur + (h^q) * fcur; % 用预测点重新计算右端函数 fpred = fun(tpred, xpred); % 校正公式：处理分数阶积分提高精度 ratio = h^q / gamma(q + 2); xnext = xcur + ratio * (fpred + (q + 1) * fcur); end

这个子函数是简化版，适合理解算法逻辑。如果做高精度的正式研究，建议把完整历史记忆项加入，实现完整的PECE分数阶Adams方法。下面留了扩展位置。

function vnext = fos_pce_variation_step(fun, q, xcur, vcur, tcur, h) % 变分方程的预测-校正推进 % 在中心轨道xcur附近求Jacobian并推进扰动向量vcur % 这里用中心差分近似Jacobian，保证与系统右端函数解耦 n = numel(xcur); epsilon = 1e-6; % 构建当前时刻Jacobian矩阵 J = zeros(n, n); for j = 1:n xp = xcur; xm = xcur; xp(j) = xp(j) + epsilon; xm(j) = xm(j) - epsilon; fp = fun(tcur, xp); fm = fun(tcur, xm); J(:, j) = (fp - fm) / (2 * epsilon); end % 对扰动向量做一阶预测-校正 fv = J * vcur; vpred = vcur + (h^q) * fv; % 用预测点的Jacobian再做一次校正 fpv = J * vpred; ratio = h^q / gamma(q + 2); vnext = vcur + ratio * (fpv + (q + 1) * fv); end

这段代码把Jacobian用中心差分近似了，好处是模型函数fun不用单独提供解析Jacobian，换系统参数时工作量小。代价是计算量增加，对四维系统影响不大。

4.3 三个系统的模型函数

模型函数格式必须满足统一签名xdot = fun(t, x)。

分数阶Lorenz：

function dx = lorenz_fos(t, x) sigma = 10; rho = 28; beta = 8/3; dx = zeros(3, 1); dx(1) = sigma * (x(2) - x(1)); dx(2) = x(1) * (rho - x(3)) - x(2); dx(3) = x(1) * x(2) - beta * x(3); end

4D分数阶Chen：

function dx = chen4d_fos(t, x) a = 35; b = 3; c = 12; d = 7; r = 0.2; dx = zeros(4, 1); dx(1) = a * (x(2) - x(1)); dx(2) = d * x(1) - x(1) * x(3) + c * x(2) - x(4); dx(3) = x(1) * x(2) - b * x(3); dx(4) = x(2) * x(3) + r * x(4); end

分数阶Duffing：

function dx = duffing_fos(t, x) delta = 0.05; gamma_amp = 0.3; omega = 1.2; dx = zeros(2, 1); dx(1) = x(2); dx(2) = -delta * x(2) + x(1) - x(1)^3 + gamma_amp * cos(omega * t); end

4.4 一键运行脚本

% 示例：计算三个系统的Lyapunov指数谱 clc; clear; close all; % 系统1：分数阶Lorenz opts.h = 0.005; opts.T = 400; opts.tran = 30; opts.teval = 1.0; q = 0.995; x0 = [20; 20; 20]; [lambda_lorenz, hist_lorenz] = fos_lyapunov_spectrum(@lorenz_fos, q, x0, opts); fprintf('分数阶Lorenz系统 Lyapunov指数谱：\n'); disp(lambda_lorenz); figure; plot(linspace(opts.tran, opts.T, size(hist_lorenz, 2)), hist_lorenz', 'LineWidth', 1.2); title('分数阶Lorenz系统Lyapunov指数谱收敛过程'); xlabel('t'); ylabel('Lyapunov指数'); grid on;

其余两个系统只需要替换函数句柄和维度参数。四维系统改为：

q4 = 0.98; x0_4 = [1; 1; 1; 1]; [lambda_chen, hist_chen] = fos_lyapunov_spectrum(@chen4d_fos, q4, x0_4, opts);

Duffing系统因为是非自治的，积分参数里需要额外关注时间步长不要设置太大，否则高频外激励(cos(\omega t))的分辨率不够。

5. 典型结果分析与参数敏感性

5.1 参考输出与谱值解释

以我常用的参考参数为例，三套系统的典型输出如下，具体数值会因为求解器精度和记忆窗口长度略有浮动：

对这个表格需要强调几点：

第一，三维系统的Lyapunov指数谱有三个值，非自治系统通过时间扩展后也可以看到第二个零指数，对应时间方向。若代码里没有显式扩展时间状态，Duffing的谱中会出现一个接近0但不严格为0的值，这是数值现象，不用紧张。

第二，4D Chen系统出现两个正Lyapunov指数是超混沌的关键证据。如果r选得不对，第二个指数会变成负值，说明系统只是普通混沌，没有形成超混沌。

第三，分数阶Lorenz在q逐渐下降时，第一个Lyapunov指数从接近整数阶系统的0.9逐渐下降，并在某个临界阶数附近变为负值或徘徊在零附近，对应混沌消失、进入周期或稳定状态。这个“临界阶数”正是很多分数阶系统同步控制论文里反复讨论的对象。

5.2 为什么谱随阶数变化如此敏感

分数阶阶数q直接影响记忆核的衰减速度。q越接近1，记忆效应越弱，系统行为越接近整数阶；q越小，记忆效应越强，系统耗散和“惯性”特性越明显，原来的混沌吸引子可能收缩成极限环甚至稳定不动点。

这也解释了为什么在分数阶混沌系统里，不能像整数阶系统那样仅凭一组参数就断言“系统是混沌的”。必须说明是在哪个阶数下得到的结论。做谱估计时同样要记录q、h、积分总时长，否则论文里的数据无法复现。

6. 常见问题与排查技巧实录

6.1 计算出的谱跳来跳去不收敛

这是最常见的问题。Lyapunov指数是时间平均量，理论上需要足够长的积分时间才能收敛。如果瞬态烧得不够，或者归一化间隔太小，谱值会在正负之间来回振荡。

我的排查步骤：

1. 先看状态变量时间序列，确认已经进入稳态吸引子，再用肉眼观察相图是否稳定在吸引子上。
2. 把tran设大一些，比如从瞬态50增加到100。
3. 记录谱收敛曲线，如果曲线缓慢但单调趋于水平，说明积分总时长不够；如果曲线一直在抖动，优先怀疑步长太大。

6.2 四维系统出现了两个相同的正指数

这通常不是真正的超混沌，而是QR分解符号处理或正交化异常导致的假象。检查一下qr(Qnew, 0)得到的R矩阵对角线是否有接近零甚至负值。严格数学上R的对角线应为正，但在数值实现中符号需要归一化，否则对数累积会出现连续加减号的混乱。

我习惯在每次QR分解后加一步符号修正，也就是代码里的：

sg = sign(diag(Ru)); sg(sg == 0) = 1; Qu = Qu .* sg'; Ru = Ru .* sg';

这一步很容易被忽略，但作用很大。

6.3 分数阶Lorenz在q接近1时与整数阶结果完全对不上

如果差距在5%以内，可能是数值误差；如果差距大到完全不在一个量级，要检查PECE子函数是否正确实现了历史权重。

分数阶Adams类的预测-校正公式不像整数阶RK那样只依赖当前状态，它需要保留所有历史点的右端函数值。完整实现里会包含这样的累加项：

[ x(t_{n+1})=x(0)+\frac{h^q}{\Gamma(q)}\sum_{j=0}^{n} b_{j,n+1} f(t_j,x_j)+\frac{h^q}{\Gamma(q+1)}f(t_{n+1},x^p_{n+1}) ]

如果不加历史项，只靠当前点做欧拉型预测，那本质上还是整数阶结构，只是把步长换成了(h^q)，大规模计算时误差会累积。

6.4 步长h如何确定

分数阶系统的积分步长对结果的影响比整数阶更敏感，因为历史记忆点数量与(t/h)成正比。h太大，历史分辨率不足，记忆效应被抹平；h太小，总步数暴涨，计算量和舍入误差同步上升。

经验上分三步做收敛性验证：先用(h=0.01)算一组谱，再把h缩小到0.005，如果两次结果差异小于工程可接受范围，说明步长够了；如果不一致，继续缩小。实际计算中，Lorenz这类三维系统用0.005基本够用，但分数阶Duffing因为外激励频率较高，最好从0.002开始试。

6.5 关于“零指数”的判定

理论上，连续自治系统的极限集若存在连续时间流，其Lyapunov指数谱中至少有一个零指数对应流方向。但在分数阶系统里，零指数的数值表现经常不是0，而是(10^{-3})量级的正数或负数。不要急着把这当成混沌判据，要结合相图结构、功率谱或Poincare截面一起判断。

我个人的判断标准是：如果算出来的最小正指数与最大指数相比小于千分之一量级，并且它随着积分时间增加而缓慢趋近于零，那么基本可以认为是零指数。

这套代码框架后续扩展空间很大，最简单的一个方向是把系统的阶数q从标量改成向量，让Lorenz的三个状态分别取不同阶数，这样就能分析“不同状态不同记忆强度”下的混沌临界现象。另一个方向是把模型函数里的参数也设计成可变输入，通过循环扫描参数，直接画出“参数-Lyapunov指数谱”分岔图，那会比只看相图判断混沌要可靠得多。我自己现在做分数阶混沌同步的控制器设计时，也是拿这套代码先算开环系统的谱，确认当前参数确实处于混沌状态，再去做控制器参数整定。