做非线性系统仿真的人，几乎都绕不过相图这个工具。我第一次画相图是研究一个带阻尼的摆，当时用ode45算了一堆时域曲线，曲线确实不少，但系统最后到底是趋于静止还是持续振荡，光看那些波形很难一眼得出结论。后来把状态变量放到同一个平面里画相图，事情一下子清楚了：所有轨迹都在向同一个点收缩，系统的长期行为一目了然。从那以后，我养成了一个习惯——拿到非线性系统，先把相图画出来，再回去看时域曲线。这篇文章就把“利用Matlab绘制非线性系统相图”这件事，从理论原理讲到能直接运行的代码，适合正在学非线性动力学、做控制仿真、或者被课程作业逼着画相图的读者。
1. 非线性系统相图：从“为什么画”到“画的是什么”
1.1 相图到底在表达什么
相图的本质很简单：把系统的状态变量作为坐标轴，把系统随时间演化的轨迹画在这个状态空间里。以二阶自治系统为例，系统可以写成：
dx1/dt = f(x1, x2) dx2/dt = g(x1, x2)这里的核心关键词是“自治”，也就是方程右边不显含时间t。对于这类系统，给定一个初始状态，就有一条确定的轨迹线，这条轨迹线在状态空间中扫出的曲线，就是相轨迹。
相图就是大量相轨迹的组合。它不是某一条具体解曲线，而是系统所有可能运动的“地形图”。看相图就像看一张地图：哪里是盆地（吸引子）、哪里是山顶（不稳定点）、哪里是山口（鞍点），全都能直观看出来。对非线性系统来说，这价值太大了——因为绝大多数非线性微分方程没有解析解，相图成了少数能直接洞察系统定性行为的工具。
我常用一个比喻：相图是“风向图加旅行轨迹图”的结合体。每个点上都标着系统在这一点附近“下一步往哪走”的方向，这就是向量场；而每一条轨迹，就是一个小球顺着风在图上走出来的路径。
1.2 相图和时域波形之间的关系
很多初学者容易把相图和时域波形搞混。时域波形是x(t)随时间t变化的曲线，横轴是时间，纵轴是某个状态变量。相图则是把两个状态变量分别放到横轴和纵轴上，时间在这里不直接出现，而是隐含在轨迹的行进方向里。
举个例子，二阶线性系统：
dx1/dt = x2 dx2/dt = -x1时域波形是正弦和余弦，看起来在“振荡”；相图则是一个圆（或者椭圆），系统状态在圆周上匀速转圈。圆的半径由初始条件决定。
再看一个耗散系统：
dx1/dt = x2 dx2/dt = -x1 - 0.2*x2时域波形是衰减振荡，振幅不断缩小；相图则是一条向内螺旋的曲线，最终收敛到原点。这里，螺旋的方向和收敛速度，包含了时域波形中不容易直接看出来的信息。
我自己的实操习惯是：先用相图判断“系统的运动结构是什么”，再用时域波形看“演化速度有多快”。两者配合，信息量比单看任何一种曲线都大得多。
1.3 理论先行：平衡点、向量场与零倾线
画相图之前，我建议至少手算一遍以下三个概念，这会让你对结果有预判，而不是画出来一脸懵。
第一个是平衡点。平衡点是满足f(x1,x2)=0且g(x1,x2)=0的点。系统一旦处于平衡点，就静止不动。这就像地图上的“盆地底部”或“山顶”。
第二个是向量场。在状态空间每个点上，计算系统在该点的速度(f, g)，画成一个个小箭头，就得到向量场。箭头指向是运动方向，箭头长度是运动速度大小。
第三个是零倾线。零倾线是满足f(x1,x2)=0或g(x1,x2)=0的曲线。在前者上，轨迹方向是竖直的；在后者上，轨迹方向是水平的。零倾线把相空间划分成不同的区域，在每个区域内，轨迹的大致走向是确定的。它是手工画相图时代最重要的辅助工具，在Matlab里用contour或fimplicit可以直接画。
理论分析的核心任务，是在平衡点附近做线性化。计算雅可比（Jacobian）矩阵：
J = [df/dx1, df/dx2 dg/dx1, dg/dx2]然后把平衡点坐标代进去，求特征值。特征值的实部符号决定了这个平衡点的局部稳定性：
| 特征值情况 | 平衡点类型 | 局部稳定性 |
|---|---|---|
| 实部均为负 | 稳定结点或焦点 | 稳定 |
| 实部有正有负 | 鞍点 | 不稳定（但有稳定流形） |
| 实部均为正 | 不稳定结点或焦点 | 不稳定 |
| 实部为零，虚部非零 | 中心 | 临界稳定 |
这套理论的价值在于：画图之前你就知道哪些区域是关键区域，比如平衡点附近、鞍点附近，这些地方要加密网格或者多取初值。
2. Matlab绘制相图的前期准备与工具选型
2.1 为什么用Matlab画相图
画相图不是只有Matlab一种工具，Python加SciPy也能做，但Matlab在几个方面确实顺手。
第一，数值积分函数成熟。ode45、ode15s这些求解器经过大量工程验证，误差控制和刚性检测都做得不错，拿来就用。第二，二维可视化函数齐全。quiver、contour、fimplicit、plot这些函数组合起来，几乎覆盖了相图所需的所有元素。第三，交互式体验好。数据在变量管理器里随时查看，图形可以缩放旋转，对调试很有帮助。
如果你要做的只是一次性分析，用Matlab脚本几十分钟就能搞定。如果你是想把这套流程做成可复用的工具，Matlab的函数封装和脚本机制也很方便。
2.2 需要用到的核心函数与环境检查
画相图常用的函数并不多，新手上手以下这几个就够了：
ode45：求解非刚性常微分方程组，生成相轨迹的核心工具quiver：绘制二维向量场，表现系统各点的运动方向contour或fimplicit：绘制零倾线、隐式曲线plot、hold on、axis equal：绘制轨迹并控制图形坐标系meshgrid：生成向量场网格点odeset：设置求解器精度、步长等参数fsolve：数值求解平衡点
fimplicit是R2016b版本之后才有的函数。如果你用的是老版本，建议用contour替代，我会在后面给出替代写法。
建议你打开Matlab，先跑一下这个最简单的例子，测试环境是否正常：
% 环境自检：线性中心系统 f = @(t, x) [x(2); -x(1)]; [t, x] = ode45(f, [0 20], [1; 0]); figure; plot(x(:,1), x(:,2)); axis equal; xlabel('x_1'); ylabel('x_2'); title('线性中心：x1''=x2, x2''=-x1'); grid on;如果这段代码能画出一个圆形轨迹，说明你的Matlab基本环境没有问题，可以继续往下写。
2.3 把系统方程写成Matlab能吃的格式
ode45只能处理一阶常微分方程组，所以高阶方程必须先降阶。这个步骤新手很容易忽略。比如摆的方程：
θ'' + (b/m)θ' + (g/L)sin θ = 0令x1 = θ，x2 = θ'，写成：
dx1/dt = x2 dx2/dt = -(g/L)*sin(x1) - (b/m)*x2在Matlab里最直接的写法是匿名函数：
g = 9.81; L = 1.0; b = 0.1; m = 1.0; f = @(t, x) [x(2); -(g/L)*sin(x(1)) - (b/m)*x(2)];这里x是列向量，x(1)就是x1，x(2)就是x2，返回的列向量第一个元素是dx1/dt，第二个是dx2/dt。匿名函数的好处是定义在脚本里，参数可以自由捕获，修改参数时不用到处找函数文件。
如果系统比较复杂，或者你打算反复使用，建议写成独立函数文件：
function dx = pendulumODE(t, x, g, L, b, m) dx = zeros(2,1); dx(1) = x(2); dx(2) = -(g/L)*sin(x(1)) - (b/m)*x(2); end调用时用函数句柄传参：
f = @(t, x) pendulumODE(t, x, g, L, b, m);关键点就一句话：所有高阶微分方程，进ode45之前先降阶为状态空间形式。
3. 核心实操：向量场、零倾线与轨迹的组合绘制
3.1 第一步：用quiver画向量场
向量场是相图的地基。用meshgrid在状态空间生成网格点，然后在每个网格点计算(f, g)，用quiver绘制。
我以Van der Pol振荡器为例。系统方程为：
dx/dt = y dy/dt = μ(1 - x^2)y - x先画向量场：
% 参数与网格设置 mu = 1.0; xrange = -3:0.4:3; yrange = -3:0.4:3; [X, Y] = meshgrid(xrange, yrange); % 计算向量场 U = Y; V = mu * (1 - X.^2) .* Y - X; % 绘制向量场 figure; quiver(X, Y, U, V, 'Color', [0.6 0.6 0.6]); axis equal; xlim([-3.5 3.5]); ylim([-3.5 3.5]); xlabel('x'); ylabel('y'); grid on;quiver的第五个参数是比例因子，默认情况下会自动缩放箭头长度。如果箭头太长、叠成一团黑，可以把比例因子调小，比如quiver(X, Y, U, V, 0.6)。
更常用的做法是归一化方向，让每个箭头等长，这样能清晰表达方向，不会因为某一点速度太大而覆盖其他信息：
L = sqrt(U.^2 + V.^2); quiver(X, Y, U./L, V./L, 0.6, 'Color', [0.5 0.5 0.5]);归一化之后，速度大小信息暂时丢了。想要保留速度大小，可以用背景色表达。pcolor可以先画速度大小的底色，再叠加归一化箭头：
figure; hold on; pcolor(X, Y, L); colormap(parula); shading interp; quiver(X, Y, U./L, V./L, 0.6, 'Color', [0.2 0.2 0.2]); axis equal; colorbar;不过背景色如果控制不好会显得很花哨，我一般只在分析时这么做，正式出图时还是以干净的箭头图为主。
3.2 第二步：用ode45生成相轨迹
向量场给出了“路标”，相轨迹就是真正“走出来的路”。用ode45从某个初值积分系统，把得到的(x, y)点画在相平面上。
hold on; % 初值列表：多个初值才能看出全局结构 x0_list = [-2.5 -0.5 0.5 2.5]; for x0 = x0_list [t, x] = ode45(f, [0 30], [x0; 0]); plot(x(:,1), x(:,2), 'LineWidth', 1.5); end这段代码会从x轴上的四个不同初值出发，各自算出一条轨迹。你会看到，有的轨迹从外面向内转，有的从里面向外转，最后都会逼近同一个闭合曲线，这就是极限环。
这里有几个经验点需要强调一下。第一个是积分时间TSPAN的选择。[0 30]表示从t=0积分到t=30。太短的话轨迹还没跑完，看不出长期行为；太长的话轨迹会在极限环上转几十圈，线条重叠又密又乱，白白增加计算量。我一般先用较短时间试跑一次，观察轨迹快收敛时的时间点，然后再定最终时间。第二个是初值选取。建议在平衡点附近、远离平衡点的地方都取几个初值，一组轨迹同时覆盖局部和全局行为，信息量最大。第三个是同一张图上叠加多条轨迹时，用不同的颜色区分。Matlab默认颜色循环已经够用，如果轨迹多，可以用lines或parula色图手动指定。
3.3 第三步：叠加零倾线与平衡点标记
光有向量场和轨迹，还不够“理论”。把零倾线加上去，分析会更清晰。
零倾线是f=0和g=0的等值线。在Matlab里可以直接画：
% 零倾线：dx/dt = 0 和 dy/dt = 0 figure; hold on; % 方法一：新版Matlab用fimplicit fimplicit(@(x, y) y, [-3.5 3.5 -3.5 3.5], 'k', 'LineWidth', 1.5); % dx/dt = 0 fimplicit(@(x, y) mu*(1-x.^2).*y - x, [-3.5 3.5 -3.5 3.5], 'k--', 'LineWidth', 1.5); % dy/dt = 0如果fimplicit在你的版本里不可用，用contour代替：
[X, Y] = meshgrid(-3.5:0.05:3.5, -3.5:0.05:3.5); U = Y; V = mu * (1 - X.^2) .* Y - X; contour(X, Y, U, [0 0], 'k', 'LineWidth', 1.5); contour(X, Y, V, [0 0], 'k--', 'LineWidth', 1.5);注意contour画等值线时，[0 0]表示只画数值为0那一层。
平衡点用fsolve求。Van der Pol系统只有一个平衡点(0,0)，但复杂系统通常有多个，需要从不同初始猜测出发多次求解：
% 求平衡点 fun = @(x) [x(2); mu*(1-x(1)^2)*x(2) - x(1)]; x_eq = fsolve(fun, [0; 0]); % 在图上标记 plot(x_eq(1), x_eq(2), 'ro', 'MarkerFaceColor', 'r', 'MarkerSize', 8);求到平衡点之后，用雅可比矩阵做局部线性化：
syms xs ys J = jacobian([ys; mu*(1-xs^2)*ys - xs], [xs, ys]); J_eq = double(subs(J, [xs, ys], [x_eq(1), x_eq(2)])); eig(J_eq)Van der Pol在原点处的雅可比矩阵是：
J = [0, 1 -1, μ]特征值为(μ ± sqrt(μ^2 - 4)) / 2。当μ>0时，特征值实部为正，原点不稳定。这个结论跟相图完全吻合：轨迹从原点附近出发，会螺旋向外，最终被极限环捕获。
3.4 第四步：调整视觉参数，让相图真正可读
基础图形画出来之后，视觉调整决定了这张图是能直接放进论文，还是只能自己看个大概。我总结了几条高频调整项。
axis equal必须加。如果不加，Matlab会自动按数据范围缩放横纵轴，一个本来圆形的极限环可能被拉伸成椭圆，这是新手最容易踩的坑。
范围控制也很重要。xlim和ylim要结合系统状态范围手动设置，避免轨迹画出视野，也避免空白太多。网格密度要适当，网格太密，箭头挤成一团；太稀，看不出方向变化。经验值是整个绘制区间内网格数控制在15×15到25×25之间，效果比较均衡。
多轨迹时，建议用颜色循环和线宽区分：
ax = gca; ax.ColorOrder = lines(7);轨迹末端加箭头标注方向：
hold on; for x0 = x0_list [t, x] = ode45(f, [0 20], [x0; 0]); plot(x(:,1), x(:,2), 'LineWidth', 1.5); % 在轨迹末端加箭头 quiver(x(end-1,1), x(end-1,2), x(end,1)-x(end-1,1), x(end,2)-x(end-1,2), 0, 'Color', [0.3 0.3 0.3]); endquiver的第三个参数是比例因子，这里设为0，表示不缩放，箭头长度就是最后两步的实际位移。
出图格式建议用高分辨率：
exportgraphics(gcf, 'phase_portrait.png', 'Resolution', 300);老版本没有exportgraphics，用print(gcf, '-dpng', '-r300', 'phase_portrait.png')。
4. 经典案例拆解：Van der Pol振荡器的极限环
4.1 系统模型与理论预期
Van der Pol振荡器是非线性动力学里最经典的例子之一，方程为：
dx/dt = y dy/dt = μ(1 - x^2)y - x当μ=0时，系统退化为线性中心，相图是一圈圈同心圆。当μ>0时，非线性项μ(1 - x^2)y起作用：在|x|<1范围内，1-x^2>0，阻尼是负的，系统从环境中吸收能量，振幅增大；在|x|>1范围内，1-x^2<0，阻尼是正的，系统耗散能量，振幅减小。这两种趋势平衡的结果，就是出现一个稳定的极限环。
理论上，平衡点在原点。当0<μ<2时，原点是不稳定焦点，特征值是一对实部为正的共轭复根，轨迹螺旋向外；当μ>2时，原点变成不稳定结点。无论哪种情况，所有轨迹最终都收敛到同一个极限环上。这个“从任意初始状态都收敛到同一闭合曲线”的现象，是线性系统永远不会有的，也是相图最能体现非线性魅力的地方。
4.2 完整可运行代码
下面这套代码可以直接复制运行，覆盖向量场、零倾线、平衡点和多条轨迹，生成一张标准的Van der Pol相图。
% Van der Pol 振荡器相图 % 系统：dx/dt = y, dy/dt = mu*(1-x^2)*y - x clear; close all; clc; % 参数 mu = 1.0; % 系统方程 f = @(t, x) [x(2); mu*(1 - x(1)^2)*x(2) - x(1)]; % 图1：向量场 + 零倾线 + 轨迹 figure('Position', [100 100 600 500]); hold on; % 向量场 [X, Y] = meshgrid(-3:0.4:3, -3:0.4:3); U = Y; V = mu * (1 - X.^2) .* Y - X; L = sqrt(U.^2 + V.^2); quiver(X, Y, U./L, V./L, 0.6, 'Color', [0.8 0.8 0.8], 'LineWidth', 0.8); % 零倾线 fimplicit(@(x, y) y, [-3.5 3.5 -3.5 3.5], 'k', 'LineWidth', 1.2); fimplicit(@(x, y) mu*(1 - x.^2).*y - x, [-3.5 3.5 -3.5 3.5], 'k--', 'LineWidth', 1.2); % 平衡点 plot(0, 0, 'ro', 'MarkerFaceColor', 'r', 'MarkerSize', 8); % 轨迹（多个初值） x0_list = [-2.5 -1.5 -0.5 0.5 1.5 2.5]; colors = lines(length(x0_list)); for i = 1:length(x0_list) [t, x] = ode45(f, [0 30], [x0_list(i); 0]); plot(x(:,1), x(:,2), 'Color', colors(i,:), 'LineWidth', 1.5); quiver(x(end-1,1), x(end-1,2), x(end,1)-x(end-1,1), x(end,2)-x(end-1,2), ... 0, 'Color', colors(i,:), 'LineWidth', 1); end % 图例和格式 xlim([-3.5 3.5]); ylim([-3.5 3.5]); axis equal; grid on; xlabel('x'); ylabel('y'); title(['Van der Pol 相图, \mu = ', num2str(mu)]); set(gca, 'FontSize', 12); % 导出 exportgraphics(gcf, 'vanderpol_phase.png', 'Resolution', 300);运行之后，你会看到灰色箭头构成向量场，黑色实线和虚线是零倾线，红色圆点是平衡点，六条彩色轨迹从不同位置出发，最终都缠绕到一个闭合曲线上。这个闭合曲线就是极限环。
如果轨迹画得太密、线条堆叠，可以把积分时间从30缩短到15，或者减少初值数量。如果希望轨迹更平滑，用odeset提高精度：
opts = odeset('RelTol', 1e-6, 'AbsTol', 1e-8); [t, x] = ode45(f, [0 30], [x0; 0], opts);4.3 参数变化下的相图演化
Van der Pol系统最有意思的地方，是参数μ对相图结构的影响。把不同μ的相图并排画出来，可以看到极限环从圆变成张弛振荡的整个过程。
mu_list = [0 0.5 1 2 5]; figure('Position', [100 100 900 600]); for k = 1:length(mu_list) mu = mu_list(k); f = @(t, x) [x(2); mu*(1 - x(1)^2)*x(2) - x(1)]; subplot(2, 3, k); hold on; % 向量场 [X, Y] = meshgrid(-4:0.5:4, -4:0.5:4); U = Y; V = mu * (1 - X.^2) .* Y - X; L = sqrt(U.^2 + V.^2); quiver(X, Y, U./L, V./L, 0.5, 'Color', [0.8 0.8 0.8]); % 轨迹 for x0 = [0.5 2.5] [t, x] = ode45(f, [0 50], [x0; 0]); plot(x(:,1), x(:,2), 'b', 'LineWidth', 1.2); end axis equal; xlim([-4 4]); ylim([-4 4]); grid on; title(['\mu = ', num2str(mu)]); end观察这组图你会发现：μ=0时相图是同心圆，没有极限环；μ=0.5时极限环接近圆形；μ越大，极限环越“扁”，轨迹在x接近±1时出现急剧的方向转折，这就是张弛振荡的特征。从应用角度看，这种振荡行为在电子振荡器、生物节律模型里都有体现。用Matlab把参数扫一遍，比读理论推导直观得多。
5. 常见问题与排查技巧实录
5.1 轨迹直接发散或出现NaN
这是画相图遇到最多的一个问题。表现是ode45算出来的轨迹直接飞向无穷，或者图形上出现断线、NaN警告。
最常见的原因是系统本身不稳定，比如非线性项有正反馈，初始状态又在吸引域之外。这时候缩短积分时间，比如从[0 100]改成[0 10]，先看轨迹早期走向，再确认是否是数值问题。
另一个原因是积分步长失控。ode45是变步长算法，遇到陡峭变化时会自适应缩小步长，但遇到刚性系统或者参数突变时，可能判断失误。处理办法是增加求解器精度设置：
opts = odeset('RelTol', 1e-6, 'AbsTol', 1e-8, 'MaxStep', 0.1); [t, x] = ode45(f, [0 30], x0, opts);如果还不行，检查系统是否刚性。Van der Pol在μ很大时就是典型的刚性系统，ode45会算得很慢甚至失败，这时候改用ode15s：
[t, x] = ode15s(f, [0 30], x0, opts);5.2 向量场箭头糊成一团或稀疏不均
网格太密是箭头糊成一团的主因。网格步长从0.4改成0.8，箭头数量立刻下降，画面清爽很多。网格太稀则看不出方向变化，这个自己调一两次就有感觉。
quiver的缩放参数也值得多试。默认自动缩放有时很激进，箭头会交叉重叠。可以试试关闭自动缩放，直接用单位方向向量：
quiver(X, Y, U./L, V./L, 0.5);这样所有箭头长度相同，只表达方向，视觉上最干净。想要同时表达速度大小，就把scale设为0.5~0.8，让速度大的点箭头略长、速度小的点箭头略短。
5.3 轨迹线看不出时间信息
相图上轨迹是曲线，时间信息隐藏在行进方向里。想看系统在某个时刻大致走到哪里，有几种办法。
最简单的是在轨迹上叠加等时间间隔的圆点：
[t, x] = ode45(f, [0 30], x0); step = 200; plot(x(1:step:end,1), x(1:step:end,2), 'o', 'MarkerSize', 4, 'Color', [0.2 0.2 0.8]);点之间隔越远，说明系统移动越快；隔越近，说明系统接近停滞。对极限环来说，你会看到系统在环的一部分快速移动、另一部分缓慢爬行。
另一种方法是对轨迹分段着色。把时间轴切成若干段，每段用不同颜色画，能直观看到轨迹演化的先后顺序。这个方法在展示慢快系统时特别有效。
5.4 同一个图叠加多条轨迹后信息混乱
初值取太多、颜色区分度不够、线型全一样，是轨迹图变乱的三个原因。解决思路是克制。
初值数量控制在5到8个，均匀分布在平衡点附近和远处。颜色用lines或hsv色图循环，不要所有轨迹都用同一种颜色。线型可以区分：从内部出发的轨迹用实线，从外部出发的轨迹用虚线，图例会清晰很多。
如果还是乱，可以把图拆成子图，或者只保留代表性轨迹。画相图的目的是看清结构，不是把所有可能轨迹都画上去。
5.5 版本差异导致函数不可用
fimplicit在R2016b之前不存在，exportgraphics在R2020a附近才加入。遇到版本问题，优先用更通用的替代方案。
fimplicit可以用contour替代，前面已经写过。exportgraphics可以用print替代。pcolor、quiver、ode45这些老牌函数，在近二十年的版本里都稳定可用。
如果你在旧版本里遇到ode45的Name, Value传参方式不一致的问题，统一改用odeset结构体传参，兼容性最好。
最后分享一个我自己的使用习惯：画相图之前，先手算一次平衡点和线性化特征值，心里有数之后再开Matlab。这样图还没画出来，你已经知道大概该看到什么。如果图的结果和理论预期对不上，要么是代码有bug，要么是系统本身存在数值解之外的复杂行为，比如分岔或者混沌。这时候相图的价值就真正体现出来了——它不是点缀，而是探索非线性系统行为的第一把钥匙。
