A system of non-linear equations with no real-valued solution

Question

I am trying to solve the following system of non-linear equations with mathematica. I have tried NSolve, Solve and findRoot. All the solutions they return are complex-valued. I am still a little doubtful that there might be a mistake in my coding. I am sharing with you snippets I wrote with NSolve and Solve. Would some body please hint me if there is any mistake in my coding which prohibits finding real solutions:

fr = {
 1/ x - 1/x^6 - 1/(y - x) - 1/(z - x) == 0, 
 1/ y - 1/y^6 + 1/(y - x) - 1/(z - y) == 0, 
 1/ z - 1/z^6 + 1/(z - x) + 1/(z - y) == 0
 };
sol=Solve[fr,{x,y,z}]// Simplify;
fr /. sol[[1]] // Simplify
tab60 = Table[sol[[i]] // N[#, 60] &, {i, Length[sol]}];
tab60[[1]]
tab60[[2]]
tab60[[3]]
tab60[[4]]


fr = {(1/x)-1/(y-x) - 1/(z-x)-1/(x^6)==10^(-15), (1/y)+1/(y-x) - 1/(z-y)-1/(y^6)==10^(-5),(1/z)+1/(z-x) + 1/(z-y)-1/(z^6)==0};
fnprec=SetPrecision[fr,48]
sol=NSolve[fnprec,{x,y,z},Reals,WorkingPrecision->3 $MachinePrecision]
{x1,y1,z1}=sol[[1]]
{(1/x)-1/(y-x) - 1/(z-x)-1/(x^6)}/.sol[[1]]

Regards

Answer 1

All of the solutions to your equation are complex-valued. But you can see how close you can get to fulfilling the equation:

fr = {1/x - 1/x^6 - 1/(y - x) - 1/(z - x), 
      1/y - 1/y^6 + 1/(y - x) - 1/(z - y), 
      1/z - 1/z^6 + 1/(z - x) + 1/(z - y)};
NMinimize[fr.fr, {x, y, z}, Reals]

{2.00819*10^-57, {x -> 1.51436*10^29, y -> 8.39089*10^28, z -> -4.3816*10^28}}

As you can see, you can get arbitrarily close to a real solution by letting x, y, and z diverge towards infinity.

Answer 2

The NSolve results might not be deemed entirely trustworthy, especially if there are relatively small imaginary parts. So here is a slightly elaborate method to show that all solutions have nonzero imaginary parts.

First rewrite as rational expressions, then pull out numerators (which we will implicitly set to zero) and denominators. For the latter we make new polynomials that force them not to vanish. (Idea: if d is a denom factor, then create a new variable, rd, and a new polynomial d*rd-1.)

fr = {1/x - 1/x^6 - 1/(y - x) - 1/(z - x) == 0, 
   1/y - 1/y^6 + 1/(y - x) - 1/(z - y) == 0, 
   1/z - 1/z^6 + 1/(z - x) + 1/(z - y) == 0};
rats = Together[Apply[Subtract, fr, {1}]];
polys = Numerator[rats]
dens = Rest[FactorList[Apply[Times, Denominator[rats]]]][[All, 1]]
rvars = Array[r, Length[dens]];
rpolys = rvars*dens - 1

(* Out[212]= {-x^2 + 3 x^7 + x y - 2 x^6 y + x z - 2 x^6 z - y z + 
  x^5 y z, x y - y^2 - 2 x y^6 + 3 y^7 - x z + y z + x y^5 z - 
  2 y^6 z, -x y + x z + y z - z^2 + x y z^5 - 2 x z^6 - 2 y z^6 + 
  3 z^7}

Out[213]= {x, x - y, y, x - z, y - z, z}

Out[215]= {-1 + x r[1], -1 + (x - y) r[2], -1 + 
  y r[3], -1 + (x - z) r[4], -1 + (y - z) r[5], -1 + z r[6]} *)

Now we will eliminate all variables except x and show that most solutions in that variable are (nontrivially) complex. We handle the remaining few separately.

Timing[
 xpoly = GroebnerBasis[Join[polys, rpolys], x, Join[{y, z}, rvars], 
   MonomialOrder -> EliminationOrder]]

(* Out[216]= {9.21875, {-32 + 1344 x^5 - 28640 x^10 + 409120 x^15 - 
   4366640 x^20 + 36807216 x^25 - 252986936 x^30 + 1445929000 x^35 - 
   6955983380 x^40 + 28370424720 x^45 - 98450253550 x^50 + 
   290839419062 x^55 - 729814727767 x^60 + 1547728545736 x^65 - 
   2751022291683 x^70 + 4048160321502 x^75 - 4844726590707 x^80 + 
   4595458598910 x^85 - 3323344363479 x^90 + 1720732564224 x^95 - 
   567829946250 x^100 + 89680668750 x^105}} *)

We can make this easier by noticing that x^5 divides every power of x. So we reduce from degree 195 to 21.

x5poly = xpoly /. x -> x5^(1/5)

(* Out[217]= {-32 + 1344 x5 - 28640 x5^2 + 409120 x5^3 - 4366640 x5^4 + 
  36807216 x5^5 - 252986936 x5^6 + 1445929000 x5^7 - 
  6955983380 x5^8 + 28370424720 x5^9 - 98450253550 x5^10 + 
  290839419062 x5^11 - 729814727767 x5^12 + 1547728545736 x5^13 - 
  2751022291683 x5^14 + 4048160321502 x5^15 - 4844726590707 x5^16 + 
  4595458598910 x5^17 - 3323344363479 x5^18 + 1720732564224 x5^19 - 
  567829946250 x5^20 + 89680668750 x5^21} *)

Solutions in x for the original will be the fifth roots of solutions for x5 in the reduced degree polynomial above.

x5solns = x5 /. Solve[x5poly == 0];
N[Im[x5solns], 20]

(* Out[223]= {0, 0, 0, -0.18323622574130903849, 0.18323622574130903849, \
-0.24531128434883849207, 0.24531128434883849207, \
-0.24506750536413966115, 0.24506750536413966115, \
-0.24013452763998733202, 0.24013452763998733202, \
-0.22069093929289484601, 0.22069093929289484601, \
-0.20983222428714758263, 0.20983222428714758263, \
-0.16725390167994680894, 0.16725390167994680894, \
-0.10776268855283983645, 0.10776268855283983645, \
-0.071144606061285339602, 0.071144606061285339602} *)

So the first three are real valued, and each has a fifth root that is real valued and appears in a solution of the non-reduced polynomial. To handle these three separately we first create a Groebner basis in `{x,y,z}, eliminating the reciprocal variables.

Timing[
 gbxyz = GroebnerBasis[Join[polys, rpolys], {x, y, z}, rvars, 
    MonomialOrder -> EliminationOrder];]

(* Out[229]= {5.4375, Null} *)

Now we plug in the real valued fifth roots (the ones obtained by taking values to the (1/5) power), and solvefor {y,z}. This can be done exactly using Solve but I show with NSolve so as not to take too long. As the solution set sizes are modes, and solutions are not logarithmically large relative to the precision used, I am fairly confident that the values are accurate.

Table[NSolve[gbxyz /. x -> x5solns[[j]]^(1/5)], {j, 3}]

(* Out[233]= {{{y -> 0.689058 - 0.192267 I, 
   z -> 0.689058 + 0.192267 I}, {y -> 0.689058 + 0.192267 I, 
   z -> 0.689058 - 0.192267 I}}, {{y -> 0.205295 - 0.834545 I, 
   z -> 0.205295 + 0.834545 I}, {y -> 0.205295 + 0.834545 I, 
   z -> 0.205295 - 0.834545 I}}, {{y -> -0.677503 - 0.539908 I, 
   z -> -0.677503 + 0.539908 I}, {y -> -0.677503 + 0.539908 I, 
   z -> -0.677503 - 0.539908 I}}} *)

All are seen to have nonzero imaginary parts.