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).

Indulge

I'm not familiar with sympy, but it appears to me that you need to somehow declare that 2WAa - WAA - Waa != 0 in order for a computational algebra system to give that answer. – Naamana

Although operator precedence is well-defined, consider adding parentheses to p**2*WAA. That would make code more readable. – Keverne

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).

Christoper

I did this, and the equilibria became simpler: (-B/(-2*B - 2*a) + (B + 2*a)/(2*(B + a))). It turns out that solver( ) wasn't simplifying by default, so using sm.solve( , simplify=True) worked! Also, I had to set sm.symbols('p', negative=False) because p is allowed to be zero. Thanks for all of the help this worked perfectly! – Indulge

What you really want is nonnegative=True, which is the same as negative=False and real=True. – Christoper

I like that sympy is sensitive to these details – Salsbury