Edited to provide general analytical answer.
With,
f = -16 x^2 (t - 2 z - 2 z y) + (t + 2 z + 2 z y) ((-4 + t - 2 z) (4 + t - 2 z)
- 4 z^2 y^2);
the Discriminant
of this cubic in t
is
d = Discriminant[f, t] // Factor
(* 4096 (4 + 12 x^2 + 12 x^4 + 4 x^6 - 8 z^2 + 20 x^2 z^2 + x^4 z^2 - 8 y z^2 +
38 x^2 y z^2 - 8 x^4 y z^2 + y^2 z^2 + 20 x^2 y^2 z^2 - 8 x^4 y^2 z^2 + 4 z^4 +
8 y z^4 - 2 x^2 y z^4 + 2 y^2 z^4 + 2 x^2 y^2 z^4 - 2 y^3 z^4 + 8 x^2 y^3 z^4 +
4 x^2 y^4 z^4 + y^2 z^6 + 2 y^3 z^6 + y^4 z^6) *)
FindInstance
suggests that d
is non-negative everywhere,
FindInstance[d < 0, {x, y, z}, Reals]
(* {} *)
as does random sampling of points {x, y, z}
. On this basis, it seems safe to conclude the roots of f
, too lengthy to reproduce here,
rt = t /. Solve[f == 0, t]
are real, despite the fact that I
appears in the rt[[2]]
and rt[[3]]
, just as observed in the question. As a further test, Im[rt[[2]]]
for several slices in, for instance, x
.
Plot3D[Table[Im[rt[[2]]] // Chop, {x, -5, 5}], {y, -5, 5}, {z, -5, 5},
PlotRange -> {-10^-8, 10^-8}]
is zero everywhere, again substantiating that all rt
are real.
The general solution, based on 3.8.2 of Abramowitz and Stegun, is derived as follows.
a = CoefficientList[f, t]
(* {-32 z + 32 x^2 z - 32 y z + 32 x^2 y z + 8 z^3 + 8 y z^3 -
8 y^2 z^3 - 8 y^3 z^3, -16 - 16 x^2 - 4 z^2 - 8 y z^2 - 4 y^2 z^2, -2 z + 2 y z, 1} *)
q = Simplify[a[[2]]/3 - a[[3]]^2/9]
(* -(16/9) (3 + 3 x^2 + (1 + y + y^2) z^2) *)
r = Simplify[(a[[2]] a[[3]] - 3 a[[1]])/6 - a[[3]]^3/27]
(* 32/27 z (-9 x^2 (1 + 2 y) + (2 + y) (9 + (-1 - y + 2 y^2) z^2)) *)
Factor[q^3 + r^2];
The last expression is proportional to the negative of d
and, therefore, is non-positive. We can then write the key expression,
s = ComplexExpand[(r + I Sqrt[-q^3 - r^2])^(1/3), TargetFunctions -> {Re, Im}]
which, after a bit of additional manipulation becomes
(* (c1 + c2^2)^(1/6) Cos[1/3 ArcTan[c2, Sqrt[c1]]] +
I (c1 + c2^2)^(1/6) Sin[1/3 ArcTan[c2, Sqrt[c1]]] *)
where
rul = {c1 -> (4096/729 (3 + 3 x^2 + (1 + y + y^2) z^2)^3 - 1024/729 z^2 (-9 x^2 (1 + 2 y)
+ (2 + y) (9 + (-1 - y + 2 y^2) z^2))^2),
c2 -> 32/27 z (-9 x^2 (1 + 2 y) + (2 + y) (9 + (-1 - y + 2 y^2) z^2))}
The roots of f
are, then,
re = {2 s[[1]] - a[[3]]/3, -s[[1]] + I s[[2]] Sqrt[3] - a[[3]]/3,
-s[[1]] - I s[[2]] Sqrt[3] - a[[3]]/3} // Simplify
(* {2/3 (z - y z + 3 (c1 + c2^2)^(1/6) Cos[1/3 ArcTan[c2, Sqrt[c1]]]),
1/3 (-2 (-1 + y) z - 3 (c1 + c2^2)^(1/6) Cos[1/3 ArcTan[c2, Sqrt[c1]]] -
3 Sqrt[3] (c1 + c2^2)^(1/6) Sin[1/3 ArcTan[c2, Sqrt[c1]]]),
-(2/3) (-1 + y) z - (c1 + c2^2)^(1/6) Cos[1/3 ArcTan[c2, Sqrt[c1]]] +
Sqrt[3] (c1 + c2^2)^(1/6) Sin[1/3 ArcTan[c2, Sqrt[c1]]]} *)
which explicitly answers the question posed by the OP. That this is equivalent to the solutions rt
, two of which do contain I
explicitly, can be seen from, for instance,
re /. rul /. {x -> 1, y -> 0.5, z -> 1}
(* {6.81903, -6.03763, 0.218601} *)
rt /. {x -> 1, y -> 0.5, z -> 1}
(* {6.81903 - 4.44089*10^-16 I, -6.03763 - 6.66134*10^-16 I, 0.218601 + 1.33227*10^-15 I}*)
or, if exact values are desired,
re /. rul /. {x -> 1, y -> 1/2, z -> 1}
(* {1/3 (1 + 4 Sqrt[31] Cos[1/3 ArcTan[(9 Sqrt[367])/8]]),
1/3 (1 - 2 Sqrt[31] Cos[1/3 ArcTan[(9 Sqrt[367])/8]] -
2 Sqrt[93] Sin[1/3 ArcTan[(9 Sqrt[367])/8]]),
1/3 (1 - 2 Sqrt[31] Cos[1/3 ArcTan[(9 Sqrt[367])/8]] +
2 Sqrt[93] Sin[1/3 ArcTan[(9 Sqrt[367])/8]])} *)
ComplexExpand[Simplify[rt /. {x -> 1, y -> 1/2, z -> 1}]] // FullSimplify
produces the same answer.
{x,y,z}
for which the real root count switches between 1 and 3, that's a whole different matter. You'd want to read up on discriminants, I think. $\endgroup$