企业官网建设流程全解析

刚接触非线性动力学的人，拿到的往往是一个看起来非常简单的三阶微分方程组。可真到了要讲清楚系统性质的时候，光盯着公式什么都看不出来。要把 Lorenz、Rossler、Chen 这类三维混沌系统的行为讲明白，庞加莱截面图、二维与三维相图、分岔图这三类图基本是绕不开的。这篇文章把我实际调试通过的一套 MATLAB 程序拆开讲，从 RK4 数值积分的写法，到截面怎么切、分岔图怎么扫，全都按可以直接复现的方式给出来。不管你是刚学混沌理论的本科生，还是需要用可视化支撑论文结果的研究生，照着敲一遍就能跑通，跑通之后再换成自己的系统也不难。

1. 为什么说相图、庞加莱截面和分岔图是混沌分析的“三件套”

1.1 相图回答“轨迹长什么样”

相图，或者说相空间轨迹图，是把系统的三个状态变量当成三维空间里的坐标，然后把时间演化的轨迹画出来。对 Lorenz 系统来说，经典参数下看到的就是那只著名的“蝴蝶”：轨道绕着一个中心转几圈，然后突然跳到另一个中心，来回往复。这张图回答的是最直观的问题——系统在状态空间里到底是怎么运动的。

相图的价值在于你一眼就能分辨出系统处于什么状态。稳定不动点在相图里就是一个点，周期解是一条闭合曲线，混沌则是一条永不闭合、也不自我重复的扭曲轨道。很多人第一次跑出蝴蝶形状时都会有点兴奋，因为这个形状本身就是混沌最直观的证据。但相图也有明显的短板：轨道在三维空间里缠绕得非常密，线条叠着线条，看不出多少精细结构。想进一步分析，就必须引入庞加莱截面。

1.2 庞加莱截面回答“轨迹的骨架结构”

庞加莱截面的思路很朴素：与其盯着整条连续轨迹看，不如在相空间里放一个“切面”，只记录轨道穿过这个切面时的点。连续的三维流就这样被降成了二维平面上的点集。

这么一降维，信息密度会高很多。周期解穿过截面时每次都落在同一个位置，在庞加莱图上就只会留下一个点；倍周期解会留下两个点；混沌系统的截面点则是一团具有自相似结构的密集点集，很多情况下呈现出分形特征。这个“离散化”过程才是真正能用来做定量分析的，后续很多 Lyapunov 指数估计、分岔分析、第一返回映射，都是以庞加莱截面为基础的。

1.3 分岔图回答“参数如何决定动力学”

相图和庞加莱截面都是在固定参数下观察系统行为，分岔图则把参数也作为变量，一次性回答“系统随参数变化会经历哪些状态转换”。做法是让某个参数连续扫描，对每个参数值都做一次数值积分，然后记录下某个状态变量在稳定轨迹上的极值（或者庞加莱截面交点），最后把这些点画在“参数-状态”平面上。

分岔图里最常见的图景就是倍周期分岔通向混沌：随着参数增大，原来的单周期变成双周期，再变成四周期、八周期，然后在某个临界参数附近突然进入混沌区。这张图能极好地展示混沌不是凭空出现的，而是动力学系统沿着参数轴演化的自然结果。

2. Lorenz系统数值求解的程序骨架：固定步长RK4的写法与参数封装

2.1 方程形式与经典参数

为了让讲解落在实处，这里统一用 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。sigma 和 rho 在气象情境里对应普朗特数和瑞利数，beta 是几何因子，但对数值分析来说，你只需要知道这组参数能让系统进入混沌状态就行。

选用 Lorenz 系统演示还有一层原因：这个系统虽然只有三个方程，但它的动力学行为非常丰富，从稳定不动点到周期解，再到混沌，只需要改变 rho 一个参数就能全部看到。这为后面画分岔图提供了绝佳素材。

2.2 RK4单步函数与主循环

MATLAB 里求解常微分方程首选往往是 ode45，但我在做这类混沌可视化时更喜欢自己写固定步长的四阶龙格-库塔（RK4）。原因后面展开说，先看代码。

下面这段是 RK4 主循环的完整写法：

