Solve simplest algebraic equations with dummy parameter

[nickme68]

Very strange behaviour of solve command in SymPy, when I try to solve quadratic equation with parameter a (which is dummy here, because solution doesn't depends on it):

from sympy import *
x, a = symbols("x a")
eq = Eq(0, (4 - 4*x + x**2)/(4*a**2))
print(solve(eq, x))
print(solve(simplify(eq), x))

Output:

[2 - sqrt(a**2 - 1)/a, 2 + sqrt(a**2 - 1)/a]
[2]

Just solve gives two (!) solutions, which depend on a. After symplifying it gives only solution x=2, which is correct. What happens? Command solveset works correctly, but I am interested in using solve command. My SymPy version is 1.8.

[nickme68]

Simpler example:

from sympy import *
x, a = symbols("x a")
eq = Eq(0, x / (2 * a))
print(solve(eq, x))
print(solve(simplify(eq), x))

Output:

[(a - 1)/a]
[0]

But if we swap sides of the equation eq = Eq(x / (2 * a), 0), then solve solves it correctly!

[oscargus]

Not sure what happens, but if you do

eq = (4 - 4*x + x**2)/(4*a**2)
print(solve(eq, x))

You get the expected result.

I guess that there should be something like

if isinstance(expr, Eq):
   expr = expr.rhs - expr.lhs

in solve, but haven't checked (and it looks here like not).

Also

eq = Eq((4 - 4*x + x**2)/(4*a**2), 0)
print(solve(eq, x))

works, so looks like a bug.

[nickme68]

Thank you, Oscar! I obtained this quadratic equation (with much more parameters) after substitution in the solution of ordinary differential equation, so it's not very convenient for me to do these manual reorderings and simplifications, I just want to solve it as is. The problem is somewhere in the order of parts of the equation: x/(2*a)=0 solved correctly, but 0=x/(2*a) - incorrectly. Another hint: 0=x/(1*a) gives correct answer, but using any other constant instead of 1 and -1 in denominator leads to wrong answer.

[ThePauliPrinciple]

>>> from sympy import *
>>> x, a = symbols("x a")
>>> eq = Eq(0, (4 - 4*x + x**2)/(4*a**2))
>>> print(solveset(eq, x))
>>> print(solveset(simplify(eq), x))
FiniteSet(2)
FiniteSet(2)

solveset seems to give the right answer in the given example.

[oscarbenjamin]

The solve function has lots of weird code paths depending on exactly how the input is provided e.g. you can get different results for solve(eq, x) vs solve([eq], x) etc or different results for Eq vs Expr:

In [2]: solve(eq, x)
Out[2]: 
⎡       ________         ________⎤
⎢      ╱  2             ╱  2     ⎥
⎢    ╲╱  a  - 1       ╲╱  a  - 1 ⎥
⎢2 - ───────────, 2 + ───────────⎥
⎣         a                a     ⎦

In [3]: solve([eq], x)
Out[3]: [(2,)]

The main branch is here:

sympy/sympy/solvers/solvers.py

Lines 1106 to 1109 in d822fcb

	if bare_f:
	solution = _solve(f[0], *symbols, **flags)
	else:
	solution = _solve_system(f, symbols, **flags)

Try working backwards though and see if you can figure out when bare_f is true or false or what the point of it is.

[smichr]

It's a problem with as_numer_denom()