多方法求李指数：从Lorenz系统到实测数据的MATLAB实现

简介：面向混沌系统研究中常用的Lyapunov指数计算需求，压缩包内整理了多种求解思路与配套程序，适合物理、生物、经济等领域的非线性系统分析者，以及希望上手混沌判据的研究生和工程师。压缩包共106个文件、约8.83MB，以Matlab脚本和Simulink模型为主体，辅以figure图形、文本说明和doc文档，覆盖从系统建模、参数设置到指数求解与结果可视化的完整链路。已有359人学习下载，可作为理解混沌理论的入门参考。内容重点讨论直接法、近似法和矩阵幂法三种Lyapunov指数计算思路，并给出洛伦兹系统、Logistic映射等经典混沌实例的演示与代码。通过运行这些程序，读者可以逐步掌握最大Lyapunov指数和平均Lyapunov指数的求解细节，学会利用正指数判别系统混沌行为，并迁移应用到自身研究或工程项目的稳定性分析中。

1. 多方法求李指数：为什么直接按定义算必然翻车

拿到一段实测的混沌电路电压，或者一组 Lorenz 系统的数值解，很多人的第一反应是直接套定义算 Lyapunov 指数。结果算出来的要么是 NaN，要么是正无穷，折腾半天也不知道问题出在哪。这里面的核心矛盾是：相邻轨道的发散速度会随着时间呈指数增长，直接积分数值根本无法稳定。所以工程上才要区分多方法求李指数——已知微分方程用 Jacobian 法，只有离散时间序列则用 Wolf 法或 Rosenstein 法。这篇文章就沿着这条主线，把每种方法的原理、参数设置和坑位讲透。如果你正在跑 lorenz 混沌 matlab 分岔图代码，或者手头有一段待验证的实验数据，这篇笔记正好对症。

2. 基于已知方程的 Jacobian 法：从 Lorenz 系统算出 Lyapunov 指数谱

Lorenz 系统是验证多方法求李指数的最佳对象。它有明确的三维常微分方程，Jacobian 矩阵可以用解析式写出来。所谓“多方法”在已知方程时的常见做法，是对切空间的变分方程做数值积分，再配合 Gram-Schmidt 正交化。

2.1 为什么 Jacobian 方法需要 Gram-Schmidt 正交化

直接对 Jacobian 矩阵连乘，数值上很快会出问题。微分方程解的切向量会指数增长，几个基向量会迅速坍缩到同一个方向。如果不做任何处理，切空间里的三个方向最后只剩下一个有效方向，算出来的结果要么是 NaN，要么是虚假的正无穷。Gram-Schmidt 正交化在这里的作用，是每积分一小步就把基向量重新掰正交，让三个方向保持独立。这样R矩阵对角线的对数才真正对应三个不同方向上的平均发散速率。这也是“多方法求李指数”里最基础、最精确的一条路径。

从数学上看，Lyapunov 指数定义为切空间流映射的奇异值对角元的对数时间平均。难点在于切向量会随 t 指数增长，超出浮点数表示范围。Gram-Schmidt 正交化相当于每隔一小段时间就对切向量做一次归一化，把累计增长量记录在对数里，而不是等它溢出再处理。工程实现上，这比直接求整个状态转移矩阵稳定得多。

2.2 用 Matlab 写 Lorenz 系统的 Jacobian 积分循环

