$x,y,z \geqslant 0$, $x+y^2+z^3=1$, prove $x^2y+y^2z+z^2x < \frac12$

The standard way of solving the problem on a conditional extremum is the method of Lagrange multipliers, which reduces it to a system of equations.

The greatest value of function $$f(x,y,z,\lambda) = x^2y+y^2z+z^2x+\lambda(x+y^2+z^3-1)$$ on the interval $$x,y,z\in[0,1]$$ is reached or at its edges, or in the inner stationary point.

$\color{brown}{\textbf{Inner stationary points.}}$

The inner stationary points has zero partial derivatives $$\begin{cases} f'_\lambda = x + y^2 + z^3 - 1 = 0\\ f'_x = z^2 + 2xy + \lambda = 0\\ f'_y = x^2 + 2yz + 2\lambda y = 0\\ f'_z = y^2 + 2zx + 3\lambda z^2 = 0. \end{cases}$$ After the excluding of parameter $\lambda$ get the system $$\begin{cases} x + y^2 + z^3 - 1 = 0\\ x^2 + 2yz = 2y(z^2 + 2xy)\\ y^2 + 2zx = 3z^2(z^2 + 2xy), \end{cases}$$ or $$\begin{cases} x + y^2 + z^3 - 1 = 0\\ (1-3y^2-z^3)^2-4y^4+2yz(1-z)=0\\ 2z(1-3yz)(1-y^2-z^3)+y^2-3z^4=0. \end{cases}$$

