Sympy solve looking for complex solutions, even if real=True?

Question

I mean to solve a symbolic equation via sympy.solve. It was taking surprisingly long, so I tried a few variations, until I found a posible cause. It looks for complex solutions, while I am only interested in real solutions.

Code below...

import sympy as sp
import numpy as np
x, y, z = sp.symbols('x, y, z', real=True, positive=True)

# Option X1: works
expr3 = 0.6*x**(9/8)*y**(7/5)*z**2.2
expr4 = 0.9*x**(1/2)*y**(4/5)*z**1.2
s3 = sp.solve(expr3 - expr4, x)
print('Exponents of x: 9/8, 1/2. Solution =', s3)

# Option X3: works
expr3 = 0.6*x**(11/9)*y**(7/5)
expr4 = 0.9*x**(1/2)*y**(4/5)
s3 = sp.solve(expr3 - expr4, x, check=True, minimal=True, rational=True, force=True)
print('Exponents of x: 11/9, 1/2. Solution =', s3)

# Option X2: takes "forever"
expr3 = 0.6*x**(11/9)*y**(7/5)
expr4 = 0.9*x**(1/2)*y**(4/5)
s3 = sp.solve(expr3 - expr4, x)
print('Exponents of x: 11/9, 1/2. Solution =', s3)

... produces this output

Exponents of x: 9/8, 1/2. Solution = [1.91313675093869/(y**(24/25)*z**(8/5)), 0.351080874270527*(-y**(-0.12)*z**(-0.2) - 0.726542528005361*I*y**(-0.12)*z**(-0.2))**8, 0.351080874270527*(-y**(-0.12)*z**(-0.2) + 0.726542528005361*I*y**(-0.12)*z**(-0.2))**8, 1.28055023108529*(0.324919696232906*y**(-0.12)*z**(-0.2) - I*y**(-0.12)*z**(-0.2))**8, 1.28055023108529*(0.324919696232906*y**(-0.12)*z**(-0.2) + I*y**(-0.12)*z**(-0.2))**8]
Exponents of x: 11/9, 1/2. Solution = [3*2**(8/13)*3**(5/13)/(4*y**(54/65)), (-2**(12/13)*3**(1/13)*cos(pi/13)/(2*y**(3/65)) - 2**(12/13)*3**(1/13)*I*sin(pi/13)/(2*y**(3/65)))**18, (-2**(12/13)*3**(1/13)*cos(pi/13)/(2*y**(3/65)) + 2**(12/13)*3**(1/13)*I*sin(pi/13)/(2*y**(3/65)))**18, (2**(12/13)*3**(1/13)*cos(2*pi/13)/(2*y**(3/65)) - 2**(12/13)*3**(1/13)*I*sin(2*pi/13)/(2*y**(3/65)))**18, (2**(12/13)*3**(1/13)*cos(2*pi/13)/(2*y**(3/65)) + 2**(12/13)*3**(1/13)*I*sin(2*pi/13)/(2*y**(3/65)))**18, (-2**(12/13)*3**(1/13)*cos(3*pi/13)/(2*y**(3/65)) - 2**(12/13)*3**(1/13)*I*sin(3*pi/13)/(2*y**(3/65)))**18, (-2**(12/13)*3**(1/13)*cos(3*pi/13)/(2*y**(3/65)) + 2**(12/13)*3**(1/13)*I*sin(3*pi/13)/(2*y**(3/65)))**18, (2**(12/13)*3**(1/13)*cos(4*pi/13)/(2*y**(3/65)) - 2**(12/13)*3**(1/13)*I*sin(4*pi/13)/(2*y**(3/65)))**18, (2**(12/13)*3**(1/13)*cos(4*pi/13)/(2*y**(3/65)) + 2**(12/13)*3**(1/13)*I*sin(4*pi/13)/(2*y**(3/65)))**18, (-2**(12/13)*3**(1/13)*cos(5*pi/13)/(2*y**(3/65)) - 2**(12/13)*3**(1/13)*I*sin(5*pi/13)/(2*y**(3/65)))**18, (-2**(12/13)*3**(1/13)*cos(5*pi/13)/(2*y**(3/65)) + 2**(12/13)*3**(1/13)*I*sin(5*pi/13)/(2*y**(3/65)))**18, (2**(12/13)*3**(1/13)*cos(6*pi/13)/(2*y**(3/65)) - 2**(12/13)*3**(1/13)*I*sin(6*pi/13)/(2*y**(3/65)))**18, (2**(12/13)*3**(1/13)*cos(6*pi/13)/(2*y**(3/65)) + 2**(12/13)*3**(1/13)*I*sin(6*pi/13)/(2*y**(3/65)))**18]

from options X1, and X3, and then it is stuck at option X2. I have to manually interrupt calculation.
X3 is one of many attempts at tinkering with options of solve to get a solution for X2. From its result, I would understand why X2 takes so long.

