pythonwolfram-mathematicasympy
Can sympy solve these polynomial square root equations?
I've attempted solving an equation in SymPy in 2 different forms with an error and an indefinite loop. This is on Python 3.8.8 and SymPy 1.7.1 with the same results on two devices. In contrast, Mathematica solves both in under a second.
Equations
Here are equivalent forms of the same equation
from sympy import *
x,xh,xhh,p,q,r = symbols('x xh xhh p q r')
eq1 = ((2*q**2-2*r**2)/(-q - r + xh) - 2*q - 2*r + x + xhh
+ sqrt(4*p*x + 4*p*xh + 4*q**2 + x**2 + 2*x*xh + xh**2)
+ sqrt(4*p*xh + 4*q**2 + xh**2 + xhh**2 + 2*xhh*(2*p + xh))
)
eq2 = (4*q**2 + 4*q*r - q*x - 2*q*xh - r*x - 2*r*xh + x*xh + xhh*(-q - r + xh)
+ (-q - r + xh)*sqrt(4*p*xh + 4*q**2 + xh**2 + xhh**2 + 2*xhh*(2*p + xh))
+ (-q - r + xh)*sqrt(4*p*x + 4*p*xh + 4*q**2 + x**2 + 2*x*xh + xh**2)
)
A,C,D,F,G,H = symbols('A C D F G H')
eq3 = A + xhh + C/D + sqrt(F) + sqrt(G + H*xhh + xhh**2)
Solving Attempts
Solving the first form solve(eq1,xhh) returns the error
Traceback (most recent call last):
File "<ipython-input-24-9923c6f02abd>", line 1, in <module>
solve(eq1,xhh)
File "/Users/jacobrichardson/anaconda/envs/env_sympy/lib/python3.8/site-packages/sympy/solvers/solvers.py", line 1097, in solve
solution = _solve(f[0], *symbols, **flags)
File "/Users/jacobrichardson/anaconda/envs/env_sympy/lib/python3.8/site-packages/sympy/solvers/solvers.py", line 1460, in _solve
gen = f_num.match(D.xreplace({d: w}))[w]
TypeError: 'NoneType' object is not subscriptable
Solving the second solve(eq2,xhh) runs for over an hour without completing on both my devices!
Indeed with my manual assistance, the 3rd form solves solve(eq3,xhh) giving
[(-A**2*D**2 - 2*A*C*D - 2*A*D**2*sqrt(F) - C**2 - 2*C*D*sqrt(F) - D**2*F + D**2*G)/(D*(2*A*D + 2*C + 2*D*sqrt(F) - D*H))]
Is there any way to get SymPy to solve equations like eq1 and eq2 as easily and conveniently as Mathematica Solve[eq1,xhh] and Solve[eq2,xhh]?
Equialent Equations
print(simplify(eq1 * (-q - r + xh) - eq2)) # "0" equations 1 and 2 equivalent up to one singularity
print(simplify(eq1 -
eq3.subs([(A,-2*q - 2*r + x),(C,2*q**2 - 2*r**2),(D,-q - r + xh),
(F,4*p*x + 4*p*xh + 4*q**2 + x**2 + 2*x*xh + xh**2),
(G,4*p*xh + 4*q**2 + xh**2),(H,2*(2*p + xh))]))) # "0" equations 1 and 3 equivalent
Solution
I tried solve(eq1, xhh) and I see the error you showed when testing with sympy 1.7.1 (the latest release). I think that's a bug that has already been fixed though because I don't see that with the latest "master" version of sympy. In fact here is the pull request that fixed it:
https://github.com/sympy/sympy/pull/20842
Instead with master I find that it is just slow as you report for eq2. Interrupting it I can see that it is slow in the checking code which can be disabled. The check=False it takes 4 seconds for eq1 and 5 seconds for eq2:
In [2]: %time solve(eq1, xhh, check=False)
CPU times: user 4.21 s, sys: 45.3 ms, total: 4.25 s
Wall time: 4.33 s
Out[2]:
⎡ ______________________
⎢ 2 2 2 4 3 3 3 3 ╱ 2
⎢- 2⋅p⋅q ⋅x - 4⋅p⋅q⋅r⋅x + 4⋅p⋅q⋅x⋅xh - 2⋅p⋅r ⋅x + 4⋅p⋅r⋅x⋅xh - 2⋅p⋅x⋅xh - 8⋅q - 16⋅q ⋅r + 4⋅q ⋅x + 8⋅q ⋅xh + 4⋅q ⋅╲╱ 4⋅p⋅x + 4⋅p⋅xh + 4⋅q
⎢────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
⎢
⎢
⎣
_____________________ ___________________________________________
2 2 2 2 2 2 2 ╱ 2 2 2 2 2 2 2
+ x + 2⋅x⋅xh + xh - 8⋅q ⋅r + 8⋅q ⋅r⋅x + 16⋅q ⋅r⋅xh + 8⋅q ⋅r⋅╲╱ 4⋅p⋅x + 4⋅p⋅xh + 4⋅q + x + 2⋅x⋅xh + xh - q ⋅x - 7⋅q ⋅x⋅xh - q ⋅x⋅╲╱
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
___________________________________________ ___________________________________________
╱ 2 2 2 2 2 2 ╱ 2 2 2 2 2
4⋅p⋅x + 4⋅p⋅xh + 4⋅q + x + 2⋅x⋅xh + xh - 2⋅q ⋅xh - 6⋅q ⋅xh⋅╲╱ 4⋅p⋅x + 4⋅p⋅xh + 4⋅q + x + 2⋅x⋅xh + xh + 4⋅q⋅r ⋅x + 8⋅q⋅r ⋅xh + 4⋅q
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
____________
2 2 2 3 2 2 2 2 ╱
- 2⋅p⋅q - 4⋅p⋅q⋅r + 4⋅p⋅q⋅xh - 2⋅p⋅r + 4⋅p⋅r⋅xh - 2⋅p⋅xh - 4⋅q - 8⋅q ⋅r + q ⋅x + 5⋅q ⋅xh + q ⋅╲╱ 4⋅p⋅x + 4⋅p
___________________________________________ ___________________________________________
2 ╱ 2 2 2 2 ╱ 2 2 2
⋅r ⋅╲╱ 4⋅p⋅x + 4⋅p⋅xh + 4⋅q + x + 2⋅x⋅xh + xh - 2⋅q⋅r⋅x - 10⋅q⋅r⋅x⋅xh - 2⋅q⋅r⋅x⋅╲╱ 4⋅p⋅x + 4⋅p⋅xh + 4⋅q + x + 2⋅x⋅xh + xh - 4⋅q⋅r⋅
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
_______________________________ ___________________________________________ _
2 2 2 2 ╱ 2 2 2 ╱
⋅xh + 4⋅q + x + 2⋅x⋅xh + xh - 4⋅q⋅r + 2⋅q⋅r⋅x + 6⋅q⋅r⋅xh + 2⋅q⋅r⋅╲╱ 4⋅p⋅x + 4⋅p⋅xh + 4⋅q + x + 2⋅x⋅xh + xh - 2⋅q⋅x⋅xh - 2⋅q⋅xh⋅╲╱
___________________________________________ _________________________________________
2 ╱ 2 2 2 2 2 ╱ 2 2
xh - 8⋅q⋅r⋅xh⋅╲╱ 4⋅p⋅x + 4⋅p⋅xh + 4⋅q + x + 2⋅x⋅xh + xh + 2⋅q⋅x ⋅xh + 4⋅q⋅x⋅xh + 2⋅q⋅x⋅xh⋅╲╱ 4⋅p⋅x + 4⋅p⋅xh + 4⋅q + x + 2⋅x⋅xh + xh
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
__________________________________________ ___________________________________________ ________
2 2 2 2 2 2 ╱ 2 2 2 ╱
4⋅p⋅x + 4⋅p⋅xh + 4⋅q + x + 2⋅x⋅xh + xh + r ⋅x + r ⋅xh + r ⋅╲╱ 4⋅p⋅x + 4⋅p⋅xh + 4⋅q + x + 2⋅x⋅xh + xh - 2⋅r⋅x⋅xh - 2⋅r⋅xh⋅╲╱ 4⋅p⋅x +
__ ___________________________________________ ___________________________________________
2 2 ╱ 2 2 2 2 2 2 2 ╱ 2 2 2 2
+ 2⋅q⋅xh ⋅╲╱ 4⋅p⋅x + 4⋅p⋅xh + 4⋅q + x + 2⋅x⋅xh + xh - r ⋅x - 3⋅r ⋅x⋅xh - r ⋅x⋅╲╱ 4⋅p⋅x + 4⋅p⋅xh + 4⋅q + x + 2⋅x⋅xh + xh - 2⋅r ⋅
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
___________________________________ ___________________________________________
2 2 2 2 3 2 ╱ 2 2 2
4⋅p⋅xh + 4⋅q + x + 2⋅x⋅xh + xh + x⋅xh - xh + xh ⋅╲╱ 4⋅p⋅x + 4⋅p⋅xh + 4⋅q + x + 2⋅x⋅xh + xh
___________________________________________ __________________________________________
2 2 ╱ 2 2 2 2 2 ╱ 2 2 2
xh - 2⋅r ⋅xh⋅╲╱ 4⋅p⋅x + 4⋅p⋅xh + 4⋅q + x + 2⋅x⋅xh + xh + 2⋅r⋅x ⋅xh + 4⋅r⋅x⋅xh + 2⋅r⋅x⋅xh⋅╲╱ 4⋅p⋅x + 4⋅p⋅xh + 4⋅q + x + 2⋅x⋅xh + xh
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
_ ___________________________________________ ___________________________________________⎤
2 ╱ 2 2 2 2 2 3 2 ╱ 2 2 2 ⎥
+ 2⋅r⋅xh ⋅╲╱ 4⋅p⋅x + 4⋅p⋅xh + 4⋅q + x + 2⋅x⋅xh + xh - x ⋅xh - x⋅xh - x⋅xh ⋅╲╱ 4⋅p⋅x + 4⋅p⋅xh + 4⋅q + x + 2⋅x⋅xh + xh ⎥
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────⎥
⎥
⎥
⎦