clear; clc; close all; % 系统参数：经典Lorenz sigma = 10; rho = 28; beta = 8/3; % 右端函数，用匿名函数封装，方便后续替换成Rossler等其他系统 f = @(t, y) [sigma*(y(2) - y(1)); y(1)*(rho - y(3)) - y(2); y(1)*y(2) - beta*y(3)]; h = 0.01; % 固定步长 tspan = 0:h:200; % 总时长 y0 = [1; 1; 1]; Y = zeros(length(tspan), 3); Y(1, :) = y0(:)'; for n = 1:length(tspan) - 1 t_n = tspan(n); y_n = Y(n, :)'; k1 = f(t_n, y_n); k2 = f(t_n + h/2, y_n + h/2*k1); k3 = f(t_n + h/2, y_n + h/2*k2); k4 = f(t_n + h, y_n + h*k3); Y(n + 1, :) = (y_n + h/6*(k1 + 2*k2 + 2*k3 + k4))'; end % 丢弃前20%作为瞬态预热部分 warmup = round(length(tspan) * 0.2); Y = Y(warmup + 1:end, :);

这段代码没有调用任何工具箱，纯基础 MATLAB 就能跑。有几个位置要特别说明一下。

第一，右端函数f用匿名函数写，参数直接捕获取值。这样换系统时只需要改这一行，后面所有积分、画图、扫描代码都能复用，省下大量重复劳动。

第二，RK4 的四个斜率 k1 到 k4 的计算顺序严格按照标准格式，这里没有做自适应步长，但精度对 Lorenz 系统来说完全够。h = 0.01是我试过的比较稳的取值，再大一点比如 0.05，在分岔点附近会出现偏差；再小也不是不行，就是计算量翻几倍。

第三，也是最容易翻车的一点：积分完成后必须丢弃前面一段轨迹。初始点[1; 1; 1]并不在吸引子上，系统需要一段时间“预热”才能跑到吸引子上。如果不丢，相图上会多出一条从初始点直奔吸引子的拖尾线，庞加莱截面上也会混入大量瞬态点，整个图都会变脏。这里取的是丢弃前 20%，实际使用时可以根据系统复杂的程度调整到 30%~50%。

2.3 手动RK4和ode45的取舍

很多人会有疑问：MATLAB 自带 ode45 那么方便，为什么自己写 RK4？我来说说实际经验。

ode45 是自适应步长的 Runge-Kutta 方法，精度高，也能处理刚性不强的大多数系统。它的问题在于两点。第一，批量参数扫描时分岔图要对每个 rho 值都跑一次完整积分，ode45 的自适应步长机制会带来大量额外开销，而且步长序列在参数变化时也不稳定，偶尔会在某些临界点突然把步长压得非常小，导致一个参数点要跑很久。第二，庞加莱截面的穿越检测依赖等间隔采样，虽然 ode45 的输出也能插值处理，但自己控制固定步长显然更方便判断和调试。

当然这不是说 ode45 不能用。如果只是单次画一张相图，ode45 完全没问题，我早期也这么干。但如果目标是批量画分岔图、需要精细控制采样密度，手写 RK4 的固定步长方案反而更省心，整个循环可预测、可并行、可向量化优化。

3. 庞加莱截面：穿越判据、插值投影与“切出来”的离散动力学

3.1 截面选择的基本原则

庞加莱截面的第一步是选一个合适的“切面”。原则主要有两条：截面必须横截轨迹流，也就是说轨道不能贴着截面走，要真正“穿”过去；截面最好取在动力学比较有代表性的位置，通常会选择某个坐标等于常数的平面。

对 Lorenz 系统来说，最常用的截面是y = 0平面。为什么选它？因为 Lorenz 系统两只“蝴蝶翅膀”的对称轴就在 y 方向，轨道在两条翅膀之间切换时必然穿越y = 0这个平面，而且穿越方向清晰，适合做方向滤波。

这里还要引入一个细节：只有从负 y 穿越到正 y 的点才记录。因为 Lorenz 轨道穿越y = 0时是一来一回两个方向各穿一次，如果两个方向都记录，截面上会混入两组不同方向的点，分析价值会下降。加一个方向判据，只保留单侧穿越，图会更加干净。

3.2 穿越检测与线性插值

实际的数值积分结果是离散点列，轨道“穿越平面”在程序里就表现为相邻两个采样点分别位于平面的两侧。实现如下：

