[Python] Issue with sympy solver, giving different results based off order

So, I have noticed when I switch the order of the symbols to be solved in the sympy.solve I get slightly different results. My code:

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) #this has same results as the next line
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)


Thus, sol and sol_2 produce the exact same results for the first two solutions for the system of equations,which are (p_f,p_b,p_t)=(0,0,0) and (-(-2*c - 2*d + r)/(2*gamma), 0, (-d + r)/gamma). However, in the third solution all but p_b values are the same. For sol the solution for p_b ends up being an A4 page in length in latex (thus I won't put it into this question), while sol_2 the solution is 'merely' (-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)). The values for the other parameters in the third solution are p_t=0 and p_f=(2**n_b*(c*n_b + d - r) + r)/(2**n_b*gamma) (for p_f they appear slightly different but a bit of arithmetic shows they are equivalent).

What is going on? I have tried simplifying the large polynomial in sympy to no avail. A simple inequality test indicates the two are not equal. I have also put the symbols in an "denested" list of symbols as indicated in the docs (i.e. solve([dbdt,dtdt,dfdt],p_t,p_b,p_f,rational=True)), and I get the same result depending on the order of the symbols as when I place them in []. And I have tried to solve it by hand (inserting the values for p_f and p_t), and I seem to make an arithmetic error somewhere due to my answer being different to both sol and sol_2 (and if the solution is truly as large as the solution in sol, it would be impossible for me to do by hand).

Continue reading...