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