[SOLVED] Sympy's solve() function doesn't yield all the solutions

These are your equations:

In [45]: from sympy import  solve
    ...: 
    ...: from sympy import symbols
    ...: x,y,a,b,c,d,g,h = symbols('x y a b c d g h')
    ...: a = 0.12
    ...: b = -2.031529100521498e-30
    ...: c = b
    ...: d = 1
    ...: g = 11
    ...: h = 0
    ...: F = (y + a/2)*x*(2*x**2 - 2*y*a + 2*d**2 + b + c) - (x*y*g**2)
    ...: G = (x**2 - y*a + b + d**2)*(x**2 - y*a + c+d**2) - 4*((y+a/2)**2)*x**2 - (h**2 + (x**2)*(g*
    ...: *2))

I'm going to convert the floats to rationals because floats and polynomials don't mix well in general. I will just note here that your value for b is so small as to be effectively zero (in standard 64-bit floating point). If you want b to be nonzero then you need to avoid floats from the outset because it has already disappeared after executing the code shown above.

In [46]: F, G = nsimplify(F), nsimplify(G) # don't use floats with polynomials

In [47]: F
Out[47]: 
                        ⎛   2   6⋅y    ⎞
-121⋅x⋅y + x⋅(y + 3/50)⋅⎜2⋅x  - ─── + 2⎟
                        ⎝        25    ⎠

In [48]: G
Out[48]: 
                                            2
     2           2        2   ⎛ 2   3⋅y    ⎞ 
- 4⋅x ⋅(y + 3/50)  - 121⋅x  + ⎜x  - ─── + 1⎟ 
                              ⎝      25    ⎠

The equation F factorises so let's take the two factors in turn:

In [56]: F.factor()
Out[56]: 
  ⎛      2         2        2               ⎞
x⋅⎝1250⋅x ⋅y + 75⋅x  - 150⋅y  - 74384⋅y + 75⎠
─────────────────────────────────────────────
                     625                     

In [57]: _, F1, F2 = F.factor().args

In [58]: F1
Out[58]: x

In [59]: solve([F1, G], [x, y])
Out[59]: [(0, 25/3)]

I believe that this is the solution that you referred to. You describe this as incorrect but it is definitely correct for the equations as shown in the code.

That was the easy one. Now for F2 we compute a Groebner basis:

In [61]: gb = groebner([F2, G], [x, y], method='f5b')

In [62]: print(gb[0])
x**2 + 16*y**4/121 + 595216*y**3/9075 + 892608*y**2/75625 + 1844017554*y/1890625 - 75643/75625

In [63]: print(gb[1])
y**5 + 74411*y**4/150 + 148777*y**3/1250 + 1845583341*y**2/250000 + 691424307*y/781250 - 113451/125000

We can see that gb[1] is a quintic in y. This quintic is irreducible and likely there are no radical formulae for its roots (Abel-Ruffini) but when the coefficients are explicit rational numbers SymPy can represent the roots as RootOf:

In [66]: r1, r2, r3, r4, r5 = Poly(gb[-1]).all_roots()

In [67]: r1
Out[67]: 
       ⎛          5               4               3                 2                              
CRootOf⎝18750000⋅x  + 9301375000⋅x  + 2231655000⋅x  + 138418750575⋅x  + 16594183368⋅x - 17017650, 0

⎞
⎠

Each of these roots can be used in the first equation to solve for x e.g.:

In [69]: print(solve([gb[0],y-r1], [x, y]))
[(-sqrt(-2250000*CRootOf(18750000*x**5 + 9301375000*x**4 + 2231655000*x**3 + 138418750575*x**2 + 16594183368*x - 17017650, 0)**4 - 200836800*CRootOf(18750000*x**5 + 9301375000*x**4 + 2231655000*x**3 + 138418750575*x**2 + 16594183368*x - 17017650, 0)**2 + 17019675 - 16596157986*CRootOf(18750000*x**5 + 9301375000*x**4 + 2231655000*x**3 + 138418750575*x**2 + 16594183368*x - 17017650, 0) - 1116030000*CRootOf(18750000*x**5 + 9301375000*x**4 + 2231655000*x**3 + 138418750575*x**2 + 16594183368*x - 17017650, 0)**3)/4125, CRootOf(18750000*x**5 + 9301375000*x**4 + 2231655000*x**3 + 138418750575*x**2 + 16594183368*x - 17017650, 0)), (sqrt(-2250000*CRootOf(18750000*x**5 + 9301375000*x**4 + 2231655000*x**3 + 138418750575*x**2 + 16594183368*x - 17017650, 0)**4 - 200836800*CRootOf(18750000*x**5 + 9301375000*x**4 + 2231655000*x**3 + 138418750575*x**2 + 16594183368*x - 17017650, 0)**2 + 17019675 - 16596157986*CRootOf(18750000*x**5 + 9301375000*x**4 + 2231655000*x**3 + 138418750575*x**2 + 16594183368*x - 17017650, 0) - 1116030000*CRootOf(18750000*x**5 + 9301375000*x**4 + 2231655000*x**3 + 138418750575*x**2 + 16594183368*x - 17017650, 0)**3)/4125, CRootOf(18750000*x**5 + 9301375000*x**4 + 2231655000*x**3 + 138418750575*x**2 + 16594183368*x - 17017650, 0))]

There are other ways to do this but that gives all the solutions in terms of RootOf. I don't know if that is what you wanted but it is possible that you were hoping for SymPy to return something that cannot possibly be returned. Here is the full set of solutions:

In [71]: sx1, sx2 = solve(gb[0], x)

In [72]: sols = solve([F1, G], [x, y], dict=True)

In [73]: sols.extend([{x:sx1.subs(y, r), y:r} for r in [r1,r2,r3,r4,r5]])

In [74]: sols.extend([{x:sx2.subs(y, r), y:r} for r in [r1,r2,r3,r4,r5]])

This is what they look like numerically:

In [76]: for s in sols:
    ...:     print(s[x].n(3), s[y].n(3))
    ...: 
0 8.33
-0.0610 -496.
-10.9 -0.121
-0.0917 0.00102
-7.71 + 0.091*I -0.045 - 3.86*I
-7.71 - 0.091*I -0.045 + 3.86*I
0.0610 -496.
10.9 -0.121
0.0917 0.00102
7.71 - 0.091*I -0.045 - 3.86*I
7.71 + 0.091*I -0.045 + 3.86*I

The question is ambiguous as to whether you are interested in getting a solution in terms of a, b, etc as symbols rather than explicit numbers. The fact that gb[-1] is an irreducible quintic makes it very unlikely that any meaningful explicit solution exists for the case of symbolic coefficients:

https://en.wikipedia.org/wiki/Abel%E2%80%93Ruffini_theorem

In general when faced with this kind of problem it is better to consider a different approach. Explicit solutions can be impossible (as here) or just too complicated to be of much use (e.g. when using Cardano's quartic formula). Few people seem to heed this advice when I give it but actually it is better not to seek explicit solutions at all in cases like this. Your original polynomial equations are already a nicer representation of the solutions than anything that any radical formula could give: for symbolic work it is usually better to use those polynomials as an implicit representation of the solutions to the system. I can't advise how to do that because I don't know what your intended next steps are (ask a new question if you want an answer to that).