How to solve the following system of inequalities

edit

A bit of setup :

Vars = var('x y z')
eqns = [
    x>=0,
    0 <= y,
    0 <= z,
    x*y*z + x^2 + y^2 + z^2 <= 2*(x*y + x*z + y*z),
    2*(x^2 + y^2 + z^2) < x^2*y + x*z^2 + y^2*z - 27
]
S0=solve(eqns, [x, y, z])

The original system uses four "comprehensive" inequalities (i. e. $\leq$ or $\geq$).

TL; DR : Sage's answer separates the answers for strict inequalities and thos for equations. Hence a list of 16 implicit solutions :

sage: print("[%s]"%",\n".join(map(repr, S0)))
[[0 < x, 0 < y, 0 < z, x^2*y + y^2*z + x*z^2 - 2*x^2 - 2*y^2 - 2*z^2 - 27 > 0, -x*y*z - x^2 + 2*x*y - y^2 + 2*x*z + 2*y*z - z^2 > 0, [y == z], [y == z]],
[0 < x, 0 < y, 0 < z, x^2*y + y^2*z + x*z^2 - 2*x^2 - 2*y^2 - 2*z^2 - 27 > 0, x*y*z + x^2 - 2*x*y + y^2 - 2*x*z - 2*y*z + z^2 == 0],
[x == 0, 0 < y, 0 < z, y^2*z - 2*y^2 - 2*z^2 - 27 > 0, -y^2 + 2*y*z - z^2 > 0],
[x == 0, 0 < y, 0 < z, y^2*z - 2*y^2 - 2*z^2 - 27 > 0, y^2 - 2*y*z + z^2 == 0],
[x == 0, y == 0, 0 < z, -2*z^2 - 27 > 0, -z^2 > 0],
[x == 0, y == 0, 0 < z, -2*z^2 - 27 > 0, z^2 == 0],
[x == 0, y == 0, z == 0, -27 > 0, 0 == 0],
[x == 0, y == 0, z == 0, -27 > 0, 0 > 0],
[x == 0, z == 0, 0 < y, -2*y^2 - 27 > 0, -y^2 > 0],
[x == 0, z == 0, 0 < y, -2*y^2 - 27 > 0, y^2 == 0],
[y == 0, 0 < x, 0 < z, x*z^2 - 2*x^2 - 2*z^2 - 27 > 0, -x^2 + 2*x*z - z^2 > 0],
[y == 0, 0 < x, 0 < z, x*z^2 - 2*x^2 - 2*z^2 - 27 > 0, x^2 - 2*x*z + z^2 == 0],
[y == 0, z == 0, 0 < x, -2*x^2 - 27 > 0, -x^2 > 0],
[y == 0, z == 0, 0 < x, -2*x^2 - 27 > 0, x^2 == 0],
[z == 0, 0 < x, 0 < y, x^2*y - 2*x^2 - 2*y^2 - 27 > 0, -x^2 + 2*x*y - y^2 > 0],
[z == 0, 0 < x, 0 < y, x^2*y - 2*x^2 - 2*y^2 - 27 > 0, x^2 - 2*x*y + y^2 == 0]]

which have to be reworked to get explicit solutions.

Let's examine them :

• $ \left[0 < x, 0 < y, 0 < z, x^{2} y + y^{2} z + x z^{2} - 2 \, x^{2} - 2 \, y^{2} - 2 \, z^{2} - 27 > 0, -x y z - x^{2} + 2 \, x y - y^{2} + 2 \, x z + 2 \, y z - z^{2} > 0\right] $

That's a restatement of the original system in the special case of all strict inequalities ; this means that Sage does not find a solution in this case.

• $ \left[0 < x, 0 < y, 0 < z, x^{2} y + y^{2} z + x z^{2} - 2 \, x^{2} - 2 \, y^{2} - 2 \, z^{2} - 27 > 0, x y z + x^{2} - 2 \, x y + y^{2} - 2 \, x z - 2 \, y z + z^{2} = 0\right] $

This is a restatement oif the original system where the last inequality is an equation. Again Sage finds no solution.

• $ \left[x = 0, 0 < y, 0 < z, y^{2} z - 2 \, y^{2} - 2 \, z^{2} - 27 > 0, -y^{2} + 2 \, y z - z^{2} > 0\right] $

In this case, $x$ is fixed at 0, the other inequalities are strict. $x$ is eliminated from the other inequalities. It is easy to see (by factorization) that the last inequality, is incompatible wit real values of $y$ and $z$ :

sage: S2=(S1:=S0[2][4]).operator()(*map(factor, S1.operands())) ; S2
-(y - z)^2 > 0
sage: solve(S2, y)
[]

which Sage can also state directly :

sage: S0[2][4].solve(y)
[]

• $\left[x = 0, 0 < y, 0 < z, y^{2} z - 2 \, y^{2} - 2 \, z^{2} - 27 > 0, y^{2} - 2 \, y z + z^{2} = 0\right]$

