# Can't obtain Real results using NSolve

I am trying to find the critical properties of the Carnahan-Starling equation of state,

$$p=\rho R T\frac{1+b\rho /4+(b\rho/4)^2-(b\rho/4)^3}{(1-b\rho /4)^3}-a\rho^2$$ using this system of equations,

$$\frac{\partial p}{\partial \rho}=0$$

$$\frac{\partial^2p}{\partial\rho^2}=0.$$

And the following simple code was prepared for Mathematica,

P[T_, rho_] := rho R T (1 + (b rho)/4 + ((b rho)/4)^2 - ((b rho)/4)^3)/(1 - (b rho)/4)^3 - a rho^2

NSolve[
{
D[P[T, rho], {rho, 1}] == 0,
D[P[T, rho], {rho, 2}] == 0
},
{T, rho}
] // MatrixForm


Once evaluated, these results obtained (converted to Latex format for the sake of readability),

$$\left\{\left\{T\to -\frac{20.4876 a}{b R},\rho \to -\frac{5.13672}{b}\right\},\left\{T\to -\frac{38.7085 a}{b R},\rho \to -\frac{1.69831}{b}\right\},\left\{T\to \frac{0.377315 a}{b R},\rho \to \frac{0.521776}{b}\right\},\left\{T\to \frac{(33.6447\, -16.1071 i) a}{b R},\rho \to \frac{13.1566\, -7.20194 i}{b}\right\},\left\{T\to \frac{(33.6447\, +16.1071 i) a}{b R},\rho \to \frac{13.1566\, +7.20194 i}{b}\right\}\right\}$$

But the above output contains complex numbers. Using Reals in NSolve, I tried to remove unphysical results,

NSolve[
{
D[P[T, ρ], {ρ, 1}] == 0,
D[P[T, ρ], {ρ, 2}] == 0
},
{T, ρ},
Reals
] // TableForm


Unfortunately that was not helpful; because, for example one of the answers for $\rho$ was,

ρ ->
Root[-1024. + 1280. b #1 + 1280. b^2 #1^2 + 64. b^3 #1^3 -
20. b^4 #1^4 + b^5 #1^5 &, 1]


and the answers for $T$ are much longer.

I am totally unable to interpret this. How can I tell Mathematica to only return Real results?

• If output = NSolve[ ... ], then Select[output, FreeQ[#, HoldPattern[Complex[__]]]&]. – march Jan 14 '16 at 4:22
• Sorry, would you please elaborate on more? I am new to Mathematica. – Shaqpad Jan 14 '16 at 4:25
• Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – Michael E2 Jan 14 '16 at 4:27
• You just need to copy and paste. Modify your code so that you have output = NDSolve[ <code that you already have written goes here> ], and then evaluate Select[output, FreeQ[#, HoldPattern[Complex[__]]]&]. This selects from your list those pairs of T and rho that don't have the complex number I in them. (Although you will have to get rid of theMatrixForm part.) – march Jan 14 '16 at 4:30
• @MichaelE2 Of course; But I'm new to Mathematica itself too! – Shaqpad Jan 14 '16 at 4:36

You get a root-expression, because the exponent is 5!

{D[P[T, rho], {rho, 1}] == 0, D[P[T, rho], {rho, 2}] == 0} // Simplify


You can only solve the system numerical.

P[T_, rho_] := rho R T (1 + (b rho)/4 + ((b rho)/4)^2 - ((b rho)/4)^3)/(1 - (b rho)/4)^3 - a rho^2
Solve[{D[P[T, rho], {rho, 1}] == 0, D[P[T, rho], {rho, 2}] == 0}, {T, rho}, Reals];


b and R just as an example

% /. {b -> 1, R -> 1} // N
{{T -> -20.4876 a, rho -> -5.13672}, {T -> -38.7085 a, rho -> -1.69831}, {T -> 0.377315 a, rho -> 0.521776}}