Sympy求解ODE平衡点方程组在Google Colab中耗时过长，求排查

你的问题核心是三元高次多项式方程组的解析求解复杂度太高，导致SymPy的solve函数陷入长时间计算。以下是几个实用的优化方案：

1. 先简化方程再求解 #

先对每个方程展开并合并同类项，减少SymPy的计算负担：

from sympy import solve, expand, simplify
from sympy.abc import c, d, e, q, m, p, u, C, D, M

# 定义原方程
eq1 = c*C*(1-C-D-M) - e*C - (d*D+ m+p*d*M)*C + (1-p)*(p*q+(1-p)*c)*(1-C-D-M)*M
eq2 = p*(p*q + (1-p)*c)*(1-C-D-M)*M-(e + u)*D+q*D*(1-C-D-M)
eq3 = (d*D+m)*C - (e+p*u)*M

# 展开并简化每个方程
simplified_eq1 = simplify(expand(eq1))
simplified_eq2 = simplify(expand(eq2))
simplified_eq3 = simplify(expand(eq3))

# 用简化后的方程求解
solutions = solve([simplified_eq1, simplified_eq2, simplified_eq3], C, D, M, dict=True)

2. 添加变量假设缩小搜索范围 #

平衡点通常对应非负实数（比如种群模型的变量C、D、M），给变量添加实数/非负假设，让SymPy聚焦合理解空间：

from sympy import symbols, solve, simplify, expand

# 定义带约束的变量
C = symbols('C', nonnegative=True, real=True)
D = symbols('D', nonnegative=True, real=True)
M = symbols('M', nonnegative=True, real=True)
c, d, e, q, m, p, u = symbols('c d e q m p u', positive=True, real=True)

# 重新定义并简化方程
eq1 = c*C*(1-C-D-M) - e*C - (d*D+ m+p*d*M)*C + (1-p)*(p*q+(1-p)*c)*(1-C-D-M)*M
eq2 = p*(p*q + (1-p)*c)*(1-C-D-M)*M-(e + u)*D+q*D*(1-C-D-M)
eq3 = (d*D+m)*C - (e+p*u)*M

simplified_eq1 = simplify(expand(eq1))
simplified_eq2 = simplify(expand(eq2))
simplified_eq3 = simplify(expand(eq3))

solutions = solve([simplified_eq1, simplified_eq2, simplified_eq3], C, D, M, dict=True)

3. 分步求解降低维度 #

从简单方程解出一个变量，代入其余方程减少变量数量：

from sympy import solve, expand, simplify
from sympy.abc import c, d, e, q, m, p, u, C, D, M

eq1 = c*C*(1-C-D-M) - e*C - (d*D+ m+p*d*M)*C + (1-p)*(p*q+(1-p)*c)*(1-C-D-M)*M
eq2 = p*(p*q + (1-p)*c)*(1-C-D-M)*M-(e + u)*D+q*D*(1-C-D-M)
eq3 = (d*D+m)*C - (e+p*u)*M

# 从第三个方程解出M的表达式
M_expr = solve(eq3, M)[0]

# 代入前两个方程，转化为二元方程组
eq1_sub = simplify(expand(eq1.subs(M, M_expr)))
eq2_sub = simplify(expand(eq2.subs(M, M_expr)))

# 求解二元方程组
solutions_CD = solve([eq1_sub, eq2_sub], C, D, dict=True)

# 补全M的解
solutions = []
for sol in solutions_CD:
    sol_full = sol.copy()
    sol_full[M] = M_expr.subs(sol)
    solutions.append(sol_full)

4. 转数值解（若无需解析解） #

如果只需要具体数值平衡点，给参数赋值后用nsolve快速求解：

from sympy import nsolve, symbols

# 给参数赋值（根据你的模型实际取值调整）
c, d, e, q, m, p, u = 1, 0.5, 0.2, 0.3, 0.1, 0.4, 0.1
C, D, M = symbols('C D M')

# 定义方程
eq1 = c*C*(1-C-D-M) - e*C - (d*D+ m+p*d*M)*C + (1-p)*(p*q+(1-p)*c)*(1-C-D-M)*M
eq2 = p*(p*q + (1-p)*c)*(1-C-D-M)*M-(e + u)*D+q*D*(1-C-D-M)
eq3 = (d*D+m)*C - (e+p*u)*M

# 传入初始猜测值求解
sol = nsolve([eq1, eq2, eq3], (C, D, M), (0.1, 0.1, 0.1))
print(f"C={sol[0]}, D={sol[1]}, M={sol[2]}")

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