简介这是一份面向非线性动力学与混沌理论学习者的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 而不是直接看 LorenzLogistic 映射的形式是 x_{n1} 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 系统仿真sigma10, rho28, beta8/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 开始使用默认色序parulaR2023b 又调整了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 之间的几个点在截面图上会看到从几个孤立点向密集曲线的转变过程。这个验证比单纯肉眼看轨迹形状可靠得多。我自己做这套流程的固定习惯是先跑最小代码确认不出错再调参数看变化最后才批量扫描参数范围。顺序反了你会花大量时间在排错上而且很难分清是数值问题还是混沌本身。希望帮到你。本文还有配套的精品资源点击获取