Maximize finds solution outside the constraint

Question

I'm trying to maximize a (quite simple) polynomial inside a sphere. The command is simply:

f[x1_, x2_, x3_] := -(x1/36) + x1^3/67 + x2/25 - (x1^2  x2)/12 + (x1  x2^2)/15 - x2^3/30 + 
   x3/27 - (2  x1^2  x3)/11 + (x1  x2  x3)/35 - (x2^2  x3)/11 + 
   x1^2  x2^2  x3 + (x1  x3^2)/22 - (x2  x3^2)/18 - x3^3/49;
Maximize[{f[x1, x2, x3], x1^2 + x2^2 + x3^2 <= 1}, {x1, x2, x3}] // N

The result I get is:

(* {-149457., {x1 -> -0.873639, x2 -> 64.5, x3 -> -48.}} *)

which grossly violates the constraint x1^2 + x2^2 + x3^2 <= 1, and even the function's value is incoherent, since it has positive values both within the constraint and outside.

To me, this seems like a textbook usage of Maximize. Is this a Mathematica bug?

Answer 1

If you drop N, you will see that the result is a very long and complex expression, full of high-degree Root objects. I would say that the result isn't really a bug, but shows an accumulation of numerical errors, which you can avoid by using arbitrary-precision calculations (set by the second argument of N), for example,

Maximize[{f[x1, x2, x3], x1^2 + x2^2 + x3^2 <= 1}, {x1, x2, x3}] // N[#, 5]&
(* {0.057092, {x1 -> -0.87364, x2 -> -0.026513, x3 -> -0.48585}} *)

Alternatively, you can use numerical maximization with NMaximize:

NMaximize[{f[x1, x2, x3], x1^2 + x2^2 + x3^2 <= 1}, {x1, x2, x3}]
(* {0.0570927, {x1 -> -0.873678, x2 -> -0.0266629, x3 -> -0.485778}} *)