%% 用 Jacobian 法计算 Lorenz 系统 Lyapunov 指数谱 clc; clear; close all; % 系统参数（经典混沌取值） sigma = 10; rho = 28; beta = 8/3; % 数值参数 dt = 0.01; % 积分步长：RK4 要求小于系统最小时间尺度 T = 100; % 总积分时间：足够长才能让指数收敛 N = round(T / dt); % 初始值 x = [1; 1; 1]; % 状态向量 Q = eye(3); % 切空间正交基 lyap_sum = zeros(1, 3); for k = 1:N % 原系统 RK4 推进 k1 = f(x, sigma, rho, beta); k2 = f(x + 0.5*dt*k1, sigma, rho, beta); k3 = f(x + 0.5*dt*k2, sigma, rho, beta); k4 = f(x + dt*k3, sigma, rho, beta); x = x + (dt/6)*(k1 + 2*k2 + 2*k3 + k4); % 变分方程 RK4 推进 J = jac(x, sigma, rho, beta); K1 = J * Q; K2 = jac(x + 0.5*dt*k1, sigma, rho, beta) * (Q + 0.5*dt*K1); K3 = jac(x + 0.5*dt*k2, sigma, rho, beta) * (Q + 0.5*dt*K2); K4 = jac(x + dt*k3, sigma, rho, beta) * (Q + dt*K3); Q = Q + (dt/6)*(K1 + 2*K2 + 2*K3 + K4); % Gram-Schmidt 正交化，Q 会被拆成正交阵 q 和上三角阵 r [Q, R] = qr(Q); lyap_sum = lyap_sum + log(abs(diag(R))); % 把 Q 的对角线符号吸收进来，保证后续迭代的数值稳定性 Q = Q * diag(sign(diag(R))); end lyap_exp = lyap_sum / (N * dt); fprintf('Lyapunov 指数谱： %.4f, %.4f, %.4f\n', lyap_exp); fprintf('最大 Lyapunov 指数 LLE = %.4f\n', lyap_exp(1)); function dx = f(x, sigma, rho, beta) dx = [sigma*(x(2)-x(1)); x(1)*(rho-x(3))-x(2); x(1)*x(2)-beta*x(3)]; end function J = jac(x, sigma, rho, beta) J = [-sigma, sigma, 0; rho - x(3), -1, -x(1); x(2), x(1), -beta]; end

这段代码里最关键的是每次 RK4 推进一步之后调用qr(Q)。R的对角线元素表示当前正交基各个方向上模长的增长倍数，取对数后累加，再除以总时间N*dt，就得到完整的 Lyapunov 指数谱。Lorenz 系统在rho=28时，理论值大约是(0.9, 0, -14.5)，三个指数之和约等于系统 Jacobian 的迹，也就是-(sigma+1+beta)。这是一个非常有效的自检手段。

参数上，dt不能取太大，否则 RK4 会引入虚假的数值耗散；T至少要覆盖 20 个以上的典型轨道周期。对于 Lorenz 系统，dt=0.01、T=100是比较稳的起点。如果算的时候发现结果在0.9附近反复横跳，先别怀疑算法，检查一下T够不够长。

2.3 参数怎么设：步长 dt 与总积分时间 T 对计算结果的影响

我一般会先跑一个快速扫描，把dt分别设为0.001、0.01、0.1，观察 LLE 的变化。结果通常如下表：

如果dt太大，状态会在混沌吸引子上跳跃，切空间基向量跟着失真；如果dt太小，累计步数过多，舍入误差会慢慢累积。实际工程中，我习惯先用较大的dt快速估算出大概量级，再逐步加密，直到 LLE 的变化小于 1%。

还需要注意一个容易混淆的点：如果你的系统模型是离散映射（比如 Logistic 映射或 Henon 映射），那就不存在dt的积分步长概念，每一步迭代本身就是一步。连续系统和离散系统的 Lyapunov 指数定义差一个“时间”量纲，算出来的数值不能直接比。多方法求李指数时，先分清对象是 ODE 还是迭代映射。

3. 基于时间序列的 Wolf 法：重构相空间与核心 Matlab 实现

当我们手里没有方程，只有一组离散采样数据时，多方法求李指数的思路就得换成时间延迟嵌入。Wolf 法是这类方法里历史最久、也最容易踩坑的一种。

3.1 重构相空间：嵌入维数 (m) 和时间延迟 (τ) 怎么定

