# How to extract only real solutions from the result of Solve

For Example, I have the following polynomial equations:

$$\left\{ \begin{array}{ll} \dfrac{2025 (208 x+5 (y-3446))}{52 (y+90)^2}+\dfrac{300 (8 x-21 y+1920) (y-80)}{7 (x+30)^3}+\frac{300 (80 x-21 z+4500) (z-100)}{7 (x+30)^3}-\dfrac{1521 (-539 x+50 z+15425)}{539 (z+39)^2} =0\\[30pt] \frac{2025 (445 y-208 z+49410) (z-45)}{52 (y+90)^3}-\dfrac{2025 (x-85) (208 x+5 (y-3446))}{52 (y+90)^3}-\frac{300 (8 x-21 y+1920)}{7 (x+30)^2}-\frac{1521 (-539 y+90 z+19680)}{539 (z+39)^2} =0\\[30pt] \frac{2025 (445 y-208 z+49410)}{52 (y+90)^2}+\frac{300 (80 x-21 z+4500)}{7 (x+30)^2}+\frac{1521 (x-25) (539 x-25 (2 z+617))}{539 (z+39)^3}+\frac{1521 (y-30) (539 y-30 (3 z+656))}{539 (z+39)^3}=0 \\ \end{array} \right.$$

The mathematica expression:

{{(300 (1920 + 8 x - 21 y) (-80 + y))/(7 (30 + x)^3) + (
2025 (208 x + 5 (-3446 + y)))/(52 (90 + y)^2) + (
300 (4500 + 80 x - 21 z) (-100 + z))/(7 (30 + x)^3) - (
1521 (15425 - 539 x + 50 z))/(
539 (39 + z)^2)}, {-((300 (1920 + 8 x - 21 y))/(7 (30 + x)^2)) - (
2025 (-85 + x) (208 x + 5 (-3446 + y)))/(52 (90 + y)^3) + (
2025 (49410 + 445 y - 208 z) (-45 + z))/(52 (90 + y)^3) - (
1521 (19680 - 539 y + 90 z))/(
539 (39 + z)^2)}, {-((2025 (49410 + 445 y - 208 z))/(
52 (90 + y)^2)) - (300 (4500 + 80 x - 21 z))/(7 (30 + x)^2) - (
1521 (-25 + x) (539 x - 25 (617 + 2 z)))/(539 (39 + z)^3) - (
1521 (-30 + y) (539 y - 30 (656 + 3 z)))/(539 (39 + z)^3)}}


How to extract only those real solutions as new Rule? I mean a general approach to such solutions with mixed real and complex solutions in a rule.

• For all Solve, Reduce, NSolve you can specify a domain for your unknowns, Reals in your case. Commented Aug 4, 2014 at 11:23

eqns = {(300 (1920 + 8 x - 21 y) (-80 + y))/(7 (30 + x)^3) + (2025 (208 x +
5 (-3446 + y)))/(52 (90 + y)^2) + (300 (4500 + 80 x - 21 z) (-100 +
z))/(7 (30 + x)^3) - (1521 (15425 - 539 x +
50 z))/(539 (39 + z)^2) ==
0, -((300 (1920 + 8 x - 21 y))/(7 (30 + x)^2)) - (2025 (-85 +
x) (208 x + 5 (-3446 + y)))/(52 (90 + y)^3) + (2025 (49410 +
445 y - 208 z) (-45 + z))/(52 (90 + y)^3) - (1521 (19680 - 539 y +
90 z))/(539 (39 + z)^2) ==
0, -((2025 (49410 + 445 y - 208 z))/(52 (90 + y)^2)) - (300 (4500 +
80 x - 21 z))/(7 (30 + x)^2) - (1521 (-25 + x) (539 x -
25 (617 + 2 z)))/(539 (39 + z)^3) - (1521 (-30 + y) (539 y -
30 (656 + 3 z)))/(539 (39 + z)^3) == 0} // Simplify;

vars = Variables[Level[eqns, {-1}]];

allSolns = Solve[eqns, vars] // Simplify;

Length[allSolns]


47

Since the exact results contain Root objects they are neither explicitly Real nor explicitly Complex.

ro = {z /. allSolns[[47]]};

{Cases[ro, _Real], Cases[ro, _Complex]}


{{}, {}}

The root objects need to be converted to numbers to be explicitly Complex or Real. Machine precision calculations are too inaccurate when converting the root objects. Specifying a precision will require Mathematica to maintain sufficient precision in the calculations.

realSolns = Cases[allSolns,
_?(FreeQ[N[#, 16], Complex] &)];

Length[realSolns]


7

realSolns // N[#, 16] & // Grid[#, Alignment -> Left] &


Version 9. Specifying a domain may speed things up or may greatly slow things down.

Cases[N[Solve[{...}, {x,y,z}], 30], {x->_Real, y->_Real, z->_Real}] Result in 28 seconds.

N[Solve[{...}, {x,y,z}, Reals], 30] Stopped it after 6 hours with no result.


Setting expr equal to the OP's expression, the equations are given by Flatten@expr == 0.

expr = {{(300 (1920+8 x-21 y) (-80+y))/(7 (30+x)^3)+(2025 (208 x+5 (-3446+y)))/(52 (90+y)^2)+(300 (4500+80 x-21 z) (-100+z))/(7 (30+x)^3)-(1521 (15425-539 x+50 z))/(539 (39+z)^2)},{-((300 (1920+8 x-21 y))/(7 (30+x)^2))-(2025 (-85+x) (208 x+5 (-3446+y)))/(52 (90+y)^3)+(2025 (49410+445 y-208 z) (-45+z))/(52 (90+y)^3)-(1521 (19680-539 y+90 z))/(539 (39+z)^2)},{-((2025 (49410+445 y-208 z))/(52 (90+y)^2))-(300 (4500+80 x-21 z))/(7 (30+x)^2)-(1521 (-25+x) (539 x-25 (617+2 z)))/(539 (39+z)^3)-(1521 (-30+y) (539 y-30 (656+3 z)))/(539 (39+z)^3)}};
eqns = Flatten@expr == 0;


Somewhat unexpectedly the quickest way to get solutions was with NSolve and setting WorkingPrecision. Using the default MachinePrecision was slower.

nrealsol = NSolve[
eqns, {x, y, z}, Reals,
WorkingPrecision -> 16]; // AbsoluteTiming
NSolve[
eqns, {x, y, z}, Reals]; // AbsoluteTiming
(*
{2.780016, Null}
{5.659495, Null}
*)


Both were faster than Solve without a domain specification. Solve[eqns, {x, y, z}, Reals] took longer than I was willing to wait, which was much less than the six hours of Bill. We can get the real solution by picking the ones for which {x, y, z} is an element of the Reals. Beware, just using N on the exact solutions without specifying a desired precision results in too much error.

AbsoluteTiming[
sol = Solve[Thread[Flatten@expr == 0], {x, y, z}];
realsol = Pick[sol, {x, y, z} \[Element] Reals /. # & /@ sol];
N@N[realsol, \$MachinePrecision] // Sort
]
(*
{14.680846, {{x -> -264.933, y -> -516.956, z -> 176.326},
{x -> -104.422, y -> 219.102,  z -> -146.78},
{x -> -78.0687, y -> 74.0697,  z -> 131.459},
{x -> 28.8908,  y -> 31.8585,  z -> -43.8798},
{x -> 77.1941,  y -> 125.209,  z -> 506.349},
{x -> 184.305,  y -> -229.324, z -> -238.648},
{x -> 618.239,  y -> -1340.34, z -> 616.787}}}
*)


Check:

Sort[{x, y, z} /. N[realsol, 16]] == Sort[{x, y, z} /. nrealsol]
(* True *)


Given the equation :

eqns = {(300 (1920 + 8 x - 21 y) (-80 + y))/(7 (30 + x)^3) + (2025 (208 x +
5 (-3446 + y)))/(52 (90 + y)^2) + (300 (4500 + 80 x - 21 z) (-100 +
z))/(7 (30 + x)^3) - (1521 (15425 - 539 x + 50 z))/(539 (39 + z)^2) == 0,
-((300 (1920 + 8 x - 21 y))/(7 (30 + x)^2)) - (2025 (-85 + x) (208 x + 5
(-3446 + y)))/(52 (90 + y)^3) + (2025 (49410 + 445 y - 208 z) (-45 + z))/
(52 (90 + y)^3) - (1521 (19680 - 539 y + 90 z))/(539 (39 + z)^2) == 0,
-((2025 (49410 + 445 y - 208 z))/(52 (90 + y)^2)) - (300 (4500 + 80 x - 21 z))/
(7 (30 + x)^2) - (1521 (-25 + x) (539 x - 25 (617 + 2 z)))/(539 (39 + z)^3) -
(1521 (-30 + y) (539 y - 30 (656 + 3 z)))/(539 (39 + z)^3) == 0} // Simplify;


Using NSolve to obain Real values and Complex having Real parts, thus making filtering simpler:

allSolns = NSolve[eqns, WorkingPrecision -> 8] // Simplify;


Extracting out the Real roots:

Realroots = Cases[allSolns, {_ -> _Real, _ -> _Real, _ -> _Real}]

{{x -> 618.23891, y -> -1340.3394, z -> 616.78706}, {x -> 184.30453,
y -> -229.32371, z -> -238.64759}, {x -> 77.194088, y -> 125.20903,
z -> 506.34885}, {x -> -78.068738, y -> 74.069706,
z -> 131.45885}, {x -> -104.42229, y -> 219.10230,
z -> -146.77983}, {x -> -264.93286, y -> -516.95612,
z -> 176.32593}, {x -> 28.890843, y -> 31.858477, z -> -43.879780}}

• Hi ! Please, read the help centre to get an insight on how to properly format your code Commented Aug 6, 2014 at 8:29

You can extract the roots and then select real ones:

x^3 - x + 1 == 0 // Solve[#, {x}]& // (x /. #)& // Select[# ∈ Reals &]