SageMath求解9阶非线性微分方程组平衡态报错求助

核心问题排查 #

你的代码存在两个关键问题导致报错：

1. 缺失dS_hdt的定义：在solve调用中包含了dS_hdt==0，但代码里从未定义这个导数表达式，未定义的符号表达式传入求解器会引发接口异常。
2. 非线性方程组复杂度过高：方程组包含大量变量乘积项（如alpha*I_hr*(I_ms +I_msr)），直接用Sage默认的Maxima求解器处理9维非线性符号方程组，容易因计算量过大导致进程崩溃（报错中的Aborted就是进程被终止的表现）。

修复步骤 #

1. 补全缺失的方程 #

根据传染病模型的守恒关系（N_h = S_h + I_hr + I_hs + I_hsr + R_h），补全dS_hdt的表达式（需根据你的模型实际转移逻辑调整）：

dS_hdt = mu_h*N_h - beta*S_h*(I_mr + I_ms + I_msr) - mu_h*S_h

2. 分步简化求解 #

避免一次性求解全部9个方程，先解线性或低复杂度的方程，代入后减少变量数量：

# 先解dR_hdt=0，得到R_h的表达式
sol_R = solve(dR_hdt == 0, R_h)[0]
R_h_expr = sol_R.rhs()

# 解dS_mdt=0，得到S_m的表达式
sol_Sm = solve(dS_mdt == 0, S_m)[0]
S_m_expr = sol_Sm.rhs()

# 将R_h和S_m的表达式代入其他方程
subs_dict = {R_h: R_h_expr, S_m: S_m_expr}
eqs_subs = [
    dS_hdt.subs(subs_dict) == 0,
    dI_hrdt.subs(subs_dict) == 0,
    dI_hsdt.subs(subs_dict) == 0,
    dI_hsrdt.subs(subs_dict) == 0,
    dI_mrdt.subs(subs_dict) == 0,
    dI_msdt.subs(subs_dict) == 0,
    dI_msrdt.subs(subs_dict) == 0
]
vars_subs = [S_h, I_hr, I_hs, I_hsr, I_mr, I_ms, I_msr]

# 求解简化后的方程组
soln = solve(eqs_subs, vars_subs)
show(soln)

3. 切换求解器 #

尝试使用SymPy作为求解器，替代默认的Maxima，可能更稳定处理复杂非线性方程组：

soln = solve(eqs_subs, vars_subs, algorithm='sympy')
show(soln)

4. 数值求解备选（若符号解不可行） #

如果符号求解仍失败，可先给参数赋值，求数值平衡态：

# 给参数赋值示例
params = {
    mu_h: 0.01, beta: 0.5, alpha: 0.1, gamma_r: 0.2, gamma_s: 0.3,
    omega: 0.05, mu_m: 0.1, epsilon: 0.4, theta: 0.02, N_h: 1000, N_m: 500
}
eqs_num = [eq.subs(params) for eq in eqs_subs]
# 使用scipy数值求解
from scipy.optimize import fsolve
import numpy as np

def system(x):
    S_h, I_hr, I_hs, I_hsr, I_mr, I_ms, I_msr = x
    return [eq.subs(dict(zip(vars_subs, x))).n() for eq in eqs_num]

# 初始猜测值
guess = [900, 50, 30, 10, 20, 15, 5]
sol_num = fsolve(system, guess)
print("数值平衡态：", sol_num)

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