Wolf 法的第一步是把一维时间序列ts重构到高维相空间。重构参数只有两个：时间延迟tau和嵌入维数dim。tau不能太小，否则ts(i)和ts(i+tau)几乎完全相关，重构出的轨迹会被压扁在对角线附近；tau也不能太大，否则前后状态完全失去关联，噪声主导了整个相空间。我一般先用互信息（Mutual Information）函数取第一个极小值点当作tau，这个方法比单纯看自相关函数更稳。

互信息的计算本身不复杂：把x的取值范围分成若干区间，统计联合分布。当tau从 0 开始增大时，互信息通常会有一个明显的第一个极小值，这就是最佳延迟。自相关函数只能捕捉线性相关，对非线性系统无能为力，所以混沌系统里优先用互信息。

dim的选择通常用虚假最近邻点（FNN）算法或 Cao 方法。简单说，就是从dim=1开始逐步增加维数，观察相空间里哪些原本相距很近的点的距离比例是否急剧变化。对于 Lorenz 这类三维系统，理论上dim=3就够了，但实测含噪数据往往需要取到dim=5甚至更高。

这里有个新手常犯的错误：dim和tau是互相耦合的。先定好tau，再找最小dim才是正序。如果你先用错误的tau去跑 FNN，得到的最小dim也会是错，后面 Wolf 法算出来的 LLE 自然不靠谱。

3.2 Wolf 算法核心步骤与 Matlab 代码

Wolf 算法的核心是“演化 + 替换”：找一个初始最近邻点，让两者在相空间里共同演化一小段距离；然后把基准点前移，在角度限制下重新找替换点，防止基准点跑到吸引子外头。

function LLE = lyap_wolf(ts, tau, dim, evo_step, fs) % ts : 输入时间序列 % tau : 时间延迟（单位：采样点数） % dim : 嵌入维数 % evo_step : 演化步数（通常取 1~3 个 tau） % fs : 采样率，用于计算物理时间 dt = 1 / fs; N = length(ts); M = N - (dim - 1) * tau; % 相空间重构，X 的维度是 M x dim X = zeros(M, dim); for i = 1:dim X(:, i) = ts(1 + (i-1)*tau : M + (i-1)*tau); end % 寻找初始最近邻点 idx_cur = 1; idx_neighbor = 0; min_dist = inf; for j = 2:M if abs(j - idx_cur) > 10 % 移除时间上过近的点 dist = norm(X(idx_cur, :) - X(j, :)); if dist < min_dist && dist > 1e-8 min_dist = dist; idx_neighbor = j; end end end sum_lyap = 0; count = 0; while idx_cur + evo_step <= M && idx_neighbor + evo_step <= M d0 = norm(X(idx_cur, :) - X(idx_neighbor, :)); d1 = norm(X(idx_cur + evo_step, :) - X(idx_neighbor + evo_step, :)); if d0 > 1e-12 && d1 > 1e-12 sum_lyap = sum_lyap + log(d1 / d0); count = count + 1; end idx_cur = idx_cur + evo_step; % 寻找替换点：Wolf 算法的灵魂是角度限制 angle_limit = pi / 6; % 30 度 idx_new_neighbor = 0; best_ratio = inf; for j = 1:M if abs(j - idx_cur) <= evo_step continue; end vec1 = X(idx_cur + evo_step, :) - X(idx_neighbor + evo_step, :); vec2 = X(idx_cur + evo_step, :) - X(j, :); n1 = norm(vec1); n2 = norm(vec2); if n1 < 1e-12 || n2 < 1e-12 continue; end cos_theta = sum(vec1 .* vec2) / (n1 * n2); if cos_theta > cos(angle_limit) && n2 < best_ratio best_ratio = n2; idx_new_neighbor = j; end end if idx_new_neighbor == 0 break; % 找不到合适替换点，提前终止 end idx_neighbor = idx_new_neighbor; end if count > 0 LLE = sum_lyap / (count * evo_step * tau * dt); else LLE = NaN; end end

