# NSolve erroneously gives no solution to a polynomial system

I have a polynomial system with three equations in three unknowns, the maximum degree is 26. Two equations are symmetric, i.e. eq1(x,y,z)=eq2(y,x,z). If I search for a symmetric solution with y=x, NSolve returns a number of solutions. However, if I apply NSolve to the original system, {} is returned after a while, without any errors. What could be going on here?

Edit: below are the polynomials in question.

eq1 = -400 z^2 (12500 y^3 z^5 (90 (-1+z) z+y (-130-1967 z+1701 z^2))+x^6 (25+420 z+164 z^2) (100 z^2+20 y z (10+189 z)+y^2 (100+3780 z+3321 z^2))^2+6250 x y^2 z^4 (180 (-1+z) z^2+4 y z (-350-3655 z+3213 z^2)+y^2 (-1300-35800 z-148177 z^2+120285 z^3))+5000 x^2 y^2 z^3 (60 z^2 (-25-200 z+126 z^2)+2 y z (-2375-54375 z-242440 z^2+176112 z^3)+y^2 (-3250-129775 z-1201160 z^2-2076368 z^3+1507086 z^4))+2 x^5 (25+420 z+164 z^2) (-10000 z^4-8000 y z^3 (5+92 z)+300 y^2 z^2 (-275-8545 z-45294 z^2+492 z^3)+20 y^3 z (-4250-188475 z-2030715 z^2-57492 z^3+278964 z^4)+y^4 (-32500-1942000 z-27773325 z^2-13103595 z^3+14755203 z^4+4901796 z^5))+500 x^3 y z^2 (-400 z^3 (10+89 z)+60 y z^2 (-525-10795 z-41594 z^2+11076 z^3)+4 y^2 z (-15000-535375 z-5230875 z^2-10618835 z^3+6404076 z^4)+y^3 (-32500-1699500 z-24450075 z^2-98053445 z^3-5080482 z^4+31443228 z^5))+25 x^4 z (-2000 z^4 (10+89 z)-800 y z^3 (200+4460 z+17477 z^2)+120 y^2 z^2 (-4500-153900 z-1409110 z^2-2323793 z^3+82656 z^4)+8 y^3 z (-90625-4401875 z-66277500 z^2-300577200 z^3-20833965 z^4+47108016 z^5)+y^4 (-325000-20987500 z-408812500 z^2-2528313000 z^3-1459453050 z^4+653500107 z^5+366000768 z^6)))
eq2 = -400 z^2 (10000 y^4 z^4 (-5 z (10+89 z)-2 y (25+420 z+164 z^2)+y^2 (25+420 z+164 z^2))+4000 x y^3 z^3 (-50 z^2 (10+89 z)-5 y z (200+4460 z+17477 z^2)-4 y^2 (125+4400 z+39460 z^2+15088 z^3)+y^3 (250+8925 z+81020 z^2+30996 z^3))+600 x^2 y z^2 (1875 (-1+z) z^4+500 y z^3 (-25-200 z+126 z^2)+50 y^2 z^2 (-525-10795 z-41594 z^2+11076 z^3)+5 y^3 z (-4500-153900 z-1409110 z^2-2323793 z^3+82656 z^4)+y^5 (2500+136500 z+2227025 z^2+11086740 z^3+4087044 z^4)+y^4 (-6875-329125 z-4766350 z^2-20412560 z^3-7221576 z^4+80688 z^5))+40 x^3 z (28125 (-1+z) z^5+625 y z^4 (-350-3655 z+3213 z^2)+250 y^2 z^3 (-2375-54375 z-242440 z^2+176112 z^3)+50 y^3 z^2 (-15000-535375 z-5230875 z^2-10618835 z^3+6404076 z^4)+5 y^4 z (-90625-4401875 z-66277500 z^2-300577200 z^3-20833965 z^4+47108016 z^5)+y^6 (25000+1837500 z+42668750 z^2+338995125 z^3+386232300 z^4+102937716 z^5)+y^5 (-106250-6496875 z-130624375 z^2-885247500 z^3-350209800 z^4+107736192 z^5+45750096 z^6))+x^4 (12500 z^5 (-130-1967 z+1701 z^2)+y^6 (25+420 z+164 z^2) (100+3780 z+3321 z^2)^2+6250 y z^4 (-1300-35800 z-148177 z^2+120285 z^3)+5000 y^2 z^3 (-3250-129775 z-1201160 z^2-2076368 z^3+1507086 z^4)+500 y^3 z^2 (-32500-1699500 z-24450075 z^2-98053445 z^3-5080482 z^4+31443228 z^5)+25 y^4 z (-325000-20987500 z-408812500 z^2-2528313000 z^3-1459453050 z^4+653500107 z^5+366000768 z^6)+2 y^5 (-812500-62200000 z-1515303125 z^2-12310874375 z^3-9689455125 z^4+4170740580 z^5+4478607612 z^6+803894544 z^7)))
eq3 = -400 (12500000 (-1+y) y^6 z^8 (25 z+2 y^2 (105+82 z)+5 y (5+84 z))+5000000 x (-1+y) y^5 z^7 (250 z^2+75 y z (10+147 z)+25 y^2 (20+693 z+4652 z^2)+y^3 (4200+60610 z+44772 z^2))+4 x^9 (105+82 z) (2500 z^4+500 y z^3 (20+273 z)+25 y^2 z^2 (600+16380 z+67481 z^2)+10 y^3 z (1000+40950 z+337405 z^2+201474 z^3)+y^4 (2500+136500 z+1687025 z^2+2014740 z^3+544644 z^4))^2+62500 x^2 y^4 z^5 (625 z^4 (-74+9 z)+500 y z^3 (-470-5607 z+378 z^2)+50 y^2 z^2 (-7550-193325 z-1082282 z^2+33228 z^3)+8 y^5 z (147000+4127900 z+25870635 z^2+17756198 z^3)+20 y^3 z (-10250-434275 z-5137625 z^2-16185017 z^3+61992 z^4)+y^4 (-16250-1576375 z-30283800 z^2-106933140 z^3+234123840 z^4+242064 z^5))+12500 x^3 y^3 z^4 (6250 z^5 (-50+9 z)+625 y z^4 (-3610-40630 z+6237 z^2)+500 y^2 z^3 (-11800-278675 z-1440017 z^2+168174 z^3)+50 y^3 z^2 (-135750-5012550 z-54425905 z^2-161386282 z^3+12110364 z^4)+8 y^6 z (1470000+61344500 z+729108450 z^2+2509338345 z^3+1543671402 z^4)+10 y^4 z (-323750-16777875 z-284394650 z^2-1725057270 z^3-2682939630 z^4+44504352 z^5)+y^5 (-406250-30198750 z-718189875 z^2-5601198800 z^3-5235243980 z^4+21387326112 z^5+86416848 z^6))+3125 x^4 y^2 z^3 (12500 z^6 (-74+9 z)+2500 y z^5 (-3610-40630 z+6237 z^2)+125 y^2 z^4 (-269600-6271800 z-31535234 z^2+4879809 z^3)+100 y^3 z^3 (-610500-22016600 z-231267240 z^2-676924767 z^3+95566338 z^4)+10 y^4 z^2 (-5582500-277080750 z-4514172950 z^2-27429959045 z^3-47636823957 z^4+5572844748 z^5)+4 y^7 z (14700000+814100000 z+14268765000 z^2+89719815000 z^3+175946555925 z^4+89184527498 z^5)+10 y^5 z (-2362500-151852500 z-3404834250 z^2-31477458950 z^3-107192270485 z^4-70258903867 z^5+4014973872 z^6)+y^6 (-3250000-266250000 z-7734616250 z^2-91274467875 z^3-353314512600 z^4+81822358510 z^5+682796135744 z^6+7739512272 z^7))+500 x^5 y z^2 (-2500000 z^7+62500 y z^6 (-470-5607 z+378 z^2)+12500 y^2 z^5 (-11800-278675 z-1440017 z^2+168174 z^3)+625 y^3 z^4 (-610500-22016600 z-231267240 z^2-676924767 z^3+95566338 z^4)+500 y^4 z^3 (-1077500-53309750 z-859760550 z^2-5217504375 z^3-9155373012 z^4+1373964516 z^5)+50 y^5 z^2 (-8225000-525798750 z-11613295500 z^2-108006433425 z^3-394956201355 z^4-373424805542 z^5+57710474136 z^6)+4 y^8 z (73500000+5073775000 z+117600525000 z^2+1077107192500 z^3+3626828197125 z^4+3688382662340 z^5+1139229494676 z^6)+10 y^6 z (-15375000-1227393750 z-34936781250 z^2-440165818750 z^3-2389380496125 z^4-4259671880765 z^5-543757517062 z^6+199369046880 z^7)+y^7 (-20312500-2011500000 z-71374225000 z^2-1119790209375 z^3-7315791538750 z^4-12681549935250 z^5+12305201074940 z^6+11922339291856 z^7+378251627040 z^8))+50 x^6 z (-6250000 z^8-2500000 y z^7 (30+431 z)+62500 y^2 z^6 (-7550-193325 z-1082282 z^2+33228 z^3)+12500 y^3 z^5 (-135750-5012550 z-54425905 z^2-161386282 z^3+12110364 z^4)+625 y^4 z^4 (-5582500-277080750 z-4514172950 z^2-27429959045 z^3-47636823957 z^4+5572844748 z^5)+500 y^5 z^3 (-8225000-525798750 z-11613295500 z^2-108006433425 z^3-394956201355 z^4-373424805542 z^5+57710474136 z^6)+50 y^6 z^2 (-53875000-4317337500 z-122350656250 z^2-1545553044375 z^3-8658001464750 z^4-17447002625280 z^5-4281527765016 z^6+1116768103056 z^7)+4 y^9 z (367500000+30385250000 z+876127612500 z^2+10549119450000 z^3+50543244451875 z^4+77207282058450 z^5+46560234887460 z^6+9670806703512 z^7)+50 y^7 z (-17562500-1748150000 z-61641381250 z^2-992616036250 z^3-7418598629625 z^4-21656294896665 z^5-10779809973576 z^6+837951801292 z^7+633852330048 z^8)+y^8 (-101562500-12916093750 z-561037093750 z^2-11032005028125 z^3-98895162378125 z^4-320995927116250 z^5+11541934933200 z^6+321942723499160 z^7+136511538445872 z^8+5502013053120 z^9))+x^8 (12500000 z^8 (-185+256 z)+5000000 y z^7 (-3700-43285 z+71528 z^2)+4 y^9 (105+82 z) (2500+136500 z+1687025 z^2+2014740 z^3+544644 z^4)^2+62500 y^2 z^5 (-16250-1576375 z-30283800 z^2-106933140 z^3+234123840 z^4+242064 z^5)+12500 y^3 z^4 (-406250-30198750 z-718189875 z^2-5601198800 z^3-5235243980 z^4+21387326112 z^5+86416848 z^6)+y^8 (25+420 z+164 z^2)^2 (-29225000-1286281250 z-16042309375 z^2+8122707000 z^3+54598605300 z^4+30273492096 z^5+803894544 z^6)+3125 y^4 z^3 (-3250000-266250000 z-7734616250 z^2-91274467875 z^3-353314512600 z^4+81822358510 z^5+682796135744 z^6+7739512272 z^7)+500 y^5 z^2 (-20312500-2011500000 z-71374225000 z^2-1119790209375 z^3-7315791538750 z^4-12681549935250 z^5+12305201074940 z^6+11922339291856 z^7+378251627040 z^8)+50 y^6 z (-101562500-12916093750 z-561037093750 z^2-11032005028125 z^3-98895162378125 z^4-320995927116250 z^5+11541934933200 z^6+321942723499160 z^7+136511538445872 z^8+5502013053120 z^9)+10 y^7 (-101562500-19520312500 z-1065143281250 z^2-25668581968750 z^3-289577442890625 z^4-1286535678011875 z^5-471989907683000 z^6+1653988149325800 z^7+1479996563635440 z^8+361065214964880 z^9+13535440402176 z^10))+10 x^7 (-6250000 z^8 (5+79 z)-12500000 y z^7 (20+663 z+4211 z^2)+125000 y^2 z^6 (-10250-434275 z-5137625 z^2-16185017 z^3+61992 z^4)+12500 y^3 z^5 (-323750-16777875 z-284394650 z^2-1725057270 z^3-2682939630 z^4+44504352 z^5)+3125 y^4 z^4 (-2362500-151852500 z-3404834250 z^2-31477458950 z^3-107192270485 z^4-70258903867 z^5+4014973872 z^6)+500 y^5 z^3 (-15375000-1227393750 z-34936781250 z^2-440165818750 z^3-2389380496125 z^4-4259671880765 z^5-543757517062 z^6+199369046880 z^7)+250 y^6 z^2 (-17562500-1748150000 z-61641381250 z^2-992616036250 z^3-7418598629625 z^4-21656294896665 z^5-10779809973576 z^6+837951801292 z^7+633852330048 z^8)+8 y^9 z (262500000+25286875000 z+872210062500 z^2+13019934956250 z^3+81022758665625 z^4+165148040303000 z^5+147358962336600 z^6+59875852069440 z^7+8997991630992 z^8)+20 y^7 z (-60156250-7708359375 z-334196875000 z^2-6638447356250 z^3-62930398093750 z^4-245935395612250 z^5-188063456267100 z^6+5474243584200 z^7+28342274761872 z^8+4103157149568 z^9)+y^8 (-101562500-19520312500 z-1065143281250 z^2-25668581968750 z^3-289577442890625 z^4-1286535678011875 z^5-471989907683000 z^6+1653988149325800 z^7+1479996563635440 z^8+361065214964880 z^9+13535440402176 z^10)))