When the last inequality is replaced by an equation, the latter has a solution :

sage: S0[3][4].solve(y)
[y == z]

which can be added to the constraint $x=0$, which is an element of solution :

sage: S3=[S0[3][0]]
sage: S3.extend(S0[3][4].solve(y)) ; S3
[x == 0, y == z]

This partial solution can be substituted in S0[3], an we can eliminate the elements known true. The resultant system is :

sage: S31=[v for v in [u.subs(S3) for u in S0[3]] if not v] ; S31
[0 < z, 0 < z, z^3 - 4*z^2 - 27 > 0]

which, alas, Sage does not solve. But we can solve the second inequation and check the first two:

sage: S32=S31[2].solve(y)[0] ; S32
[z > 5.056147144240078]
sage: bool(S31[0].subs(z=S32[0].rhs()))
True

We can add this inequality to our solution :

sage: S3.append(S32[0]) ; S3 [x == 0, y == z, z > 5.056147144240078]

One note that the lower bound of the $z$ part of the soution is a numerical approximation. An exact value (and its radical expression) can be obtained as :

 sage: (MagicRoot:=S0[3][3].subs(S3[:2]).lhs().roots(ring=AA, multiplicities=False)[0], MagicRoot.radical_expression())

$ \left(5.056147018090691?, \frac{1}{3} \, \left(\frac{1}{2}\right)^{\frac{1}{3}} {\left(27 \, \sqrt{985} + 857\right)}^{\frac{1}{3}} + \frac{32 \, \left(\frac{1}{2}\right)^{\frac{2}{3}}}{3 \, {\left(27 \, \sqrt{985} + 857\right)}^{\frac{1}{3}}} + \frac{4}{3}\right) $

Further implicit solutions can be treated in similar fashion :

• Solution 4 has no solution, since no real has a negative square (last inequation).

• Similarly, Solution 5 includes both O<z and="" z^2="=0," which="" are="" mutually="" incompatible="" ;="" no="" solution.<="" p="">

• Solution 6 reduces to -27>0, which is false. No solution. Ditto for Solution 7, which has also other antilogies...

• Similarly, the last two elements of both Solution 8 and Solution 8 are antilogies for $y\in\mathbb{R}$.

• Solution 10 has no solutions for the same reasons as for Solution 9 up to a permutation of variables.

• Solution 11 is identical to Solution 3 up to a variable permutation, and admits a similar solution :

$ \left[y = 0, z = x, x > 5.056147018090691?\right] $

• Solutions 12 and 13 have the same antilogies as Solutions 6, 7, 8 and 9.

• Solution 14 has no solution, for the same reasons as Solutions 2 and 10 up to a permutatin of variables

• Solution 15 is analogous to Solutions 3 and 10 up to a permutation of variables, and admits a similar solution :

$ \left[z = 0, y = x, x > 5.056147018090691?\right] $

Overall, Sage finds three (infinite) solutions :

$$ \left[\left[x = 0, y > 5.056147018090691?, z = y\right], \left[y = 0, z = x, x > 5.056147018090691?\right], \left[z = 0, y = x, x > 5.056147018090691?\right]\right] $$

This happens to appear similar to the solution found by the Wolfram engine :

sage: try:
....:     MSols=mathematica(eqns).Solve(Vars)
....: except:
....:     try:
....:         mathematica("Sol=Solve[%s, %s]"%tuple(map(lambda u:repr(mathematica(u)), (eqns, (x, y, z)))))
....:     except:
....:         pass
....:     MSols=mathematica("Sol")
....: 
sage: MSols
{{x -> ConditionalExpression[0, y > Root[-27 - 4*#1^2 + #1^3 & , 1, 0]], 
  z -> ConditionalExpression[y, y > Root[-27 - 4*#1^2 + #1^3 & , 1, 0]]}, 
 {y -> ConditionalExpression[0, x > Root[-27 - 4*#1^2 + #1^3 & , 1, 0]], 
  z -> ConditionalExpression[x, x > Root[-27 - 4*#1^2 + #1^3 & , 1, 0]]}, 
 {y -> ConditionalExpression[x, x > Root[-27 - 4*#1^2 + #1^3 & , 1, 0]], 
  z -> ConditionalExpression[0, x > Root[-27 - 4*#1^2 + #1^3 & , 1, 0]]}}

[ One notes that Sage intercepts as an error a warning emitted by the Wolfram engine, hence the workaround...)

Fricas gives answers similar to those of Sage ; Giac doesn't solve the system at all ; Sympy can solve only for one variable at a time.

It must be noted that the three infinite solutions, unidimensional varieties, are totally missed by implicit_plot_3d when plotting the envelopes of the inequalities... I'm not (yet) sure of the (in-)correctness of this behaviour.

I'd be interested by Maple's answer...

EDIT : I forgot to add that Sage's solutions check.

HTH,