MATLAB混沌系统可视化工具包：相轨迹、庞卡莱截面与多模式分岔图一键生成

本文还有配套的精品资源，点击获取

简介：一套开箱即用的MATLAB混沌分析工具，无需额外工具箱，直接运行即可绘制三维相图（y_phase.m）、庞卡莱截面（y_poincare.m，输出Poincare.eps）、单参数正向/反向分岔图（y1_bifurcation_C1_Forward_and_Backward.m）、双参数联合分岔图、以及初值敏感性引发的共存吸引子分岔（y1_bifurcation_C2_CoexistingBifurcation.m）。支持对不同初始条件（如y1_bifurcation_InitialValue_x.m）进行扫描，识别单涡卷与双涡卷吸引子的共存与演化；通过getmax_y1.m~getmax_y4.m系列函数自动提取状态变量极值序列，构建高精度分岔图；LTF_Bifur1_y1y2y3_text等脚本辅助多变量同步可视化。system.m统一定义动力学方程，模块化设计便于替换洛伦兹、Chen、Lü等常见混沌模型。所有脚本含中文注释，结构清晰，适配R2018a及以上标准MATLAB环境。

1. 这不是“画图工具”，而是一套混沌系统诊断工作台

你有没有试过，在MATLAB里敲完一个洛伦兹方程，跑出一堆乱七八糟的曲线，却说不清它到底处在周期态、准周期态还是混沌态？或者明明改了参数，相图看起来差不多，但分岔图上却突然冒出几条“幽灵分支”，反复检查代码也没错——最后发现是初值选在了两个共存吸引子的边界上？我干过太多次这种事。这套工具包，就是我在带本科生做混沌实验、指导研究生跑数值仿真、甚至自己调试新提出的四维超混沌系统时，从零开始攒出来的“混沌诊断工作台”。它不叫“可视化插件”，也不叫“绘图脚本集”，它更像一把带刻度、有校准、还能自动识别异常的示波器：你把系统扔进去，它告诉你此刻的运动本质是什么，参数往哪调会失稳，初值差0.001会导致轨迹跳到另一个完全不同的世界里去。

核心关键词——混沌相图、庞卡莱截面、分岔图绘制、初值敏感性、共存吸引子——这五个词不是并列关系，而是层层递进的诊断链条。相图是“望闻问切”的“望”，让你一眼看出轨迹是否缠绕、是否发散；庞卡莱截面是“切”，用一个超平面去“切”连续流，把三维混沌压缩成二维离散点集，看它是单点（不动点）、两点（周期2）、一团密密麻麻的云（混沌）；分岔图是“问”，系统性地追问：“当这个参数从2.8变到3.2，系统行为究竟经历了哪些质变？”；初值敏感性是“验”，验证李雅普诺夫指数是否为正——不是靠算指数，而是直接让两组仅差1e-12的初值跑一遍，看它们的轨迹何时彻底分离；而共存吸引子，则是整个链条的“压轴难题”，它意味着同一个参数下，系统可能有多个稳定归宿，你的初值落在哪片“盆地”，就决定了你最终看到的是单涡卷还是双涡卷。这套工具包的全部设计，都围绕着如何把这五个抽象概念，变成你鼠标一点、命令一敲就能看见、能比对、能存档、能写进论文图注里的具体图像。它不需要Symbolic Math Toolbox去解符号方程，不需要Parallel Computing Toolbox去加速——所有计算都在标准ODE求解器ode45框架内完成，所有图形输出都控制在figure句柄可导出EPS/PNG的范围内。我把它部署在实验室三台不同年代的电脑上（R2018a、R2020b、R2023a），没装任何额外工具箱，运行零报错。这不是炫技，而是为了确保你今晚十点赶论文图时，不会因为某个工具箱版本不兼容而抓狂。

2. 工具包整体架构与设计逻辑拆解

2.1 模块化分层：为什么必须把system.m单独拎出来？

很多人第一次接触混沌仿真，习惯把微分方程直接写死在主脚本里。比如洛伦兹系统，直接在y_phase.m开头写：

dx = sigma*(y-x); dy = x*(rho-z) - y; dz = x*y - beta*z;

