做非线性系统仿真的人几乎都绕不过相图这个工具。我第一次画相图是研究一个带阻尼的摆当时用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(线性中心x1x2, 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)就是x1x(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]表示从t0积分到t30。太短的话轨迹还没跑完看不出长期行为太长的话轨迹会在极限环上转几十圈线条重叠又密又乱白白增加计算量。我一般先用较短时间试跑一次观察轨迹快收敛时的时间点然后再定最终时间。第二个是初值选取。建议在平衡点附近、远离平衡点的地方都取几个初值一组轨迹同时覆盖局部和全局行为信息量最大。第三个是同一张图上叠加多条轨迹时用不同的颜色区分。Matlab默认颜色循环已经够用如果轨迹多可以用lines或parula色图手动指定。3.3 第三步叠加零倾线与平衡点标记光有向量场和轨迹还不够“理论”。把零倾线加上去分析会更清晰。零倾线是f0和g0的等值线。在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^20阻尼是负的系统从环境中吸收能量振幅增大在|x|1范围内1-x^20阻尼是正的系统耗散能量振幅减小。这两种趋势平衡的结果就是出现一个稳定的极限环。理论上平衡点在原点。当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要么是系统本身存在数值解之外的复杂行为比如分岔或者混沌。这时候相图的价值就真正体现出来了——它不是点缀而是探索非线性系统行为的第一把钥匙。