how to solve this equations with sympy?

[zdlspace0528]

import sympy

x, y, z, e1, e2, c1, c2, c3, c5, s1, s2, s3, g1, g2 = sympy.symbols("x y z e1 e2  c1 c2 c3  c5 s1 s2 s3 g1 g2")

fx = x*(x - 1)*(c1 + c2*y*z - c2*y - c2*z - e1 - e2*y*z + e2*y + 

e2*z - g1*y - g2*z - s1*y)

fy = y*(y - 1)*(c3 + g1*x - g1 + s1*x)

fz = -z*(z - 1)*(-c5 + s2 + s3)

sympy.solve([fx, fy, fz], [x, y, z])

Why can not the solve finish? I run for long long time, but it can not finish.

[oscarbenjamin]

I think that this is a valid performance issue. It should not be this slow.

[oscarbenjamin]

One problem is that the solvers don't just factorise the polynomials to consider different cases: instead everything is passed to groebner. Since computing Groebner bases has double-exponential complexity every effort should be made to reduce the problem first. Here the different factors of each equation can be considered separately reducing the complexity substantially. In this case there are 27 combinations of factors of which most give very simple linear systems. Only the 3rd factor in fx leads to anything nontrivial. Solving the second two equations for (x,z) and substituting the first gves an equation for y that is easily solvable:

In [21]: sols = solve([fy, fz], [x, z], dict=True)

In [22]: sols
Out[22]: 
⎡⎧   -c₃ + g₁      ⎫  ⎧   -c₃ + g₁      ⎫⎤
⎢⎨x: ────────, z: 0⎬, ⎨x: ────────, z: 1⎬⎥
⎣⎩   g₁ + s₁       ⎭  ⎩   g₁ + s₁       ⎭⎦

In [23]: fx.subs(sols[0])
Out[23]: 
           ⎛-c₃ + g₁    ⎞                                      
(-c₃ + g₁)⋅⎜──────── - 1⎟⋅(c₁ - c₂⋅y - e₁ + e₂⋅y - g₁⋅y - s₁⋅y)
           ⎝g₁ + s₁     ⎠                                      
───────────────────────────────────────────────────────────────
                            g₁ + s₁                            

In [24]: solve(fx.subs(sols[1]), y)
Out[24]: 
⎡c₁ - c₂ - e₁ + e₂ - g₂⎤
⎢──────────────────────⎥
⎣       g₁ + s₁        ⎦

In [25]: solve(fx.subs(sols[0]), y)
Out[25]: 
⎡     c₁ - e₁     ⎤
⎢─────────────────⎥
⎣c₂ - e₂ + g₁ + s₁⎦

That gives 2 of 10 possible solutions. We can also solve the second two equations for y and z and use those in fx to find x giving:

In [26]: sols = solve([fy, fz], [y, z], dict=True)

In [27]: sols
Out[27]: [{y: 0, z: 0}, {y: 0, z: 1}, {y: 1, z: 0}, {y: 1, z: 1}]

In [28]: solve(fx.subs(sols[0]), x)
Out[28]: [0, 1]

In [29]: solve(fx.subs(sols[1]), x)
Out[29]: [0, 1]

In [30]: solve(fx.subs(sols[2]), x)
Out[30]: [0, 1]

In [31]: solve(fx.subs(sols[3]), x)
Out[31]: [0, 1]

All combinations gives another 8 solutions.

Ways that solve etc can be improved:

• Factorise polynomials before calculating Groebner bases.
• Linear equations (factors) should be solved before any remaining polynomial equations.