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

Question

I need to solve the following nonlinear system of equations :

I am using sympy's solve() function to obtain the solutions to the system of equations. Normally, solve() yields all the solutions. I have specified the values of the constants a,b,c,d,g,h and expecting to get all the solutions.

But solve() yeilds only one solution and that too incorrectly. My question is:

How can I get all the symbolic solutions using a single python utility ? I don't want numerical solution.

Here, is the portion of the code that implements solve() :

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))
sol = solve([F,G],(x,y),dict=True,simplify=False)
print(sol)

Could you include the code too, with desired and actual outputs? — doneforaiur, Sep 02 '23 at 16:52
you can use the `solveset` function from the `sympy` library,It returns a set of all the solutions to the system of equations. — Sauron, Sep 02 '23 at 17:01
@bajju does solveset() solve a system of equations, or it solves only single equation ?? Actually I tried , but it showed ValueError : invalid sympy expression — bubucodex, Sep 02 '23 at 17:04
Can you include complete code that defines all of a, b, x, y etc? — Oscar Benjamin, Sep 02 '23 at 17:04
Oh, my bad but `4*((y+a/2)**2)*x**2` should be `4*(y+a/2)*x**2`. There's no square in the inner part. — doneforaiur, Sep 02 '23 at 17:06
Oh you've edited your image, I thought I'm losing my mind lol. Just wanted to make sure the code is correct. :^) — doneforaiur, Sep 02 '23 at 17:11
@jared Because by independent means I could verify it. When I plotted the two functions they crosses at points none of which is same as the obtained solution. — bubucodex, Sep 02 '23 at 17:37
While there are definitely other solutions that Sympy is not returning, the solution I get when running your code (0,8.33...) is valid. — jared, Sep 02 '23 at 17:51
agree with @jared I am getting the output `[{x: 0.0, y: 8.33333333333333}]`, it satisfies the equations when you plug in this values. — Sauron, Sep 02 '23 at 17:58
Plug in x = 0 and y = 1/0.12 and you'll see that F and G both equal 0. — jared, Sep 02 '23 at 17:58

score 0 · Accepted Answer · answered Sep 02 '23 at 17:57

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).

I was also working on this question to find solutions, and this is quite a approach but the values of `sols` are not valid, other than first values( 0, 8.33). — Sauron, Sep 02 '23 at 18:04
No, all of the solutions are valid. Try e.g. `F.subs(sols[3]).n()` — Oscar Benjamin, Sep 02 '23 at 18:05
You just need to use more than 3 digits if you want to test the solutions numerically by substituting them in. — Oscar Benjamin, Sep 02 '23 at 18:09
@OscarBenjamin Actually initially i wanted to get the solutions in terms of a, b, c... etc using solve() . But then the program keeps on running. So, i then set numerical values to the constants. Can you please explain how to obtain implicit representations ? — bubucodex, Sep 02 '23 at 18:17
@OscarBenjamin I countour plotted the two polynomials. But they didn't intersect at (0, 25/3). This made me declare that the solution is incorrect . — bubucodex, Sep 02 '23 at 18:21
It makes no sense to declare a,b,c,... as free variables and then shortly after redefine them as numerical constants. It's not wrong, just useless. — Lutz Lehmann, Sep 02 '23 at 22:06
The original equations already are an implicit representation. The question is how to use them to do what you want to do but no one can answer that without knowing what you want to do with the solutions (and as I said before: ask a new question if you are interested in that). — Oscar Benjamin, Sep 02 '23 at 23:14

score -1 · Answer 2 · answered Sep 02 '23 at 18:01

Try this:

from sympy import symbols, Eq, solve


x, y, a, b, c, d, g, h = symbols('x y a b c d g h')


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))


solution = solve([Eq(F, 0), Eq(G, 0)], (x, y))

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

2 Answers2