简介一套面向混沌理论学习的MATLAB源码包围绕洛伦兹系统与Logistic映射展开适合希望直观理解非线性动力学、混沌现象及“蝴蝶效应”的本科生、研究生和科研人员。压缩包内包含6个文件主体为5个.m脚本另附1个.asv自动保存副本整体仅2KB代码精简、无冗余依赖可直接在MATLAB中运行。这些脚本分别实现洛伦兹系统的庞加莱截面图、最大李雅普诺夫指数求解以及Logistic映射分叉图绘制等功能覆盖了从连续系统到离散映射的典型混沌分析路径。已有812人学习下载对于混沌理论入门与演示场景具有较实用的参考价值。通过运行这些源码读者可以直观观察系统轨迹在截面上的复杂无规则分布借助李雅普诺夫指数的正负判断系统是否处于混沌状态并对比Logistic映射随参数变化产生的周期分叉与混沌区域从而获得从模型搭建、数值计算到结果解读的完整实践体验。1. 洛伦兹系统与 Logistic 映射为什么这两张分叉图常被放在同一个包里从网上下载一个叫“洛伦兹系统.rar”的压缩包解压后大概率会看到两个关键词Logistic 和 Lorenz。一个是一维离散迭代一个是三维连续常微分方程但它们都是混沌理论课上的“入门必画”。很多人拿到这个包第一个问题是哪个文件是主脚本改了参数为什么图就散了第二个问题是Logistic 分叉图我懂洛伦兹系统也有分叉图吗怎么画实际上这两张图放在同一个包里很合理——Logistic 映射是研究倍周期分岔的最简模型Lorenz 系统则是第一个被发现混沌解的三维连续系统二者在混沌的“普适性”上遥相呼应。这篇笔记不讲空理论直接带你用 Matlab 把两张图从零画出来参数怎么设、瞬态怎么抛、分叉点怎么看以及那个让无数人翻车的“局部极大值法”到底怎么避坑。2. 用 Matlab 把 Logistic 分叉图画出来向量化写法与两个关键参数2.1 Logistic 迭代公式与分叉图原理Logistic 映射的迭代式是x_{n1} r * x_n * (1 - x_n)其中r是控制参数x在 0 到 1 之间。当r从 2.8 增长到 4.0 时系统的稳态行为从固定点变成周期分叉最终进入混沌。分叉图的本质是横轴取r纵轴取迭代替掉瞬态后的x的长期值。如果你把每个r下最后 100 个x都画出来就能看到周期分裂成两条线、四条线最终变成一片混沌带。Matlab 里画分叉图最忌讳用for循环嵌套for循环因为r可能取上千个点每个点又要迭代上千次写成纯循环会卡到怀疑人生。常见做法是向量化让所有r同时迭代用矩阵运算一次性推进。这也是网上那些 20 行的分叉图脚本能跑得飞快的原因。2.2 向量化绘制 Logistic 分叉图从 r2.5 到 4.0下面是最常见的一种向量化实现我一般会直接用到 2.5 到 4.0 的完整区间这样能看到清晰的周期窗口。% logistic_bifurcation.m clear; clf; r 2.5:0.001:4.0; % 参数范围步长 0.001 约 1500 个点 N 800; % 每个 r 迭代总次数 M 200; % 抛弃前 M 次瞬态 x 0.2 * ones(size(r)); % 所有 r 共享同一个初始值也可以随机 % 先迭代 M 次抛弃瞬态不绘图 for i 1:M x r .* x .* (1 - x); end % 继续迭代 N-M 次每次绘图 hold on; for i 1:(N-M) x r .* x .* (1 - x); plot(r, x, ., MarkerSize, 1); end hold off; xlabel(r); ylabel(x); title(Logistic Map Bifurcation);这段代码的逻辑是先把x初始化为 0.2然后用for循环迭代M200次但不画图目的是把系统从任意初始值“拉”到吸引子上。之后每迭代一次就把当前所有r对应的x值画成散点。由于r是一个行向量x也是行向量r .* x .* (1-x)就是逐元素运算一次迭代就是整个向量的一步更新。参数说明r的步长决定了分叉图的横向分辨率0.001 已经能清晰看到分叉点。N越大画出来的吸引子越完整但点也越多MarkerSize调到 1 避免图片发黑。M的取值要足够大否则在周期窗口会出现细线或噪声点尤其是r在 3.83 附近的周期 3 窗口瞬态可长达数百步。2.3 提高绘制速度预分配矩阵与一次性绘图上面的脚本每次迭代都调用一次plot当N800时要画几百次在旧版 Matlab 上会有点慢。更高效的做法是把长期值存入一个大矩阵最后用scatter一次性绘图。% logistic_bifurcation_fast.m r 2.5:0.001:4.0; N 500; % 采样点数 M 300; % 抛弃瞬态数 x 0.2 * ones(size(r)); for i 1:M x r .* x .* (1 - x); end X zeros(length(r), N); % 预分配矩阵每行一个 r每列一次采样 for i 1:N x r .* x .* (1 - x); X(:, i) x; end % 将矩阵展开成向量 R repmat(r(:), N, 1); Xvec X(:); plot(R, Xvec, ., MarkerSize, 1); xlabel(r); ylabel(x);预分配矩阵的好处是避免反复调用绘图函数内存占用也很小1500 个 r 乘 500 次采样约 75 万个点单次plot完全扛得住。repmat是把r纵向复制N次让每一个采样点都能找到对应的横坐标。如果你希望初始值不固定可以把x rand(size(r)) * 0.9 0.05但要注意避免初始值恰好取到 0 或 1那样迭代会卡死在边界。3. Lorenz 系统三维相图与分叉图ode45 和 Poincare 截面3.1 Lorenz 微分方程与经典参数Lorenz 系统是 1963 年从大气对流模型中抽象出的三维常微分方程组dx/dt sigma*(y - x) dy/dt x*(rho - z) - y dz/dt x*y - beta*z经典参数是sigma 10、beta 8/3rho 28。在rho 28时系统处于混沌状态相轨线在三维空间里像一只蝴蝶的翅膀。画相图很简单用ode45直接数值积分即可但如果要画洛伦兹分叉图就需要把rho当作变化参数对每个rho求解一次微分方程再提取长期行为的某个度量。这里有一个关键认知Lorenz 系统是连续系统它的“分叉图”不像 Logistic 那样有明确的迭代序列。常见做法是把连续轨迹的局部极大值比如x坐标的峰值按顺序提出来观察它们随rho的变化。这也是标题里“洛伦兹分叉图”最常指代的画法。3.2 用 ode45 求解 Lorenz 系统并绘制三维相图% lorenz_phase.m sigma 10; beta 8/3; rho 28; 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]); plot3(Y(:,1), Y(:,2), Y(:,3), LineWidth, 0.5); xlabel(x); ylabel(y); zlabel(z); title([Lorenz attractor, rho , num2str(rho)]);ode45是 Matlab 自带的变步长四阶-五阶 Runge-Kutta 求解器对非刚性系统足够。这里的匿名函数f接收时间t和状态向量Y返回[dx/dt; dy/dt; dz/dt]。初始条件[1;1;1]不是定式只要不是原点附近的不稳定平衡点就行。时间跨度[0 50]在rho28时能画出一圈完整的蝴蝶形状大约 5000 个输出点绘图很流畅。需要注意的是ode45的输出点不均匀默认会按内部步长加密。如果你只想看长期吸引子可以把前 10 秒的轨迹丢弃比如改成[t,Y] ode45(f, [10 50], [1;1;1])这样Y里就是稳态后的轨迹相图不会包含从初始点到吸引子的过渡段。3.3 绘制 Lorenz 系统关于 rho 的分叉图局部极大值法这是重头戏也是网络上各种 rar 包里最容易“画错”的部分。我推荐的做法是对每个rho从rho20到rho100用ode45积分足够长时间抛弃瞬态后找到x(t)的局部极大值然后把这些极大值作为纵轴散点画在rho对应的横坐标上。% lorenz_bifurcation.m sigma 10; beta 8/3; rho_list 20:0.5:100; hold on; for k 1:length(rho_list) rho rho_list(k); 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]); % 抛弃前 50 秒瞬态只保留 t50 的数据 idx t 50; x Y(idx, 1); % 用差分找局部极大值 dx diff(x); peaks x(1:end-1) 0 dx(1:end-1) 0 dx(2:end) 0; peak_x x(peaks); if ~isempty(peak_x) plot(rho * ones(size(peak_x)), peak_x, ., MarkerSize, 2); end end hold off; xlabel(rho); ylabel(local maxima of x); title(Lorenz bifurcation diagram (local maxima));逻辑说明diff(x)得到每一步的变化量若某点的左边变化量为正、右边变化量为负那么该点就是局部极大值。由于x是列向量dx(1:end-1)和dx(2:end)的长度都正好与x(1:end-1)对齐。这里的peaks是逻辑索引筛选出满足条件的点。参数说明rho_list的步长取 0.5能看出分叉结构但细节可能不够。如果你想把从周期到混沌的分岔点看得更准可以把步长改成 0.1但计算量会成倍增加因为每个rho都要调一次ode45时间跨度越长越慢。积分时间[0 100]里t50相当于丢弃前 50 秒这取决于系统初始瞬态的衰减速度。rho接近 24.74 时瞬态极长可能需要积分到 200 秒才稳定。4. 洛伦兹分叉图的三种画法局部极值、截面、峰间间隔4.1 方法一局部极大值法的优缺点上面写的局部极大值法是最直观的因为x(t)的峰值反映了振荡幅度。对rho从 24 到 28 之间系统会经历周期激增如周期 1、周期 2、混沌局部极大值能清楚揭示这种分岔。缺点是当系统处于混沌时局部极大值的数量会很多且不规则散布成一条带子当系统处于周期状态时极大值的数量等于周期数。因此画出来的图在周期区是几条线在混沌区是一片点非常好读。但要注意rho太小时比如小于 24.74系统收敛到固定点此时迭代很久x是常数没有局部极大值peaks全为假画不出点。这是正常的不要以为是代码 bug。4.2 方法二Poincare 截面法Poincare 截面法更适合混沌研究在相空间里取一个截面例如z rho - 1平面记录轨迹每次穿过该截面时的(x, y)坐标然后观察这些点的分布。对 Lorenz 系统常见截面是y 0且dy/dt 0。% lorenz_poincare.m sigma 10; beta 8/3; rho 28; 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 200], [1; 1; 1], odeset(MaxStep, 0.01)); % 找 y0 且 dy/dt 0 的点 dy sigma*(Y(:,2) - Y(:,1)); % 实际上 dy/dt 就是方程第一项不对应该用 Y(:,2) 的导数 % 更简单用 y0 前后符号变化 y Y(:,2); sign_y sign(y); idx find(diff(sign_y) 2); % 从负变正且是上升 x_p Y(idx, 1); z_p Y(idx, 3); plot(x_p, z_p, .);上面的代码需要仔细说明sign_y取y的符号diff后如果从 -1 变为 1差值为 2表示轨迹从负半轴穿到正半轴也就是截面y0并正向穿过。MaxStep限制最大步长避免漏掉截面。不过由于ode45是变步长直接取Y(idx,1)可能不精确因为交点不在采样点上。更精确的做法是用插值但作为快速查看够用了。4.3 方法三峰间间隔法峰间间隔法把相邻两个局部极大值的时间差取出来观察间隔的分布。这个方法适合研究混沌系统的“相位”特性比如 Lorenz 系统在rho较大时会出现峰值交替现象。代码跟局部极大值法类似但要把峰值对应的时刻也记录下来然后画间隔随峰值序号的变化。这个方法画出的不是经典分叉图而是一种动力学诊断图适合进阶分析。三种方法的对比如下方法横轴纵轴适合场景计算成本局部极大值法rhox 的峰值观察分岔和混沌带中等Poincare截面rho截面点坐标研究轨道拓扑结构较高峰间间隔法峰值序号时间间隔分析间歇混沌中等我一般会先用局部极大值法因为它最容易被非专业读者接受。如果你在复现别人的 rar 包时看到“分叉图”是一片密密麻麻的散点八成就是局部极值法。而如果你的包里有三个子图可能就是把x、y、z三个分量分别提峰值后画在一起。5. 避坑与常见问题从 NaN 到假分叉五个必查点5.1 现象一Logistic 分叉图出现横线或竖线跟染料一样原因r的步长太大导致每个r的迭代值在绘图时错位或者M抛得不够瞬态点也被画上去了。更常见的是你在迭代后没有把x限制在 0 到 1 之间当r超过 4 时x会迅速溢出到无穷大。解决检查r是否最大到 4.0不要超过 4.0。把M设到 500 以上。如果图上出现孤立点打印几个r下的x值看看是不是数值发散。5.2 现象二Lorenz 分叉图在 rho24 附近出现“毛刺”甚至有几条断线原因瞬态未完全消除。当rho接近 24.74霍普夫分岔点时系统需要极长的时间才能收敛到极限环你只积分到 100 秒后 50 秒仍然包含衰减振荡。解决把积分时间从[0 100]加长到[0 300]并把丢弃段提高到t 200。如果你用ode45且MaxStep默认 0.5 左右步长太大也可能导致局部极值检测出错建议设置odeset(MaxStep, 0.05)让微分方程更平滑。5.3 现象三ode45 跑了半天不出图CPU 占用 100%原因rho_list太密集比如步长 0.01范围 0 到 200要解两万个微分方程每段积分 100 秒这在一台普通笔记本上可能要跑几小时。解决先降低分辨率把rho_list步长改为 1快速预览大致形态。然后按需细化到 0.1 或 0.05。还有一种做法是用parfor并行循环但需要先准备parpool。如果你只是画图完全没有必要对每个rho都积分那么长时间先试[0 30]加瞬态丢弃看形状对不对再加大。5.4 现象四分叉图上的点阵出现“带状空隙”或者缺失某些 rho原因局部极大值检测条件有误。比如我用dx 0和dx(2:end) 0来判峰值但dx里的负数可能是接近零的浮点误差导致把平台误判为峰值或者漏掉极值。解决给判定条件加一个小阈值例如dx(1:end-1) 1e-8且dx(2:end) -1e-8同时要求x本身大于一个阈值比如 0.5避免把噪声当作峰值。这是我在实际中踩过最大的坑——不加阈值时周期性波形上每个数字噪声都会被当成“局部极大值”导致分叉图看起来像两条粗线。5.5 现象五别人解压你的 rar 后点运行报错找不到函数或工具箱原因代码里用了findpeaks、diff等工具箱或者路径里有中文名。findpeaks属于 Signal Processing Toolbox并不是每个人都有。解决尽量不要依赖工具箱自写差分极大值检测。另外把主脚本和函数文件放在同一层目录不要用addpath加上层目录因为路径移动后容易失效。打包成 rar 时文件名不要用中文和空格比如lorenz_logistic_bifurcation.m这种命名更稳。发布前用一台干净的 Matlab 环境跑一遍确认不是只有你的电脑能运行。6. 把脚本整理成可复用的 .m 并打包目录结构与导出高清图最后分享一个我习惯的收尾动作。开发完分叉图代码后不要只扔一个孤零零的脚本至少要有一个main.m入口一个functions子目录放公共函数一个figs目录放导出的图。例如lorenz_logistic/ ├── main.m ├── functions/ │ ├── logistic_bifurcation.m │ ├── lorenz_phase.m │ └── lorenz_bifurcation.m └── figs/main.m里可以写一个简单的菜单让使用者选择画哪张图。导出图片用exportgraphics或print指定 300 dpi 以上避免在论文里模糊。我通常用print(gcf, figs/logistic_bifurcation, -dpng, -r300)。如果你希望别人直接复现可以在代码开头加rng(42)固定随机种子这样初始值如果用了随机数每次运行结果都一样。另外给脚本加注释时别偷懒。每个函数的第一行写上“输入参数、输出参数、依赖项、参考来源”不然三个月后你自己都看不懂。这就是我从翻车中总结出来的最后一条血泪经验代码能跑只是开始能被人拿来就跑才是发布的意义。希望这些细枝末节能帮到你尤其是那些在rho和r之间反复横跳的朋友。愿你画出的每一条分叉线都清晰如初。本文还有配套的精品资源点击获取