% 截面：y = 0，只记录由负到正穿越 cross_idx = find(Y(1:end-1, 2) < 0 & Y(2:end, 2) > 0); px = zeros(length(cross_idx), 1); pz = zeros(length(cross_idx), 1); for k = 1:length(cross_idx) i = cross_idx(k); % 线性插值，找到穿越点的精确位置 frac = (0 - Y(i, 2)) / (Y(i+1, 2) - Y(i, 2)); px(k) = Y(i, 1) + frac * (Y(i+1, 1) - Y(i, 1)); pz(k) = Y(i, 3) + frac * (Y(i+1, 3) - Y(i, 3)); end figure; plot(px, pz, 'k.', 'MarkerSize', 4); xlabel('$x$', 'Interpreter', 'latex'); ylabel('$z$', 'Interpreter', 'latex'); title('Poincare section: y = 0, \dot{y} > 0', 'Interpreter', 'latex');

穿越点为什么需要线性插值，而不是直接用采样点本身？因为步长h = 0.01时，相邻两个采样点在时间上差了 0.01，位置差通常在 0.1 量级。如果直接把跨过截面后的下一个点当作穿越点，误差最大可能达到这个量级。线性插值虽然简单，但可以把误差降到二阶小量，对绘制庞加莱截面图来说足够了。

如果对精度还有更高要求，可以用二分法在穿越处做局部细化，不过实际画图时一般用不上。对混沌吸引子来说，截面点的分形结构本质上对初始条件敏感，过分追求单个穿越点的精度没有太大意义。

3.3 第一返回映射的思路与极值法

庞加莱截面画出来之后，还可以继续做“第一返回映射”：把截面上每个点映射到轨道下一次穿回截面的那个点，得到一组p_neighbor到p_next的对应关系。对一维映射来说，这种图能非常直观地看出动力学是周期性的，还是混沌的，还是拟周期的。实现上就是在cross_idx基础上取相邻穿越点的序号做对应，代码非常短。

另一个经常用的替代方案是“局部极值法”。以 Lorenz 的 z 变量为例，记录每个局部极大值z_max，得到一组一维数据。这个方法在分岔图里特别有用，因为每个周期对应一个极大值，周期翻倍时极大值的个数也会翻倍，比截面上点的统计更直观。后面画分岔图时，主要就是靠这个极值法。

4. 二维与三维相图的绘制：轨迹结构、视角与时间维度的信息取舍

4.1 三维轨迹：不要一上来就画整条线

三维相图通常直接用 plot3 画轨迹线，但这里有个很实际的坑：积分 200 秒、步长 0.01，得到的是两万多个点，plot3 把这些点全部连起来，画出的线条会密到完全看不见结构，只剩一个黑乎乎的线团。

我常用的做法是降低可视化密度。计算时保持细步长，绘制时每隔 5 到 10 个点取一个样再连线，线条反而更清晰。代码如下：

figure('Position', [100 100 1200 400]); Y_plot = Y(1:5:end, :); % 抽样显示 subplot(1, 3, 1); plot3(Y_plot(:,1), Y_plot(:,2), Y_plot(:,3), '-', 'LineWidth', 0.3); xlabel('$x$', 'Interpreter', 'latex'); ylabel('$y$', 'Interpreter', 'latex'); zlabel('$z$', 'Interpreter', 'latex'); grid on; view(120, 25);

view(120, 25)这个视角是我反复试出来比较好看的角度，能看到蝴蝶结构的两个“翅膀”同时展开。如果你用默认视角，经常会觉得图形扁成一片，稍微转一下视角，立体感立刻就出来了。

如果想进一步看出轨迹演化的时间顺序，可以用颜色渐变代替纯色线条：

% 按时间着色，早期轨迹用冷色，后期用暖色 c = (1:size(Y_plot, 1))'; scatter3(Y_plot(:,1), Y_plot(:,2), Y_plot(:,3), 1, c, '.'); colormap(parula); colorbar;

这个处理特别适合展示“瞬态如何进入吸引子”的过程。我写论文插图时经常用这种着色方式，审稿人会比较喜欢，因为它比单纯的黑线多了一层信息。

4.2 二维投影怎么看才不白画

三维相图虽然立体，但有些结构细节会被遮挡，这时需要二维投影。对 Lorenz 系统，最常用的投影面是 XY 平面，因为它能直接呈现蝴蝶左右两叶的对称结构；XZ 平面更适合观察变量间的非线性耦合关系；YZ 平面用的相对少，但在检查折叠结构时很管用。

subplot(1, 3, 2); plot(Y_plot(:,1), Y_plot(:,2), '.', 'MarkerSize', 1); xlabel('$x$', 'Interpreter', 'latex'); ylabel('$y$', 'Interpreter', 'latex'); grid on; subplot(1, 3, 3); plot(Y_plot(:,1), Y_plot(:,3), '.', 'MarkerSize', 1); xlabel('$x$', 'Interpreter', 'latex'); ylabel('$z$', 'Interpreter', 'latex'); grid on;

