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).
This is another interesting case where aggressive symbol management leads to a nice solution for this set of equations. Using something to group symbols not of interest (like the condense function here) your equations can be written as ce, R where ce are your equations with symbols you are not solving for recursively lumped together into replacements, R:
... your code then
from sympy import var,solve,nsimplify
eqs = [nsimplify(i,rational=True) for i in [dbdt,dtdt,dfdt]]
C0, C1, C10, C11, C12, C13, C2, C4, C5, C6, C7=var('C0, C1, C10, C11, C12, C13, C2, C4, C5, C6, C7')
ce = [C2*(C6*p_b + C7*p_b + p_b*p_f), C2*p_t*(C4 + p_f), C13*(C10*p_f + C11*p_f + C12*p_f*(p_b + p_t) + p_b)]
R = [(C13, C5), (C12, -2**n_b*gamma/r), (C11, -2**n_b*d/r), (C10, 2**n_b), (C7, C1*r/gamma), (C6, -C0/gamma), (C5, r/2**n_b), (C4, -c/gamma - d/gamma + r/gamma/2), (C2, gamma), (C1, 1 - 1/2**n_b), (C0, c*n_b + d)]
In this form any permutations of symbols and equations gives the same result (and it takes only about 3 seconds to do the 36 permutations):
from itertools import permutations
v=p_b, p_f, p_t
p=[solve(e,vi) for e in permutations(ce) for vi in permutations(v)]
>>> {frozenset(k) for pi in p for k in pi}
{frozenset({0, -(C10 + C11)*(C6 + C7)/(C12*C6 + C12*C7 - 1), -C6 - C7}), frozenset({0, -C4, -(C10 + C11)/C12}), frozenset({0})}
Your solution set is the following (which satisfies all equations as is shown).
>>> [[i.simplify() for i in Tuple(*si).subs(R)] for si in p[0]]
[[0, 0, 0], [0, (c + d - r/2)/gamma, (-d + r)/gamma], [-(d - r)*(2**n_b*(c*n_b + d) - r*(2**n_b - 1))/(2**n_b*gamma*(c*n_b + d - r)), c*n_b/gamma + d/gamma - r/gamma + r/(2**n_b*gamma), 0]]
>>> {ei.subs(dict(zip(v,si))).cancel().simplify() for ei in eqs for si in _}
{0}
As a side note, it is interesting to see that groebner can handle this system of equations though the resulting complexity varies with the permutation order of the variables:
from sympy import groebner, count_ops
>>> [count_ops(groebner(ce,p)) for p in permutations(v)]
[63, 59, 63, 42, 32, 32]
>>> [count_ops(groebner(eqs,p)) for p in permutations(v)]
[946, 344, 946, 213, 156, 156]
If you search for Groebner here in [sympy] you will find lots of applications shown by Oscar Benjamin.