如何获取SIR模型ODE系统的所有稳态解？SymPy求解疑问

首先明确：SIR模型的稳态要求所有导数均为0，即dSdt=0、dIdt=0、dRdt=0。从dRdt=g*I=0推导，由于恢复率g通常为正数，必然得出I=0。将I=0代入另外两个方程，dSdt=-b*S*0=0、dIdt=b*S*0 -g*0=0均自动成立，因此所有满足I=0的(S, I, R)都是稳态。

SymPy返回的[(S, 0, R)]是稳态的通解，表示S和R可取任意值（若模型隐含总人数守恒约束S+I+R=N，则R=N-S）。你提到的[(g/b, 0, R)]并非独立的第二个稳态解，只是通解中S=g/b的一个特例。

若需显式获取包含S=g/b的解 #

如果你想针对性得到S=g/b的情况（本质仍属于I=0的稳态），可以通过添加额外条件求解：

import sympy as sp
S, I, R, b, g = sp.symbols('S I R b g', positive=True)
dSdt = -b * S * I
dIdt = b * S * I - g * I
dRdt = g * I
# 加入条件b*S - g = 0，锁定S的取值
steady_states = sp.solve([dSdt, dIdt, dRdt, sp.Eq(b*S, g)], [S, I, R])
steady_states

运行结果为[(g/b, 0, R)]，符合你的预期。

关于地方病平衡点的说明 #

如果你误将地方病平衡点（I≠0时dIdt=0的点）当作稳态，需要注意：此时dRdt=g*I≠0，不满足稳态的所有导数为0的要求。若仅需寻找dIdt=0的点，可忽略dRdt=0的约束求解：

import sympy as sp
S, I, R, b, g = sp.symbols('S I R b g', positive=True)
dSdt = -b * S * I
dIdt = b * S * I - g * I
equilibria = sp.solve([dSdt, dIdt], [S, I])
equilibria

结果会返回[(S, 0), (g/b, 0)]，但需注意这两个点中只有I=0的情况是真正的稳态。

内容的提问来源于stack exchange，提问作者YoshiroVilchez