1. 演化博弈与Lotka-Volterra模型的理论根基1.1 为什么演化博弈和LV模型会出现在同一个标题里先说清楚一个容易让新手绕晕的点演化博弈和Lotka-Volterra简称LV模型看起来是两个东西一个来自博弈论一个来自生态学但它们之间有一条非常深的数学纽带。早期研究里Taylor和Jonker在1978年提出的复制动态方程本质上就可以通过变量替换转化成Lotka-Volterra的形式。换句话说一个演化博弈系统完全可以写成一个捕食-被捕食或竞争共存的生态动力学系统。反过来很多生态学中的多物种竞争模型也能映射到演化博弈的分析框架里。所以你在MATLAB里做这两类模型的稳定点分析和相位图绘制底层用的是一套数学工具微分方程的平衡点求解、雅可比矩阵的特征值判定、相空间轨迹的可视化。这也是为什么把这俩放在一篇博文里讲而不是拆成两篇——学会了一套流程两个模型都能吃透。这篇内容面向的读者我默认是已经在读研或者做科研的大概率是做管理科学、系统工程、环境经济、生物数学这类方向的。你可能已经看过几篇演化博弈的论文知道复制动态方程长什么样也可能只是听说过LV模型但不太清楚怎么在MATLAB里把相位图画得既漂亮又可信。这篇就是把这两类模型的建模、分析、画图流程完整地过一遍代码直接可以改成你自己的参数用。1.2 复制动态方程与LV方程的数学共性先写一下双方演化博弈的标准形式。假设有两个种群或两个参与方各自有两个策略用x表示种群1选择策略A的比例y表示种群2选择策略B的比例那么复制动态方程一般写成dx/dt x * (1 - x) * (f1(x,y) - g1(x,y)) dy/dt y * (1 - y) * (f2(x,y) - g2(x,y))其中f1和g1分别是策略A和策略B的期望收益f2和g2同理。这里的关键是有x(1-x)这个项它保证了x永远在[0,1]区间内因为只要x不在0或1x(1-x)就大于0而符号由收益差决定。这个结构天然适合做演化博弈。再看LV模型的标准形式dx/dt x * (r1 a11*x a12*y) dy/dt y * (r2 a21*x a22*y)注意这里的x和y不再代表比例而是种群密度或规模。但数学结构上两者都是dx/dt x * 某个多项式的形式。如果把复制动态方程展开x(1-x)乘上收益差之后因为收益函数本身通常是x和y的线性函数展开后就会得到一个关于x和y的二次多项式形式上和LV方程非常接近。这意味着什么意味着你在MATLAB里写求解函数的时候可以用同一种模式去写——先定义dx x .* 多项式再把两边用ode45去积分。这个共性会让你写三方博弈的时候省很多事后面我会专门讲到怎么利用这个结构。1.3 双方演化的稳定点判定从雅可比矩阵到特征值稳定点分析是这类模型绕不开的核心步骤。道理很简单演化博弈关心的最终问题是系统会收敛到哪个稳定状态而稳定状态就是让所有dx/dt 0同时成立的解也就是平衡点。但平衡点不一定是稳定的得算雅可比矩阵把每个平衡点代进去求特征值。雅可比矩阵的形式是这样的J [d(dx/dt)/dx, d(dx/dt)/dy; d(dy/dt)/dx, d(dy/dt)/dy]算出特征值后判定规则如下特征值情况稳定性判定所有特征值实部均小于0渐近稳定点ESS所有特征值实部均大于0不稳定点源点特征值实部异号鞍点存在实部等于0的特征值临界情况需要更高阶判定竞品很多传统做法是手推雅可比矩阵再代值但MATLAB里完全可以自动完成。我的习惯是用符号计算工具箱先定义符号变量x和y把微分方程写成符号表达式然后调用jacobian函数求雅可比矩阵再subs代入数值平衡点最后用eig算特征值。这个流程快且不易出错尤其适合三方博弈——三方博弈的雅可比是3×3矩阵手推极其容易漏项。2. MATLAB实操以双方演化博弈为例的完整流程2.1 从收益矩阵到复制动态方程的编程实现先给一个具体的双方演化博弈例子。假设有两个参与方参与方1有策略U和V参与方2有策略C和D收益矩阵如下CDU(a, e)(b, f)V(c, g)(d, h)其中a,b,c,d是参与方1在不同策略组合下的收益e,f,g,h是参与方2的收益。对于参与方1选择策略U的期望收益是a*y b*(1-y)选择策略V的期望收益是c*y d*(1-y)。参与方2类似选择策略C的期望收益是e*x g*(1-x)选择策略D的期望收益是f*x h*(1-x)。把这个写成MATLAB的微分方程函数长这样function dydt twoPlayerReplicator(t, y, payoff1, payoff2) % y [x; y] x y(1); y_val y(2); % 参与方1的期望收益 fU payoff1(1)*y_val payoff1(2)*(1-y_val); % 策略U的收益 fV payoff1(3)*y_val payoff1(4)*(1-y_val); % 策略V的收益 f1_avg x*fU (1-x)*fV; % 参与方2的期望收益 fC payoff2(1)*x payoff2(2)*(1-x); % 策略C的收益 fD payoff2(3)*x payoff2(4)*(1-x); % 策略D的收益 f2_avg y_val*fC (1-y_val)*fD; dx x * (1-x) * (fU - fV); dy y_val * (1-y_val) * (fC - fD); dydt [dx; dy]; end这个函数是核心后面的数值求解和相位图都靠它。这里有个细节我把参与方2的策略C定义为对应x方向的策略所以当参与方1选择U、参与方2选择C时收益矩阵的(1,1)位置是payoff1(1)和payoff2(1)。如果你的矩阵定义不同注意对应关系别搞反。2.2 用ode45求解复制动态方程并画时间演化图回到命令行用ode45对这个系统做数值积分。先设置参数然后从不同的初始点出发看系统收敛到哪% 收益矩阵参数 a 3; b 0; c 0; d 1; % 参与方1 e 3; f 0; g 0; h 1; % 参与方2 payoff1 [a, b; c, d]; payoff2 [e, f; g, h]; % 时间范围 tspan [0 50]; % 多个初始点 init1 [0.2 0.2]; init2 [0.8 0.8]; init3 [0.1 0.9]; init4 [0.9 0.1]; allInit [init1; init2; init3; init4]; figure; hold on; for i 1:size(allInit, 1) [t, y] ode45((t,y) twoPlayerReplicator(t, y, payoff1, payoff2), tspan, allInit(i,:)); plot(t, y(:,1), LineWidth, 1.5); plot(t, y(:,2), --, LineWidth, 1.5); end xlabel(时间 t); ylabel(策略比例); title(双方演化博弈的时间演化); legend(x (参与方1策略U比例), y (参与方2策略C比例)); grid on;这里有个重要的实操经验ode45默认的误差容限是1e-3但演化博弈在接近鞍点或边界时会非常敏感我建议把误差容限调严一点比如options odeset(RelTol, 1e-6, AbsTol, 1e-8); [t, y] ode45(..., tspan, init, options);不调容限直接跑有时候会发现轨迹在稳定点附近振荡不收敛或者出现微小的负值。尤其是后面三方博弈的时候这个坑更明显。万一算出了负值别慌在相位图绘制时做了裁剪就行因为演化博弈的比例变量理论上必须在[0,1]数值误差可能导致越界这属于正常现象。2.3 双方博弈的相位图绘制向量场与轨迹叠加相位图是这个分析里可视化效果最直观的图。横轴是x纵轴是y把每个点上的运动方向用小箭头画出来再把几条代表性的轨迹叠上去。画向量场的常见做法是用meshgrid生成网格然后计算每个网格点的dx/dt和dy/dt用quiver函数画箭头[x_grid, y_grid] meshgrid(0:0.05:1, 0:0.05:1); dx zeros(size(x_grid)); dy zeros(size(y_grid)); for i 1:numel(x_grid) dydt_temp twoPlayerReplicator(0, [x_grid(i), y_grid(i)], payoff1, payoff2); dx(i) dydt_temp(1); dy(i) dydt_temp(2); end figure; quiver(x_grid, y_grid, dx, dy, 1.5, Color, [0.6 0.6 0.6], LineWidth, 0.8); hold on; % 叠加从不同初始点出发的轨迹 for i 1:size(allInit, 1) [t, y] ode45((t,y) twoPlayerReplicator(t, y, payoff1, payoff2), tspan, allInit(i,:), options); plot(y(:,1), y(:,2), LineWidth, 1.8); plot(y(1,1), y(1,2), o, MarkerSize, 8, MarkerFaceColor, r); plot(y(end,1), y(end,2), s, MarkerSize, 10, MarkerFaceColor, g); end xlabel(x (参与方1策略U比例)); ylabel(y (参与方2策略C比例)); title(双方演化博弈相位图); xlim([0 1]); ylim([0 1]); axis square; grid on;这一步有很关键的点同时也最容易被忽略网格密度和箭头缩放的配合。quiver的缩放系数1.5是我常用的经验值但如果你的收益参数很大比如超过10箭头的长度会急剧增加导致图面混乱。我的做法是先把dx和dy做归一化再乘一个固定的箭头长度这样不管参数多大图面永远是清爽的norm_dx dx ./ sqrt(dx.^2 dy.^2); norm_dy dy ./ sqrt(dx.^2 dy.^2); quiver(x_grid, y_grid, norm_dx, norm_dy, 0.5, Color, [0.6 0.6 0.6], LineWidth, 0.8);归一化之后箭头只表示方向不表示大小但这反而更适合看流场的整体走势。2.4 雅可比矩阵的符号计算与自动判定收益参数一变稳定点就可能变相位图也完全不一样。为了快速判定稳定点我建了个脚本用符号计算完成整个过程syms xs ys % 收益矩阵参数 a 3; b 0; c 0; d 1; e 3; f 0; g 0; h 1; fU a*ys b*(1-ys); fV c*ys d*(1-ys); fC e*xs g*(1-xs); fD f*xs h*(1-xs); dxs xs*(1-xs)*(fU - fV); dys ys*(1-ys)*(fC - fD); % 求雅可比矩阵 J jacobian([dxs; dys], [xs, ys]); % 候选平衡点边界点和内点 eq_points [0 0; 0 1; 1 0; 1 1]; % 内点通过数值求解 s vpasolve([dxs0, dys0], [xs, ys], [0.2 0.8; 0.2 0.8]); eq_points [eq_points; double([s.xs, s.ys])]; % 逐个判定 for k 1:size(eq_points, 1) Jk double(subs(J, [xs, ys], eq_points(k,:))); eigvals eig(Jk); fprintf(平衡点 (%.4f, %.4f): 特征值 %.4f, %.4f\n, eq_points(k,1), eq_points(k,2), real(eigvals(1)), real(eigvals(2))); if all(real(eigvals) 0) fprintf( 判定渐近稳定点 (ESS)\n); elseif all(real(eigvals) 0) fprintf( 判定不稳定点\n); else fprintf( 判定鞍点\n); end end这个脚本我几乎每个演化博弈课题都会用到改改收益参数就能跑。有几个需要注意的地方vpasolve必须给初始搜索区间不然可能漏掉内点解边界上的平衡点虽然在数学上不稳定或鞍点但实际仿真中轨迹可能会在边界上停留很长时间分析时要结合时间演化图一起看不能只依赖特征值判定。3. Lotka-Volterra模型的稳定点分析与相位图3.1 经典捕食-被捕食模型在MATLAB中的求解LV模型最经典的形态是捕食-被捕食模型。设x是猎物密度y是捕食者密度方程为dx/dt x * (alpha - beta * y) dy/dt y * (-gamma delta * x)这组方程有个著名的性质它的平衡点是中心型稳定轨迹绕着平衡点形成闭合轨道没有渐进稳定性。这说明经典LV模型过于理想化现实中绝大多数捕食系统是有阻尼或密度制约的。所以做稳定点分析的时候我一般会建议用带logistic项的修改版LV模型dx/dt x * (r1 - a11*x a12*y) dy/dt y * (r2 a21*x - a22*y)这样模型就有了内部密度制约项平衡点会变成螺旋焦点或结点轨迹会收敛或发散分析起来更有实际意义。MATLAB求解很简单同样是ode45function dydt LV_model(t, y, params) % params [r1, a11, a12, r2, a21, a22] r1 params(1); a11 params(2); a12 params(3); r2 params(4); a21 params(5); a22 params(6); x y(1); z y(2); % 用z表示捕食者避免与函数名y冲突 dx x * (r1 - a11*x a12*z); dz z * (r2 a21*x - a22*z); dydt [dx; dz]; end3.2 LV模型平衡点稳定性与相图绘制的核心差异LV模型的相图和演化博弈的相图有个重要的差异演化博弈的变量始终限制在[0,1]×[0,1]方块内而LV模型的变量是种群密度范围可以很大而且必须保证非负。所以画向量场的时候网格范围和密度要重新考虑params [1.5, 0.8, 0.6, 0.3, 0.5, 0.7]; x_range 0:0.1:3; z_range 0:0.1:3; [x_grid, z_grid] meshgrid(x_range, z_range); dx zeros(size(x_grid)); dz zeros(size(z_grid)); for i 1:numel(x_grid) dydt_temp LV_model(0, [x_grid(i), z_grid(i)], params); dx(i) dydt_temp(1); dz(i) dydt_temp(2); end figure; quiver(x_grid, z_grid, dx, dz, 1.2); hold on; % 数值求解多条轨迹 tspan [0 30]; for init_x [0.5 1.5 2.5] for init_z [0.5 1.5 2.5] [t, y_traj] ode45((t,y) LV_model(t, y, params), tspan, [init_x init_z], odeset(RelTol,1e-6)); plot(y_traj(:,1), y_traj(:,2), LineWidth, 1.5); end end % 标记平衡点 eq_pts [0 0; 3/0.8 0; 0 3/0.7; 0 0]; % 需要根据参数具体算 xlabel(猎物密度 x); ylabel(捕食者密度 z); title(Lotka-Volterra模型相位图); axis tight; grid on;LV模型的平衡点求解也是解代数方程组。零平衡点(0,0)总是存在非零平衡点可能有一个或多个。符号计算同样派得上用场syms xs zs eq1 xs*(r1 - a11*xs a12*zs); eq2 zs*(r2 a21*xs - a22*zs); solutions solve([eq10, eq20], [xs, zs]);这里特别注意零平衡点的稳定性取决于r的符号。如果r1 0且r2 0原点稳定两个种群都灭绝如果r1 0且r2 0原点沿猎物方向不稳定。这个结论在手算和代码里都要验证一下别想当然。我用过不少次这套流程处理带修正的LV模型发现一个很实用的技巧把平衡点画在相位图上时用scatter加不同颜色区分稳定和不稳定一眼就能看出系统的归宿。代码就是stable [1 0 0]; saddle [1 1 0]; for k 1:size(eq_pts,1) if isstable(eq_pts(k,:)) scatter(eq_pts(k,1), eq_pts(k,2), 100, filled, MarkerFaceColor, r); else scatter(eq_pts(k,1), eq_pts(k,2), 100, filled, MarkerFaceColor, y); end end3.3 从演化博弈复制动态到LV模型的转换方法这个转换是理解两类模型关系的关键。假设我们有一个2×2对称博弈复制动态方程是dx/dt x(1-x)((a-c) (bd-a-c)y)做变量替换u x/(1-x)把比例变量映射到正实数轴上数学推导可以得到du/dt u * ((a-c) (bd-a-c) * v/(1v))这个形式再经过整理就能改写成LV型的方程。虽然实际科研中很少真的去做这个变量替换去转换方程但理解这个等价关系能帮你做一件事把生态学中关于物种共存、竞争排斥的结论直接类比到演化博弈中或者反过来把演化博弈的均衡概念用到生态建模里。我在实际课题里跨界用这个关系的场景是从博弈模型推导出合作者比例的动态后用LV模型的既有理论判断参数在什么范围内会出现周期震荡。这对于管理科学中的囚徒困境惩罚机制这类模型特别有用——惩罚强度过大时系统可能出现不收敛的周期解此时用演化博弈的框架看是没有ESS用LV的框架看是极限环。4. 三方演化博弈的建模与相位图绘制4.1 三方博弈的复制动态方程结构三方演化博弈是当前论文里的热门方向因为大多数现实问题确实涉及三个参与方——比如政府、企业、公众三方协同治理或者供应商、制造商、零售商三方供应链。三方博弈和双方博弈的核心区别在于每个参与方的复制动态方程多了一个变量系统从二维变成三维。假设三个参与方的策略选择比例分别是x、y、z各自有收益函数f_i那么方程组是这个样子dx/dt x*(1-x)*(f1(x,y,z) - g1(x,y,z)) dy/dt y*(1-y)*(f2(x,y,z) - g2(x,y,z)) dz/dt z*(1-z)*(f3(x,y,z) - g3(x,y,z))其中f_i和g_i是各参与方两个策略的期望收益函数它们不再只是x和y的线性组合而是包含z的三元线性函数——因为参与方3的策略选择会直接影响参与方1和参与方2的收益。我在处理三方博弈时发现了一个让很多人崩溃的问题收益矩阵的维度。双方的收益矩阵是2×2×2每个参与方2个策略而三方的收益矩阵是2×2×2×3三个参与方每个2个策略加上参与方索引写起来极其容易出错。我的解决方法是把收益参数命名为结构体数组比如params.A_UC表示参与方1选择U、参与方2选择C、参与方3选择对应的某种策略时的收益这样写函数的时候能减少很多索引错误。4.2 三方系统平衡点搜索与高维雅可比矩阵判定三方博弈的平衡点数量比双方多得多。每个变量的取值可以是0、1或者内部解所以边界平衡点最多有2^3 8个加上可能存在的内点均衡实际要判定的平衡点可能超过10个。代码上我分成两步走。第一步先算所有边界组合——也就是x,y,z分别取0或1的所有组合combinations [0 0 0; 0 0 1; 0 1 0; 0 1 1; 1 0 0; 1 0 1; 1 1 0; 1 1 1];第二步对每一组边界组合代入复制动态方程检查该点是否满足所有dx/dt 0。注意如果某个变量是0或1它对应的方程自动为0所以只需要检查不是0或1的那个变量的方程是否为0。用这个筛选条件可以快速排除大量无效组合。筛选出平衡点后接着算3×3的雅可比矩阵。和双方博弈一样我还是用符号计算来求syms xs ys zs % 定义三个复制动态方程 dxs, dys, dzs % ... J jacobian([dxs; dys; dzs], [xs, ys, zs]); % 代入每个平衡点求特征值这里有一个实打实的教训3×3矩阵的特征值符号判定位数高而且经常出现一对共轭复根加一个实根的情况。判定规则变成了条件判定所有实部 0渐近稳定点所有实部 0不稳定点实部有正有负鞍点实部都 0 且存在虚部不为0稳定焦点螺旋收敛我在论文里见过很多人只写特征值均为负故稳定完全不提虚部的情况。如果特征值是复数系统在收敛过程中会呈现螺旋轨迹相位图上是转着圈进去的这个在图上非常明显如果不说明会显得分析不到位。另外三方博弈的相位图绘制是个硬骨头。三维相空间没法直接画完整的相位图业界通用做法是画三维投影图或切面图。我用得比较顺手的方案是固定第三个变量z在若干取值比如0.1, 0.5, 0.9在每个z平面上画x和y的二维向量场。这样既能看到系统在不同z层的行为差异又不用面对杂乱无章的三维箭头。z_fixed [0.1 0.5 0.9]; for z_val z_fixed [x2d, y2d] meshgrid(0:0.05:1, 0:0.05:1); dx2d zeros(size(x2d)); dy2d zeros(size(y2d)); for i 1:numel(x2d) dydt_temp threePlayerReplicator(0, [x2d(i), y2d(i), z_val], params); dx2d(i) dydt_temp(1); dy2d(i) dydt_temp(2); end subplot(1, 3, find(z_fixed z_val)); quiver(x2d, y2d, dx2d, dy2d, 1.5); title([z , num2str(z_val)]); end4.3 三方博弈中常见的周期解与收敛问题三方演化博弈最让人头疼的问题我实战下来觉得有两个。第一个是参数敏感导致的非收敛解——轨迹在一个闭曲线附近不断绕圈不趋于任何平衡点。这时候你的特征值判定会发现平衡点全都不稳定或者存在一对纯虚根。这个结果从理论上是合理的但在论文里不好呈现因为系统一直震荡意味着政策建议不好写。我遇到这种情况的习惯做法是试着调整参数找物理含义更合理的区域或者在论文里如实报告并解释震荡在现实中的含义。第二个是数值积分误差导致的伪周期解。三方系统比双方更容易出现误差积累所以我一般会提高ode45的容差甚至改用ode15s处理刚性问题。怎么判断是真实的周期还是数值误差把时间步长减半再算一次如果轨迹明显变样就是数值问题如果轨迹基本重合说明是系统的固有行为。相位图上还有个小技巧为了看清螺旋收敛可以在时间演化图的末尾只画最后一段时间比如t从20到50这样能更清楚地看到轨迹是否真的收敛到平衡点。5. MATLAB实现中的细节打磨与问题排查5.1 相位图的美化与可视化技巧相位图画得好不好直接影响审稿人对你工作的第一印象。我总结几条实战经验箭头分层处理。先用浅灰色画方向场再用粗线条画代表性的轨迹最后用醒目的标记标注平衡点。三层叠加层次感立刻出来。有一个常见的坑是箭头太密导致轨迹线条全被遮挡这时候可以把网格间距从0.05改成0.1或者把箭头缩放系数改小。轨迹线用渐变色。三条从不同起点出发的轨迹用同一个色系的不同深浅来画这样可以避免用红色、蓝色、绿色这种默认配色图的观感会提升一个档次。我常用的方案是colormap配Color属性colors [0.2 0.2 0.8; 0.4 0.6 0.9; 0.8 0.4 0.2]; for i 1:size(allInit, 1) [t, y] ode45(..., allInit(i,:), options); plot(y(:,1), y(:,2), Color, colors(i,:), LineWidth, 2); end字体和导出。xlabel、ylabel、title一定要指定字体大小不然默认字体在期刊排版时显得很小气。导出图片用exportgraphics输出300dpi的PNG或者矢量PDF千万不要截图或者用print的默认设置。5.2 常见报错为什么ode45算出来全是NaN这是我最常被问的问题。NaN的来源几乎都是一样的数值积分的初始步长太大导致某个变量在某个中间时刻越过零变成负值然后负值被传入收益函数后出现了log或sqrt的NaN。等一下——我们的复制动态方程按理说没有log和sqrt怎么会出现NaN问题出在收益函数的实现上。如果你在收益函数里用y^0.5这种幂运算或者用到log(1-y)那么在y略微超出[0,1]区间时就会出错。另外x*(1-x)看起来没问题但如果x是负的比如-0.01x*(1-x)是负的可能把轨迹推得更远最终发散。排查思路是固定的调小初始步长odeset(InitialStep, 1e-3)。检查收益函数里有没有非多项式运算比如log、sqrt、除零风险。在积分循环中对变量做钳位x min(max(x, 0), 1)虽然这改变了原方程但对演化博弈的比例变量来说是合理的数值处理。如果以上都没用打印中间变量看是哪一步开始出问题。5.3 符号计算与数值计算并用的加速策略三方博弈频繁求平衡点和特征值如果用纯数值方法找所有平衡点容易漏如果用纯符号方法算可能卡死在符号表达式的化简上。我的实践是混合策略先用符号计算推导复制动态方程的表达式jacobian求雅可比矩阵这个步骤非常快然后把符号表达式转成MATLAB函数句柄用matlabFunction数值积分和向量场绘制全部用函数句柄速度比直接喂符号表达式给ode45快一个数量级。% 符号推导后转成函数句柄 J_func matlabFunction(J, Vars, [xs, ys, zs]); % 在某个平衡点(peq)处快速求雅可比 J_num J_func(peq(1), peq(2), peq(3));这个技巧我强烈建议用特别是当你需要对多个平衡点反复判定或者要做参数敏感性分析扫几百组参数的时候。5.4 常见问题速查表现象可能原因解决办法ode45输出NaN变量越界或收益函数有非法运算调小InitialStep对变量钳位到[0,1]轨迹不收敛始终震荡系统存在极限环或数值容差过大调小RelTol和AbsTol检查特征值有无纯虚根平衡点漏检vpasolve搜索区间设置不当扩大搜索区间或用fsolve多初值扫描相位图箭头方向混乱网格间距过大或参数量级差异悬殊加密网格对箭头做归一化符号计算卡死高维符号表达式太复杂用assume简化或改用数值方法求雅可比三方博弈轨迹穿出[0,1]数值误差积累用钳位函数改用更稳定的ode15s6. 从论文到复现我的一些经验与提醒做演化博弈加LV模型的课题最大的优势在于图表好看、结论直观最大的坑在于看起来简单实际参数和代码细节一堆。收益矩阵的参数设定不能随手拍脑袋。很多论文的模型之所以能得出想要的均衡是因为参数选得很巧。我的习惯是先画一组大范围参数下的相位图观察系统行为随参数变化的趋势再去文献里找支撑参数取值的依据最后在论文里明确写出参数的现实含义。这方面审稿人很敏感上来就令a5, b3, c2却不解释容易被要求补说明。复现别人的模型时先别急着跑相位图。我建议第一件事是打印平衡点和特征值对照原文的判定结果确认完全一致后再画图。这一步能省掉你后面大量为什么我的图和论文长得不一样的困惑。因为很多时候差异不在求解代码而在矩阵的转置、收益函数的正负号这些细微之处。还有一个常被忽略的细节时间范围的设置。复制动态方程收敛的速度取决于特征值的实部大小如果参数接近临界值收敛时间可能非常长。我习惯先跑一次长时间(比如tspan[0 500])的模拟看轨迹在大约什么时间稳定下来再调整时间范围避免画出来的轨迹因为时间不够长而歪在半路上。最后再说说代码组织。这类项目我强烈建议按参数定义-方程函数-求解脚本-绘图脚本四段式的结构组织并且把参数独立封装成结构体。因为你会反复修改参数跑实验如果参数散落在各个脚本里改起来非常痛苦。这也是我做课题以来踩坑最深的一条经验分享出来给大家避坑。