这段代码执行完，LLE就是预估的最大 Lyapunov 指数。其中evo_step * tau * dt是把“演化步数”折算成真实物理时间，这样不同采样率的数据算出的结果才有可比性。角度限制30 度是为了避免替换点跳到轨道的侧面甚至反方向。

代码里找最近邻用的双层循环，时间复杂度是 O(M^2)。如果你的数据量超过几万点，建议换 k-d 树来搜邻域，不然跑一次要等很久，尤其多方法求李指数需要在多个参数间来回比较时，计算效率会直接拖累调试进度。

3.3 演化步长、角度限制和数据长度的经验值

实际使用中，evo_step取tau的 1 到 2 倍效果较好。如果取得太大，发散过程已经饱和，对数距离曲线进入平台期；如果取得太小，每一步的增长率受局部几何结构影响太大，噪声会被放大。角度限制一般取 15 到 30 度，太松会导致轨迹跳跃，太紧则可能找不到替换点而提前退出。

数据长度是 Wolf 法的硬约束。经验公式是M > 10^(dim+1)。也就是说，dim=3时至少需要几千个点，dim=5时就需要几十万个点。如果你的实验数据只有几百个点，别硬上 Wolf 法，下面这章的小数据量方法才是正解。

4. 小数据量 Rosenstein 法：处理实验噪声数据的稳妥选择

多方法求李指数在实测场景里，最常遇到的情况是数据短、噪声大、系统参数未知。Wolf 法在数据量不足时极容易中断，这时候 Rosenstein 法更合适。它的思路不是追踪单条轨道，而是统计所有邻近点对的距离发散规律。

4.1 Rosenstein 算法：为什么它能抑制噪声

Rosenstein 法也叫小数据量法，核心是让统计规律来对抗噪声。具体来说，对相空间里每一个参考点，都去找它的最近邻点；然后让所有点对共同演化k步，计算平均对数距离S(k)。由于混沌吸引子是遍历的，大量点对平均之后，局部噪声的影响会被压制下去，剩下的趋势主要反映系统的整体发散率。最后对S(k)的线性区段做最小二乘拟合，斜率就是 LLE。

为什么它特别适合实验数据？因为实验数据通常混有高斯白噪声，噪声在相空间里会造成小尺度上的随机抖动。Wolf 法只追踪一对轨道，一旦在这个小尺度上被噪声误导，结果就会断掉或者偏高。而 Rosenstein 法对成千上万个点取平均，噪声在平均过程中被大量抵消，留下的是系统动态的净趋势。

Rosenstein 法本身也不需要找替换点，这就省掉了最麻烦的角度限制参数。对新手来说，它的参数设定直觉得多：只要给出tau、dim和max_iter，剩下的交给统计平均。

4.2 Matlab 实现要领：找最近邻点与平均发散

function LLE = lyap_rosenstein(ts, tau, dim, max_iter) % ts : 输入时间序列 % tau : 时间延迟 % dim : 嵌入维数 % max_iter : 最大演化步数（建议 200~500） N = length(ts); M = N - (dim - 1) * tau; % 相空间重构 X = zeros(M, dim); for i = 1:dim X(:, i) = ts(1 + (i-1)*tau : M + (i-1)*tau); end % 找每个点的最近邻点（排除时间窗太近的点） min_dist = inf * ones(M, 1); idx_nn = zeros(M, 1); for i = 1:M for j = 1:M if abs(i - j) > 10 dist = norm(X(i, :) - X(j, :)); if dist < min_dist(i) min_dist(i) = dist; idx_nn(i) = j; end end end end % 计算平均对数距离 S(k) S = zeros(max_iter, 1); cnt = zeros(max_iter, 1); for i = 1:M j = idx_nn(i); if j == 0 continue; end for k = 1:max_iter if i + k <= M && j + k <= M dk = norm(X(i+k, :) - X(j+k, :)); if dk > 1e-12 S(k) = S(k) + log(dk); cnt(k) = cnt(k) + 1; end end end end % 线性拟合：只取统计量足够的区段 valid = find(cnt > 0.5 * max(cnt)); y = S(valid) ./ cnt(valid); t = valid * tau; % 实际时间 p = polyfit(t, y, 1); LLE = p(1); end

