error in using solve_poly_system and solve_triangulated
36 views
Skip to first unread message
Junwei Huang
Nov 23, 2014, 4:13:02 PM11/23/14
to sy...@googlegroups.com
Hi
I am quite new to sympy. I found sympy as I was searching a way to solve a system of 3 6-order polynomials for 3 unknowns. I tried to solve this system using either solve_poly_system or solve_triangulated but got the same error. Here is the part of code:
"
eq1 = 0.
eq2 = 0.
eq3 = 0.
m=0
for k in range(0,7):
for j in range(0,7-k):
for i in range(0,7-k-j):
if mod(i+j+k,2)==0:
eq1 = eq1 + e1C[m]*p**i*q**j*r**k
eq2 = eq2 + e2C[m]*p**i*q**j*r**k
eq3 = eq3 + e3C[m]*p**i*q**j*r**k
m=m+1
#rr = solve_poly_system([eq1, eq2, eq3], p, q, r)
rr= solve_triangulated([eq1, eq2, eq3], p, q, r)
"
e1C, e2C, and e3C are constant coefficients and eq1, eq2, and eq3 are the three polynomial equations. I got the following error:
--------------------------------------------------------------------------
KeyError Traceback (most recent call last
<ipython-input-12-9cbae56e6534> in <module>()
----> 1 rr = sympy.solve_triangulated([eq1, eq2, eq3], p, q, r)
/usr/local/lib/python2.7/dist-packages/sympy/solvers/polysys.pyc in solve_trngulated(polys, *gens, **args)
263
264 """
--> 265 G = groebner(polys, gens, polys=True)
266 G = list(reversed(G))
267
/usr/local/lib/python2.7/dist-packages/sympy/polys/polytools.pyc in groebner, *gens, **args)
6380
6381 """
-> 6382 return GroebnerBasis(F, *gens, **args)
6383
6384
/usr/local/lib/python2.7/dist-packages/sympy/polys/polytools.pyc in __new__(s, F, *gens, **args)
6420 polys[i] = ring.from_dict(poly.rep.to_dict())
6421
-> 6422 G = _groebner(polys, ring, method=opt.method)
6423 G = [Poly._from_dict(g, opt) for g in G]
6424
/usr/local/lib/python2.7/dist-packages/sympy/polys/groebnertools.pyc in groeer(seq, ring, method)
43 seq = [ s.set_ring(ring) for s in seq ]
44
---> 45 G = _groebner(seq, ring)
46
47 if orig is not None:
/usr/local/lib/python2.7/dist-packages/sympy/polys/groebnertools.pyc in _buchberger(f, ring)
236 # ordering divisors is on average more efficient [Cox] page 111
237 G1 = sorted(G, key=lambda g: order(f[g].LM))
--> 238 ht = normal(h, G1)
239
240 if ht:
/usr/local/lib/python2.7/dist-packages/sympy/polys/groebnertools.pyc in normal(g, J)
102
103 def normal(g, J):
--> 104 h = g.rem([ f[j] for j in J ])
105
106 if not h:
/usr/local/lib/python2.7/dist-packages/sympy/polys/rings.pyc in rem(f, G)
1419 c1 = get(m1, zero) - c*cg
1420 if not c1:
-> 1421 del f[m1]
1422 else:
1423 f[m1] = c1
KeyError: (0, 0, 7)
-------------------------------------------------------
eq1, eq2, and eq3 are like this:
In[13]: eq1.simplify()
Out[13]: 5105.00458755661*p**6 - 108.473959633689*p**5*q + 2402.41285238498*p**5*r + 9008.47267975219*p**4*q**2 - 1255.39607773283*p**4*q*r + 11277.0701891541*p**4*r**2 - 1307.26969182159*p**4 + 2011.04932840783*p**3*q**3 + 4868.72300813206*p**3*q**2*r + 2350.26623001839*p**3*q*r**2 + 83.2810004274752*p**3*q + 5306.9841358631*p**3*r**3 - 410.133684217865*p**3*r + 3581.05378966718*p**2*q**4 - 1494.32987338524*p**2*q**3*r + 9663.32292693404*p**2*q**2*r**2 - 1603.980521043*p**2*q**2 - 1554.39827080096*p**2*q*r**3 + 146.877127581136*p**2*q*r + 6189.03791479042*p**2*r**4 - 1951.0971844993*p**2*r**2 + 104.459609097912*p**2 + 2059.68147982275*p*q**5 + 2434.36150406603*p*q**4*r + 4573.37297182227*p*q**3*r**2 - 129.155224299263*p*q**3 + 5326.78736764999*p*q**2*r**3 - 414.231578475927*p*q**2*r + 2541.58989447865*p*q*r**4 - 163.881993246563*p*q*r**2 - 7.32935355549153*p*q + 2912.55844639838*p*r**5 - 459.09272674491*p*r**3 + 16.3862149483427*p*r - 267.46300018096*q**6 - 251.079215546565*q**5*r - 551.942706049301*q**4*r**2 - 383.602775058541*q**4 - 518.132756933654*q**3*r**3 + 48.9590425270454*q**3*r - 281.245218262232*q**2*r**4 - 921.212155372886*q**2*r**2 + 67.1547741681612*q**2 - 264.01718641355*q*r**5 + 51.0079896560767*q*r**3 - 2.33274553975296*q*r - 541.968982299977*r**4 + 79.5536696127415*r**2 - 2.47202741879275
In [14]: eq2.simplify()
Out[14]: -533.530507020419*p**6 + 4905.43084014754*p**5*q - 251.079215546566*p**5*r - 582.162658435345*p**4*q**2 + 2434.36150406603*p**4*q*r - 1101.00564042678*p**4*r**2 - 114.412696339151*p**4 + 8753.94762999568*p**3*q**3 - 1494.32987338524*p**3*q**2*r + 10767.2083390326*p**3*q*r**2 - 723.759577198517*p**3*q - 518.132756933655*p**3*r**3 + 48.9590425270454*p**3*r + 2518.77143858425*p**2*q**4 + 4868.72300813206*p**2*q**3*r + 2371.36169096872*p**2*q**2*r**2 - 570.416841704303*p**2*q**2 + 5326.78736764999*p**2*q*r**3 - 414.231578475927*p**2*q*r - 561.023034195364*p**2*r**4 - 380.66073794717*p**2*r**2 + 29.9539788526177*p**2 + 3767.68958665181*p*q**5 - 1255.39607773283*p*q**4*r + 9621.24207100615*p*q**3*r**2 - 610.744183471118*p*q**3 - 1554.39827080096*p*q**2*r**3 + 146.877127581136*p*q**2*r + 5907.79269652819*p*q*r**4 - 812.539281619356*p*q*r**2 + 22.3950435470699*p*q - 264.01718641355*p*r**5 + 51.0079896560767*p*r**3 - 2.33274553975296*p*r + 2559.17857546841*q**6 + 2402.41285238498*q**5*r + 5653.28315129873*q**4*r**2 - 655.344482123229*q**4 + 5306.9841358631*q**3*r**3 - 410.133684217865*q**3*r + 3102.61292867402*q**2*r**4 - 978.100220594949*q**2*r**2 + 52.3664159395267*q**2 + 2912.55844639838*q*r**5 - 459.09272674491*q*r**3 + 16.3862149483428*q*r - 271.69327358712*r**4 + 39.8808744205781*r**2 - 1.23924660588265
In [15]: eq3.simplify()
Out[15]: 5638.53509457703*p**5*r + 624.630294795803*p**4*q*r + 2653.49206793155*p**4*r**2 - 205.066842108932*p**4 + 10215.2656329833*p**3*q**2*r - 1036.26551386731*p**3*q*r**2 + 97.9180850540909*p**3*q + 12378.0758295808*p**3*r**3 - 975.548592249651*p**3*r + 3472.3673313955*p**2*q**3*r + 5326.78736764999*p**2*q**2*r**2 - 414.231578475927*p**2*q**2 + 3961.13372056658*p**2*q*r**3 - 272.271365596867*p**2*q*r + 5825.11689279675*p**2*r**4 - 918.18545348982*p**2*r**2 + 32.7724298966855*p**2 + 4534.64968247843*p*q**4*r - 1036.26551386731*p*q**3*r**2 + 97.9180850540908*p*q**3 + 11253.0949565319*p*q**2*r**3 - 866.875718496121*p*q**2*r - 1056.0687456542*p*q*r**4 + 204.031958624307*p*q*r**2 - 9.33098215901183*p*q + 6750.06094898578*p*r**5 - 1083.93796459995*p*r**3 + 39.7768348063708*p*r + 2826.64157564937*q**5*r + 2653.49206793155*q**4*r**2 - 205.066842108932*q**4 + 6205.22585734804*q**3*r**3 - 489.050110297474*q**3*r + 5825.11689279675*q**2*r**4 - 918.18545348982*q**2*r**2 + 32.7724298966855*q**2 + 3383.85814693625*q*r**5 - 543.386547174239*q*r**3 + 19.9404372102891*q*r + 3176.57563281193*r**6 - 765.15107460148*r**4 + 56.1568814642871*r**2 - 1.16333498638407
Could anyone please suggest what is causing this error? Is there any other way of solving this system of polynomial equations in Python? Thanks very much.
I am quite new to sympy. I found sympy as I was searching a way to solve a system of 3 6-order polynomials for 3 unknowns. I tried to solve this system using either solve_poly_system or solve_triangulated but got the same error. Here is the part of code:
"
eq1 = 0.
eq2 = 0.
eq3 = 0.
m=0
for k in range(0,7):
for j in range(0,7-k):
for i in range(0,7-k-j):
if mod(i+j+k,2)==0:
eq1 = eq1 + e1C[m]*p**i*q**j*r**k
eq2 = eq2 + e2C[m]*p**i*q**j*r**k
eq3 = eq3 + e3C[m]*p**i*q**j*r**k
m=m+1
#rr = solve_poly_system([eq1, eq2, eq3], p, q, r)
rr= solve_triangulated([eq1, eq2, eq3], p, q, r)
"
e1C, e2C, and e3C are constant coefficients and eq1, eq2, and eq3 are the three polynomial equations. I got the following error:
--------------------------------------------------------------------------
KeyError Traceback (most recent call last
<ipython-input-12-9cbae56e6534> in <module>()
----> 1 rr = sympy.solve_triangulated([eq1, eq2, eq3], p, q, r)
/usr/local/lib/python2.7/dist-packages/sympy/solvers/polysys.pyc in solve_trngulated(polys, *gens, **args)
263
264 """
--> 265 G = groebner(polys, gens, polys=True)
266 G = list(reversed(G))
267
/usr/local/lib/python2.7/dist-packages/sympy/polys/polytools.pyc in groebner, *gens, **args)
6380
6381 """
-> 6382 return GroebnerBasis(F, *gens, **args)
6383
6384
/usr/local/lib/python2.7/dist-packages/sympy/polys/polytools.pyc in __new__(s, F, *gens, **args)
6420 polys[i] = ring.from_dict(poly.rep.to_dict())
6421
-> 6422 G = _groebner(polys, ring, method=opt.method)
6423 G = [Poly._from_dict(g, opt) for g in G]
6424
/usr/local/lib/python2.7/dist-packages/sympy/polys/groebnertools.pyc in groeer(seq, ring, method)
43 seq = [ s.set_ring(ring) for s in seq ]
44
---> 45 G = _groebner(seq, ring)
46
47 if orig is not None:
/usr/local/lib/python2.7/dist-packages/sympy/polys/groebnertools.pyc in _buchberger(f, ring)
236 # ordering divisors is on average more efficient [Cox] page 111
237 G1 = sorted(G, key=lambda g: order(f[g].LM))
--> 238 ht = normal(h, G1)
239
240 if ht:
/usr/local/lib/python2.7/dist-packages/sympy/polys/groebnertools.pyc in normal(g, J)
102
103 def normal(g, J):
--> 104 h = g.rem([ f[j] for j in J ])
105
106 if not h:
/usr/local/lib/python2.7/dist-packages/sympy/polys/rings.pyc in rem(f, G)
1419 c1 = get(m1, zero) - c*cg
1420 if not c1:
-> 1421 del f[m1]
1422 else:
1423 f[m1] = c1
KeyError: (0, 0, 7)
-------------------------------------------------------
eq1, eq2, and eq3 are like this:
In[13]: eq1.simplify()
Out[13]: 5105.00458755661*p**6 - 108.473959633689*p**5*q + 2402.41285238498*p**5*r + 9008.47267975219*p**4*q**2 - 1255.39607773283*p**4*q*r + 11277.0701891541*p**4*r**2 - 1307.26969182159*p**4 + 2011.04932840783*p**3*q**3 + 4868.72300813206*p**3*q**2*r + 2350.26623001839*p**3*q*r**2 + 83.2810004274752*p**3*q + 5306.9841358631*p**3*r**3 - 410.133684217865*p**3*r + 3581.05378966718*p**2*q**4 - 1494.32987338524*p**2*q**3*r + 9663.32292693404*p**2*q**2*r**2 - 1603.980521043*p**2*q**2 - 1554.39827080096*p**2*q*r**3 + 146.877127581136*p**2*q*r + 6189.03791479042*p**2*r**4 - 1951.0971844993*p**2*r**2 + 104.459609097912*p**2 + 2059.68147982275*p*q**5 + 2434.36150406603*p*q**4*r + 4573.37297182227*p*q**3*r**2 - 129.155224299263*p*q**3 + 5326.78736764999*p*q**2*r**3 - 414.231578475927*p*q**2*r + 2541.58989447865*p*q*r**4 - 163.881993246563*p*q*r**2 - 7.32935355549153*p*q + 2912.55844639838*p*r**5 - 459.09272674491*p*r**3 + 16.3862149483427*p*r - 267.46300018096*q**6 - 251.079215546565*q**5*r - 551.942706049301*q**4*r**2 - 383.602775058541*q**4 - 518.132756933654*q**3*r**3 + 48.9590425270454*q**3*r - 281.245218262232*q**2*r**4 - 921.212155372886*q**2*r**2 + 67.1547741681612*q**2 - 264.01718641355*q*r**5 + 51.0079896560767*q*r**3 - 2.33274553975296*q*r - 541.968982299977*r**4 + 79.5536696127415*r**2 - 2.47202741879275
In [14]: eq2.simplify()
Out[14]: -533.530507020419*p**6 + 4905.43084014754*p**5*q - 251.079215546566*p**5*r - 582.162658435345*p**4*q**2 + 2434.36150406603*p**4*q*r - 1101.00564042678*p**4*r**2 - 114.412696339151*p**4 + 8753.94762999568*p**3*q**3 - 1494.32987338524*p**3*q**2*r + 10767.2083390326*p**3*q*r**2 - 723.759577198517*p**3*q - 518.132756933655*p**3*r**3 + 48.9590425270454*p**3*r + 2518.77143858425*p**2*q**4 + 4868.72300813206*p**2*q**3*r + 2371.36169096872*p**2*q**2*r**2 - 570.416841704303*p**2*q**2 + 5326.78736764999*p**2*q*r**3 - 414.231578475927*p**2*q*r - 561.023034195364*p**2*r**4 - 380.66073794717*p**2*r**2 + 29.9539788526177*p**2 + 3767.68958665181*p*q**5 - 1255.39607773283*p*q**4*r + 9621.24207100615*p*q**3*r**2 - 610.744183471118*p*q**3 - 1554.39827080096*p*q**2*r**3 + 146.877127581136*p*q**2*r + 5907.79269652819*p*q*r**4 - 812.539281619356*p*q*r**2 + 22.3950435470699*p*q - 264.01718641355*p*r**5 + 51.0079896560767*p*r**3 - 2.33274553975296*p*r + 2559.17857546841*q**6 + 2402.41285238498*q**5*r + 5653.28315129873*q**4*r**2 - 655.344482123229*q**4 + 5306.9841358631*q**3*r**3 - 410.133684217865*q**3*r + 3102.61292867402*q**2*r**4 - 978.100220594949*q**2*r**2 + 52.3664159395267*q**2 + 2912.55844639838*q*r**5 - 459.09272674491*q*r**3 + 16.3862149483428*q*r - 271.69327358712*r**4 + 39.8808744205781*r**2 - 1.23924660588265
In [15]: eq3.simplify()
Out[15]: 5638.53509457703*p**5*r + 624.630294795803*p**4*q*r + 2653.49206793155*p**4*r**2 - 205.066842108932*p**4 + 10215.2656329833*p**3*q**2*r - 1036.26551386731*p**3*q*r**2 + 97.9180850540909*p**3*q + 12378.0758295808*p**3*r**3 - 975.548592249651*p**3*r + 3472.3673313955*p**2*q**3*r + 5326.78736764999*p**2*q**2*r**2 - 414.231578475927*p**2*q**2 + 3961.13372056658*p**2*q*r**3 - 272.271365596867*p**2*q*r + 5825.11689279675*p**2*r**4 - 918.18545348982*p**2*r**2 + 32.7724298966855*p**2 + 4534.64968247843*p*q**4*r - 1036.26551386731*p*q**3*r**2 + 97.9180850540908*p*q**3 + 11253.0949565319*p*q**2*r**3 - 866.875718496121*p*q**2*r - 1056.0687456542*p*q*r**4 + 204.031958624307*p*q*r**2 - 9.33098215901183*p*q + 6750.06094898578*p*r**5 - 1083.93796459995*p*r**3 + 39.7768348063708*p*r + 2826.64157564937*q**5*r + 2653.49206793155*q**4*r**2 - 205.066842108932*q**4 + 6205.22585734804*q**3*r**3 - 489.050110297474*q**3*r + 5825.11689279675*q**2*r**4 - 918.18545348982*q**2*r**2 + 32.7724298966855*q**2 + 3383.85814693625*q*r**5 - 543.386547174239*q*r**3 + 19.9404372102891*q*r + 3176.57563281193*r**6 - 765.15107460148*r**4 + 56.1568814642871*r**2 - 1.16333498638407
Could anyone please suggest what is causing this error? Is there any other way of solving this system of polynomial equations in Python? Thanks very much.
Mateusz Paprocki
Nov 24, 2014, 3:12:48 AM11/24/14
to sympy
Hi,
SymPy is unable to detect (floating point) zero, so intermediate
polynomials have erroneous terms, causing this error. This is a bug
and has to be fixed, thought at this point I'm not sure where exactly
the issue occurs.
Mateusz
> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/sympy/f2516e57-6bf9-4d41-b82a-8b788c6024a5%40googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
polynomials have erroneous terms, causing this error. This is a bug
and has to be fixed, thought at this point I'm not sure where exactly
the issue occurs.
Mateusz
> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/sympy/f2516e57-6bf9-4d41-b82a-8b788c6024a5%40googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
Junwei Huang
Nov 24, 2014, 1:45:08 PM11/24/14
to sy...@googlegroups.com
Thanks Mateusz
for your help. Since I have copied the three polynomial equations in the original post, could you or any developer use them to track down the bug?
btw, Mathematica(TM) NSolve can handle it and gave 6 solutions. But I would really like to see sympy can do the same work, as I really like python coding environment.
Thanks,
Junwei
btw, Mathematica(TM) NSolve can handle it and gave 6 solutions. But I would really like to see sympy can do the same work, as I really like python coding environment.
Thanks,
Junwei
Aaron Meurer
Nov 24, 2014, 7:01:26 PM11/24/14
to sy...@googlegroups.com
Mathematica's NSolve likely uses numerical algorithms. SymPy also has
nsolve(), which uses numerical algorithms (I couldn't get it to work
with your system but you may have better luck if you have a better
guess at the root than I do). solve() works symbolically, which has
advantages and disadvantages. A disadvantage is that it can be slow,
and basically overkill if you only care about a numerical solution. It
also may not find a solution if it isn't available in "closed form".
An advantage is that it might find more solutions than the numeric
algorithm because you don't have to guess at a root and then hope to
converge to it.
In short, however, SymPy's symbolic routines do not perform so well
with floating point numbers.
Aaron Meurer
> https://groups.google.com/d/msgid/sympy/d78920db-6490-40f4-a436-95971df98b73%40googlegroups.com.
nsolve(), which uses numerical algorithms (I couldn't get it to work
with your system but you may have better luck if you have a better
guess at the root than I do). solve() works symbolically, which has
advantages and disadvantages. A disadvantage is that it can be slow,
and basically overkill if you only care about a numerical solution. It
also may not find a solution if it isn't available in "closed form".
An advantage is that it might find more solutions than the numeric
algorithm because you don't have to guess at a root and then hope to
converge to it.
In short, however, SymPy's symbolic routines do not perform so well
with floating point numbers.
Aaron Meurer
Junwei Huang
Nov 24, 2014, 9:10:42 PM11/24/14
to sy...@googlegroups.com
Thank you, Aaron Meurer.
Sympy solve took so long that I lost my patient and killed it.
Sympy nsolve does give me one solution, i.e., the solution close to the initial guess, when the guess is good. Although Mathematica gave all 6 solutions, Sympy might work for me at this moment. Thanks very much.
Sympy solve took so long that I lost my patient and killed it.
Sympy nsolve does give me one solution, i.e., the solution close to the initial guess, when the guess is good. Although Mathematica gave all 6 solutions, Sympy might work for me at this moment. Thanks very much.
You received this message because you are subscribed to a topic in the Google Groups "sympy" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/sympy/_nONKfongyk/unsubscribe.
To unsubscribe from this group and all its topics, send an email to sympy+un...@googlegroups.com.
To post to this group, send email to sy...@googlegroups.com.
Visit this group at http://groups.google.com/group/sympy.
To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAKgW%3D6K%2Bb0o%3Do0jmsa60hDstt_LrUHqaNRkr%2Bbd2dtk-ygwhsQ%40mail.gmail.com.
Reply all
Reply to author
Forward
0 new messages