企业官网建设流程全解析

第一次在论文里看到 Lorenz 吸引子那双蝴蝶翅膀时，大多数人都会被它的几何美感抓住。但真到了自己动手用 Matlab 画图的时候会发现，画出一条好看的相图只是起点。做非线性动力学研究，或者自己搭了一个混沌系统想验证它的行为，通常需要把三维相图、二维相图、庞加莱截面图和分岔图放在一起看：相图告诉你吸引子长什么样，庞加莱截面告诉你轨迹在某个截面上的“离散节奏”，分岔图则直接扫描参数告诉你在哪个区间进入混沌、哪个区间又回到周期。这篇文章就是我整理好的一套 Matlab 程序合集，以 Lorenz 这个经典的三阶微分方程系统为例，把四种图形全部跑通，代码可以直接复制改成你自己的系统。适合刚接触混沌、需要给论文补图，以及想验证自己写的微分方程是否混沌的同学。

1. 为什么画了相图还要画分岔图：四种图形各回答什么问题

很多人第一次接触混沌系统时，以为只要能画出那条蝴蝶曲线就算搞定了。实际上相图只是“看起来像混沌”，它不能精确告诉你系统到底是周期、拟周期还是混沌。我个人的习惯是：先画三维相图和二维投影图建立直觉，再用庞加莱截面判断运动形态，最后用分岔图扫描参数摸清全局。这四种图各有分工，谁也不能替代谁。

1.1 四种图形工具的分工

三维相图和二维相图回答的是“吸引子长什么样”：轨迹在相空间里怎么走，是被压缩到某个平面附近，还是真的在三维空间里铺开。庞加莱截面图回答的是“轨迹穿过某个平面时留下的点有什么规律”：如果只有有限个点，对应周期运动；如果形成一条闭合曲线，对应拟周期；如果是一团带有分形结构的点云，才是混沌。分岔图回答的是“某个参数变化时，系统行为怎么切换”：稳定点、极限环、周期倍化、混沌窗口，全部能在一张图上看到。

这四种图是互补的。我见过不少同学只画了三维相图，就断定自己的系统是混沌，结果拿去算 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(i+1, 3) - Y(i, 3)); % 归一化权重 xCross(k) = Y(i,1) + w * (Y(i+1,1) - Y(i,1)); yCross(k) = Y(i,2) + w * (Y(i+1,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算一遍，如果吸引子形态差异明显，立刻怀疑是求解器不合适。

我自己在跑完 Lorenz 的整套图之后，形成了一套固定的验证流程：先画三维相图确认吸引子几何形态，再利用庞加莱截面判断运动类型，最后用分岔图扫描参数全局。这三个工具组合起来大约能覆盖 80% 的定性判断需求。如果你想进一步定量证明系统是混沌，那就要上最大 Lyapunov 指数了，但图形工具依然是定位参数区间的第一站。上面的代码全部可以直接复制到你的脚本里，把lorenz_sys替换成你自己的三阶方程，图形部分基本不用改。

最后再说一个小技巧：所有图完成后，记得用exportgraphics(gcf, 'filename.png', 'Resolution', 300)导出高分辨率图片，论文排版时比截图清晰得多。混沌系统图形最忌讳的就是细节糊成一片，高分辨率输出能保住那些细如发丝的拉伸折叠结构。