做非线性动力学研究或者上过混沌相关课程的人大概率都逃不过亲手画一张 Lorenz 系统相图的命运。早期我用 Matlab 画这些图的时候最头疼的事情是全网搜到的代码要么只给一段孤零零的脚本要么是截图里才看得到结果根本不知道参数怎么调、截面怎么取、暂态怎么去掉。这次我把自己的实现整理成一套完整的程序分享出来覆盖三阶微分方程混沌系统的庞加莱截面图、二维相图、三维相图以及分岔图基于最常见的 Lorenz 系统代码可以直接跑关键位置都有注释和参数解释。这篇文章适合正在做混沌系统数值分析的硕博生、做非线性信号处理的工程师以及准备用 Matlab 完成相关课程设计的同学。我会把每个图形的数学含义、Matlab 实现原理、以及我在实际调试中踩过的坑一起讲清楚确保你不只是把代码复制走而是能根据自己的系统改参数、换截面、调出符合预期的结果。1. 混沌系统与四类图形先搞清楚要画什么1.1 为什么是三阶微分方程在开始写代码之前有必要先把一个基本问题说透为什么混沌系统的入门研究对象几乎都是三阶自治微分方程从动力系统理论看连续时间系统要产生混沌相空间维度至少是三维。二维连续系统的轨迹受限于平面Poincaré-Bendixson 定理直接排除了混沌的可能性——平面上要么收敛到平衡点要么形成极限环不存在“既不自洽又不发散”的奇怪吸引子。因此无论是 Lorenz 系统、Rossler 系统还是 Chen 系统都是三个状态变量、三个一阶微分方程构成的耦合系统本质上是一个三阶系统。用 Matlab 求解这类问题的核心思路非常直接把三阶微分方程组写成状态向量的形式交给积分器去数值求解。我们以最经典的 Lorenz 系统为例方程形式如下dx/dt sigma * (y - x) dy/dt x * (rho - z) - y dz/dt x * y - beta * z三个参数中sigma 是 Prandtl 数rho 是 Rayleigh 数beta 是几何参数。经典的混沌参数组合是 sigma10beta8/3rho28。这个参数下系统表现出典型的蝴蝶形奇怪吸引子四个图形都能产生非常有辨识度的结果。我建议所有初学者先用这套参数跑通流程再去尝试其他参数组合。1.2 四类图分别回答什么问题很多同学拿到题目时会产生一个困惑这几种图到底有什么区别为什么要画这么多张我的理解是这样的它们分别从不同尺度回答“系统在干什么”这个问题。二维相图回答的是“两个状态变量之间的关系是什么”它把三维轨迹投影到某个坐标平面上适合做初步观察。三维相图展示系统在相空间中的完整轨迹能够直观看到奇怪吸引子的几何结构。庞加莱截面回答的是“轨迹如何穿过某个横截面”它把连续流转化为离散映射通过截面上的点集分布判断系统是周期、拟周期还是混沌。分岔图回答的是“系统行为如何随参数变化”它是参数扫描后的全局视图能看出从稳定到周期再到混沌的完整演化路径。这四类图的组合关系可以理解成相图是“单张照片”庞加莱截面是“某个特定角度的切片扫描”分岔图则是“一段视频的摘要”。后面我会按照这个逻辑把每张图的实现思路串起来。2. 环境准备与微分方程系统搭建2.1 Matlab 环境与基础设置本文的代码完全不依赖额外工具箱只需要一个基础 Matlab 环境就能运行。我自己的测试环境是 Matlab R2023b但 R2020 之后的版本应该都能直接跑通。需要提醒的是如果你的电脑上安装的是较新版本遇到启动闪退或者 license 激活异常这类问题优先检查系统环境变量和许可证文件路径这属于安装问题和绘图代码本身无关。网上搜到的所谓“2026b 下载安装教程”很多带有捆绑内容我不建议从非官方渠道获取安装包。建议在跑代码前先执行一次干净的clear; close all; clc;避免工作区残留变量干扰。数值求解混沌系统对初始条件极其敏感如果之前实验留下了某些同名变量很可能导致轨迹走向完全不同的分支。这是我在实践中踩过的真实坑有一次分岔图跑出来的结果与理论完全不符排查了半天发现是工作区里残留了一个旧版本的 rho 标量导致循环内参数覆盖出错。2.2 定义 Lorenz 系统的微分方程函数我们先把 Lorenz 系统写成一个独立的函数文件。这个函数接收时间和状态向量返回导数向量。注意这里的顺序必须是(t, y)即使方程中不显式含时间 t也要保留这个占位参数否则 ode45 会报错。function dydt lorenz_system(t, y, sigma, rho, beta) % Lorenz 系统微分方程 % y [x; y; z] x y(1); y_state y(2); z y(3); dydt zeros(3, 1); dydt(1) sigma * (y_state - x); dydt(2) x * (rho - z) - y_state; dydt(3) x * y_state - beta * z; end这里我用了y_state作为中间变量是因为主脚本里通常还需要用y作为整个状态向量的变量名两者容易混淆。函数本身很简单但有一个细节必须强调三个方程之间的耦合关系不能写错特别是第二个方程中的- y_state这一项它是系统耗散性的关键组成部分。如果这里漏掉或者写成加号系统行为会完全改变甚至可能发散。主脚本中通过匿名函数把参数传递给微分方程函数sigma 10; beta 8/3; rho 28; f (t, y) lorenz_system(t, y, sigma, rho, beta);这种做法的好处是后续做分岔图扫描 rho 时只需要在循环体内重新定义匿名函数不需要反复修改子函数文件。3. 相图的绘制二维与三维轨迹的实现3.1 数值积分的参数选择调用 ode45 求解时最影响结果质量和计算速度的是时间跨度、初始条件和求解器容差。不同参考书给出的初始条件略有差异但只要在吸引子盆地里最终轨迹经过暂态后都会收敛到同一个吸引子上。tspan [0 100]; y0 [1; 1; 1]; options odeset(RelTol, 1e-6, AbsTol, 1e-8); [t, y] ode45(f, tspan, y0, options);我习惯把RelTol设为 1e-6AbsTol设为 1e-8。这两个参数控制自适应步长的精度设置太松会导致轨迹发散或出现明显误差累积设置太紧则计算时间显著增加。对 Lorenz 系统来说上面的设置已经能保证轨迹在 100 个时间单位内保持数值稳定。需要特别注意的是tspan的前面一部分是暂态过程。无论初始点选在哪里轨迹都需要一段时间才能落到吸引子上。如果初始点距离吸引子很远这段暂态可能会很长。画相图时我通常会丢弃前 10% 的数据点只绘制收敛后的轨迹n_skip round(length(t) * 0.1); x y(n_skip:end, 1); y_state y(n_skip:end, 2); z y(n_skip:end, 3);这不是强迫症而是为了避免初始点远离吸引子时产生的“飞行路径”污染图形。你可以在自己机器上试一下把初始点改成[100; 100; 100]如果不丢弃暂态相图边缘会出现一条很突兀的长线严重影响可读性。3.2 二维相图三个平面的投影二维相图本质上就是把三维轨迹投影到坐标平面。对 Lorenz 系统最常见的三个投影平面是 x-y、x-z 和 y-z。其中 x-z 平面的蝴蝶形最经典x-y 平面看起来像一个“猫头鹰脸”y-z 平面的结构相对模糊。绘制代码非常简单figure(Color, white); subplot(1, 3, 1); plot(x, y_state, LineWidth, 0.5); xlabel(x); ylabel(y); title(x-y 相图); axis equal; subplot(1, 3, 2); plot(x, z, LineWidth, 0.5); xlabel(x); ylabel(z); title(x-z 相图); subplot(1, 3, 3); plot(y_state, z, LineWidth, 0.5); xlabel(y); ylabel(z); title(y-z 相图);有两点值得说明。第一axis equal不是必须的但对某些平面加上后能保持几何比例不失真。第二LineWidth不宜设得太粗因为混沌轨迹会反复折叠线宽过粗会把结构细节糊在一起。实际输出如果线条密集程度太高可以适当增加 tspan 结束时间轨迹更充分覆盖吸引子出图效果会更饱满。3.3 三维相图视觉层次优化三维相图是展示奇怪吸引子几何结构的最佳方式。Lorenz 系统的吸引子由两个围绕不稳定平衡点的“翅膀”构成轨迹在两侧之间非周期切换。直接plot3画出来的图虽然能看但整条线颜色单一看不出轨迹的运动方向和疏密分布。我的做法是引入颜色渐变用颜色反映时间演化figure(Color, white); t_plot t(n_skip:end); c t_plot; % 用时间作为颜色映射的依据 scatter3(x, y_state, z, 2, c, filled); colormap(jet); colorbar; xlabel(x); ylabel(y); zlabel(z); title(Lorenz 系统三维相图 (颜色表示时间演化)); grid on; view([30, 20]);这里我用了scatter3而不是plot3区别在于散点图可以精确控制每个点的颜色。点到数通常有几万个scatter3的性能完全可以接受但如果你用很旧的 Matlab 版本比较大的数组绘制时可能出现卡顿。这时候可以把点抽稀例如每 3 个点取 1 个再画。另一个提升观感的手段是旋转视角。view函数的两个参数分别是方位角和仰角我常用的角度组合是[30, 20]能看到蝴蝶形状的正面结构改成[0, 90]会得到俯视图这时候吸引子的两瓣结构非常清晰。多保存几个视角的图在论文中使用时可以根据需要选择。4. 庞加莱截面图的实现与判断4.1 截面的数学定义与选择思路庞加莱截面的思想比较直观在一个周期轨道上取一个横截面记录轨迹每一次穿过该截面的位置用这些离散点来研究系统的长期行为。对连续系统而言这是一个降维分析工具把连续流转化为离散映射。对 Lorenz 系统来说最常用的截面选取方式有两种。第一种是选取某个固定平面比如 z rho - 1 平面约等于 27这是系统的两个不稳定平衡点所在的高度附近轨迹穿过该平面时的 (x, y) 坐标可以揭示系统的折叠结构。第二种是选取局部极值面即记录 z 取局部极大值时对应的 (x, y) 坐标这种方法的优点是截面点分布均匀常用于绘制分岔图。两种方法本质上都是构造一个映射区别在于截面在相空间中的位置和方向。我在实际中更喜欢用“z 取局部极大值”的方式原因有两点一是它不需要额外求解超平面穿越实现简单二是后续画分岔图时正好需要记录局部极值点一套代码可以复用。4.2 代码实现检测局部极值点使用findpeaks函数检测 z 序列的峰值是 Matlab 中最直接的做法。findpeaks来自 Signal Processing Toolbox如果没有这个工具箱也可以用diff手动判别当diff(z_traj)由正变负时即出现局部极大值。% 方法一使用 findpeaks [~, locs] findpeaks(z); x_poincare x(locs); y_poincare y_state(locs); % 方法二手动差分判别无工具箱版本 dz diff(z); peak_idx find(dz(1:end-1) 0 dz(2:end) 0) 1; x_poincare2 x(peak_idx); y_poincare2 y_state(peak_idx);使用findpeaks时要注意一个细节它会检测所有局部峰值但如果轨迹在某个区域内出现高频小幅抖动可能产生大量意义不大的点。解决办法是设置MinPeakHeight或者MinPeakDistance作为阈值。我没有在代码里刻意加阈值因为标准 Lorenz 系统的 z 轨迹比较平滑峰值点每个周期一个稀疏度和信噪比都较好。绘制庞加莱截面图时我习惯把三个视角都展示出来x-y 平面、x-z 平面以及 y-z 平面上的投影点。figure(Color, white); subplot(1, 3, 1); plot(x_poincare, y_poincare, ., MarkerSize, 6); xlabel(x); ylabel(y); title(庞加莱截面投影 (x-y)); grid on; subplot(1, 3, 2); plot(x_poincare, z(locs), ., MarkerSize, 6); xlabel(x); ylabel(z); title(庞加莱截面投影 (x-z)); subplot(1, 3, 3); plot(y_poincare, z(locs), ., MarkerSize, 6); xlabel(y); ylabel(z); title(庞加莱截面投影 (y-z));从结果上看Lorenz 系统在 rho28 时的庞加莱截面图呈现类似“两条弧线”的结构而不是有限的几个点也不是封闭连续的曲线。这正是奇怪吸引子的典型特征截面点集看起来像一条被无限折叠的一维曲线具有分形结构。如果你的结果是一堆杂乱无章的稀疏点通常是 tspan 太短导致轨迹还没有充分覆盖吸引子如果你的结果是几个孤立的点那说明系统很可能处于周期状态需要检查参数是否正确。4.3 如何用截面图判断系统状态庞加莱截面图的价值不在于好看而在于它能提供一个清晰的状态判据。如果截面图上只有 1 个点对应的是周期 1 轨道2 个点对应周期 2 轨道有限个点对应周期轨道一条封闭曲线对应拟周期运动而具有分形结构的密集点集则对应混沌。这个经验法则适用于绝大多数连续动力系统不仅限于 Lorenz。我强烈建议在课程报告或论文中放一组对比图rho24 时的庞加莱截面单点、rho27.4 时的截面可能随初值呈现层叠结构、rho28 时的截面分形点集。这组对比比单张图有说服力得多也能直观展示系统状态的转变。5. 分岔图的绘制参数扫描与可视化5.1 分岔图的生成逻辑分岔图是由“局部极值采样”和“参数扫描”两个环节组成的。以 Lorenz 系统为例固定 sigma 和 beta让 rho 在指定范围内连续变化对每个 rho 做一次数值积分舍弃暂态后记录轨迹的局部极大值点然后把所有 rho 对应的极值点画在同一张图上就得到分岔图。它可以清晰展示系统从收敛到周期、再到混沌、再进入周期窗口的完整演化路径。这里有一个容易翻车的地方“对每个 rho 做一次积分”看起来简单但积分时间长度、暂态丢弃长度、极值采样策略都会对图形质量产生巨大影响。我调整参数后得到的核心代码段如下sigma 10; beta 8/3; rho_range 24:0.1:40; n_rho length(rho_range); tspan [0 200]; y0 [1; 1; 1]; options odeset(RelTol, 1e-6, AbsTol, 1e-8); rho_plot []; extreme_plot []; for k 1:n_rho rho rho_range(k); f (t, y) lorenz_system(t, y, sigma, rho, beta); [t, y] ode45(f, tspan, y0, options); z_traj y(:, 3); n_skip round(length(t) * 0.3); % 丢弃前 30% 暂态 z_traj z_traj(n_skip:end); [~, locs] findpeaks(z_traj); if isempty(locs) continue; end x_extremes y(locs n_skip, 1); % 注意索引补偿 n_pts length(x_extremes); rho_plot [rho_plot; rho * ones(n_pts, 1)]; extreme_plot [extreme_plot; x_extremes]; end figure(Color, white); plot(rho_plot, extreme_plot, ., MarkerSize, 3); xlabel(\rho); ylabel(z 局部极值对应的 x 值); title(Lorenz 系统分岔图 (24 \rho 40));这段代码中有三个细节需要重点说明。第一个细节是索引补偿。findpeaks返回的locs是在z_traj中的位置而z_traj是从y(:,3)中截断了一部分后得到的因此真正对应的原始数据索引是locs n_skip。我在第一次写这段代码时忽略了这一点导致取到的 x 值整体偏移画出来的分岔图完全是错的。第二个细节是暂态丢弃比例。这里我丢弃了前 30% 而不是前 10%原因在于分岔图对暂态极其敏感。如果暂态没有被充分去除在那些本应呈现单周期或双周期的参数区间图像上会产生大量不属于吸引子的杂散点图形看起来“脏”。实际调试时我发现20% 到 30% 的丢弃比例是比较稳妥的选择丢弃太少会脏太多则会不必要地浪费计算时间。第三个细节是逐点拼接 vs 预分配。上面代码用了动态拼接的方式优点是直观且能在 rho 点数较少时保持代码简洁。但如果你要提高效率把rho_plot和extreme_plot预先分配成一个大矩阵然后按索引填充速度可以提升不少。下面会专门讲性能优化。5.2 计算成本与性能优化分岔图最让人头疼的问题是计算时间。rho 从 24 扫到 40步长取 0.1总共 161 次积分每次积分 tspan200。在我自己的机器上这段代码大约需要 5 到 10 分钟取决于 Matlab 版本和 CPU 性能。如果你只是临时看一眼趋势这个时间可以接受但如果要生成高分辨率分岔图比如步长 0.01就必须考虑优化。我试过几种可行方案按推荐顺序排列使用parfor并行循环。前提是安装了 Parallel Computing Toolbox并且用parpool开启并行池。由于每个 rho 的积分是互相独立的这个场景是并行计算的标准适用场景。并行化之后耗时可以降到原来的三分之一到四分之一。缩短 tspan 和暂态比例。tspan 从 200 降到 100暂态比例从 30% 降到 20%大部分情况下图形结构不会发生肉眼可见的变化但时间直接减半。缺点是某些接近临界点的参数区间可能出现暂态残留需要肉眼检查。使用ode45时放宽容差。把RelTol从 1e-6 改为 1e-4计算速度会有明显提升。对绘制分岔图而言精度足够但如果你还要拿同一批数据计算 Lyapunov 指数或做定量分析这个做法就不合适了。把连续系统降维为离散映射。对 Lorenz 系统理论上可以预先算出每个极值点到下一个极值点的映射关系用迭代代替数值积分。这种方法速度快但不通用换其他系统就得重新推导。我不建议初学者在这个方向投入太多时间。5.3 如何阅读分岔图绘制完成后分岔图的信息量是很大的。以 Lorenz 系统为例rho 从 24 开始系统处于一个稳定的平衡点分岔图上表现为单个值。随着 rho 超过约 24.74出现 Hopf 分岔系统进入周期振荡图像分裂为两条分支。继续增大 rho会出现倍周期分岔分支不断分裂最终在 rho28 附近进入混沌这就是图上看起来像“云团”的区域。在混沌区域中间还能看到若干周期窗口比如 rho30.1 附近和 rho37.8 附近这些窗口对应着系统短暂回到周期状态的行为。看懂这张图你的非线性动力学直觉会上升一个台阶。比如你在做实验信号分析时遇到奇怪的频谱就能联想到这背后可能对应着某个参数接近了分岔点。这也是我觉得花时间把分岔图画对、画好特别值得的原因。6. 常见问题与排查技巧实录6.1 图形出现异常发散或 NaN这是最多人遇到的第一个坑。表现是轨迹数值很快变成无穷大或 NaN相图上什么都看不到。排查顺序如下第一检查微分方程的公式是否正确。Lorenz 系统的三个方程互相耦合任何一个符号出错都会导致发散。第二检查参数输入顺序匿名函数(t, y) lorenz_system(t, y, sigma, rho, beta)和函数声明的参数顺序必须保持一致。第三检查初始条件是否过大。混沌系统对初始条件敏感但并不意味着初始条件大就发散只是某些特殊初始点可能落在吸引子盆地之外。标准做法是使用[1; 1; 1]或其他靠近原点的点。第四检查options容差是否设得太紧或太松。容差太松可能导致积分误差掩盖真实动力学太紧则可能导致求解器步长过小计算时间异常长。如果上面都检查过仍然 NaN尝试用显式 Euler 法做一个简单验证取一个极小的步长看看前几步轨迹是否大致合理。这能帮你快速把问题缩小到“方程写错”还是“积分器设置不对”。6.2 庞加莱截面图的点数稀疏或不连续如果你的截面图点非常少比如只有几十个点首先要检查 tspan 是否足够长。Lorenz 系统在 rho28 时轨迹在两个翅膀之间不规则跳跃需要相当长的时间才能充分覆盖吸引子。建议 tspan 至少取 200甚至可以取到 500。其次检查是否用了“局部极大值”作为截面如果是确认使用的是 z 或 x 变量而不是固定频率采样。用等时间间隔采样点并直接画出来得到的不是庞加莱截面而是轨迹的原始采样两者概念完全不同。另一个细节是findpeaks默认忽略两端不完整的峰值如果你发现截面点在首尾处明显缺失或分布不对称通常不是代码问题而是数据集不够长。多跑一段时间即可。6.3 分岔图在周期区间出现“毛刺”效应分岔图在周期区间应该是几条清晰的分支但有时候画出来分支旁边有淡淡的杂散点。这几乎都是暂态没有被完全去除导致的。解决方法是增加暂态丢弃比例或者直接增加 tspan。特别要注意的是在倍周期分岔临界点附近系统收敛到周期轨道的过程会非常缓慢称为临界慢化。比如 rho27.8 附近的参数区间暂态可能持续相当长的时间如果 tspan 不够分岔图上会出现类似混沌但实际上只是暂态的假结构。遇到这种情况我建议对特定参数局部加密扫描比如单独对 rho 从 27 到 28 步长取 0.01把每条的轨迹增长到 500 个时间单位仔细观察临界区域的变化。这样做不仅能确认毛刺是否为暂态还能帮助定位精确的临界参数值。6.4 常见问题速查表现象可能原因排查方法轨迹全部发散方程写错、参数顺序不对、初始点异常单步验证前三步的值相图上有突兀的直线暂态未丢弃增加n_skip比例庞加莱截面点太少tspan 太短延长积分时间到 200 以上庞加莱截面看起来像连续的线用等间隔采样画了原始轨迹改用局部极值或超平面穿越检测分岔图杂点很多暂态去除不充分丢弃比例提高到 30%分岔图靠近临界点处模糊临界慢化导致收敛缓慢对该区域单独加长时间跨度三维图颜色条看不出时间方向tspan 前段暂态主导颜色映射先截断暂态再设颜色变量代码运行太慢循环内动态拼接数组预分配数组或改用parfor6.5 几个容易被忽略的通用建议我建议把代码封装成函数而不是只写脚本。脚本中所有变量都暴露在工作区一旦来回调试很容易出现变量名覆盖问题。封装成函数后输入参数清晰输出明确后续换系统、换参数也更方便。函数化之后你还可以顺手写一个简单的配置文件把参数和绘图开关统一管理。另外保存图形时尽量用矢量格式。用exportgraphics(gcf, lorenz_3d.pdf, ContentType, vector)导出的矢量图在论文缩放时完全不会糊。位图在这种情况下很容易出现锯齿特别是在线上传或评审打印时观感差异很大。7. 经验总结与扩展方向7.1 我踩过的几个印象深刻的坑第一次画分岔图时我偷懒没有丢弃暂态结果在 rho24 到 25 附近看到一团密密麻麻的散点我还一度以为混沌在更小的参数下就已经出现了。后来查阅文献才发现理论预测的最小混沌参数应该在 28 附近那些点全是暂态响应。那一次让我彻底记住了暂态处理的重要性。还有一次我在整理庞加莱截面图时误用了等时间间隔采样画出来的“截面”是一条连续曲线怎么调整参数都没有碎片化特征。把findpeaks换上去之后结果立刻变成了分形点集。这件事让我意识到对图形概念的数学理解不扎实工具再好也白搭。我也曾经在画三维图时试图把整个吸引子轨迹用plot3加超细线宽展示图确实精致但文件大小接近 50MB导出论文插图时几乎无法处理。后来改用scatter3并适当抽稀文件体积和视觉效果都更理想。7.2 这套代码可以怎么扩展当前代码的核心框架非常通用你只需要替换微分方程函数就能画其他混沌系统。最简单的练手对象是 Rossler 系统dx/dt -y - z dy/dt x a*y dz/dt b z*(x - c)参数选取a0.2b0.2c5.7时系统处于混沌状态。把lorenz_system替换成 Rossler 的版本其他代码几乎不用改。你还可以尝试在三维相图上叠加两个系统的轨迹做对比视觉冲击力很强。更进一步你可以在庞加莱截面的基础上计算 Lyapunov 指数谱判断系统混沌程度也可以把分岔图与时间序列功率谱对比形成更完整的非线性分析链条。如果你做的是工程应用还可以把混沌系统作为激励源替代传统噪声信号研究系统在混沌输入下的响应特性。我目前正在把这段代码移植成更通用的“系统定义 分析绘图”框架只需要在配置文件中改方程和参数就能同时输出上述所有图形。计划加入对离散时间系统的支持这样就能把 Lorenz 映射等离散混沌系统也纳入分析流程。可以说围绕混沌系统数值可视化的这套方法到今天依然是分析和解释非线性行为的可靠起点值得耐心打磨。