这看似简单，但一旦你要对比Chen系统和Lü系统，就得复制粘贴三份几乎一样的绘图脚本，只改中间几行方程——这是灾难性的维护成本。这套工具包强制采用“动力学模型与可视化逻辑解耦”的设计，核心就在system.m这个文件。它不是一个简单的函数，而是一个参数驱动的状态方程工厂。打开它，你会看到这样的结构：

function dxdt = system(t, x, params) % SYSTEM 混沌系统动力学方程统一接口 % 输入: t - 时间; x - 状态向量 [x1,x2,x3,...]; params - 结构体参数 % 输出: dxdt - 状态导数向量 % 注意：此函数必须返回列向量！ switch params.model_name case 'lorenz' dxdt = [params.sigma*(x(2)-x(1)); x(1)*(params.rho-x(3)) - x(2); x(1)*x(2) - params.beta*x(3)]; case 'chen' dxdt = [params.a*(x(2)-x(1)); params.c*x(1) - x(1)*x(3) + x(2); x(1)*x(2) - params.b*x(3)]; case 'lu' dxdt = [params.a*(x(2)-x(1)); params.c*x(1) - x(1)*x(3); x(1)*x(2) - params.b*x(3)]; end dxdt = dxdt(:); % 强制转为列向量，避免ode45报错

这个设计背后有三个硬性理由：
第一，可追溯性。你在y_phase.m里调用ode45(@system, tspan, x0, opts, params)，所有参数（sigma、rho、beta）都来自外部传入的params结构体。这意味着你生成的每一张相图，其底层动力学模型、参数取值、甚至模型名称（’lorenz’/’chen’），都完整记录在调用它的脚本中，而不是藏在某个函数内部。写论文时，图注可以直接写“基于system.m中定义的Chen系统（a=35, b=3, c=28）”，无需再翻源码确认。
第二，可替换性。要新增一个Liu系统？只需在system.m的switch里加一个case 'liu'分支，写好它的方程，然后在调用脚本里把params.model_name = 'liu'，其余所有绘图、分岔、庞卡莱脚本全部无缝兼容。我去年帮一位做忆阻混沌电路的同学接入他们自研的四维模型，只改了system.m的12行代码，其他37个脚本一行未动。
第三，可复现性。params结构体可以轻松保存为.mat文件。你跑完一组双参数分岔，把完整的params连同结果一起存下来，半年后导师问“当时那个红点对应的参数是多少？”，你load进来就能立刻复现，而不是对着几十个变量名猜来猜去。

2.2 极值提取引擎：getmax_y*.m系列为何要分四个文件？

分岔图的本质，是观察系统长期行为对参数的依赖关系。而“长期行为”在数值仿真中，最稳健的代理指标就是状态变量的局部极值序列。比如对洛伦兹系统的z变量，我们不画整个z(t)曲线，而是提取它所有极大值点构成的序列{z_max1, z_max2, …}，当系统进入混沌态时，这个序列会呈现无规则分布；当处于周期态时，它会收敛到有限个点上。getmax_y1.m到getmax_y4.m正是为此而生，但它们不是简单的“找极大值”函数。

以getmax_y1.m为例，它接收ODE求解器输出的完整时间序列[t, y]（其中y是N×3矩阵，对应x,y,z），然后执行一套五步清洗流程：
1.预热剔除：自动识别并丢弃前20%的轨迹点。这是关键！混沌系统初始阶段常有暂态过程（transient），比如从任意初值出发，需要绕几十圈才能真正落到吸引子上。若不剔除，分岔图底部会出现大量杂乱“毛刺”，误导你认为系统在该参数下不稳定。
2.滑动窗口极值检测：不用findpeaks那种全局搜索，而是设定一个宽度为window_len（默认50个点）的滑动窗口，在每个窗口内找最大值。这能有效抑制数值噪声导致的伪峰，尤其对高采样率数据至关重要。
3.幅值阈值过滤：剔除所有小于min_amp（默认0.1）的极值。混沌吸引子有明确的几何边界（如洛伦兹z∈[0,45]），低于此阈值的“峰”大概率是数值误差。
4.收敛性判据：检查最后100个极值的标准差是否小于tolerance（默认0.05）。若不满足，说明系统尚未进入稳态，函数会报错提示“需延长仿真时间”，而不是强行出图。
5.去重与排序：对剩余极值进行unique(round(...,3))处理，消除因浮点精度导致的重复点，并升序排列，为后续绘图做准备。