Edit 2: As is noted in the comments, the system above is huge, and it takes too long to try different solution methods. Interestingly, while working on the initial problem, I came across a much simpler system, but with the same NSolve issue. Here are the equations:

foc = {-x^4 y^2+2 x^3 y^3+x^2 y^4-2 x^4 y z+4 x^3 y^2 z-x^4 y^2 z+10 x^2 y^3 z-2 x^3 y^3 z+4 x y^4 z-x^2 y^4 z-2 x^4 z^2+8 x^2 y^2 z^2+8 x y^3 z^2+2 y^4 z^2,x^4 y^4+2 x^3 y^5-x^2 y^6+8 x^4 y^3 z+18 x^3 y^4 z-2 x^4 y^4 z+4 x^2 y^5 z-4 x^3 y^5 z-2 x y^6 z-2 x^2 y^6 z+16 x^4 y^2 z^2+44 x^3 y^3 z^2-4 x^4 y^3 z^2+28 x^2 y^4 z^2-10 x^3 y^4 z^2-8 x^2 y^5 z^2-2 y^6 z^2-2 x y^6 z^2+12 x^4 y z^3+40 x^3 y^2 z^3-2 x^4 y^2 z^3+38 x^2 y^3 z^3-6 x^3 y^3 z^3+8 x y^4 z^3-7 x^2 y^4 z^3-2 y^5 z^3-4 x y^5 z^3-y^6 z^3+3 x^4 z^4+12 x^3 y z^4+14 x^2 y^2 z^4+4 x y^3 z^4-y^4 z^4,y (2 x^2 y^3+8 x^2 y^2 z+4 x y^3 z+10 x^2 y z^2+10 x y^2 z^2-2 x^2 y^2 z^2+y^3 z^2-x y^3 z^2+4 x^2 z^3+6 x y z^3-4 x^2 y z^3+2 y^2 z^3-4 x y^2 z^3-2 x^2 z^4-3 x y z^4-y^2 z^4)}


