简介：这是一份面向非线性动力学与混沌理论学习者的Matlab源码包，围绕洛伦兹系统与Logistic映射，提供分叉图、庞加莱截面图和李雅普诺夫指数图的完整绘制代码。洛伦兹系统由三个非线性微分方程构成，是研究蝴蝶效应与确定性系统不可长期预测的经典模型；Logistic映射则直观展示从周期点到混沌的参数演化过程。包内共6个文件，以5个.m脚本为主，另有1个.asv自动保存备份，整体仅2KB，代码精炼，适合学生、研究人员及数学建模爱好者直接运行和二次修改。通过运行这些脚本，可以直观观察系统轨迹在庞加莱截面上的无规则分布，计算李雅普诺夫指数并判断系统稳定性，从而深入理解混沌现象的数学本质。目前已有812人学习下载，是入门混沌理论、开展数值实验的高性价比参考资料。
1. 洛伦兹系统、Logistic 与洛伦兹分叉图：一份 Matlab 混沌仿真到底在跑什么
打开压缩包看到“Logistic_Lorenz_matlab_洛伦兹分叉图”这组关键词的人，多半不是奔着理论来的，而是要在 Matlab 里把混沌的两条经典路径跑通：一条是离散迭代的 Logistic 映射，一条是连续微分方程的 Lorenz 系统。前者用一行迭代就能画出教科书级的分叉图，后者用 ode45 积分 50 秒就能看到那对著名的“蝴蝶翅膀”。而标题里的“洛伦兹分叉图”，则是把连续系统也用分叉图的语言讲清楚的关键一步。这篇文章会按“先立住理论、再动手复现、最后避开坑”的顺序，把这套流程完整拆开。适合刚装好 Matlab、照着网盘教程把 2023b 或更早版本跑起来，却不知道参数该往哪填的人；也适合已经能画图、但总觉得自己的蝴蝶和分叉图跟论文里长得不太一样的熟手。先说明一点：这个包的核心不是某个现成函数，而是你手里那份脚本背后的三个技术动作——迭代、积分、采样画分叉图。把这三个动作吃透，比收集任何 .m 文件都值钱。
2. Logistic 映射与分叉图：离散混沌的最小模型
2.1 周期倍增路径：为什么先讲 Logistic 而不是直接看 Lorenz
Logistic 映射的形式是 x_{n+1} = r·x_n·(1 - x_n)，它把一个区间内的实数映射回自己。r 是控制参数，x 是状态。当 r 在 2.4 到 4 之间变化时，系统从稳定不动点走向周期 2、周期 4，再走向混沌，这条路径就是周期倍增分叉路径。Matlab 社区里几乎所有混沌仿真教程都把 Logistic 放在最前面，原因很实际：它不需要求解微分方程，不需要考虑积分步长，只需要一层 for 循环迭代足够多次，就能看到完整的分叉结构。这个特征让它成为检验“你是否理解了分叉图到底在画什么”的最佳试金石。
有一个关键点需要先强调：分叉图不是“把迭代序列全部画出来”，而是“在丢弃瞬态之后，把稳态行为画在横轴 r 上”。所谓瞬态，是指从初值出发到系统收敛到吸引子之前的过渡段。如果把这个过渡段也画进去，图上会出现一堆杂乱的过渡点，淹没真正的分叉结构。这也是后面避坑章里最常见的问题之一。
2.2 用 Matlab 画出第一张 Logistic 分叉图：向量化写法
常见做法是让 r 作为行向量，一次性更新所有参数对应的 x 值，而不是对每个 r 单独循环。向量化能显著缩短运行时间，而且代码更接近迭代公式的原始形态。
% Logistic 分叉图：向量化写法 r = 2.4:0.001:4; % 控制参数范围，步长 0.001 nr = length(r); % 参数点数量 x = 0.5 * ones(1, nr); % 所有参数共用一个初值 x0 = 0.5 drop = 200; % 丢弃前 200 次迭代（瞬态） N = 600; % 总迭代次数 figure('Color', 'w'); hold on; for k = 1:N x = r .* x .* (1 - x); % 迭代公式，x 是向量 if k > drop plot(r, x, '.', 'MarkerSize', 1, 'Color', [0.2 0.4 0.8]); end end xlabel('r'); ylabel('x'); title('Logistic 分叉图');这段代码的逻辑是：每次迭代里，所有 r 对应的 x 同时更新，更新完判断当前迭代次数是否已经超过丢弃阈值。超过才画点，保证图上只显示稳态部分。注意plot放在循环内部、每次都画全量 r 的点，这对 1601 个参数点来说完全够快；如果 r 步长取到 0.0001 级别，我一般会改成分批累积到矩阵再一次性绘图，避免图形句柄刷新占用时间。
2.3 三个必调参数：r 范围、瞬态丢弃数、每 r 采样点数
第一个参数是 r 的范围。常见做法是 2.4 到 4，因为 r < 2.4 时只有稳定不动点，画出来就是一条平线，信息量不大。如果你想把分叉图做得更“满”，可以改成 2.8 到 4；但如果要完整展示周期倍增路径，从 2.4 起步更标准。第二个参数是瞬态丢弃数 drop，它直接决定图的干净程度。drop 取 100 到 300 之间一般够用，200 是稳妥值。第三个参数是每个 r 对应的稳态采样点数，也就是 N - drop。这个值决定分叉图在混沌区间（r 接近 4 时）的“密度”——混沌区间里 x 会取遍几乎整个区间，采样点太少会显得稀疏，太多则图整体变黑。对 1601 个参数点，N = 600 已经能画出清晰的轮廓；如果你把 r 步长加密到 0.0005，建议把 N 提到 800 以上。
这里有个容易混淆的细节：分叉图中的“每 r 采样点数”和你绘图时的点数不是一回事。绘图时每个 r 只画一个点也行，但稳态会在周期点之间跳动，所以你最终会看到多条分支线。混沌区间则是无数个点覆盖成一条带。有些教程用scatter代替plot，效果类似，但scatter对大规模点更吃内存，1 个点 1 个点的画法在参数点过万时会明显卡顿。我一般只在大范围参数扫描时才改用矩阵存储 + 单次绘图。
3. Lorenz 系统：从微分方程到蝴蝶翅膀
3.1 经典参数与系统结构：连续系统的混沌为什么难调
Lorenz 系统是三个一阶常微分方程：
dx/dt = σ(y - x) dy/dt = x(ρ - z) - y dz/dt = xy - βz
其中 σ 是普朗特数，ρ 是瑞利数相关的参数，β 是几何参数。教科书经典取值是 σ = 10、ρ = 28、β = 8/3，这个组合下系统进入混沌状态，相空间轨迹形成双叶吸引子，也就是俗称的蝴蝶翅膀。物理背景是大气对流简化模型，但对仿真来说，你只需要关心它作为连续混沌系统的两个核心性质：初值极敏感、轨迹在吸引子内永不重复。
连续系统的仿真比离散迭代多一层复杂度：你需要选积分方法、积分区间、输出步长，还要处理 ode45 的误差控制。很多第一次跑 Lorenz 的人会发现：直接按默认设置画出来的蝴蝶只有半边翅膀，或者轨迹很快就飞出画面。这通常不是方程写错，而是积分容差太松或积分区间太短。
3.2 ode45 的最小实现：一个能直接跑的 Lorenz 仿真
% Lorenz 系统仿真：sigma=10, rho=28, beta=8/3 sigma = 10; rho = 28; beta = 8/3; f = @(t, y) [sigma*(y(2) - y(1)); ... y(1)*(rho - y(3)) - y(2); ... y(1)*y(2) - beta*y(3)]; [t, y] = ode45(f, [0 50], [1; 1; 1], odeset('RelTol', 1e-6)); figure('Color', 'w'); plot3(y(:,1), y(:,2), y(:,3), 'LineWidth', 0.8); grid on; xlabel('x'); ylabel('y'); zlabel('z'); title('Lorenz 混沌吸引子');这里f是匿名函数，输入是时间 t 和状态向量 y，输出是三个导数值。ode45 的四个参数分别是：方程函数、时间区间 [0 50]、初值 [1;1;1]、以及通过odeset设置的积分选项。RelTol是相对误差容差，默认值是 1e-3，对 Lorenz 这种对误差极敏感的系统来说偏松，我通常调到 1e-6 甚至 1e-8。注意这里我用了[1; 1; 1]作为初值，它足够偏离不动点，能很快进入吸引子。
积分区间 [0 50] 对应的是无量纲时间。对经典参数来说，这足以让轨迹在吸引子上绕几十圈，画出完整的双叶结构。如果你把区间改成 [0 10]，可能只看到轨迹在其中一个叶附近绕圈，看起来像半个蝴蝶——这不是 bug，是观测窗口太短的体现。
3.3 观察窗口与积分选项：为什么默认设置画不出“标准图”
ode45 的步长是自适应变化的，它会根据误差容差自动加密或放宽步长。你看到的三维轨迹曲线，实际上是把内部计算点插值到输出点之后的结果。这里有一个常见误区：t向量的长度不是由你决定的，而是由 ode45 根据误差控制自动生成的。如果你希望轨迹更平滑、点数更均匀，可以额外设置'MaxStep'。对 Lorenz 系统，我一般把 MaxStep 设在 0.01 到 0.05 之间。太小会让积分变得很慢，太大则轨迹的细节可能被跳过，尤其是在两个翅膀切换的瞬间。
另外一个值得说明的点是：plot3默认的线宽 0.5 在转角处会显得纤细，蝴蝶轮廓不够清晰。我习惯把LineWidth调到 0.8 到 1.2，颜色用深蓝或深红系，而不要用亮黄色——亮色在白色背景上反光严重，分叉结构的细节会被吞掉。如果你觉得轨迹太密、看不出翼型轮廓，可以每 N 个点采样一次再画线。这个技巧放到第 6 章展开。
4. 洛伦兹分叉图的两种画法：极值采样与庞加莱截面
4.1 连续系统分叉图的难点：一个点怎么变成一条线
这里要先厘清一个概念：Lorenz 系统本身是连续的，轨迹在相空间里是一条连续的曲线，它不像 Logistic 那样天然有“迭代点序列”。所谓洛伦兹分叉图，本质上是把连续系统降到离散——通过某种采样规则，从轨迹中提取出一个离散序列，再以这个序列为纵轴、以某个系统参数（通常是 ρ）为横轴，绘制分叉图。这个做法在物理和工程文献里非常常见，核心目的是观察当 ρ 变化时，系统的动力学行为如何从周期走向混沌、再在混沌中穿插周期窗口。
采样规则有两大流派。第一种是取轨迹在某个方向上的局部极值，比如每次轨迹到达蝴蝶翅膀折返点时，记录当时的 z 值；第二种是取庞加莱截面，也就是让轨迹穿过一个指定平面时记录穿过点的坐标。两种方法各有侧重：极值法简单直观，适合快速确认混沌区间；庞加莱截面法信息量更大，能同时看到多个变量的状态，但实现时对选面方向有讲究。
4.2 局部极值法：以 ρ 为参数扫描 z 的转折点
% Lorenz 分叉图：取 z 方向的局部极大值 rho_list = 20:0.2:60; % 控制参数扫描范围 sigma = 10; beta = 8/3; figure('Color', 'w'); hold on; for rho = rho_list f = @(t, y) [sigma*(y(2)-y(1)); ... y(1)*(rho-y(3)) - y(2); ... y(1)*y(2) - beta*y(3)]; [t, y] = ode45(f, [0 100], [1; 1; 1], ... odeset('RelTol', 1e-6, 'MaxStep', 0.05)); z = y(:, 3); dz = diff(z); % 局部极大值：前一点上升、后一点下降 idx = find(dz(1:end-1) > 0 & dz(2:end) < 0) + 1; plot(rho * ones(size(idx)), z(idx), '.', 'MarkerSize', 1, ... 'Color', [0.8 0.2 0.2]); end xlabel('rho'); ylabel('z 局部极大值'); title('洛伦兹分叉图（局部极值法）');这段代码的扫描节奏很慢是正常的，它在 ρ = 20 到 60 之间以 0.2 为步长，每个参数点都要完整积分到 t = 100。MaxStep设为 0.05 是为了保证极值检测的精度——如果步长太大，轨迹在折返点的采样点过少，diff检测到的极值位置会偏离真实值，导致分叉图在周期区域出现毛刺。局部极大值的检测逻辑是：dz表示相邻采样点的差，dz(1:end-1) > 0意味着前一段在上升，dz(2:end) < 0意味着后一段在下降，两者同时满足的位置就是峰值。
这个代码跑出来的图，你会在 ρ 从 24 附近开始看到清晰的周期倍增结构，到 ρ = 28 附近进入混沌，之后混沌区间里穿插着几个周期窗口，最明显的是 ρ 在 35 到 40 之间的某几条“缝隙”。如果你把局部极大值换成局部极小值，图的形状会反过来，但周期倍增的特征位置不变。
4.3 庞加莱截面法：用一个平面截出系统的“指纹”
庞加莱截面的做法是选一个合适的平面（通常是 z = ρ - 1 附近，这是系统不动点的 z 坐标之一），当轨迹穿过这个平面且方向一致时记录交点。代码上不需要专门做三维几何求交，可以用条件判断来捕获“从平面一侧穿越到另一侧”的时刻：
% 庞加莱截面：z = z0 平面，取正向穿越 z0 = rho - 1; % 截面位置 [~, y] = ode45(f, [0 120], [1; 1; 1], ... odeset('RelTol', 1e-6, 'MaxStep', 0.02)); % 正向穿越：当前点 > z0，前一个点 < z0 cross_idx = find(y(2:end, 3) > z0 & y(1:end-1, 3) < z0); x_section = y(cross_idx, 1); y_section = y(cross_idx, 2);这里返回的x_section和y_section就是截面点集。把它画出来，周期运动对应有限个离散点，混沌运动对应一条类似曲线的密集点集。庞加莱截面相对极值法的好处是：当系统在某个参数附近出现周期窗口时，截面图上能明显看到点数突然变少、形成闭合环或几个孤立点，这是极值法容易忽略的细节。缺点是对截面位置敏感——选在轨迹中间位置时，正负穿越都可能发生，如果不去区分方向，图会变得杂乱。我一般只取正向穿越，并且MaxStep压到 0.02 以下，否则截面点的位置误差会很大。
5. 洛伦兹与 Logistic 仿真的 5 个高频坑：从玄学变成可定位的报错
5.1 分叉图右侧一片黑雾，但完全没有分叉结构
这是最多人问的现象。现象：r 取到 3.8 以上时，图右侧应该是密集的“雾状带”，但你画出来只有一条粗线，或者干脆一团黑。原因：瞬态丢弃数太少，或者总迭代次数不够，使得非稳态点混入图中。解法：把drop提到 300 以上，N提到 1000 以上再试。如果还是黑，检查plot的MarkerSize是不是取得太大——超过 2 就会让点之间互相覆盖，混沌带变成实心黑块。我自己的血泪经验是：先用MarkerSize = 1出图，确认结构对了再慢慢放大，而不是反过来。
5.2 Lorenz 轨迹只绕一边翅膀转，另一只翅膀始终不出现
现象：按经典参数 ρ = 28 跑出来的轨迹，在三维图里只在 x 正半轴或负半轴一侧绕圈。原因：积分区间太短，轨迹还没完成从一侧翅膀跳跃到另一侧的过程。Lorenz 吸引子的两个翅膀之间切换时间是不规则的，短区间里可能恰好没有切换。解法：把时间区间从 [0 50] 拉到 [0 100]，或者把初值从 [1;1;1] 改成 [10; 10; 10]——离原点更远的初值会更快进入充分的绕行状态。如果你已经积分到 100 还是只看到一只翅膀，那就要怀疑RelTol是不是太松，导致数值误差把轨迹牢牢钉在了一侧。把RelTol从默认 1e-3 改成 1e-6，这类问题九成能解决。
5.3 换了 Matlab 版本（2023b 或更新版）后，图线颜色和线型“莫名其妙”变了
现象：同样的 plot 代码，在旧版跑出来是彩色分叉图，换到新版本后所有点变成同一种颜色，或者三维轨迹的宽度变细了。原因：Matlab 从 R2014b 开始使用默认色序（parula），R2023b 又调整了plot3的默认线宽和颜色分配逻辑；用hold on叠加多次plot时，新版不再自动旋转颜色，而是固定用第一条线的颜色。解法：不要依赖默认颜色，显式给每组plot指定'Color'参数。分叉图尤其要注意：如果循环里不写颜色，所有 r 的点会共用同一个颜色，混沌带和周期区几乎无法区分。显式指定色值虽然啰嗦，但换来的是任何版本下输出一致。
5.4 洛伦兹分叉图画出来全是噪声，找不到周期窗口
现象：用 4.2 节的极值法扫描 ρ，得到的图从左到右全是密密麻麻的点，完全看不出周期倍增结构。原因：三成是MaxStep太大导致极值点偏移，七成是积分时间不够长，稳态没有建立。Lorenz 系统在某个参数下进入混沌后，吸引子上一些区域是被“排斥”的，轨迹需要时间“沉降”到吸引子附近。解法：把积分区间加长到 150 以上，并且只取后半段的极值。代码层面就是[t, y] = ode45(...)之后加一句y = y(t > 50, :);把前 50 个时间单位作为瞬态丢弃。这个动作对应的是离散系统里drop的概念，但连续系统的瞬态长度通常要比离散迭代多得多。
5.5 把 Logistic 的分叉图和 Lorenz 的轨迹画在同一张图里，结果乱了套
现象：用同一个坐标系想同时展示 Logistic 迭代点和 Lorenz 轨迹，得到的图缠在一起看不出任何结构。原因：Logistic 是 1D 迭代，Lorenz 是 3D 轨迹，两者状态空间维数不同，横纵轴含义完全不同，强行共坐标系没有意义。解法：用subplot或tiledlayout分开放置，不要共用坐标轴。更合理的做法是：先用 Logistic 分叉图理解“分叉图是什么”，再用 Lorenz 的庞加莱截面建立“连续系统的分叉视图”，这两者本来就是互补关系，不是替代关系。这个错误在网上下载的整合包里特别常见——脚本作者为了省事把所有图堆在一个 figure 里，结果每个图都失去了可读性。
6. 从画得出到读得懂：三个验证混沌的必要技巧
画图只是第一步。这条线的最终价值在于你能确认自己得到的图“真的是混沌”，而不是数值误差制造的人工产物。这里给出三个在 Matlab 里低成本就能实现的验证手段。
第一个手段是初值敏感性检验。对同一个 Lorenz 系统跑两组初值，让它们的初始差只有 1e-8 量级，然后把两条轨迹的 x 分量差值随时间画出来。如果差值随时间近似指数增长，说明系统是混沌的；如果差值始终很小或线性增长，说明你选参数落在了周期区。如果两条轨迹跑出来完全一致，那大概率是RelTol设得太松，od45 在误差控制下把轨迹“压”到了同一条近似解上。
第二个手段是 Lyapunov 指数估计。对 Logistic 映射，可以跳过复杂算法，用简化估计直接算：
% Logistic 映射的 Lyapunov 指数简化估计 r = 3.9; x = 0.1; lambda_sum = 0; for k = 1:2000 lambda_sum = lambda_sum + log(abs(r - 2*r*x)); x = r * x * (1 - x); end lambda = lambda_sum / 2000; disp(lambda);逻辑很简单：对一维映射，雅可比就是导数 f'(x) = r - 2rx，每次迭代累加导数的对数值再平均，λ > 0 说明相邻轨迹指数分离，是混沌；λ < 0 说明是周期运动。这个估算方法对 Logistic 这类单变量映射足够准，但不要直接套到多变量系统上——多变量要算雅可比矩阵的乘积在李代数意义上的特征值，复杂度完全不同。
第三个手段是验证周期窗口：在庞加莱截面上，如果构成一条闭合曲线，那就是周期运动；如果形成密集点集，则是混沌。你可以用 4.3 节的截面代码扫描 ρ = 30 到 35 之间的几个点，在截面图上会看到从几个孤立点向密集曲线的转变过程。这个验证比单纯肉眼看轨迹形状可靠得多。
我自己做这套流程的固定习惯是：先跑最小代码确认不出错，再调参数看变化，最后才批量扫描参数范围。顺序反了，你会花大量时间在排错上，而且很难分清是数值问题还是混沌本身。希望帮到你。
本文还有配套的精品资源，点击获取
