Sympy求解Logistic微分方程时solveset未识别k正性问题求助

问题场景 #

在为学生编写Logistic微分方程求解的课程笔记时，针对方程 dN/dt = kN(1000−N)，初始条件 N(0)=25，并通过 N(1)=100 确定参数k，使用SymPy代码求解时，发现solveset未识别已声明为正数的符号k，导致可能出现不符合物理意义的解或无法正确筛选解。

解决方案 #

以下是三种可行的解决办法，均能确保只得到正数k的解：

方法1：使用solve替代solveset

SymPy的solve函数会更主动地利用符号的positive假设条件，直接返回符合约束的解，代码简洁适合教学展示：

import sympy as sp

# 声明正数符号k和t
k = sp.Symbol('k', positive=True) 
t = sp.Symbol('t', positive=True)
N = sp.Function('N')
N0 = 25
N1 = 100

# 定义微分方程
eqn = sp.Eq(sp.diff(N(t), t), k*N(t)*(1000 - N(t)))
# 代入初始条件N(0)=25求特解
sol_part = sp.dsolve(eqn, ics={N(0): N0})
# 代入t=1得到N(1)的表达式
sol_part_k = sol_part.subs(t, 1)
# 使用solve求解k，自动应用positive约束
k_solution = sp.solve(sol_part_k.rhs - N1, k)
print(k_solution)  # 输出: [log(13)]（化简后的结果，原表达式为log(39/3)）

方法2：在solveset中显式指定定义域

如果坚持使用solveset，可以通过domain参数明确指定k的取值范围为正实数区间，强制筛选符合条件的解：

import sympy as sp

k = sp.Symbol('k', positive=True) 
t = sp.Symbol('t', positive=True)
N = sp.Function('N')
N0 = 25
N1 = 100

eqn = sp.Eq(sp.diff(N(t), t), k*N(t)*(1000 - N(t)))
sol_part = sp.dsolve(eqn, ics={N(0): N0})
sol_part_k = sol_part.subs(t, 1)
# 显式指定定义域为(0, +∞)的正实数区间
k_solution = sp.solveset(sol_part_k.rhs - N1, k, domain=sp.Interval.open(0, sp.oo))
print(k_solution)  # 输出: {log(13)}

方法3：手动过滤所有解中的正数解

先求解所有可能的k值，再通过SymPy的is_positive属性筛选出符合正数约束的解，适合需要展示完整求解过程的教学场景：

import sympy as sp

k = sp.Symbol('k', positive=True) 
t = sp.Symbol('t', positive=True)
N = sp.Function('N')
N0 = 25
N1 = 100

eqn = sp.Eq(sp.diff(N(t), t), k*N(t)*(1000 - N(t)))
sol_part = sp.dsolve(eqn, ics={N(0): N0})
sol_part_k = sol_part.subs(t, 1)
# 先获取所有可能的解
all_solutions = sp.solveset(sol_part_k.rhs - N1, k)
# 筛选出正数解
positive_solutions = {sol for sol in all_solutions if sol.is_positive}
print(positive_solutions)  # 输出: {log(13)}

总结 #

三种方法中，方法1使用solve最适合课程笔记，代码简洁且直接利用符号约束得到正确解；方法2和3则更灵活，可根据教学需求选择展示不同的求解逻辑。

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