search

sympy / sympy

A computer algebra system written in pure Python
https://sympy.org/
Other
12.88k stars 4.42k forks source link

Nonlinsolve takes much time or hangs #22932

Open proy87 opened 2 years ago

proy87 commented 2 years ago

nonlinsolve([x**2*y**2 - 150*x**2*y - (2*x**2*y - 225*x**2 + 6750000)**2/360000, x**2*(y - 88)**2 - 5625*x**2 - (2*x**2*y - 225*x**2 + 10125000)**2/360000], x, y) takes time, but not much. nonlinsolve([x**2*y**2 - 150*x**2*y - (2*x**2*y - 225*x**2 + 6750000)**2/360000, x**2*(y - 88)**2 - 5625*x**2 - (2*x**2*y - 225*x**2 + 10125000)**2/360000], y, x) - I didn't wait for the end of the calculations.

oscarbenjamin commented 2 years ago

The time seems to be taken in this call to expand: https://github.com/sympy/sympy/blob/c4e836cdf73fc6aa7bab6a86719a0f08861ffb1d/sympy/polys/polyroots.py#L378 That call can often make the expressions for the roots more complicated than they need to be and is also often a cause of slowness. I think it would be good to try removing it to see if it causes any problems.

As an example of what that line does. For one of the roots computed during the call to nonlinsolve that is slow we have

ipdb> p alpha
sqrt(3)*sqrt(-(719416405318363*127**(1/3) + 307120199*(1707119728864301702681 + 16443911250*sqrt(64948227782682816267)*I)**(1/3) + 127**(2/3)*(1707119728864301702681 + 16443911250*sqrt(64948227782682816267)*I)**(2/3))/(1707119728864301702681 + 16443911250*sqrt(64948227782682816267)*I)**(1/3))/2418
ipdb> p alpha.expand(complex=True)
-sqrt(3)*I*(94322816633799601*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**4 + 91365883475432101*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*sin(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 78008530546*sqrt(5664696104869)*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**3*sin(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) + 91365883475432101*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 91365883475432101*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 182731766950864202*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) + 78008530546*sqrt(5664696104869)*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)*cos(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) + 78008530546*sqrt(5664696104869)*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) + 188645633267599202*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 91365883475432101*sin(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 78008530546*sqrt(5664696104869)*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)*sin(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 91365883475432101*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 91365883475432101*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 94322816633799601*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**4 + 182731766950864202*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) + 78008530546*sqrt(5664696104869)*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**3*cos(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) + 78008530546*sqrt(5664696104869)*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**3)**(1/4)*sin(-I*log(94322816633799601*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**4 + 91365883475432101*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*sin(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 78008530546*sqrt(5664696104869)*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**3*sin(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) + 91365883475432101*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 91365883475432101*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 182731766950864202*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) + 78008530546*sqrt(5664696104869)*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)*cos(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) + 78008530546*sqrt(5664696104869)*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) + 188645633267599202*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 91365883475432101*sin(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 78008530546*sqrt(5664696104869)*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)*sin(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 91365883475432101*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 91365883475432101*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 94322816633799601*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**4 + 182731766950864202*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) + 78008530546*sqrt(5664696104869)*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**3*cos(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) + 78008530546*sqrt(5664696104869)*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**3)/4 + I*log(-307120199*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 - 127*sqrt(5664696104869)*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) - 127*sqrt(5664696104869)*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)*cos(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) - 127*sqrt(5664696104869)*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)*sin(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) - 307120199*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 - 127*sqrt(5664696104869)*I*sin(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) + 127*sqrt(5664696104869)*I*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)*cos(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) + 127*sqrt(5664696104869)*I*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3))/2)/2418 + sqrt(3)*(94322816633799601*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**4 + 91365883475432101*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*sin(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 78008530546*sqrt(5664696104869)*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**3*sin(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) + 91365883475432101*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 91365883475432101*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 182731766950864202*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) + 78008530546*sqrt(5664696104869)*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)*cos(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) + 78008530546*sqrt(5664696104869)*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) + 188645633267599202*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 91365883475432101*sin(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 78008530546*sqrt(5664696104869)*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)*sin(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 91365883475432101*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 91365883475432101*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 94322816633799601*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**4 + 182731766950864202*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) + 78008530546*sqrt(5664696104869)*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**3*cos(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) + 78008530546*sqrt(5664696104869)*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**3)**(1/4)*cos(-I*log(94322816633799601*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**4 + 91365883475432101*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*sin(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 78008530546*sqrt(5664696104869)*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**3*sin(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) + 91365883475432101*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 91365883475432101*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 182731766950864202*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) + 78008530546*sqrt(5664696104869)*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)*cos(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) + 78008530546*sqrt(5664696104869)*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) + 188645633267599202*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 91365883475432101*sin(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 78008530546*sqrt(5664696104869)*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)*sin(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 91365883475432101*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 91365883475432101*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 + 94322816633799601*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**4 + 182731766950864202*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2*cos(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) + 78008530546*sqrt(5664696104869)*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**3*cos(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) + 78008530546*sqrt(5664696104869)*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**3)/4 + I*log(-307120199*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 - 127*sqrt(5664696104869)*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) - 127*sqrt(5664696104869)*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)*cos(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) - 127*sqrt(5664696104869)*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)*sin(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) - 307120199*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)**2 - 127*sqrt(5664696104869)*I*sin(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)*cos(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) + 127*sqrt(5664696104869)*I*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3)*cos(2*atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3) + 127*sqrt(5664696104869)*I*sin(atan(16443911250*sqrt(64948227782682816267)/1707119728864301702681)/3))/2)/2418