注意这里用了散点而不是连线。混沌轨迹在二维投影里如果连线，线条交叉严重，散点反而能看出空间的覆盖范围。MarkerSize = 1可能会让你觉得点太小，但随着采样点数增加，小点画出来的密集区域更均匀，不会糊成一团。

判断一张二维相图画得是否合格，关键是看它能不能让你一眼看出吸引子的拓扑结构。Lorenz 的 XY 投影应该是两个近似对称的圆叶，中间有细的连接部分；如果画出来是乱七八糟的斑块，先检查是否忘了丢弃瞬态部分，再看步长是不是太大导致积分“跳”出了吸引子。

4.3 画图时的几个通用优化技巧

实际作图时还有几个经验。第一，坐标轴范围不要自动缩放到底，稍微留一点空白，图会耐看很多，可以用xlim、ylim手动微调。第二，MATLAB 默认字体对中文支持不好，内嵌的 latex 解释器则很美观，xlabel('$x$', 'Interpreter', 'latex')这种写法推荐直接沿用。第三，保存图片时优先用矢量格式，比如exportgraphics(gcf, 'phase.eps', 'ContentType', 'vector')，论文里放大也不会糊。

5. 分岔图的参数扫描流程：瞬态丢弃、峰值记录与批量计算优化

5.1 分岔图的记录方式选择

分岔图的标准做法是固定其他参数，只让一个参数连续变化。对 Lorenz 系统来说，扫描 rho 是最经典的方案，因为它的动力学随 rho 变化经历了“稳定不动点 -> 倍周期分岔 -> 混沌”的完整过程。

记录什么数据？两种主流方式。第一种用庞加莱截面交点，画出来的是截面上的某个坐标随参数的变化；第二种用局部极大值，比如记录 z 变量的每个局部峰值。我推荐第二种，因为它实现简单、抗噪性好，而且对周期倍化过程分辨得很清楚。系统处于周期 1 状态时 z 只有一个峰值，周期 2 状态出现两个不同的峰值，周期 4 就出现四个，倍周期过程在图上会显示成清晰的分叉。

峰值检测不复杂，因为 RK4 积分出来的曲线非常光滑，直接用局部比较即可：

% 记录z的局部极大值 z = Y(ceil(nt * 0.6):end, 3); % 丢弃前60%瞬态，只分析后半段 pks = z(2:end-1); mask = z(1:end-2) < pks & pks > z(3:end); z_peaks = pks(mask);

这里的ceil(nt * 0.6)是关键。扫描参数时，每个 rho 值下系统都需要重新“稳定”到新吸引子上，如果直接从头记录，前面的瞬态峰值会混入图里，产生大量噪声点。丢弃的比例可以视收敛速度调整，一般 50% 到 70%。

5.2 参数续传和固定初值的差异

扫描参数时有一个容易忽略的选择：每个 rho 值是从固定初值开始积分，还是沿用上一个 rho 值积分结束时的状态作为初值？这两种方式得到的分岔图可能非常不同。

固定初值的好处是实现简单，但坏处也很明显：某些参数下系统可能有多个共存吸引子，从同一个初值出发可能只跳到其中一个分支上，导致分岔图缺分支；而且每次从瞬态重新开始，需要丢弃更多的预热段。

我实际更推荐“参数续传”方式：先在小参数下从任意初值跑一段，让系统落到吸引子上，然后把这个终点状态作为下一个参数点的初值。这样吸引子会沿着参数轴“连续演化”，看到的动力学分支是沿着同一条路径展开的，物理意义更明确。