那为什么需要getmax_y1到getmax_y4四个文件？因为不同混沌系统维度不同，关注变量不同。y1通常指第一个状态变量（如洛伦兹的x），y2是第二个（y），以此类推。而getmax_y4.m的存在，就是为了支持四维系统（如Qi系统、超混沌Lorenz-Haken系统）。你不必记住哪个文件对应哪个变量——所有分岔脚本（如y1_bifurcation_C1_Forward_and_Backward.m）在调用时，会根据params.plot_var = 1或params.plot_var = 3自动选择对应的getmax_y*.m。这种设计，把“用户意图”（我要看z变量的分岔）和“底层实现”（用getmax_y3提取z的极值）彻底隔离开，降低误操作风险。

2.3 分岔图策略库：正向/反向扫描与共存吸引子分离的物理意义

分岔图不是简单地“参数扫一遍，画个点”。y1_bifurcation_C1_Forward_and_Backward.m和y1_bifurcation_C2_CoexistingBifurcation.m代表了两种根本不同的物理问题，它们的算法逻辑也截然不同。

C1：单参数正向/反向扫描解决的是路径依赖性（hysteresis）问题。想象一个非线性弹簧，当外力缓慢增大时，它在某个临界点突然失稳（snap-through）；而当外力缓慢减小时，它又在另一个更低的临界点才恢复原状——这就是滞后环。混沌系统同样存在这种现象。C1脚本的流程是：
- 正向扫描：参数p从p_min到p_max，每次迭代使用上一步的终值作为下一步的初值；
- 反向扫描：参数p从p_max回到p_min，同样使用上一步终值；
- 最终将两组结果画在同一张图上，用不同颜色区分。

如果某段参数区间内正向与反向结果不重合，就出现了滞后环，这往往是系统存在多稳态（multi-stability）的铁证。我用它诊断过一个忆阻Chua电路，发现当控制参数在1.2~1.5区间时，正向扫描显示单周期，反向扫描却显示双周期——这意味着电路在此区间具有双稳态，可通过脉冲触发在两种模式间切换。

C2：共存吸引子分岔则直击混沌理论中最棘手的问题之一：初值决定论下的不可预测性。同一组参数下，不同初值可能导致轨迹收敛到完全不同的吸引子。y1_bifurcation_C2_CoexistingBifurcation.m的破解思路是“暴力枚举+聚类识别”。它会：
- 固定参数p；
- 在预设的初值空间（如x0∈[-2,2], y0∈[-2,2], z0∈[-2,2]）内，按网格生成N组初值；
- 对每组初值独立运行ODE，提取其稳态极值序列；
- 使用DBSCAN聚类算法，根据极值序列的欧氏距离，将N组结果自动分为K类；
- 每一类代表一个共存吸引子，其分岔分支被单独绘制。

我曾用它分析一个改进型Lü系统，在参数ρ=32时，发现存在三个共存吸引子：一个单涡卷、一个双涡卷、一个四涡卷。它们的z变量极值范围分别是[15,25]、[5,35]、[0,45]，在分岔图上表现为三条完全分离的、互不干扰的分支。没有C2脚本，你只会看到一片混乱的“雾”，而C2把它清晰地解构成三个独立世界。这才是真正揭示混沌系统内在复杂性的利器。

3. 核心功能实操详解与关键参数配置

3.1 相轨迹绘制（y_phase.m）：不只是画线，更是理解吸引子几何

y_phase.m是整个工具包的“门面”，但它绝非简单的plot3(x,y,z)。它的核心价值在于可控的轨迹演化呈现和吸引子结构解析。运行它之前，你需要准备一个params结构体，至少包含：