这里有个容易忽略的点：S(k)并不是全程都能线性拟合。初期噪声主导、后期轨道折叠效应显现，只有中间一段表现出近似线性的增长。我一般会先画出t与S(k)的关系图，肉眼确认线性区段，再自动截取拟合。

参数选择上，max_iter也不是越大越好。如果数据长度只有 1000 点，max_iter取 500 会导致后面几百步只有少数几对点还在统计范围内，方差急剧增大。所以我通常取min(200, M/3)作为最大迭代步数。

4.3 与 Wolf 法的边界：什么数据用 Rosenstein，什么数据用 Wolf

两者没有绝对优劣，只有适用边界的区别。我这里整理了一张对比表：

我的经验法则是：数据长度超过 5000 且信噪比高，优先用 Wolf；数据只有几百到几千点，或者明显带噪，优先用 Rosenstein。实际项目里，我习惯两种方法都跑一遍，把结果放在一起对比，这本身就是一种验证。

顺手提一下，Rosenstein 方法还有一位“近亲”，叫 Kantz 法。Kantz 法不是找一个最近邻点，而是取一个空间小邻域内所有点的平均。它的抗噪能力更强，但计算量也更大。如果你的数据量中等、噪声水平中等，可以考虑把 Kantz 作为第三种交叉验证手段。

5. 多方法求李指数的避坑指南：从分岔图对照到参数敏感性分析

不管用哪种方法，最后都要落到“结果可不可信”上。多方法求李指数的最大风险，是算出一个数但不知道对不对。这一章专门写常见的翻车现场和排查方法。

5.1 先画分岔图，再算 Lyapunov 指数：混沌判据的一致性

在正式算 LLE 之前，先把系统的分岔图画出来。比如 Lorenz 系统，以rho为横轴，画出状态变量x的局部极值，就能看出周期窗口和混沌区域。这时再叠加 LLE 的符号：分岔图上周期 1 对应的参数区，LLE 应该是负的；倍周期分岔点附近，LLE 应该逼近 0；混沌区 LLE 大于 0。如果 LLE 的符号和分岔图对不上，那一定是代码或参数有问题。这是最直观的交叉验证方法。

画分岔图本身也有讲究。对连续的 Lorenz 系统，不能直接采样x的全部点，而是记录轨线穿过某个截平面时的截面点（Poincare 截面）。常见做法是取y=0且dy/dt>0时的x值，这样得到的点云才反映出系统的稳定周期和混沌结构。如果你把连续系统每隔固定步长直接取点画分岔图，看到的多半是一团糊状的散点，什么信息都提取不出来。

5.2 常见问题 1：延迟时间 τ 选太小导致相空间轨迹被挤压

现象：Wolf 法算出来的 LLE 大得离谱，甚至接近log(2)/dt。原因：tau取太小（比如直接取 1），重构的相空间轨迹被压扁在对角线上，邻近点之间的距离被人为拉长。解决：用互信息函数找出第一个极小值点作为tau。对于 Lorenz 标准数据，tau通常取 8 到 15。如果信号是周期信号叠加噪声，tau应大于信号的自相关时间。

为什么会被压扁？当tau远小于系统的特征时间尺度，x(t)和x(t+tau)几乎是同一个变量，重构出来的三维轨迹看起来像一条线段，而不是立体的吸引子。这时最近邻点的距离会被错误地放大，导致累计对数距离虚高。处理方式也很简单：宁可把tau调大一点，也不要让轨迹去找中央对角线。

5.3 常见问题 2：嵌入维数 m 太大导致计算量激增且噪声主导