The first element in the solution list is what I am looking for, possibly with later simplification. My questions are as follows:

1. Q1: Why is solve looking for and returning complex solutions, if I set real=True, positive=True?
2. Q2: Even if not understanding Q1, is there any way to tell solve not to look for complex solutions? (Note: I also want positive=True for x, not x=0)
3. Q3: I would be ok with using options for solve that lead to calculations in reasonable time, even if I get complex solutions.
   But how can I then automatically pick the real solution?
   I would pick element 0, but I am not sure it will always be the correct one.

Related:

1. Ignore imaginary roots in sympy
2. Sympy very slow at solving equations using solve

Answer 1

By default, solve converts floating point numbers in your expression into rational numbers, which then are used to solve some polynomial equation. When the computation takes a long time, we can play with keyword arguments. For the third case:

expr3 = 0.6*x**(11/9)*y**(7/5)
expr4 = 0.9*x**(1/2)*y**(4/5)
s3 = sp.solve(expr3 - expr4, x, check=False)
print('Exponents of x: 11/9, 1/2. Solution =', s3)

The computation is immediate. Another option would be to ask solve to avoid converting floating point numbers to rational:

s3 = sp.solve(expr3 - expr4, x, rational=False)
print('Exponents of x: 11/9, 1/2. Solution =', s3)

Again, the computation is really fast.

Now to the other question. One way to work around the issue might be by telling solve to only look for specific regions. The following commands takes a few seconds to compute:

sp.solve([expr3 - expr4, x >= 0, y >= 0], x)
# [(1.75314834631487/y**(54/65),)]

However, note that the trivial solution x=0 is omitted.

Let's try on your first case:

expr3 = 0.6*x**(9/8)*y**(7/5)*z**2.2
expr4 = 0.9*x**(1/2)*y**(4/5)*z**1.2
sp.solve([expr3 - expr4, x >= 0, y >= 0, z >= 0], x)
# [(1.91313675093869/(y**(24/25)*z**(8/5)),)]

Answer 2

With rationals the last one runs faster:

In [45]: # Option X1: works
    ...: expr3 = Rational(3,5)*x**(Rational(9,8))*y**(Rational(7,5))
    ...: expr4 = Rational(9,10)*x**(Rational(1,2))*y**(Rational(4,5))
    ...: s3 = sp.solve(expr3 - expr4, x)

In [46]: print(s3)
[0, 3*2**(2/5)*3**(3/5)/(4*y**(24/25)), (-2**(4/5)*3**(1/5)/(8*y**(3/25)) + 2**(4/5)*3**(1/5)*sqrt(5)/(8*y**(3/25)) - 2**(4/5)*3**(1/5)*I*sqrt(sqrt(5)/8 + 5/8)/(2*y**(3/25)))**8, (-2**(4/5)*3**(1/5)/(8*y**(3/25)) + 2**(4/5)*3**(1/5)*sqrt(5)/(8*y**(3/25)) + 2**(4/5)*3**(1/5)*I*sqrt(sqrt(5)/8 + 5/8)/(2*y**(3/25)))**8, (-2**(4/5)*3**(1/5)*sqrt(5)/(8*y**(3/25)) - 2**(4/5)*3**(1/5)/(8*y**(3/25)) - 2**(4/5)*3**(1/5)*I*sqrt(5/8 - sqrt(5)/8)/(2*y**(3/25)))**8, (-2**(4/5)*3**(1/5)*sqrt(5)/(8*y**(3/25)) - 2**(4/5)*3**(1/5)/(8*y**(3/25)) + 2**(4/5)*3**(1/5)*I*sqrt(5/8 - sqrt(5)/8)/(2*y**(3/25)))**8]

Answer 3

I don't know the answer for Q1 and Q2. But I managed to get code that gives me the single real positive solution that one would get by hand calculation for expr3 and expr4 of the same form as in the OP (expr3 = number * product(xi**bi, i, 1, n), with n variables and exponents bi typically rational). So far code worked for all use cases, including those that were not properly solved with any of what was posted. And it works fast.

# Option X6: works
expr3 = 0.6*x**(4/9)*y**(1/3)
expr4 = 0.9*x**(2/5)*y**(2/5)
s3_allsoln = sp.solve([expr3 - expr4, x>0, y>0], x, check=False, minimal=True, rational=True, force=True, dict=True)
s3_real = [s3[x] for s3 in s3_allsoln if s3[x].is_real]
s3_realpos = [s3 for s3 in s3_real if s3.is_positive]
if(len(s3_realpos) != 1):
    print('Warning, found', len(s3_realpos), 'real positive solutions')
s3 = s3_realpos[0]
print('Exponents of x: 4/9, 2/5. Solution =', s3)

Output:

Exponents of x: 4/9, 2/5. Solution = 31381059609*sqrt(6)*y**(3/2)/8388608

I still need to simplify the numerical factor, into a decimal number. But that is likely much easier.