用Matlab绘制Lorenz吸引子：相图、庞加莱截面与分岔图全攻略

第一次在论文里看到 Lorenz 吸引子那双蝴蝶翅膀时大多数人都会被它的几何美感抓住。但真到了自己动手用 Matlab 画图的时候会发现画出一条好看的相图只是起点。做非线性动力学研究或者自己搭了一个混沌系统想验证它的行为通常需要把三维相图、二维相图、庞加莱截面图和分岔图放在一起看相图告诉你吸引子长什么样庞加莱截面告诉你轨迹在某个截面上的“离散节奏”分岔图则直接扫描参数告诉你在哪个区间进入混沌、哪个区间又回到周期。这篇文章就是我整理好的一套 Matlab 程序合集以 Lorenz 这个经典的三阶微分方程系统为例把四种图形全部跑通代码可以直接复制改成你自己的系统。适合刚接触混沌、需要给论文补图以及想验证自己写的微分方程是否混沌的同学。1. 为什么画了相图还要画分岔图四种图形各回答什么问题很多人第一次接触混沌系统时以为只要能画出那条蝴蝶曲线就算搞定了。实际上相图只是“看起来像混沌”它不能精确告诉你系统到底是周期、拟周期还是混沌。我个人的习惯是先画三维相图和二维投影图建立直觉再用庞加莱截面判断运动形态最后用分岔图扫描参数摸清全局。这四种图各有分工谁也不能替代谁。1.1 四种图形工具的分工三维相图和二维相图回答的是“吸引子长什么样”轨迹在相空间里怎么走是被压缩到某个平面附近还是真的在三维空间里铺开。庞加莱截面图回答的是“轨迹穿过某个平面时留下的点有什么规律”如果只有有限个点对应周期运动如果形成一条闭合曲线对应拟周期如果是一团带有分形结构的点云才是混沌。分岔图回答的是“某个参数变化时系统行为怎么切换”稳定点、极限环、周期倍化、混沌窗口全部能在一张图上看到。图形核心做法直观信息典型判定三维相图直接绘制 x-y-z 空间轨迹吸引子的几何形态、翅膀结构轨迹是否被限制在低维流形二维相图选择坐标面投影轨迹疏密、对称性、翼的分支投影是否出现折叠与拉伸庞加莱截面记录轨迹穿过某平面的点截点离散分布有限点/闭合曲线/分形点云分岔图扫描参数并记录稳态特征量全局分岔路径倍周期分支、混沌带、周期窗口这四种图是互补的。我见过不少同学只画了三维相图就断定自己的系统是混沌结果拿去算 Lyapunov 指数发现最大指数为负——问题就出在“看起来乱”不等于“混沌”。1.2 为什么是三阶微分方程系统先解释一个很多人刚学时容易绕进去的点标题里说的“三阶微分方程系统”指的是一个连续时间自治系统里有三个状态变量也就是三个一阶常微分方程组成的方程组。Lorenz、Rossler、Chen 这些经典混沌系统全部是这种形式。这背后有一个基本结论连续时间系统要出现混沌相空间维数至少是 3。二维平面上的连续自治系统可以被 Poincaré-Bendixson 定理限制住要么趋于平衡点、要么趋于极限环不允许出现奇怪吸引子。所以当你打算验证一个“三阶微分方程系统”是否混沌时最低要求就是在三维相空间里看轨迹。如果你手头是一个四阶或者更高阶系统那画图方法完全一样只是 Lyapunov 谱的计算会更复杂一些。本文以 Lorenz 为样例系统但整套绘图流程对所有高阶系统通用。2. Lorenz 系统建模与求解先把轨迹算对画图的前提是数值积分要准。混沌系统对初始误差和数值耗散极其敏感所以求解器参数这块值得多花两分钟。2.1 方程写进 Matlab三行就够Lorenz 系统的标准形式如下dx/dt sigma * (y - x) dy/dt x * (rho - z) - y dz/dt x * y - beta * z经典参数是sigma 10、rho 28、beta 8/3这个组合下系统处于混沌状态。写成 Matlab 的微分方程函数就是一个文件末尾的一个子函数function dydt lorenz_sys(~, y, sigma, rho, beta) % y 是长度 3 的状态向量 [x; y; z] dydt zeros(3,1); dydt(1) sigma * (y(2) - y(1)); dydt(2) y(1) * (rho - y(3)) - y(2); dydt(3) y(1) * y(2) - beta * y(3); end注意上面的入参里第一个~是时间 t因为 Lorenz 系统是自治系统方程里不显含时间但ode45调用时仍需保持这个接口。2.2 求解器参数对混沌系统不能心慈手软数值积分混沌系统最常见的问题是容差设置太宽松。ode45的默认相对容差是1e-3对大多数工程问题够用但对混沌系统这个精度会让轨迹在几秒后明显偏离真实解甚至导致吸引子形态变形。我自己习惯把RelTol和AbsTol都设到1e-8积分时长较长的时候再往下调一档到1e-10。sigma 10; rho 28; beta 8/3; x0 [1; 0; 0]; % 初始状态尽量取在吸引域内 tspan [0 100]; % 积分 100 秒 opts odeset(RelTol, 1e-8, AbsTol, 1e-8); [t, Y] ode45((t, y) lorenz_sys(t, y, sigma, rho, beta), ... tspan, x0, opts); % 看一眼波形是否合理 figure; plot(t, Y(:,1), LineWidth, 0.8); xlabel(t); ylabel(x); title(Lorenz 时间序列);这段代码跑完后Y的三列分别就是 x、y、z 的时间序列。后面所有图形都以这里的计算结果为数据源。注意混沌系统对初值敏感。如果你改了初始条件后面的庞加莱截面和相图形态不会变但轨迹的“相位”会不同这是混沌的固有性质不是代码 bug。3. 三维相图和二维相图把吸引子画出来算好轨迹之后绘图本身很简单但有几个细节直接决定图好不好看、能不能放进论文。3.1 剔除瞬态是关键一步ode45从x0 [1;0;0]出发后轨迹会先经过一段“瞬态过程”然后才被吸引到 Lorenz 吸引子上。如果直接把全部轨迹画出来瞬态部分会在图中拉出一条长长的尾巴把吸引子的精细结构盖住。所以我在画相图之前通常只保留t 5或者t 10之后的数据idx_ss t 5; % 去掉前 5 秒的瞬态 Yss Y(idx_ss, :);对于 Lorenz 这种收敛较快的系统5 秒足够如果换了别的系统收敛时间不确定可以先画出时间序列看振幅是否进入稳定波动再决定截断点。3.2 二维投影不同角度看蝴蝶Lorenz 吸引子最经典的视角是三维图但二维投影能更清楚地显示结构。x-y 投影能看到绕两个中心点的盘旋x-z 投影是那张最著名的“蝴蝶脸”y-z 投影则能看出翅膀的厚度。我的习惯是把三个投影放在同一张 figure 的子图里figure; subplot(1,3,1); plot(Yss(:,1), Yss(:,2), Color, [0.2 0.4 0.8], LineWidth, 0.3); xlabel(x); ylabel(y); title(x-y 投影); axis equal; subplot(1,3,2); plot(Yss(:,1), Yss(:,3), Color, [0.8 0.3 0.2], LineWidth, 0.3); xlabel(x); ylabel(z); title(x-z 投影); axis equal; subplot(1,3,3); plot(Yss(:,2), Yss(:,3), Color, [0.2 0.7 0.3], LineWidth, 0.3); xlabel(y); ylabel(z); title(y-z 投影); axis equal;axis equal这一步容易被忽略。不加的话Matlab 会自动拉伸坐标轴蝴蝶会被压扁或者拉长视觉上完全失真。3.3 让图形更科研的细节处理三维相图用plot3直接画连续线条就够但线条宽度不宜太大否则轨迹密集区域会糊成一团。我一般设LineWidth为 0.3 到 0.5颜色用偏暗的蓝或紫比默认的黄色好看很多。如果需要展示轨迹在吸引子上随时间推进的流向可以把时间作为颜色维度用scatter3或者patch做渐变效果更直观figure; t_ss t(idx_ss); scatter3(Yss(:,1), Yss(:,2), Yss(:,3), 1, t_ss, .); xlabel(x); ylabel(y); zlabel(z); title(Lorenz 三维相图颜色随时间变化); colormap(jet); colorbar; view(3); grid on;还有一个经验画三维相图时view角度默认是(-37.5, 30)但 Lorenz 吸引子在view(-45, 20)或者view(120, 25)下更容易看到双翼的分离结构。多转几个角度截图挑一张信息量最大的放论文里。4. 庞加莱截面图把连续流打回离散映射庞加莱截面是判断混沌最直观的工具之一。它的原理不复杂在相空间中选一个横截面每当轨迹穿过这个截面时记录下交点坐标。连续系统被这样“采样”后变成离散映射原本难以分析的连续流变成了平面上的一堆点。4.1 截面图的数学原理和选择技巧对 Lorenz 系统最常用的截面之一是z rho - 1这个平面因为这个位置大致位于两个翼盘旋中心的中间高度。也可以用固定值比如z 20效果差不多。选择截面时要注意两点一是截面不能与轨迹运动方向相切否则交点会非常稀疏甚至没法连续记录二是尽量避开系统的对称平面否则交点会大量重合在一条线上看不出结构。4.2 方法一轨迹数据后处理插值最简单的实现方式不需要额外设置 ODE 选项直接在后处理时从 (Y) 矩阵里找穿越点。思路是找到相邻两步中 z 分量跨越截面高度 (z_0) 的索引然后用线性插值估计交点坐标z0 rho - 1; crossIdx find(diff(sign(Y(:,3) - z0)) ~ 0); xCross zeros(size(crossIdx)); yCross zeros(size(crossIdx)); for k 1:length(crossIdx) i crossIdx(k); d z0 - Y(i, 3); % 截面高度与当前步的差 w d / (Y(i1, 3) - Y(i, 3)); % 归一化权重 xCross(k) Y(i,1) w * (Y(i1,1) - Y(i,1)); yCross(k) Y(i,2) w * (Y(i1,2) - Y(i,2)); end figure; plot(xCross, yCross, ., MarkerSize, 4); xlabel(x); ylabel(y); title([庞加莱截面 z num2str(z0)]);这种方法在ode45步长足够密时完全够用误差很小。如果你积分时长很长、步长很稀疏线性插值会丢失精度那就需要事件检测。4.3 方法二ode45 事件检测ode45自带Events选项能够在积分过程中精确定位轨迹穿过某个平面的时刻和状态。这样做最大的好处是不管步长多大穿越点都不会被漏掉且坐标精度更高适合需要长程庞加莱截面或者后续做统计分析的情况。z0 rho - 1; opts odeset(RelTol, 1e-10, AbsTol, 1e-10, ... Events, (t, y) poincare_events(t, y, z0)); [t, Y, te, Ye] ode45((t, y) lorenz_sys(t, y, sigma, rho, beta), ... tspan, x0, opts); figure; plot(Ye(:,1), Ye(:,2), ., MarkerSize, 4); xlabel(x); ylabel(y); title([庞加莱截面 z num2str(z0) 事件检测]); grid on; function [value, isterminal, direction] poincare_events(~, y, z0) value y(3) - z0; % 穿越 z z0 平面 isterminal 0; % 不终止积分 direction 0; % 正方向和负方向穿越都记录 end这里的Ye是事件点的状态矩阵每一行对应一次穿越截面时的 ([x, y, z])。direction 0表示双向穿越都记录如果只想记录从下往上穿越改成direction 1。4.4 庞加莱截面上怎么判断混沌我最初看庞加莱截面时最困惑的问题就是什么样子算混沌经验法则如下。有限个孤立点周期运动轨迹一遍又一遍穿过同一位置。一条光滑闭合曲线拟周期运动截点连成环。一团不可数、带有自相似结构的点云混沌。Lorenz 系统在经典参数下庞加莱截面会呈现两簇对称分布的点云每一簇都像拉伸折叠后留下的细密条纹这基本就是混沌无疑。5. 分岔图扫描参数看系统“变脸”分岔图和前面三类图完全不是一个量级的东西因为它不是画一条轨迹而是把整个参数轴上每个取值对应的稳态行为压缩到一张图里。对 Lorenz 系统来说最经典的扫参对象是 (\rho)也就是 Rayleigh 数相关的那个参数。5.1 核心思路每个参数值只保留稳态行为分岔图的基本逻辑是固定其他参数改变 (\rho)对每个 (\rho) 做一次长时间积分。初期的瞬态完全丢弃只记录稳态后的特征量。这个特征量有两种常见取法一是记录穿截面的交点坐标二是记录某个状态变量的局部极大值。对 Lorenz 系统我推荐取 x 分量的局部极大值因为 x 的峰值序列在混沌区会形成清晰的两支带便于观察周期窗口。选择稳态后的局部极值有个好处不需要额外定义截面代码更加通用。换个新系统的时候只要把微分方程函数替换掉就能直接跑分岔图。5.2 扫参策略和延拓技巧直接对每个 (\rho) 都从同一个初始条件[1;0;0]开始积分代码简单但在某些参数区间尤其是临界点附近瞬态会很长积分不够久的话图上会出现假点。以 Neumann 边界条件那种情况来对比用一个从上一个参数状态延续下来的“延拓法”会稳很多。延拓法的意思是把上一个 (\rho) 积分结束时的状态作为下一个 (\rho) 的初始状态因为吸引子随参数连续变化这样初值已经贴近吸引子瞬态极短。下面是基于延拓的局部极值法分岔图代码% bifurcation_lorenz.m clear; clc; close all; sigma 10; beta 8/3; rhoList 10:0.1:50; figure; hold on; state [1; 0; 0]; % 初始状态逐步延拓 for rho rhoList f (t, y) [sigma * (y(2) - y(1)); y(1) * (rho - y(3)) - y(2); y(1) * y(2) - beta * y(3)]; [~, Y] ode45(f, [0 100], state, odeset(RelTol, 1e-8, AbsTol, 1e-8)); idx_ss t_ss 50; % 前 50 秒当瞬态丢弃 x_ss Y(idx_ss, 1); % 找局部极大值差分符号从 变 - 的位置 dx diff(x_ss); signChange diff(sign(dx)); peakIdx find(signChange 0) 1; peaks x_ss(peakIdx); if isempty(peaks) continue; end plot(rho * ones(size(peaks)), peaks, ., MarkerSize, 1, ... Color, [0.1 0.3 0.8]); % 更新初始状态使用当前参数下最后一次积分的末状态 state Y(end, :); end xlabel(\rho); ylabel(x); title(Lorenz 系统分岔图);这个例子里的rhoList从 10 扫到 50步长 0.1。想观察更细的结构比如周期窗口内的倍周期分岔可以把步长缩小到0.01但计算量会成倍增加。parfor并行可以缓解但延拓法本身是串行的并行化需要改成每个 (\rho) 独立积分更适合参数网格很大的情况。5.3 怎么判断分岔图分岔图读起来很直观横轴是参数纵轴是稳态特征量。图上每一个点表示“在这个参数下系统稳定后经过的所有峰值位置”。看图的要点是一条水平线系统收敛到稳定平衡点峰值只有一个。一个位置裂成两个点发生倍周期分岔系统出现二周期振荡。不断分裂成多个点再转入一条竖带系统沿倍周期路径进入混沌。竖带中间突然出现一条细缝里面只剩几个点混沌突变为周期窗口。Lorenz 系统在 (\rho \approx 24.74) 附近发生 Hopf 分岔从稳定点变成极限环随后经历倍周期分岔在 (\rho \approx 28) 附近进入混沌。把上面代码的扫参范围改成 0 到 50你会发现低参数区很简单高参数区则呈现出非常丰富的混沌带和周期窗口。这是我每次验证新系统都先跑一遍分岔图的原因——它用一张二维图就把系统的“性格”摸清了。6. 我踩过的坑和调参心得这部分是我最想分享的内容。这些坑不是从教科书上看来的是实打实跑代码跑出来的。6.1 容差造成的“伪混沌”最早我给自定义系统画分岔图时用默认容差1e-3跑结果在高参数区出现了一大片看起来像混沌的散点但局部放大后发现其实是一堆锯齿状的数值振荡。后来把AbsTol调低到1e-8那些假散点全部消失。记住混沌系统对误差极其敏感宽容差下的“混沌”可能只是数值噪声。画图前先固定一组参数对比不同容差下吸引子是否重合是最简单的验证。6.2 瞬态剔除长度怎么定我见过有人设置t 200只取后面 10 秒的数据结果窗口太短峰值样本不足分岔图变成稀疏的几个点。也有人只踢掉前 1 秒瞬态尾巴很长吸引子还没收敛就画进去图上多出一条“拖尾”。稳妥做法是先画该参数下的时间序列看振幅和局部结构何时稳定再决定剔除长度。或者更简单剔除掉总时长前 20%然后剩下 80% 用于统计。6.3 分岔图总花屏多半是步长和积分时长的锅分岔图最容易出的问题有三个参数步长太大、积分时间不够、特征量取法不当。参数步长太大会跳过窄的周期窗口。Lorenz 系统在 (\rho 30) 附近就有极窄的周期窗粗步长完全看不见。积分时间不够会导致点落在瞬态过程中图上有大片杂乱无章的点。特征量取法不当则表现为图上所有点挤成一条粗带无法分辨倍周期结构。我一般先用步长 0.2 快速扫一遍全局确认混沌区位置后再对感兴趣的区间用步长 0.01 加密。6.4 别把数值振荡当混沌最后一类问题来自方程组本身。有些系统在部分参数区间刚性很强ode45不适合会解出高频伪振荡。这时画相图会看到一团实心圆盘分岔图更是一整片均匀噪声。解决办法是改用ode15s或者把时间跨度缩小观察。我刚学混沌时在这一步浪费了不少时间后来养成了交叉验证的习惯同一参数下分别用ode45和ode15s算一遍如果吸引子形态差异明显立刻怀疑是求解器不合适。现象可能原因检查方法相图有拖尾长线瞬态未剔除画时间序列确认收敛点庞加莱截面全是散乱点参数在混沌区或容差过松降低容差重算分岔图边界粗糙杂乱参数步长过大或积分时长不足缩小步长、加长积分时间低参数区出现密集振荡方程组刚性数值伪振荡换 ode15s 对比我自己在跑完 Lorenz 的整套图之后形成了一套固定的验证流程先画三维相图确认吸引子几何形态再利用庞加莱截面判断运动类型最后用分岔图扫描参数全局。这三个工具组合起来大约能覆盖 80% 的定性判断需求。如果你想进一步定量证明系统是混沌那就要上最大 Lyapunov 指数了但图形工具依然是定位参数区间的第一站。上面的代码全部可以直接复制到你的脚本里把lorenz_sys替换成你自己的三阶方程图形部分基本不用改。最后再说一个小技巧所有图完成后记得用exportgraphics(gcf, filename.png, Resolution, 300)导出高分辨率图片论文排版时比截图清晰得多。混沌系统图形最忌讳的就是细节糊成一片高分辨率输出能保住那些细如发丝的拉伸折叠结构。