# find roots with specified coefficients

Please assume we have an equation which we want to obtain its roots. We can use

Roots[-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) == 0, t]


in order to access its roots. We encounter to answers as

-(2/3) (-z +
y z) + ((1 + I Sqrt) (-48 - 48 x^2 - 16 z^2 - 16 y z^2 -
16 y^2 z^2))/(3 2^(
2/3) (1152 z - 576 x^2 z + 576 y z - 1152 x^2 y z - 128 z^3 -
192 y z^3 + 192 y^2 z^3 + 128 y^3 z^3 + Sqrt[
4 (-48 - 48 x^2 - 16 z^2 - 16 y z^2 - 16 y^2 z^2)^3 + (1152 z -
576 x^2 z + 576 y z - 1152 x^2 y z - 128 z^3 - 192 y z^3 +
192 y^2 z^3 + 128 y^3 z^3)^2])^(1/3)) - (1/(
6 2^(1/3)))(1 - I Sqrt) (1152 z - 576 x^2 z + 576 y z -
1152 x^2 y z - 128 z^3 - 192 y z^3 + 192 y^2 z^3 + 128 y^3 z^3 +
Sqrt[4 (-48 - 48 x^2 - 16 z^2 - 16 y z^2 -
16 y^2 z^2)^3 + (1152 z - 576 x^2 z + 576 y z -
1152 x^2 y z - 128 z^3 - 192 y z^3 + 192 y^2 z^3 +
128 y^3 z^3)^2])^(1/3)


in which there is a term as (1 + I Sqrt) and (1 - I Sqrt). But when we use

Plot[-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) /. {x -> 1,
y -> 0.5, z -> 1}, {t, -10, 10}]


we see there are three roots.

It means the imaginary parts in roots are not important albeit {x -> 1, y -> 0.5, z -> 1}. In my case x, y and z are in [0,1]. For this reason we can we have three real roots.

With the condition which is governed on the coefficients (x, y and z) that must be in [0,1], how can we obtain the real configuration of the roots of the main equation? ( for example an equation as a root without an imaginary part)

• If you want roots extressed as radicals, that cannot be done. Check this link for further detail. If you are looking for conditions on {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. – Daniel Lichtblau Aug 16 '16 at 14:32
• You can't generally express the roots of a depressed cubic via radicals without using complex numbers. That's mathematics, sorry. – John Doty Aug 16 '16 at 14:33
• @ Daniel,No I don't want to find roots in a special interval. the coefficient of polynomial is in an intervals that roots are real and imaginary parts are eliminated. – Unbelievable Aug 16 '16 at 14:42

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[] and rt[], just as observed in the question. As a further test, Im[rt[]] for several slices in, for instance, x.

Plot3D[Table[Im[rt[]] // 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[]/3 - a[]^2/9]
(* -(16/9) (3 + 3 x^2 + (1 + y + y^2) z^2) *)
r = Simplify[(a[] a[] - 3 a[])/6 - a[]^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[] - a[]/3, -s[] + I s[] Sqrt - a[]/3,
-s[] - I s[] Sqrt - a[]/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 (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 (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 Cos[1/3 ArcTan[(9 Sqrt)/8]]),
1/3 (1 - 2 Sqrt Cos[1/3 ArcTan[(9 Sqrt)/8]] -
2 Sqrt Sin[1/3 ArcTan[(9 Sqrt)/8]]),
1/3 (1 - 2 Sqrt Cos[1/3 ArcTan[(9 Sqrt)/8]] +
2 Sqrt Sin[1/3 ArcTan[(9 Sqrt)/8]])} *)

ComplexExpand[Simplify[rt /. {x -> 1, y -> 1/2, z -> 1}]] // FullSimplify


• Your explanation clarifies more things. – Unbelievable Aug 16 '16 at 18:02
roots = t /. {Roots[-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) == 0, t] //
ToRules} // Simplify;

roots // LeafCount

(*  861  *)


To get a representation without I use ComplexExpand

roots2 = ComplexExpand[roots, TargetFunctions -> {Re, Im}] // Simplify;

roots2 // LeafCount

(*  25973  *)


The first representation using I is much less complicated.

ex1 = roots /. {x -> 1, y -> 1/2, z -> 1} // FullSimplify

(*  {Root[9 - 41 #1 - #1^2 + #1^3 &, 3],
Root[9 - 41 #1 - #1^2 + #1^3 &, 1],
Root[9 - 41 #1 - #1^2 + #1^3 &, 2]}  *)

ans1 = ex1 // N

(*  {6.81903, -6.03763, 0.218601}  *)


Converting the Root expressions

ex1R = ex1 // ToRadicals // ComplexExpand // Simplify

(*  {1/3 (1 + 4 Sqrt Cos[1/3 ArcTan[(9 Sqrt)/8]]),
1/3 (1 - 2 Sqrt Cos[1/3 ArcTan[(9 Sqrt)/8]] -
2 Sqrt Sin[1/3 ArcTan[(9 Sqrt)/8]]),
1/3 (1 - 2 Sqrt Cos[1/3 ArcTan[(9 Sqrt)/8]] +
2 Sqrt Sin[1/3 ArcTan[(9 Sqrt)/8]])}  *)


For the second form

ex2 = roots2 /. {x -> 1, y -> 1/2, z -> 1} // FullSimplify

(*  {(1/3)*(1 + 4*Sqrt*
Cos[(1/3)*ArcTan[(9*Sqrt)/
8]]), (1/3)*
(1 - 2*Sqrt*
Cos[(1/3)*ArcTan[(9*Sqrt)/
8]] - 2*Sqrt*
Sin[(1/3)*ArcTan[(9*Sqrt)/
8]]), (1/3)*
(1 - 2*Sqrt*
Cos[(1/3)*ArcTan[(9*Sqrt)/
8]] + 2*Sqrt*
Sin[(1/3)*ArcTan[(9*Sqrt)/
8]])}  *)


This is identical to earlier result

ex2 === ex1R

(*  True  *)


Demonstrating that the original Root expression is equivalent

(ex2 // N) == ans1

(*  True  *)

• It is so amazing! – Unbelievable Aug 16 '16 at 18:01