Incorrect result of Maximize

Question

Maximize is a well done command of Mathematica. However, there are spots on the Sun too. Trying in version 13.1 on Windows 10

Maximize[{Sqrt[x*y]/(z + 1) + Sqrt[y*z]/(x + 1) + Sqrt[x*z]/(y + 1), 
1/(x + 1) + 1/(y + 1) + 1/(z + 1) == 2 && x > 0 && y > 0 &&z > 0}, {x, y, z}]

, I obtain (in a dozen minutes) a warning

Maximize::wksol:Warning:there is no maximum in the region in which the objective function is defined and the constraints are satisfied;a result on the boundary will be returned.

which is OK and

Sqrt[2], {x -> -1, y -> 0, z -> -1}}

which is not correct. The change of variables {x->r^2,y->s^2,z->t^2} does not help and produces another incorrect result.

Maximize[{r*s/(t^2 + 1) + t*s/(r^2 + 1) + r*t/(s^2+1), 
1/(r^2 + 1) + 1/(s^2 + 1) + 1/(t^2 + 1) == 2 && r > 0 && s > 0 && t > 0}, {r, s, t}]

{\[Infinity], {r -> Indeterminate, s -> Indeterminate, t -> Indeterminate}}

and the same warning. The NMaximize command produces

NMaximize[{Sqrt[x*y]/(z + 1) + Sqrt[y*z]/(x + 1) + Sqrt[x*z]/(y + 1), 
1/(x + 1) + 1/(y + 1) + 1/(z + 1) == 2 && x > 0 && y > 0 &&z > 0}, {x, y, z}]

{1.41417,{x->24842.9,y->0.0000201261,z->0.0000201261}}

without any warning.

The question arises: what is a workaround?

Addition. Another example of such sort is

Minimize[{x/Sqrt[1 - x^2] + y/Sqrt[1 - y^2] + z/Sqrt[1 - z^2], 
x^2 + y^2 + z^2 == 1 && {x, y, z} >= 0 && {x, y, z} < 1}, {x, y, z}]

where Mathematica is running without any response for a long time.

Answer 1

Workaround, done in the Cloud

Maximize isn't realy good with square roots. Get rid of them.

$Version

cond1 = {a > 0, b > 0, c > 0};

f = 
 Sqrt[x*y]/(z + 1) + Sqrt[y*z]/(x + 1) + 
    Sqrt[x*z]/(y + 1) /. {x -> a^2, y -> b^2, z -> c^2} // 
  Simplify[#, cond1] &;

cond2 = {1/(x + 1) + 1/(y + 1) + 1/(z + 1) == 2} /. {x -> a^2, 
    y -> b^2, z -> c^2};

(max = Maximize[{f, Join[cond1, cond2]}, {a, b, c}]) // AbsoluteTiming

13.1.0 for Linux x86 (64-bit) (June 16, 2022)

Maximize::natt: The maximum is not attained at any point satisfying the given constraints.  >>

 {0.378328,{Sqrt[2],{a->0,b->0,c->Indeterminate}}}

Get an impression of the three equivalent maxima

ContourPlot3D[
 Sqrt[x*y]/(z + 1) + Sqrt[y*z]/(x + 1) + Sqrt[x*z]/(y + 1) == 
  Sqrt[2] - .01, {x, 0, 25}, {y, 0, 25}, {z, 0, 25}, 
 ContourStyle -> Red, FaceGrids -> All]