rho_list = 0.5:0.25:50; Nrho = length(rho_list); sigma = 10; beta = 8/3; h = 0.01; tmax = 800; nt = round(tmax / h); % 初始状态 y_current = [1; 1; 1]; results = cell(Nrho, 1); for r = 1:Nrho rho = rho_list(r); f = @(t, y) [sigma*(y(2) - y(1)); y(1)*(rho - y(3)) - y(2); y(1)*y(2) - beta*y(3)]; Y = zeros(nt, 3); Y(1, :) = y_current'; for n = 1:nt - 1 k1 = f(0, Y(n, :)'); k2 = f(0, Y(n, :)' + h/2*k1); k3 = f(0, Y(n, :)' + h/2*k2); k4 = f(0, Y(n, :)' + h*k3); Y(n+1, :) = (Y(n, :)' + h/6*(k1 + 2*k2 + 2*k3 + k4))'; end % 续传：把当前rho下的终点作为下一点的初值 y_current = Y(end, :)'; % 丢弃前60%，记录z局部极大值 z = Y(ceil(nt * 0.6):end, 3); if length(z) > 3 pks = z(2:end-1); mask = z(1:end-2) < pks & pks > z(3:end); z_peaks = pks(mask); results{r} = [rho * ones(length(z_peaks), 1), z_peaks(:)]; end end all_data = cell2mat(results); figure; plot(all_data(:,1), all_data(:,2), 'k.', 'MarkerSize', 1); xlabel('$\rho$', 'Interpreter', 'latex'); ylabel('$z_{max}$', 'Interpreter', 'latex'); grid on;

用results这个 cell 数组来收集数据，是因为不同 rho 值下峰值个数不一样，无法直接用矩阵存储。最后用cell2mat合并成一个两列矩阵，第一列是 rho，第二列是峰值，一次性画出来。

5.3 计算量估算与并行加速

不要小看这段代码的计算量。rho 从 0.5 到 50、步长 0.25，共 198 个参数点；每个点要积分 800 秒、步长 0.01，也就是 8 万次 RK4 迭代；总共约 1600 万次迭代。单线程跑大概需要几分钟到十几分钟，具体取决于机器性能。第一次调试时建议把rho_list改成0.5:2:50，先验证代码能跑通，再加密参数步长。

如果电脑有多核，可以加一行parfor替换for。不过需要注意，results数组需要换一种合并方式，因为parfor对循环变量索引的要求比较严格：

rho_collect = []; peak_collect = []; parfor r = 1:Nrho % 这段代码需要把y_current改为独立初始化，不能依赖上一轮的y_current ... end

也就是说，parfor和参数续传是冲突的，因为续传方式本身就存在循环内的数据依赖。如果用并行，就得放弃续传，改成固定初值。两种方式取舍，我的建议是：优先保证分岔图的连续性和可解释性，单线程跑几分钟完全能接受；实在嫌慢，再放弃续传换 parfor。

分岔图画完之后，可以验证一下经典结论：rho 在约 1 到 13 附近系统趋向一个稳定不动点，图上只有一条接近常数的水平线；rho 继续增大到约 24 附近，出现周期倍化分岔；到 28 附近进入混沌区；再往上可能会看到周期窗口，这些窗口通常出现在混沌区内部，表现为一条突然“变干净”的竖直线。能把这几段都看到，说明扫描程序没有问题。相比代码本身，对图的解读才是真正体现功力的地方。

从Lorenz换到Rossler：改一处就能复用的核心逻辑

整套程序最有价值的部分是它的可迁移性。换系统时，只需要改右端函数f和对应参数。比如 Rossler 系统：

dx/dt = -y - z dy/dt = x + ay dz/dt = b + z(x - c)

取经典的 a = 0.2、b = 0.2、c = 5.7，右端函数改成：

a = 0.2; b = 0.2; c = 5.7; f = @(t, y) [-y(2) - y(3); y(1) + a * y(2); b + y(3) * (y(1) - c)];

庞加莱截面的位置就需要重新想一下了。Rossler 的吸引子像一个纸环，截面取y = 0时轨道穿越方向很明确，也能用。Chen 系统、Chua 电路系统也都按同样的思路处理：先确定系统的状态变量范围，再根据轨道走向选一个横截截面。

我自己在换系统时有一个习惯：先把相图跑出来，盯着轨迹看它在哪个平面穿越最干净，再定截面位置。直接套用别的系统的截面参数往往会得到一团没有规律的点，这并不是程序问题，而是截面没有真正横截轨迹流。

最后再分享一个细节习惯：画庞加莱截面和分岔图时，散点图的MarkerSize不要一上来就调大。先用 1 到 2 画一遍，如果点的密度结构看不清，再调大。我见过很多人直接设成 6 到 8，结果吸引子的精细结构全糊在一块，落点密集区变成了一个实心黑块。用小尺寸散点才能看到混沌吸引子截面那种“内疏外密”的分布层次。这套程序跑完之后，建议你顺手在纸上把每个图对应的物理意义写一遍：相图看拓扑，截面看结构，分岔图看演化路径。能讲清这三者的区别，你对混沌系统的理解就已经超过大部分只会贴代码的人了。