params.model_name = 'lorenz'; % 必须与system.m中定义一致 params.sigma = 10; params.rho = 28; params.beta = 8/3; % 洛伦兹经典参数 params.x0 = [1,1,1]; % 初值，列向量！ params.tspan = [0, 100]; % 仿真时间，足够长以覆盖暂态 params.save_eps = true; % 是否保存为EPS矢量图（论文首选） params.show_trajectory = true; % 是否显示完整轨迹线（true）或仅显示终态点（false） params.point_size = 12; % 终态点大小，用于突出吸引子骨架

最关键的参数是params.show_trajectory。当设为true时，脚本会用comet3动画式绘制，你能亲眼看到轨迹如何从初值出发，绕行、拉伸、折叠，最终被“捕获”进混沌吸引子的螺旋结构中——这对教学演示极其有力。但若设为false，它会先运行完整ODE，然后只提取最后20%时间点的状态，用scatter3绘制为彩色点云，并自动添加坐标轴标签、标题、以及一个吸引子尺寸标注框（bounding box）。这个框的长宽高会实时计算并显示在图标题里，例如“Lorenz Attractor (x: -20→20, y: -25→25, z: 0→45)”，让你量化比较不同参数下吸引子的膨胀程度。

提示：y_phase.m内置了“轨迹稀疏化”功能。当tspan很长（如1000秒）导致点数过多时，它会自动对时间序列进行降采样，保证绘图流畅且不失真。你无需手动干预，它会根据length(t)自动选择采样间隔。

3.2 庞卡莱截面生成（y_poincare.m）：如何选对“切面”才能看见真相？

庞卡莱截面的质量，90%取决于截面的选择。y_poincare.m不提供“智能推荐”，而是给你三把精准的“手术刀”：

1. 平面截面（Plane Section）：最常用，如z=27平面。脚本要求你指定params.poincare_plane = 'z'和params.poincare_value = 27。它会遍历ODE输出的所有相邻点对(t_i,y_i)和(t_{i+1},y_{i+1})，当y_i(3)*y_{i+1}(3) < 0（即z坐标跨过27）时，用线性插值计算精确交点。注意：它默认只记录正向穿越（z从小于27变为大于27），这是为了避免同一周期内多次穿越造成的点堆积。若需双向穿越，修改params.direction = 'both'。

2. 事件驱动截面（Event-Based Section）：适用于复杂系统。例如，在Chen系统中，你想观察当x*y > 50时的截面。此时设置params.poincare_event = @(t,x) x(1)*x(2) - 50，脚本会调用odeset的Events选项，让ODE求解器在事件发生时精确停驻并记录状态。

3. 等时截面（Isochronous Section）：对研究同步现象至关重要。设置params.poincare_iso_period = 0.5，脚本会在t=0.5, 1.0, 1.5…这些固定时刻记录状态，无论系统当前处于什么相位。这能揭示系统是否具有近似周期性。

无论哪种方式，y_poincare.m都会输出一个Poincare.eps文件，并在命令行打印关键统计：

Poincare Section on z=27: 1247 points generated. Mean distance between points: 0.83 | Std: 0.41 Estimated fractal dimension (box-counting): ~2.15

这个分形维数估算，是判断混沌与否的定量依据——若维数接近整数（如1.02），很可能是周期轨道；若在2.0~2.5之间，则高度疑似混沌。

3.3 分岔图构建全流程：从参数扫描到图像精修

以最常用的单参数正向/反向分岔为例，y1_bifurcation_C1_Forward_and_Backward.m的执行流程如下：

Step 1：参数网格定义

params.p_range = linspace(2.5, 3.5, 1000); % 1000个参数点 params.p_step = 0.01; % 步长，影响分辨率

Step 2：初值继承策略
这是C1脚本的灵魂。它内部维护一个x_last变量，存储上一次ODE仿真的终值。当p从2.500增至2.501时，x0被设为x_last，而非重新随机初始化。这模拟了物理系统中参数缓慢变化的真实场景。

Step 3：极值提取与存储
对每个p_i，调用ode45得到[t,y]，再调用getmax_y1(t,y,params)提取x变量的稳态极大值序列。脚本会自动判断：若序列长度<50，视为未收敛，跳过此点；若>50，则取最后50个点，存入结果矩阵bif_data。

