SymPy solve函数求解结果随符号顺序变化的问题咨询

我在使用sympy.solve求解方程组时，发现切换待求解符号的顺序会得到不同的结果。代码如下：

from sympy.solvers import solve
import sympy as sy

r=sy.Symbol('r')
p_f=sy.Symbol("p_f")
p_b=sy.Symbol("p_b")
n_b=sy.Symbol("n_b")
n_t=sy.Symbol("n_t")
gamma=sy.Symbol("gamma")
p_t=sy.Symbol("p_t")
d=sy.Symbol("d")
c=sy.Symbol("c")
n_t = 1

dtdt= r*(1-2**-n_t)*p_t+p_f*p_t*gamma-p_t*(d+c*n_t)
dbdt= r*(1-2**-n_b)*p_b+p_f*p_b*gamma-p_b*(d+c*n_b)
dfdt = r*p_f+p_b*r*2**-n_b-gamma*(p_t+p_b)*p_f-d*p_f
dtdt=sy.simplify(dtdt)

# sol=solve([dtdt, dbdt,dfdt],[p_b,p_t,p_f],rational=True) # 与下一行结果一致
sol=solve([dbdt,dtdt,dfdt],[p_t,p_b,p_f],rational=True) 
sol_2=solve([dfdt,dbdt,dtdt],[p_f,p_b,p_t],rational=True)

现象说明 #

• sol和sol_2给出的前两个解完全一致：
  1. (p_f,p_b,p_t)=(0,0,0)
  2. (-(-2*c - 2*d + r)/(2*gamma), 0, (-d + r)/gamma)
• 第三个解中，除p_b外的其他值一致：
  • sol中的p_b表达式异常冗长，LaTeX渲染后能占满一页A4纸；
  • sol_2中的p_b表达式为：(-2**n_b*c*d*n_b + 2**n_b*c*n_b*r - 2**n_b*d**2 - 2**n_b*r**2 + 2**(n_b + 1)*d*r - d*r + r**2)/(2**n_b*gamma*(c*n_b + d - r))
• 第三个解中p_t=0，p_f的表达式看似不同，但经运算可证明等价；
• 尝试简化sol中冗长的p_b表达式无效，简单不等式测试显示两个p_b表达式不相等；
• 按文档要求使用非嵌套符号列表调用solve，结果仍随符号顺序变化；
• 手动代入p_f和p_t求解时出现运算错误，若sol中的解真如此冗长，手动计算不现实。

请问这是怎么回事？

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