Using of Groebner basis allows to get the positive solutions $$ \genfrac{[}{.}{0}{0}{x\approx 0.16367,\quad y\approx 0.761982,\quad z\approx 0.634724,\quad f\approx 0.454882} {x\approx 0.558113,\quad y\approx 0.498412,\quad z\approx 0.578371,\quad f\approx 0.485622}. $$ Alternative way is shown below.

$\color{brown}{\textbf{The edges.}}$

The edges of the field are achieved when $x = 0$, $y = 0$ or $z = 0$.

Substitution $x = 0$ in the expressions for the partial derivatives leads to the system $$ \begin{cases} x=0\\ y^2+z^3 = 1\\ 2yz+2\lambda y = 0\\ y^2+3\lambda z^2 = 0 \end{cases} $$ with solutions $$ \genfrac{[}{.}{0}{} {x=0,\quad y=\sqrt{0.75},\quad z=\sqrt[3]{0.25}\approx 0.629991,\quad f\approx 0.47247} {x=0,\quad y=0,\quad z=1,\quad f=0}$$

Substitution $y = 0$ in the expressions for the partial derivative leads to the system $$ \begin{cases} y=0\\ x+z^3=1\\ z^2+\lambda = 0\\ 2zx+3\lambda z^2 = 0 \end{cases} $$ with solution $$x=0.6,\quad y=0,\quad z=\sqrt[3]{0.4}\approx 0.736806,\quad f\approx 0.32573.$$

Substitution $z = 0$ in the expressions for the partial derivative leads to the system $$ \begin{cases} z=0\\ x+y^2=1\\ 2xy+\lambda = 0\\ x^2+2\lambda y = 0 \end{cases} $$ with solutions $$ \genfrac{[}{.}{0}{} {x=0.8,\quad y=\sqrt{0.2}\approx 0.447214,\quad z=0,\quad f\approx 0.286217} {x=0,\quad y=0,\quad z=1,\quad f=0} $$

The values of function at the vertices of the area (unit parallelepiped) equal to zero.

So the greatest value approximately equals to $0.485622$. Given accuracy of calculations provides the inequality $$\boxed{x^2y+y^2z+z^2x<1/2.}$$

$\color{brown}{\textbf{System resolving, alternative way.}}$

Set condition allows us to reduce the problem to finding unconditional extremes of $$f(y,z)=(1-y^2-z^3)^2y+y^2z+z^2(1-y^2-z^3).$$ Necessary optimality conditions in the field have the form: $$\begin{cases} f'_y = (1-y^2-z^3)^2-4y^2(1-y^2-z^3)+2yz-2yz^2 = 0\\ f'_z = -6yz^2(1-y^2-z^3)+y^2+2z(1-y^2-z^3)-3z^4 = 0, \end{cases}$$ or $$\begin{cases} 5y^4+y^2(6z^3-6)+y(2z-2z^2)+(z^3-1)^2 = 0\\ 6y^3z^2+y^2(1-2z)+6y(z^3-1)z^2-5z^4+2z = 0. \end{cases}$$

If to consider the coefficient of the highest power of $y$ as denominator in equation, and the remaining coefficients - numerators, we get the above equations. Then we can subtract the second equation factor $y$ from the first and repeat subtraction with factor $1$, obtaining the system

$$\begin{cases} C_{2,2}(z)y^2 + C_{2,1}(z)y + C_{2,0}(z) = 0\\ 6y^3z^2+y^2(1-2z)+6y(z^3-1)z^2-5z^4+2z = 0, \end{cases}$$ where $$C_{2,2}(z) = 36z^7-36z^4+20z^2-20z+5,$$ $$C_{2,1}(z) = 18z^6+102z^5-30z^2,$$ $$C_{2,0}(z) = 36z^{10}-72z^7+50z^5+11z^4-20z^2+10z.$$

Thus, the order of the first equation in $y$ reduced from fourth to second. Likewise, lowering the order of the second equation by a first equation, we obtain

$$\begin{cases} C_{2,2}(z)y^2 + C_{2,1}(z)y + C_{2,0}(z) = 0\\ D_{2,2}(z)y^2 + D_{2,1}(z)y + D_{2,0}(z) = 0,\\ \end{cases}\qquad(1)$$ where $$D_{2,2}(z) = 180z^8+576z^7-72z^5-144z^4+40z^3-60z^2+30z-5,$$ $$D_{2,1}(z) = 180z^7-30z^6-30z^5-60z^3+30z^2,$$ $$D_{2,0}(z) = 180z^{11}-252z^8+100z^6-28z^5+25z^4-40z^3+40z^2-10z.$$

The system $(1)$ is a linear in the unknowns $y^2$ and $y$, so $$y^2=\dfrac{\Delta_2(z)}{\Delta_0(z)},\quad y=\dfrac{\Delta_1(z)}{\Delta_0(z)},\qquad(2)$$ where $$\Delta_0(z) = C_{2,2}(z)D_{2,1}(z) - C_{2,1}(z)D_{2,2}(z),$$ $$\Delta_2(z) = -C_{2,0}(z)D_{2,1}(z) + C_{2,1}(z)D_{2,0}(z),$$ $$\Delta_1(z) = -C_{2,2}(z)D_{2,0}(z) + C_{2,0}(z)D_{2,2}(z),$$ or $$\Delta_0(z) = 90z^{10}-828z^9-1662z^8-144z^7+396z^6+1028z^5-200z^4+180z^3-220z^2+10z,$$ $$\Delta_2(z) = -90z^{13}+540z^{12}+30z^{11}+234z^{10}-864z^9-290z^8+256z^7+74z^6+220z^5-160z^4+100z^3-50z^2,$$ $$\Delta_1(z) = 576z^{13}-1296z^{10}+90z^9+815z^8+444z^7-5z^6-580z^5+126z^4+10z^3+80z^2-30z-5.$$

From $(2)$ for $\Delta_0(z)\not=0$ should be $$\Delta_1^2(z) - \Delta_2(z)\Delta_0(z) = 0,$$ or $$331776z^{26}-1484892z^{23}-19440z^{22}+1233720z^{21}+3079404z^{20}+195732z^{19}-3189924z^{18}-3109428z^{17}+368233z^{16}+2734116z^{15}+1243978z^{14}-741000z^{13}-805907z^{12}+75696z^{11}+164040z^{10}+82560z^9-172194z^8+53440z^7-10290z^6+32440z^5-7460z^4-4400z^3+100z^2+300z+25 = 0.$$

The coefficients are particially calculated using the Mathcad package, and $\mathcal{polyroots}()$ function is also used, which calculates all the roots of the polynomial by the "accompanying matrix" method.

Calculating values $y$ with $(2)$ and $x,f$ by the formula $$x = 1-y^2-z^3,\quad f=xy^2+yz^2+zx^2$$ and checking them by substituting in the original system, we obtain the following stationary points with the positive coordinates: $$ (x,y,z,f)\in\left[\genfrac{}{}{0pt}{0} {(0.5581125,\ 0.4984120,\ 0.5783713,\ 0.4856221)} {(0.1636702,\ 0.7619816,\ 0.6347238,\ 0.4548812)} \right. $$

My new proof (2020/12/06)

Since $y^2 + \frac{1}{4} \ge y$, it suffices to prove that $$x^2(y^2 + \tfrac{1}{4}) + y^2 z + z^2 x < \frac{1}{2}.$$ Since $y^2 = 1-x-z^3$, it suffices to prove that $$4x^3+(4z^3-5)x^2+(-4z^2+4z)x+4z^4-4z+2 > 0. \tag{1}$$

If $x = 0$, (1) is true since $4z^4-4z+2=(2z^2-1)^2+(2z-1)^2$. (due to @Yuri Negometyanov)

If $x > 0$, by noting that \begin{align} &4x^3+(4z^3-5)x^2+(-4z^2+4z)x+4z^4-4z+2\\ =\ & 4x\left(x + \frac{1}{2}z^3 - \frac{11}{16}\right)^2 + \frac{1}{2}\left(x - 2\sqrt{2z^4 - 2z + 1}\right)^2\\ &\quad + x\left[2\sqrt{2z^4 - 2z + 1} - \left(z^6 - \frac{11}{4}z^3 + 4z^2 - 4z + \frac{121}{64}\right)\right], \tag{2} \end{align} it suffices to prove that, for all $z \in [0, 1]$, $$2\sqrt{2z^4 - 2z + 1} - \left(z^6 - \frac{11}{4}z^3 + 4z^2 - 4z + \frac{121}{64}\right) > 0. \tag{3}$$ Since $2z^4 - 2z + 1 - (1 - \frac{3}{4}\sqrt[3]{2}) = \frac{1}{8}(3\sqrt[3]{4} + 4\sqrt[3]{2} z + 4z^2)(2z - \sqrt[3]{2})^2$, we have $\sqrt{2z^4 - 2z + 1} \ge \sqrt{1 - \frac{3}{4}\sqrt[3]{2}} > \frac{19}{81}$. We have $$\sqrt{2z^4 - 2z + 1} = 1 - \frac{2z - 2z^4}{1 + \sqrt{2z^4 - 2z + 1}} \ge 1 - \frac{2z - 2z^4}{1 + \frac{19}{81}}.$$ It suffices to prove that, for all $z \in [0, 1]$, $$2 - \frac{4z - 4z^4}{1 + \frac{19}{81}} - \left(z^6 - \frac{11}{4}z^3 + 4z^2 - 4z + \frac{121}{64}\right) > 0$$ or $$-1600z^6 + 5184z^4 + 4400z^3 - 6400z^2 + 1216z + 175 > 0.$$ We have \begin{align} &-1600z^6 + 5184z^4 + 4400z^3 - 6400z^2 + 1216z + 175\\ \ge\ & -1600z^6 + 5184z^4 + 4400z^3 - 6400z^2 + 1216z + 175 - 1600(1 - z^2)(z^2-\tfrac{1}{5})^2\\ =\ & 2944z^4 + 4400z^3 - 5696z^2 + 1216z + 111\\ >\ & 2\sqrt{2944z^4 \cdot 111} + 2\sqrt{4400z^3\cdot 1216z} - 5696z^2\\ =\ & (2\sqrt{2944\cdot 111} + 2\sqrt{4400\cdot 1216} - 5696)z^2\\ \ge\ & 0. \end{align} We are done.

Remark: Actually, (1) is true for all $x\ge 0, z\ge 0$.

This answer is incomplete.

Let $A=x+y^2$, $B=y^2+z^3$ and $C=z^3+x$.

Claim: $x^2y \leq \dfrac{A}{2}x^{3/2}$, $y^2z \leq \dfrac{B^{2/3}}{4^{2/3}}y^{4/3}$, $z^2x \leq \dfrac{C^{4/3}}{4^{2/3}}x^{1/3}$.

Assuming the claim, we note that $A+B+C=2$ and we need to prove that the given sum is less than $\dfrac{A+B+C}{4}$.

The method to be employed here is exactly the same as in a previous answer:

• Olympiad inequality $\sum_{cyc} \frac{x^4}{8x^3+5y^3} \geqslant \frac{x+y+z}{13}$

Let $u = x$ , $v = y^2$ , $w = z^3$ , then $u,v,w \ge 0$ , $u+v+w=1$ and the inequality to be established: $$ x^2y+y^2z+z^2x < \frac12 \quad \Longrightarrow \quad f(u,v,w) = u^2v^{1/2}+v\,w^{1/3}+w^{2/3}u < \frac12 $$ The maximum of this function inside the equilateral triangle must shown to be less than $1/2$.
There is no symmetry argument, because the latter is effectively destroyed by the "weird" condition $x+y^2+z^3=1$ , as it is called. Another proof without words is attempted by plotting a contour map of the function, as depicted. Levels (nivo) of these isolines are defined (in Delphi Pascal) as:

nivo := min + sqrt(g/grens)*(max-min); { sqrt = square root ; grens = 25 ; g = 0..grens }

The darkness of the isolines is proportional to the (positive) function values; they are almost black near the maximum and almost white near the minimum values. Maximum and minimum values of the function are observed to be:

 4.58251457205350E-0003 < f < 4.85621276951755E-0001 < 1/2

The little $\color{blue}{\mbox{blue}}$ spot is where $\left|f(u,v,w) - \mbox{max}\right|< 0.0002$ . This maximum is close to values found by other people here, but not quite. Perhaps it's interesting to know the position of the maximum as well:

(x,y,z) = ( 5.58304528246164E-0001 , 4.97693736095187E-0001 , 5.78892473099889E-0001 )

I have no idea how to convert these numerical values into something more analytical.

EDIT. What you see is what you get :-) Without doubt. Some decent error analysis reveals that the value of the maximum to be trusted in this answer is : $\;0.48562 \pm 0.00003\;$ , quite in agreement with values given elsewhere (e.g. in the comment by Steve Kass).