Step 4：图像生成与优化
最终绘图不是简单的scatter(p_vec, bif_data)。它采用抖动+半透明叠加策略：
- 对每个p_i，将其对应的50个极值点在p轴上做微小随机抖动（±0.0005），避免垂直堆叠成一条线；
- 设置MarkerFaceAlpha = 0.3，使高密度区域自然加深，低密度区域透出背景；
- 添加xlabel('Parameter \rho'),ylabel('x_{max}'),title('Lorenz Bifurcation (Forward \& Backward)')，所有LaTeX语法均被MATLAB正确渲染。

注意：bifurcation_diagram.png是该脚本的默认输出，但你可以在脚本末尾轻松添加print('-depsc2', 'MyBifurcation.eps')生成出版级矢量图。我建议永远优先保存EPS，PNG仅作快速预览。

3.4 共存吸引子深度分析（y1_bifurcation_C2_CoexistingBifurcation.m）：如何让“看不见”的分支显形？

C2脚本的威力，在于它把“初值敏感性”从一句教科书结论，变成了可量化、可绘图、可比较的实验数据。其核心是initial_grid参数：

params.initial_grid = struct(... 'x', linspace(-3,3,5), ... % x0取5个值：-3,-1.5,0,1.5,3 'y', linspace(-3,3,5), ... % y0同理 'z', linspace(-3,3,5)); % z0同理 % 总计 5x5x5 = 125组初值

脚本会遍历这125组初值，对每组独立运行ODE（使用相同的params.tspan和params.model_name），提取其稳态x变量的极值序列。然后，它启动聚类分析：

1. 将125个极值序列（每个长度约100）展平为125×100矩阵；
2. 计算每两行之间的欧氏距离，构建距离矩阵；
3. 调用clusterdata函数，设定'maxclust', 5（最多允许5类），自动确定最优聚类数；
4. 对每一类，计算其所有成员的极值序列的均值与标准差，作为该吸引子的“指纹”；
5. 最终绘图：横轴仍是参数p，纵轴是x_max，但每个点的颜色代表其所属的吸引子类别（1~K），并用不同形状标记（圆圈、方块、三角等）。

我用它分析一个双涡卷Lü系统时，在ρ=30处发现了两类吸引子：一类极值集中在[10,20]（单涡卷），另一类在[5,35]（双涡卷）。当ρ增加到31时，单涡卷类消失，所有初值都落入双涡卷类——这清晰地刻画了单涡卷向双涡卷的演化分岔点。没有C2，你只会看到ρ=30处一片模糊的“雾”，而C2把它解构成两个泾渭分明的世界。

4. 常见问题排查与独家避坑指南

4.1 “分岔图全是噪点，看不出结构！”——暂态未剔除的典型症状

现象：你运行y1_bifurcation_C1_Forward_and_Backward.m，得到的图上布满杂乱的竖线，像静电干扰，完全无法识别周期窗口或混沌带。

根源：getmax_y*.m中的预热剔除比例（warmup_ratio）设置过小。默认值0.2（20%）对大多数系统足够，但对某些慢收敛系统（如某些分数阶混沌）可能不足。

排查步骤：
1. 打开getmax_y1.m，找到warmup_ratio = 0.2;这一行；
2. 临时改为warmup_ratio = 0.5;，重新运行；
3. 若噪点显著减少，说明问题定位准确；
4. 进一步验证：用y_phase.m单独运行一个参数点（如ρ=28），开启params.show_trajectory = true，观察轨迹需要多久才稳定进入吸引子（目测约30秒），则warmup_ratio应设为30/100 = 0.3。

终极方案：在y1_bifurcation_*.m脚本中，添加一个“暂态诊断模式”。在循环体内加入：

if p == params.p_range(1) % 只对第一个参数点诊断 figure; plot(t, y(:,1)); title(['Transient Analysis at p=',num2str(p)]); line([t(end)*0.2, t(end)*0.2], ylim, 'Color','r','LineStyle','--'); legend('x(t)', 'Warmup Cutoff'); end

这张图会直观告诉你，20%的截断点是否合理。