So not only does it take a long time but the expanded expression is much more complicated.

proy87 commented 2 years ago

Why do the times of execution differ depending on the order of variables?

oscarbenjamin commented 2 years ago

Because it tries to solve the variables in a different order. I think there is probably a way to determine an order that is more likely to be optimal e.g. by checking the degrees:

In [4]: [{s: degree(eq, s) for s in [x,y]} for eq in eqs]
Out[4]: [{x: 4, y: 2}, {x: 4, y: 2}]

The order affects the computed Groebner basis:

In [8]: for p in groebner(eqs, [x, y]):
   ...:     print(p)
x**2 - 90725625*y**2/162409 + 106343434125*y/1299272 + 142694818875/2598544
y**3 - 9853*y**2/40 + 581009*y/40 + 289986841/20320

In [9]: for p in groebner(eqs, [y, x]):
   ...:     print(p)
162409*x**4/226032714843750 + 307120199*x**2/160734375000 + y - 237291/5080
x**6 + 1727551119375*x**4/649636 + 48005967041015625*x**2/649636 - 563135147094726562500/162409

Another missing optimisation here is noticing that we deflate the polynomials using the substitution z = x**2.

smichr commented 2 years ago

Notice that a common sub-expression makes the difference linear in y:

>>> eqs
(-x**4*y**2/90000 + x**4*y/400 - 9*x**4/64 + x**2*y**2 - 225*x**2*y + 16875*x**2/2 - 126562500, 
 -x**4*y**2/90000 + x**4*y/400 - 9*x**4/64 + x**2*y**2 - 577*x**2*y/2 + 59101*x**2/4 - 284765625)
>>> eqs[0] - eqs[1]
127*x**2*y/2 - 25351*x**2/4 + 158203125
>>> solve(_, y)
[25351/254 - 316406250/(127*x**2)]
>>> yy = _[0]

subtituting into either expression and solving for x**2 gives three possibilities:

>>> [i.n(chop=True) for i in solve(eqs[0].subs(y,yy),x**2)]
[-53314.1112076742, 24722.6045291583, -2630668.63188369]
>>> [i.n(chop=True) for i in solve(eqs[1].subs(y,yy),x**2)]
[-53314.1112076742, 24722.6045291583, -2630668.63188369]