# How to solve the following system of inequalities

I came across this system of inequalities in MaplePrimes, where a user complained about Maple taking too long to compute. I would like to see how SageMath performs in this case.

x, y, z = 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
]

solve(eqns, [x, y, z])


SageMath provides some solutions, but it seems that the results are not as satisfactory.

[[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], [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]]


For example, in the last solution, it appears that x and y are not independent of each other. Many solutions seem to be outputted in a similar form without much simplification.

Here's another example. When we look at the solution obtained by Mathematica, it seems that SageMath is not performing as well. I'm not sure if SageMath has more accurate commands for solving such systems of inequalities.

Reduce[{x^2 + y^2 > 1, x^4 - y^4 < 1}, {x, y}]
(*output: (x < -1 && (y < -(-1 + x^4)^(1/4) || y > (-1 + x^4)^(1/4))) || (-1 <=
x <= 1 && (y < -Sqrt[1 - x^2] || y > Sqrt[1 - x^2])) || (x >
1 && (y < -(-1 + x^4)^(1/4) || y > (-1 + x^4)^(1/4)))*)

x, y = var('x y')
inequalities = [
x^2 + y^2 > 1, x^4 - y^4 < 1
]

solution = solve(inequalities, x, y)
print(solution)

# output: [
[-x^4 + y^4 + 1 > 0, x^2 + y^2 - 1 > 0]
]

edit retag close merge delete

1

Potentially QEPCAD should be able to solve such systems.

( 2023-06-29 03:48:08 +0200 )edit

Sort by ยป oldest newest most voted

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,

more

In Maple, the best option is to use RAGlib package: https://www-polsys.lip6.fr/~safey/RAG...

( 2023-07-01 00:16:12 +0200 )edit