4.2 “庞卡莱截面空空如也，一个点都没有！”——截面未被穿越的真相

现象：y_poincare.m运行完毕，命令行显示0 points generated，Poincare.eps是空白。

根源：最常见的原因是截面位置选在了吸引子之外。例如，对洛伦兹系统，z变量的混沌吸引子范围是[0,45]，若你设params.poincare_value = 100，显然永远穿不过。

快速诊断法：
1. 先用y_phase.m绘制该参数下的相图，观察z坐标的实际范围；
2. 若z始终<50，就把poincare_value设为30（吸引子中部）；
3. 更保险的做法：用y_poincare.m的'debug'模式。在调用前加一行：params.debug_mode = true;，它会输出y矩阵中z列的最大最小值，如z_min = 0.2, z_max = 44.8，你立刻就知道可选范围。

进阶陷阱：对于高维系统，某些截面可能因轨迹过于“稀疏”而错过。此时启用params.poincare_event，定义一个更宽松的事件，如@(t,x) abs(x(3)-27) < 0.1（z在26.9~27.1之间即触发），比精确等于更鲁棒。

4.3 “共存吸引子分岔图颜色混乱，分不清哪类是哪类！”——聚类失效的应对

现象：C2脚本输出的图上，颜色斑驳，没有清晰的分支，聚类数K波动很大（有时3类，有时7类）。

根源：DBSCAN聚类对距离阈值epsilon极度敏感。默认值是根据经验设定的，但不同系统、不同参数下，极值序列的尺度差异巨大。

解决方案：手动指定epsilon。在y1_bifurcation_C2_CoexistingBifurcation.m中，找到聚类调用行：

idx = clusterdata(D, 'maxclust', max_clust, 'distance', 'euclidean');

改为：

epsilon_manual = 2.5; % 根据你的系统调整！ [idx, C] = dbscan(D, epsilon_manual, min_pts); % min_pts默认5

如何确定epsilon_manual？运行一次脚本，获取距离矩阵D，然后画它的直方图：

figure; histogram(D(:), 50); xlabel('Pairwise Distance'); ylabel('Count'); line([epsilon_manual, epsilon_manual], ylim, 'Color','r');

理想情况是，红色线左侧有大量计数（同类吸引子内距离），右侧计数趋近于0（异类吸引子间距离）。我通常把epsilon_manual设为直方图峰值右侧的“谷底”位置。

4.4 “所有脚本都报错‘Undefined function or variable ‘params’’！”——新手必踩的初始化坑

现象：刚下载工具包，双击y_phase.m，MATLAB立刻报错，提示params未定义。

真相：这不是bug，而是设计哲学——所有脚本都不自带参数，强迫你显式定义。这杜绝了“上次运行的参数污染本次结果”的隐患。

正确姿势：
1. 新建一个run_example.m脚本；
2. 在里面定义params结构体（参考前面各节的示例）；
3. 然后调用目标函数：y_phase(params);；
4. 运行run_example.m，而非直接运行y_phase.m。

我建议你建立一个params_template.m文件，里面写好所有常用参数的默认值，每次新项目时copyfile一份，再修改。这样既规范，又避免遗漏。

4.5 “为什么我的分岔图和论文里的不一样？”——采样率与精度的隐性战争

现象：你复现一篇论文的分岔图，参数完全一样，但你的图上周期窗口更窄，混沌带更“毛”。

根源：两个隐形参数在作祟——ODE求解器精度和极值提取分辨率。
-ode45的相对误差容限RelTol默认是1e-3，论文作者可能用了1e-6；
-getmax_y*.m的滑动窗口window_len默认50，若时间序列采样太密（如dt=0.001），50点只覆盖0.05秒，可能漏掉真实极值。

实操校准：
在params中显式设置求解器选项：

opts = odeset('RelTol', 1e-6, 'AbsTol', 1e-9, 'MaxStep', 0.1); params.odeset_opts = opts;

并在getmax_y1.m中，根据你的dt动态调整window_len：

dt = t(2)-t(1); window_len = round(0.5 / dt); % 确保窗口覆盖至少0.5秒物理时间