Please note that these are not symmetric anymore. NSolve[foc,Reals] gives {} without any errors. So does NSolve[Join[foc, {x >= 1, y >= 1, z >= 1}], Reals] (the constraints give the domain I am interested in). Running GroebnerBasis before NSolve still produces no results but the following error message appears from NSolve:

Infinite solution set has dimension at least 1.
Returning intersection of solutions with
(171802 x)/178835-(113492 y)/178835-(121484 z)/178835 == 1.


The best part of it, the solution exists as can be easily seen from the following contour plot (and by playing with the upper limits it seems the solution is unique):

ContourPlot3D[Evaluate@Thread@(foc == 0), {x, 1, 10}, {y, 1, 10}, {z, 1, 10}]


Solve, given the constraints on the domain, can find that unique solution. Somehow, NSolve fails in this case...

Edit 3: belisarius has pointed out that foc/.{x->0,y->0} = {0,0,0}. It is also the case that foc/.{x->0,z->0} = {0,0,0} and foc/.{y->0,z->0} = {0,0,0}. I cannot trivially delete these roots from the original equations, because the original equations do not factor so. However, these roots do seem to present the problem when using NSolve after GroebnerBasis as is evident from the error message. So, I've tried dropping all the small common factors from the Groebner basis.

gb = GroebnerBasis[foc,{x,y,z}];
NSolve[FactorList[#][[-1,1]]&/@gb,Reals]


gives the set of the remaining solutions, with the one I was looking for being among them.

So, it seems that NSolve gets caught in the positive dimensionality part of the solution, even if explicitly told the region (x>1,y>1,z>1), where that positive dimensionality is not present. When run after GroebnerBasis, NSolve at least admits it, but when run on the original equations, it seems to hide the fact.

(I accept belisarius' answer as the comment below helped me to at least have an idea of what might be going wrong.)

Edit 4: Just for fun:

addConstraints[foc_] := {Sequence @@ foc, x >= 1, y >= 1, z >= 1}
Reduce[addConstraints@Thread[foc == 0], {x, y, z}, Backsubstitution -> True] // N // Timing
{6.87911, x == 6.18083 && y == 4.58414 && z == 3.14352}

gb := GroebnerBasis[foc, {x, y, z}];
NSolve[addConstraints[FactorList[#][[-1, 1]] & /@ gb], Reals] // Timing
{2.97643, {{x -> 6.18083, y -> 4.58414, z -> 3.14352}}}

{6.84882, {{x -> 6.18083, y -> 4.58414, z -> 3.14352}}}

• Since you are looking for a numeric solution, FindRoot is sometimes more aggressive. See if it can find a solution to the original system.
– Bill
Jul 11, 2015 at 19:31
• Can you maybe include the polynomials you speak of in your post? Jul 11, 2015 at 19:38
• @J.M. I've added the polynomials. Jul 11, 2015 at 20:24
• @J.M. Economics... I am looking for a Nash equilibrium of a specific problem. I can easily find at least one equilibrium, but given that my problem has best responses whose derivatives are polynomials, I thought using NSolve on the FOCs would be a neat way to check if the solution I find is the only one. Regarding putting the constraints into NSolve - no I didn't try. Thank you for the suggestion. Trying now (but it probably will take a few hours, so results will be known later). Jul 11, 2015 at 20:38
• @J.M. It turns out I had an error in my computations, the polynomials I need are slightly different. Also, when I use NSolve to find the symmetric solution, the precision of the result is quite poor, maybe that's the core of the problem somehow. Anyways, I'll try GroebnerBasis followed by NSolve on the original equations and will let you know what comes out. Jul 11, 2015 at 21:25

One way to deal with the singularities (that lead to dimensional components of the solution space) when x, y and/or z are zero is through a change of variables that lets excessive zeros be factored out.

changevar = {y -> x s, z -> x s t};
Short /@ Factor[foc /. changevar]


We can use FactorList to get at the irreducible factors and delete constants and simple variables. We find six real solutions, only one of which has all variables positive. (It was not clear to me whether the lower bound of 1 on the solutions was imposed to try to get rid of the zeros or was part of the system.)

nfoc = Flatten@ DeleteCases[
FactorList[#][[All, 1]] & /@ (foc /. changevar),
_?NumericQ | s | x | t,
{2}];
nsol = NSolve[nfoc == {0, 0, 0}, {x, s, t}, Reals, WorkingPrecision -> 25];
sol = Thread[{x, y, z} -> #] & /@ ({x, y, z} /. changevar /. nsol)
(*
{{x ->  20.6512,   y ->  1.14137,  z -> -0.452436},
{x -> -18.1434,   y -> 28.1209,   z -> -6.50644},
{x ->   6.18083,  y ->  4.58414,  z ->  3.14352},
{x ->   1.82148,  y -> -4.7767,   z -> -1.25265},
{x ->   0.885016, y ->  1.62512,  z -> -0.298749},
{x ->  -1.20071,  y -> -0.747601, z ->  1.57008}}
*)


Verify:

foc /. sol
(*
{{0.*10^-18, 0.*10^-18, 0.*10^-20}, {0.*10^-14, 0.*10^-11, 0.*10^-15},
{0.*10^-18, 0.*10^-16, 0.*10^-18}, {0.*10^-20, 0.*10^-18, 0.*10^-20},
{0.*10^-23, 0.*10^-22, 0.*10^-23}, {-0.*10^-22, 0.*10^-21, 0.*10^-22}}
*)

• This is also the fastest solution. Removing the factors from the Groebner Basis and then using NSolve takes about 3 sec (was the fastest solution so far). This solution takes about 0.6 sec. Jul 21, 2015 at 17:35
• While this solution works for this particular instance of the problem, it is not general. E.g., changevar = {y -> z s, x -> z s t} does not solve the issue anymore, because it misses the hyperplane where y = 0, x = 0, z > 0. For the same reason, the extension of this approach to 4 or more equations can be challenging. Whereas removing simple factors from the Groebner Basis and then running NSolve should always do the trick, I guess. Jul 22, 2015 at 5:46
• @Andrei I recall checking something like that for this particular system. What happens if you fold in the plane and axis and solve each system? Something like {y -> 0, x -> z t} and {y -> 0, x -> 0}, and call NSolve on each. That's three calls total to NSolve, but each subsequent one is simpler. But if the way you suggest works, maybe it's simpler. Jul 22, 2015 at 11:08

Too much fuss :)

rf1 = Reduce[And @@ Thread[foc == 0] && x > 0 && y > 0 && z > 0, {x, y, z},
Quartics -> True, Backsubstitution -> True];

Show[ContourPlot3D[Evaluate@Thread[foc == 0], {x, 1, 7}, {y, 1, 6}, {z, 1, 6},
Mesh -> None],
Graphics3D[{Green, Sphere[{x, y, z}, .3] /. (ToRules@rf1 // N)}]]


rf1 // N
(* x == 6.18083 && y == 4.58414 && z == 3.14352 *)

• The whole problem seems caused by foc /. {x -> 0, y -> 0} being {0,0,0} Jul 13, 2015 at 10:19
• It's a good point that reduce can find the solution. Albeit, that does not answer the original question of why NSolve cannot. Also, Quartics->True is unnecessary. Jul 14, 2015 at 18:21
• @Andrei Usually questions like Why MyFun[ ] can't do ... whatever are not answerable. Mathematica inner code and algorithms are (mostly) a black box. Jul 14, 2015 at 18:33