1. 稳定状态模型在竞赛题里的真实地位先说个我自己的体验。参加建模比赛多了你会发现很多题目表面上在问某条曲线最后走到哪、某种传染病会不会爆发、种群数量会不会稳定在一个值本质上都在考同一件事——系统的长期行为。这种题你如果只会拿着 ode45 从头到尾模拟一遍画出曲线然后说它收敛了那评委顶多给你个基础分。真正得分的关键是你能不能把为什么收敛收敛到什么状态参数变化后这个状态会不会突变这套逻辑说清楚而这正是稳定状态模型steady-state model干的事情。所谓稳定状态通俗讲就是一个动态系统经历足够长时间之后变量不再随时间变化或者围绕某个基准小幅波动系统定格在某个状态上。数学表达很直白对连续时间系统就是令导数等于零$$\frac{dx}{dt} f(x) 0$$满足这个方程的点 $x^*$ 叫平衡点或平衡状态。它是系统的候选终点但并不是所有候选终点都真的能到达——有些平衡点你一靠近它就跑了这就是不稳定平衡点。稳定状态模型的核心工作就是找出所有平衡点再判断哪些是真正会被到达的稳定状态。这类模型能覆盖的题目范围非常广种群生态类捕食者与被捕食者数量是否会达到平衡引入天敌后生态会不会崩溃传染病传播类染病人数最终归零还是持续存在也就是清零阈值在哪经济运行类供需价格会不会收敛到均衡价格调整周期和阻尼怎么设计化学反应类反应物浓度会不会趋向某个平衡浓度网络传播类信息、舆情、谣言的扩散规模和消退速度。一句话只要题目里出现长期最终是否会临界条件这类词基本都可以往稳定状态模型上靠。适合的读者也不限于竞赛党写毕业论文做动力学仿真的、做控制系统的、研究生物数学的这套思路都能直接用。2. 稳定性的数学判据从一维到二维的完整推演2.1 一维系统为什么要看切线的斜率先来最简单的一维情形方程长这样$$\frac{dx}{dt} f(x)$$平衡点满足 $f(x^*) 0$。但它稳定不稳定怎么看我当年学的时候老师给过一个特别好的生活类比——山谷和山顶的区别。你往山谷里放一颗小球轻轻推它一下它会滚回谷底往山顶上放一颗碰一下它就滚没影了。平衡点从数学上就是斜率等于零的地方但斜率等于零有三种情况穿过零稳定、正切过零不稳定、还有既有正又有负的鞍点。严格判据其实就一条在平衡点附近把系统线性化看 $f(x^*)$ 的符号。$$\frac{dx}{dt} \approx f(x^)(x - x^)$$如果 $f(x^) 0$那么扰动 $\Delta x x - x^$ 会满足 $\Delta x f(x^)\Delta x$这是个指数衰减的过程扰动被压下去平衡点渐近稳定。反过来如果 $f(x^) 0$扰动指数放大不稳定。如果恰好等于零恭喜你遇上了临界情形需要看更高阶导数才能判断。这个判据用起来极其顺手。比如后面要讲的 Logistic 人口模型$f(x) rx(1 - x/K)$在 $x^0$ 处算导数得到 $r$正数所以 0 是不稳定的——意味着种群不会自然灭绝在 $x^K$ 处导数是 $-r$负数稳定——种群数量最终会稳定在环境容纳量附近。整个分析不需要任何复杂计算一个导数符号就把模型的核心行为说透了。2.2 二维与高维系统Jacobian 矩阵的特征值说了算到了二维及以上事情变复杂了但思路完全一样。系统$$\frac{dx_1}{dt} f_1(x_1, x_2), \quad \frac{dx_2}{dt} f_2(x_1, x_2)$$在平衡点 $(x_1^, x_2^)$ 附近做线性化得到的系数矩阵叫Jacobian 矩阵$$J \begin{bmatrix} \frac{\partial f_1}{\partial x_1} \frac{\partial f_1}{\partial x_2} \ \frac{\partial f_2}{\partial x_1} \frac{\partial f_2}{\partial x_2} \end{bmatrix}$$然后算 $J$ 的特征值。只要所有特征值的实部都小于零平衡点就是渐近稳定的只要有一个实部大于零就不稳定。这个结论对任意维数都成立可以说是整个稳定状态模型的地基。原理也不难理解线性化之后系统解是特征值的指数组合 $e^{\lambda t}$实部为负的项随时间衰落实部为正的项随时间增长仅此而已。竞赛场上我们经常遇到的是二维系统手算特征值往往太麻烦这时候有两个偷懒捷径。第一个是直接让 Matlab 的eig帮你算简单粗暴第二个是利用行列式和迹的判别法对于 $2\times2$ 矩阵记 $T\text{tr}(J)$$D\det(J)$则有条件稳定性结论$D 0$ 且 $T 0$渐近稳定两个特征值实部为负$D 0$ 且 $T 0$不稳定两个特征值实部为正$D 0$鞍点不稳定一正一负实部$T^2 - 4D 0$焦点螺旋收敛或发散看 $T$ 定稳定这个判别法竞赛手算、写论文推导公式时特别有用。很多参考答案为了严谨会硬算特征值你直接用 $T$、$D$ 符号说明推导路径又短又不容易错。2.3 别忽视的三种稳定定义还有一个很多教材不细讲、但实际比比皆是容易踩坑的概念。你说稳定到底指哪种稳定渐近稳定asymptotically stable扰动后系统会回到平衡点。这是大家默认的意义也是建模题里最常用的。Lyapunov 稳定stable扰动后系统不会跑远但也不一定回到平衡点可能绕着它转圈。不稳定unstable微小扰动就导致系统远离平衡点。经典例子就是无阻尼摆和理想 Lotka-Volterra 模型——它们的平衡点附近是闭合轨道扰动后围成新的轨道既不接近也不远离。这种稳定但不渐近稳定的情况如果你只做数值模拟很容易误判成一直在波动然后怀疑自己代码写错了。实际上这是系统的真实结构特征。竞赛里遇到极限环或周期振荡类描述说的就是这个。3. Matlab 实操框架求平衡点、算稳定性、画相图一整套3.1 用符号工具箱求平衡点Matlab 做稳定状态分析我的习惯是先符号后数值两板斧。第一步永远先用 Symbolic Math Toolbox 把平衡点解析求出来能算出闭式解最好可以直接看出参数影响算不出来再上vpasolve或fsolve做数值兜底。以一维 Logistic 模型为例syms x r K f r*x*(1 - x/K); % 求平衡点 eq_points solve(f 0, x) % 计算导数并代入各平衡点判断稳定性 df diff(f, x); for i 1:length(eq_points) val double(subs(df, x, eq_points(i))); fprintf(平衡点 x%sf(x*)%s%s\n, ... char(eq_points(i)), char(simplify(subs(df, x, eq_points(i)))), ... string(val 0)); end输出大概是eq_points 0 K 平衡点 x0f(x*)r稳定条件为 0 平衡点 xKf(x*)-r稳定条件为 1注意这里有个细节符号工具箱返回的r是符号变量不能直接用val 0判断。我上面的写法里double()只有在r已赋值时才有效否则会报错。竞赛实践中我更推荐下面这种写法——给参数赋具体值再判断r_val 0.5; K_val 100; eq_points_numeric double(subs(eq_points, {r, K}, {r_val, K_val})); df_numeric matlabFunction(subs(df, {r, K}, {r_val, K_val})); for i 1:length(eq_points_numeric) lambda df_numeric(eq_points_numeric(i)); fprintf(x* %.2ff(x*) %.4f稳定判定%s\n, ... eq_points_numeric(i), lambda, ... string(lambda 0)); end手工推导公式用符号法讲数值结论用赋值后的数值法别混用。3.2 二维系统的 Jacobian 自动计算二维系统的手动偏导算起来容易出错Matlab 里可以自动搞定syms x1 x2 a b c d f1 a*x1 b*x2; f2 c*x1 d*x2; % 自动求 Jacobian J jacobian([f1; f2], [x1, x2]); % 在特定平衡点处赋值 x1_star 0; x2_star 0; J_num double(subs(J, {x1, x2, a, b, c, d}, {x1_star, x2_star, -2, 1, 0.5, -1})); eig(J_num)竞赛中很多系统是非线性的比如 $f_1 x_1(1 - x_1) - 2x_1x_2$ 这种jacobian照样一把梭比自己手写 $\partial f/\partial x$ 靠谱得多。拿到 $J$ 之后算特征值、算 $T$ 和 $D$、判断类型一气呵成。这里我建议你把上面这段封装成一个函数checkStability(f, vars, point, params)决赛三天里你会反复调用它。3.3 数值模拟用 ode45 验证你的判断符号分析给出的是数学结论但它考试和论文里不好直观展示。这时候用ode45做数值模拟把时间序列画出来既验证分析结果也方便评委直观理解。继续用二维 Lotka-Volterra 捕食者-被捕食者模型举例% 参数: a1.1, b0.4, c0.1, d0.4 % dx/dt a*x - b*x*y (猎物) % dy/dt -c*y d*x*y (捕食者) a 1.1; b 0.4; c 0.1; d 0.4; lv (t, y) [a*y(1) - b*y(1)*y(2); ... -c*y(2) d*y(1)*y(2)]; % 多组初值对比 tspan [0 60]; init_conds [2 1; 3 1; 5 2; 1.5 0.8]; figure; for i 1:size(init_conds, 1) [t, y] ode45(lv, tspan, init_conds(i, :)); plot(t, y(:,1), LineWidth, 1.2); hold on; plot(t, y(:,2), --, LineWidth, 1.2); end xlabel(时间 t); ylabel(种群数量); legend(猎物 x,捕食者 y,Location,best); title(Lotka-Volterra 模型不同初值下的时间序列); grid on;画完这张图你立刻能看到种群数量呈现周期性振荡不会收敛到某一点。这就是Lyapunov 稳定但不渐近稳定的直观体现。你再用quiver画相图会更加直观——所有轨迹绕着平衡点$x^c/d$$y^a/b$转圈既不靠近也不远离。数值模拟在这个案例里不是验证收敛而是验证不收敛这本身就是一个重要结论。3.4 用 quiver 画相图和零线相图是稳定状态模型最漂亮的表达方式也是论文里最能一眼出效果的图。基本套路是先在网格点上算向量场方向再叠加几条实际轨迹。% 在状态空间画向量场 [x1g, x2g] meshgrid(0:0.2:5, 0:0.2:5); dx1 a*x1g - b*x1g.*x2g; dx2 -c*x2g d*x1g.*x2g; figure; quiver(x1g, x2g, dx1, dx2, 1.5, Color, [0.6 0.6 0.6]); hold on; % 叠加几条轨迹用刚算好的数值解 for i 1:size(init_conds, 1) [t, y] ode45(lv, tspan, init_conds(i, :)); plot(y(:,1), y(:,2), LineWidth, 1.5); end % 标记平衡点 plot(c/d, a/b, ro, MarkerSize, 8, MarkerFaceColor, r); text(c/d0.1, a/b, 平衡点, FontSize, 10); xlabel(猎物 x); ylabel(捕食者 y); title(Lotka-Volterra 相图轨迹绕平衡点做闭合运动); grid on; axis tight;如果系统是收敛型的比如 SIR 传染病模型相图里你会看到所有轨迹最终都汇入一条走廊最后指向平衡点。这时候建议把**零线nullcline**也画上去——零线就是 $\dot{x}_10$ 和 $\dot{x}_20$ 的曲线它们的交点就是平衡点。在低维系统里零线是理解轨迹走向的地图比单纯堆参数直观得多。用fimplicit可以很方便地画零线figure; fimplicit((x1, x2) a*x1 - b*x1.*x2, [0 5 0 5], b, LineWidth, 1.5); hold on; fimplicit((x1, x2) -c*x2 d*x1.*x2, [0 5 0 5], r--, LineWidth, 1.5); legend(\it dx/dt\rm0, \it dy/dt\rm0);零线的交点一眼就能看出来——这就是平衡点的几何意义。4. 两个完整案例Logistic 种群模型与 SIS 传染病模型4.1 案例一Logistic 种群模型的全流程分析这个模型几乎是我向所有新手推荐的第一道稳定状态练手题因为它的分析路径标准、结论直观、Matlab 代码短而且完全能映射到竞赛题里。模型本身一句话种群增长受资源限制增长率不是常数而是随数量增大而降低。$$\frac{dN}{dt} rN\left(1-\frac{N}{K}\right)$$其中 $r$ 是内禀增长率$K$ 是环境容纳量。第一步求平衡点$N^0$ 和 $N^K$这个在上面已经算过。第二步判断稳定性$N0$ 处不稳定$NK$ 处稳定。第三步上 Matlab 验证并出图r 0.8; K 500; f (t, N) r*N*(1 - N/K); tspan [0 30]; % 从不同初始值出发10 300 800高于K N0_list [10; 300; 800]; figure; for i 1:length(N0_list) [t, N] ode45(f, tspan, N0_list(i)); plot(t, N, LineWidth, 1.5); hold on; end yline(K, k--, LineWidth, 1, Label, K500); xlabel(时间); ylabel(种群数量 N(t)); legend(N010,N0300,N0800,Location,best); title(Logistic 种群增长无论初值如何最终都收敛到 K); grid on;注意第三个初值我故意选了 $N_0800 K$此时种群增长率是负的数量会下降直至逼近 $K$。这正是稳定状态模型的精髓——最终状态由系统结构决定和初值无关在一定范围内。这个结论写进论文配上图说服力立刻拉满。我给竞赛党一个建议你的论文里不要只写我们得到平衡点 K要写出这样一段话——无论初始种群数量低于还是高于环境容纳量系统都会渐近收敛到 $K$这意味着长期来看种群规模由资源上限决定而非初始条件。4.2 案例二SIS 传染病模型参数临界值才是得分点SIS 模型是传染病建模里最简单但最常考的一个。把人群分成易感者 S 和感染者 I 两类假设感染后能痊愈但可能再次感染。忽略出生死亡总人数 $N S I$ 守恒。$$\frac{dS}{dt} -\beta S I \gamma I, \quad \frac{dI}{dt} \beta S I - \gamma I$$因为总人数守恒其实只有一个独立变量记 $S N - I$代入得到$$\frac{dI}{dt} \beta I (N - I) - \gamma I I(\beta N - \gamma - \beta I)$$求平衡点令右边为零得到 $I^* 0$ 和 $I^* N - \gamma/\beta$。注意第二个平衡点要存在需要 $N \gamma/\beta$也就是 $\beta N / \gamma 1$。这个无量纲数就是著名的基本再生数$$R_0 \frac{\beta N}{\gamma}$$当 $R_0 1$染病人数会稳定在一个常数——疾病形成地方性流行当 $R_0 \leq 1$疫情逐步消退系统收敛到无病平衡点。关键点在于$R_0$ 不仅仅决定爆发与否还告诉你稳定状态在哪。零平衡点 $I^*0$ 在 $R_01$ 时是不稳定的——意味着疾病无法自然清零。Matlab 验证如下N 10000; beta 0.0002; % 每个感染者每天传染概率 × 接触人数 gamma 0.5; % 恢复率 1/天 R0 beta*N/gamma; sis (t, I) beta*I.*(N - I) - gamma*I; figure; I0_list [1; 100; 2000]; for i 1:length(I0_list) [t, I] ode45(sis, [0 120], I0_list(i)); plot(t, I, LineWidth, 1.5); hold on; end I_star N - gamma/beta; yline(I_star, k--, Label, sprintf(I^*%.0f, I_star)); xlabel(时间天); ylabel(感染人数 I(t)); title(sprintf(SIS 模型R_0 %.1f 1收敛到地方性流行水平, R0)); grid on;运行后你会看到三条曲线虽然起点差异巨大最终都汇到同一条水平线 $I^* N - \gamma/\beta 7500$。这就是稳定状态的吸引子属性。反过来把 $\beta$ 调小让 $R_0 1$则 $I^* 0$ 变成唯一稳定平衡点所有曲线都衰减归零。两套参数一对比一张临界相变图就出来了。事实上你还可以扫一遍 $\beta$画出稳定感染数对 $R_0$ 的曲线——在 $R_01$ 处会有一个明显的分岔。这个图几乎成了传染病建模论文的标配。我用上面代码连续扫了 20 组 $\beta$ 值得到分岔曲线后评委当场就明白了模型的全部信息。4.3 用上面两个案例提炼竞赛论文的写作模板这两个案例做完你应该能看出稳定状态模型的四段式论文写法了建立模型方程写出 $\frac{dx}{dt}f(x)$说明每个参数含义和单位求平衡点令 $f(x)0$求出所有平衡点表达式稳定性分析算 Jacobian 或 $f(x^*)$判断各平衡点的稳定性最好给出和参数的依赖关系数值验证与临界讨论ode45 模拟几条典型初值的轨迹如果系统存在分支临界值比如 $R_01$就扫参数画分岔图。这套范式可以迁移到几乎任何动态系统建模题。很多参赛队把 80% 时间花在拟合参数上模型分析只扔一个收敛结论可惜了。参数拟合固然重要但评委更想看到你对手中模型行为的理解——稳定状态分析就是展示理解深度的最佳载体。5. 实战中的高频坑和我的调试心得5.1 坑一符号法求平衡点时把参数值当符号变量混用这个错我见过太多次。符号工具箱里solve(f0, x)求出来的平衡点里带参数你如果不给参数赋值就贸然代进double()直接报错或得到一串符号表达式。正确思路是先用符号法做推导、写出公式再用subs赋参数值拿到数值平衡点。两者混在一段代码里新手很容易翻车。我的习惯是把matlabFunction用起来——符号表达式一键转函数句柄后续数值计算直接feval省去反复subs的繁琐f_func matlabFunction(f, Vars, [r, K, x]); f_star f_func(r_val, K_val, x_star);5.2 坑二ode45 遇到刚硬系统卡死或震荡有的系统时间尺度差几个数量级典型的如化学反应系统——某一成分几十毫秒就反应完另一成分要好几天才变化。此时ode45会疯狂缩小步长运行时间成倍拉长甚至报错“无法满足积分容差”。解决办法非常简单换求解器。Matlab 内置的ode15s就是干这个的一换经常立竿见影[t, y] ode15s((t,y) stiff_system(t,y), tspan, y0);我怎么判断该用ode45还是ode15s看运行时间如果ode45超过十几秒还在跑果断中断换ode15s。这不是玄学ode15s是变阶多步法专门处理特征值实部差异大的刚性问题。竞赛三天里时间就是分数别在这上面硬耗。5.3 坑三平衡点看着收敛了但数值模拟永远到不了这里面有个数学和数值的偏差要讲清楚。渐近稳定说的是 $t \to \infty$ 时才精确到达数值模拟永远只能逼近。你画出来的曲线最后一段贴着平衡线走但终值不是精确等于平衡点这很正常。不要在论文里写数值解精确达到平衡点要写在仿真时间范围内趋于平衡点。如果相图里轨迹明明应该收敛却出现小幅震荡或“卡住不动”先查绝对容差和相对容差options odeset(RelTol, 1e-6, AbsTol, 1e-8); [t, y] ode45(odefun, tspan, y0, options);特别是系统变量数量级差很大的时候比如 $S$ 上万、$I$ 只有几十AbsTol必须按变量分别设置否则小的变量会被噪声淹没曲线看起来就像没收敛。我习惯用AbsTol向量options odeset(RelTol, 1e-6, AbsTol, [1e-4 1e-6]);5.4 坑四离散模型和连续模型判据别搞混还有一个特别容易在数学建模里栽跟头的地方。上面讲的连续系统判据是 $f(x^*) 0$Jacobi 矩阵特征值实部为负但离散系统的判据完全不同。离散模型长这样$$x_{n1} g(x_n)$$它的平衡点 $x^$ 满足 $x^ g(x^)$稳定性条件是 $|g(x^)| 1$而不是 $g(x^*)0$。注意是绝对值小于 1。很多同学把连续判据套到差分方程上结果结论完全反了。最经典的例子就是离散 Logistic 映射 $x_{n1} rx_n(1 - x_n)$。$r$ 从 2.5 提到 3.2系统从稳定点走进周期 2 振荡再到周期 4、周期 8最后混沌——这全在 $|g(x^*)|$ 跨越 1 之后发生。如果你用连续系统的思维去看会完全无法理解稳定状态为什么还会长成周期震荡。遇到这类题目千万先搞清楚模型是差分方程还是微分方程。5.5 坑五相图的向量场尺度没调好quiver画出来的箭头长短取决于向量大小。在平衡点附近向量接近零所有箭头都缩成小点看不清走向远离平衡点处向量又大又长箭头互相重叠图糊成一团。解决办法是调scale参数quiver(x, y, dx, dy, 0.8); % 0.8 是箭头长度缩放因子另外建议在画相图之前把变量做无量纲化处理让状态变量的量级落在同一个区间里比如都归一化到 0~1 或 0~10否则一个变量上千、另一个变量只有个位数画出来的向量场会被大变量主导小变量的动态细节全看不出来。归一化代码如下x1_norm x1 / x1_max; x2_norm x2 / x2_max;5.6 坑六分支处的稳定性判定要小心前面提过 $f(x^*) 0$ 或特征值实部等于零的临界情况。这种中立稳定情形非常微妙比如系统刚好在分岔点上微小参数扰动就可能让稳定状态从收敛变成振荡。数值模拟在这种参数附近极其敏感稍微改下容差结果都不同。我的建议是遇到临界情况不要只凭一组参数下结论要做参数敏感度分析。扫一小段参数区间比如 $r 2.9$ 到 $3.3$把每个参数下的稳态值画出来——你立刻能看到系统在哪个点分叉。这是竞赛论文里很出彩的一步也是评委眼中这个队真的理解模型的标志。6. 我的个人经验补充从分析到论文呈现的最后一公里最后唠叨一点偏文的东西。稳定状态模型的分析做得再漂亮论文里呈现不出来也白搭。我踩过几次亏之后总结了三条呈现原则第一全文字公式和参数符号必须统一。很多队伍前面用 $N$后面写population参数一会儿叫 $\beta$ 一会儿叫传染率评委阅读体验极差。我会在正式写作前列一张符号表哪怕只是给自己看——这是稳定状态模型论文里最容易乱的地方因为涉及公式推导、代码变量、文字描述三套体系。第二图和结论必须一一对应。写了系统收敛到 K下面就要放一张展示收敛的图写了存在临界 $R_01$就要放分岔图。图是论据不是装饰每张图都得能在正文里找到对应论断。相反如果某张图你没法用一句话说明它在论证什么那就删掉。第三数值实验必须给出参数值。我看到太多人贴代码却不说 $r$、$K$、$\beta$、$\gamma$ 取了什么值。没有参数读者无法复现你的图等于白做。在附录里把每组仿真的参数列成表既专业又方便评委验证。就我个人经验来说稳定状态模型是数模竞赛性价比很高的一类工具——数学门槛不高Matlab 实现直白但对模型行为的解释力极强。不管你是第一次接触数学建模的新手还是正在准备国赛、华为杯的老手花一晚上把这套流程跑通遇到问长期趋势的题你就有了完整的分析框架。真到了赛场上你会发现最值钱的不是某个技巧而是我知道这题该往哪个方向想的笃定。