这样，你的结果才能与高精度文献对标。

5. 从工具包到研究闭环：如何把它嵌入你的科研工作流

这套工具包的价值，远不止于“画几张好看的图”。在我自己的研究中，它已深度融入从问题发现到成果固化的全链条：

阶段一：快速探索（Rapid Exploration）
当你提出一个新系统，第一反应不是埋头推导，而是用y_phase.m和y1_bifurcation_C1_Forward_and_Backward.m快速扫一遍参数空间。我通常设置p_range = linspace(0,50,200)，粗扫200个点，5分钟内就能勾勒出大致的分岔骨架：哪里有周期窗口？哪里开始混沌？是否存在滞后环？这比花一周手推雅可比矩阵快得多，且能立刻发现反直觉现象——比如某个参数区间内，混沌态竟比邻近周期态更“稳定”（Lyapunov指数更负），这往往指向新的物理机制。

阶段二：深度诊断（Deep Diagnosis）
一旦锁定关键参数区间（如ρ=27.5~28.5），立即启动C2脚本，用initial_grid精细扫描初值空间。这时，y1_bifurcation_C2_CoexistingBifurcation.m输出的不仅是分岔图，更是吸引子盆地（basin of attraction）的拓扑快照。我把它的聚类结果导入Python，用networkx构建初值点到吸引子类的映射图，再用matplotlib绘制二维初值空间着色图——一张图就清晰显示：在x-y平面上，哪些区域初值导向单涡卷，哪些导向双涡卷，它们的边界是否分形？这直接支撑了我关于“多稳态竞争机制”的论文核心论点。

阶段三：成果固化（Result Solidification）
所有图表生成后，我绝不直接截图。而是用工具包内置的save_eps功能，导出.eps矢量图，再用Adobe Illustrator进行专业排版：统一字体（Times New Roman）、调整字号（坐标轴10pt，标题12pt）、添加比例尺、插入箭头标注关键特征点。Poincare.eps和bifurcation_diagram.png（我保留PNG用于补充说明）会被整合进同一张大图，形成“相图-截面-分岔”三位一体的证据链。审稿人看到的，不是孤立的图，而是一个自洽、可追溯、可复现的混沌行为证据包。

最后分享一个小技巧：在system.m中，我预留了一个params.custom_equation字段。当你要研究一个完全自定义的系统，又不想修改system.m源码时，可以这样用：

params.model_name = 'custom'; params.custom_equation = @(t,x,p) [p.a*(x(2)-x(1)); p.b*x(1)-x(1)*x(3); x(1)*x(2)-p.c*x(3)];

然后在system.m的switch末尾加一行：

case 'custom' dxdt = params.custom_equation(t, x, params);

三行代码，无限扩展。这套工具包，本质上是一个活的、生长的混沌研究操作系统。你每一次运行，都在为它注入新的物理直觉；每一次修改，都在加固你对非线性动力学本质的理解。它不承诺给你答案，但它确保，你提出的每一个问题，都能被清晰地看见、被严谨地检验、被漂亮地呈现。

本文还有配套的精品资源，点击获取

简介：一套开箱即用的MATLAB混沌分析工具，无需额外工具箱，直接运行即可绘制三维相图（y_phase.m）、庞卡莱截面（y_poincare.m，输出Poincare.eps）、单参数正向/反向分岔图（y1_bifurcation_C1_Forward_and_Backward.m）、双参数联合分岔图、以及初值敏感性引发的共存吸引子分岔（y1_bifurcation_C2_CoexistingBifurcation.m）。支持对不同初始条件（如y1_bifurcation_InitialValue_x.m）进行扫描，识别单涡卷与双涡卷吸引子的共存与演化；通过getmax_y1.m~getmax_y4.m系列函数自动提取状态变量极值序列，构建高精度分岔图；LTF_Bifur1_y1y2y3_text等脚本辅助多变量同步可视化。system.m统一定义动力学方程，模块化设计便于替换洛伦兹、Chen、Lü等常见混沌模型。所有脚本含中文注释，结构清晰，适配R2018a及以上标准MATLAB环境。

本文还有配套的精品资源，点击获取