Sympy Solve( ) Gives Incorrect Answer

Question

I'm using sympy to work through some mathematical models, and I found that for some reason sympy.solve( ) gives me the wrong answers.

import sympy as sm
p, WAA, WAa, Waa = sm.symbols( 'p, WAA, WAa, Waa' )

num = p**2*WAA + p*(1-p)*WAa
denom = p**2*WAA + 2*p*(1-p)*WAa + (1-p)**2*Waa

dipMod = sm.Eq( num / denom , p )
eq = sm.solve(dipMod, p)
print eq

The results should be 0, 1, and (WAa - Waa)/(2WAa - WAA - Waa). I checked my formula by solving in MATLAB, which gave the correct answers.

Interestingly, sympy solves a simpler version of the equation correctly:

hapMod = sm.Eq(  WAA*p / (WAA*p + Waa*(1-p)), p )
print sm.solve(hapMod, p)

which gives a solution of 0 and 1. Am I doing something wrong or is there a bug in solve? I'd like to stick with Python if possible, rather than switching to MATLAB.

UPDATE:

I encountered this problem again:

p, WA, Wa = sm.symbols('p, WA, Wa')

hapMod = sm.Eq( p*WA / (p*WA + (1-p)*Wa) , p )
hapModEq = sm.solve( hapMod, p )

gives the correct answers. But substituting in for WA and Wa

hapMod2 = hapMod.subs( [(WA, 1+a*(1-p)), (Wa, 1+B*p)], simultaneous=True  )
hapMod3 = sm.simplify(hapMod2)
print hapMod3
hapMod3Eq = sm.solve(hapMod3, p)

again gives the incorrect answers. MATLAB gives the correct answers of 0, 1, a/(B + a). I found that if I take the polynomials out of the denominator and solve

test = sm.Eq( p*(a*p - a - 1)/(B*p - B*p + a*p - a*p - 1), p )
print sm.solve(test, p)

it works just fine. Is there something about having polynomial denominators that throws off sympy?

UPDATE UPDATE

After messing around with it, I've discovered that sympy is giving the correct answers, but leaving them in an oddly expanded form like

(B + 2*a)/(2*(B+a)) + sqrt(B**2)/(-2*B - 2*a). 

This simplifies to a/(B+a) which is correct, but sympy doesn't simplify it either when presenting the equilibrium or when explicitly asked to simplify this equation seprately. So it seems more like an issue with simplification than with solving. It solves correctly after all. It just seems weird that sympy would leave things like

sqrt(B**2) or sqrt( (WAA - WAa)**2) 

without simplifying them to B or (WAA - WAa).

Answer 1

If you want to simplify sqrt(x**2) to x, you need to set x to be positive, as this is not true otherwise. This is done by setting x = Symbol("x", positive=True).