现象：计算耗时暴增，而且结果随着dim增大不断漂移。原因：嵌入维数过高，维数灾难显现，高维空间里所有点之间的距离都变得相近，信噪比下降。解决：用 Cao 方法或虚假最近邻点法确定最小dim。一般情况下，Lorenz 系统取dim=3就有很好的结果；实测电压信号取 4 到 6 就够了，不要无脑加到 10 以上。

举个具体例子，一段 10000 点的实测信号，dim=3时最近邻搜索只需遍历 10000 个点，dim=10时不仅耗时翻倍，而且高维空间里任意两点之间的欧氏距离都差不多大，最近邻点失去了物理意义，算出来的 LLE 完全由噪声决定。

5.4 常见问题 3：Wolf 法与 Rosenstein 法结果不一致

现象：同一段数据，Wolf 算出 0.9，Rosenstein 算出 0.5。原因：两种方法对参数和数据质量的要求不同，evo_step、拟合区段、噪声水平都会影响最终数值。解决：不要只盯着最终数字。把 Rosenstein 的t与S(k)关系图打印出来，检查线性区段是否明显。如果线性区段很短，说明数据本身质量不足。这种场景下我更信任 Rosenstein 的斜率，因为它来自统计平均。

我再加一个排查步骤：把 Wolf 法的evo_step从 1 扫到 10，看 LLE 的变化曲线。如果 LLE 随evo_step单调上升，说明数据里含有大量高频噪声；如果 LLE 随evo_step大幅抖动，说明轨道折叠效应太强。真正靠谱的结果应该对evo_step呈现一段平台区，取平台区的中间值作为最终结果。

5.5 常见问题 4：数据里有趋势项或直流偏置

现象：Wolf 法直接报错或者算出负值。原因：时间序列里如果有缓慢变化的趋势项（比如传感器漂移），重构后的相空间轨迹会被趋势项拉伸成一个巨大的弧线，局部发散趋势被淹没。解决：先做去趋势和去均值预处理，再做相空间重构。差分操作要谨慎，差分会放大高频噪声，弄不好会把真实混沌信号的高频成分滤掉。

我们在处理混沌电路实测电压时，最常见的就是运放偏置电压导致的直流偏置。去掉均值之后重算，LLE 从负的变成正的，这种情况我遇到不下十次了。

6. 进阶技巧：批量扫描参数画 Lyapunov 谱，以及验证结果的三个习惯

多方法求李指数一旦跑通，正经用法是批量扫参画谱。以 Lorenz 为例，把rho从 0 扫到 300，就能得到完整的 Lyapunov 谱，用来分析倍周期分岔路径和混沌区间。批量扫描时，注意别在循环里反复分配数组。先预分配LLE_list，再循环调用函数，速度能快很多。

% 批量扫描 Lorenz 系统的 rho，计算最大 Lyapunov 指数 rho_list = 0:0.5:300; LLE_list = zeros(size(rho_list)); for idx = 1:length(rho_list) LLE_list(idx) = lyap_lorenz(rho_list(idx)); % 内部含 RK4 和 QR end plot(rho_list, LLE_list, 'LineWidth', 1.5); xlabel('rho'); ylabel('LLE'); grid on;

这个函数式封装有个好处：它天然支持多方法求李指数的横向对比。同一份rho_list，分别用 Jacobian 法和 Rosenstein 法跑一遍，把两条曲线叠在一起画，重合度越高，结果越可信。

最后分享三个验证结果的习惯。第一，先用已知结论做基准测试：Lorenz 在rho=28时 LLE 约 0.9，如果你的代码算出来是 0.5，先别急着调参数，回头检查代码逻辑。第二，检查 Lyapunov 指数谱的求和规律：所有指数之和应该等于 Jacobian 矩阵的迹。算出来若不满足，说明积分过程或正交化步骤有 bug。第三，保存好参数文件。我自己的血泪经验是，调试阶段跑了几十个参数组合，过两周再看数据，根本不记得当时用的是哪个tau和dim。把参数和脚本版本记录在文件名里，等于给自己留一颗后悔药。希望这些经验能帮你在多方法求李指数的路上少走弯路。

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