# Solving cubic equation with unknown coefficients

I am trying to solve a cubic equation with unknown coefficients (p and t), here is the code:

Rgas = 8.314; (* gas constant *)
acoef[tc_, pc_] := (27*Rgas^2*tc^2)/(64*pc);
bcoef[tc_, pc_] := (Rgas*tc)/(8*pc);
latentK = 76.9*10^3*39*1.67*10^-27*6.022*10^23;
tcK = 2223; (* Kelvin *)
pcK = 16*10^6; (*Pa *)
roots = v /.
NSolve[(p - Rgas*t)*v^3 - bcoef[tcK, pcK]*p*v^2 + acoef[tcK, pcK]*v
- acoef[tcK, pcK]*bcoef[tcK, pcK] == 0, v] // Chop;
test = roots[[1]]
(p - Rgas*t)*test^3 - bcoef[tcK, pcK]*p*test^2 + acoef[tcK, pcK]*test
- acoef[tcK, pcK]*bcoef[tcK, pcK] // Chop


1) I had to use NSolve[] because otherwise it said it did not have exact expressions. I tried to play around with Rationalize[], in vain. Is there any chance of finding an exact solution (i.e. using Solve[]) even if the coefficients (p and t) are not known?

2) In the last lines, I am plugging one of the roots (they are all the same) back into the equation, to see if it yields $0$. And it doesn't. Why? What am I doing wrong? Chop[]?

• This is the van der Waals equation, I take it? What are you solving for, critical volume? – J. M.'s technical difficulties Nov 23 '15 at 14:32
• I am solving for the volume, in terms of pressure and temperature. Critical pressure and critical temperature enter the formulae for the coeffcients,a and b (that I've called acoef and bcoef). – SuperCiocia Nov 23 '15 at 14:47
• Then, might I suggest reformulating the cubic to be in terms of the compressibility factor $Z=\frac{PV}{RT}$, and solving for $Z$ instead? – J. M.'s technical difficulties Nov 23 '15 at 14:52
• but $Z$ would not incorporate all $V$ factors...? For example, the $\frac{-bP}{RT}$ term would be written $\frac{-bZ}{V}$ – SuperCiocia Nov 23 '15 at 14:57
• $Z=\cfrac1{1-\cfrac{bP}{ZRT}}-\cfrac{aP}{(RT)^2 Z}$. I trust that you can figure out where this came from. – J. M.'s technical difficulties Nov 23 '15 at 15:03

If you clear the definition of Rgas then Solve will work

acoef[tc_, pc_] := (27*Rgas^2*tc^2)/(64*pc);
bcoef[tc_, pc_] := (Rgas*tc)/(8*pc);
latentK = 76.9*10^3*39*1.67*10^-27*6.022*10^23;
tcK = 2223;(*Kelvin*)
pcK = 16*10^6;(*Pa*)

eqn = (p - Rgas*t)*v^3 - bcoef[tcK, pcK]*p*v^2 + acoef[tcK, pcK]*v -
acoef[tcK, pcK]*bcoef[tcK, pcK];
asolns = Solve[eqn == 0, v]


And you can plug this back into your equation,

(eqn /. asolns)[[1]]


which looks pretty hairy - and maybe some combination of Simplify and ComplexExpand would show it to be zero, but I'm too lazy for that. I just test it with some numbers

% /. {Rgas -> 8.314, p -> 1, t -> 100}
(* 3.33067*10^-16 - 6.2626*10^-16 I *)


There's no reason not to use NSolve if you have Rgas defined as a real number in decimal format.

Rgas = 8.314;
solns = NSolve[eqn == 0, v];
test = (v /. solns)[[1]] // Chop
eqn /. v -> test


And then you put in values for p and t to see if it is zero,

% /. {Rgas -> 8.314, p -> 1, t -> 100}
(* -0.001301 *)


That didn't work. As you suspect, it was the Chop that messed you up,

test = (v /. solns)[[1]];
answer = eqn /. v -> test;
answer /. {Rgas -> 8.314, p -> 1, t -> 100}
(* 1.22791*10^-16 - 1.38741*10^-14 I *)


So I don't think you'll get a pretty, one-page solution to this equation, but the numerical one works as long as you don't apply Chop until you've put in values for all the variables.

• As I recall, for my thesis, I broke down and used Root[] instead to represent the roots of the van der Waals equation, since the full symbolic form was too unwieldy to be practical. – J. M.'s technical difficulties Nov 23 '15 at 15:05
• As a continuation, I now would like to use the real root (in terms of p and t) in a differnetial equation, to be solved numerically. Basically i would want dp/dt = const/realroot, where realroot = v /. asolns[[1]]. Do I somehow have to change realroot to be a function of p[t] and t? – SuperCiocia Nov 23 '15 at 19:26

Try the following: first solve exactly the cubic equation:

sl = Solve[a*v^3 - b*v^2 + c*v + d == 0, v]


The returned result is long and, therefore, I do not write it here. Get it by evaluating the code above. Then this

 Rgas = 8.314;(*gas constant*)
acoef[tc_, pc_] := (27*Rgas^2*tc^2)/(64*pc);
bcoef[tc_, pc_] := (Rgas*tc)/(8*pc);
latentK = 76.9*10^3*39*1.67*10^-27*6.022*10^23;
tcK = 2223;(*Kelvin*)pcK = 16*10^6;(*Pa*)

sl[[1, 1, 2]] /. {a -> (p - Rgas*t), b -> bcoef[tcK, pcK]*p,
c -> acoef[tcK, pcK], d -> -acoef[tcK, pcK]*bcoef[tcK, pcK]}


gives you your result, for example, for the first of the three roots, also rather long. Be aware that I did not check, if the solution is real, do it yourself.

(*
(0.0000481303 p)/(
p - 8.314 t) - (2^(
1/3) (-2.08487*10^-8 p^2 +
27.0199 (p - 8.314 t)))/(3 (6.02072*10^-12 p^3 + Sqrt[
4 (-2.08487*10^-8 p^2 +
27.0199 (p - 8.314 t))^3 + (6.02072*10^-12 p^3 -
0.0117043 p (p - 8.314 t) + 0.0351128 (p - 8.314 t)^2)^2] -
0.0117043 p (p - 8.314 t) + 0.0351128 (p - 8.314 t)^2)^(
1/3) (p - 8.314 t)) + (1/(
3 2^(1/3) (p - 8.314 t)))((6.02072*10^-12 p^3 + Sqrt[
4 (-2.08487*10^-8 p^2 +
27.0199 (p - 8.314 t))^3 + (6.02072*10^-12 p^3 -
0.0117043 p (p - 8.314 t) + 0.0351128 (p - 8.314 t)^2)^2] -
0.0117043 p (p - 8.314 t) + 0.0351128 (p - 8.314 t)^2)^(1/3))
*)


Have fun!

• thanks, why is there only one root though? Shouldn't there be 3? – SuperCiocia Nov 23 '15 at 15:01
• @Super, yes, there should be three, but in this case, you'll have either one real root (full vapor case) or three real roots (where only the two extreme roots have any physical meaning). – J. M.'s technical difficulties Nov 23 '15 at 15:11
• @SuperCiocia It is only an example to show you, how to calculate. All the rest you can do yourself. – Alexei Boulbitch Nov 23 '15 at 16:00
Clear[p, t, v];
Rgas = 8.314 // Rationalize;(*gas constant*)
acoef[tc_, pc_] := (27*Rgas^2*tc^2)/(64*pc);
bcoef[tc_, pc_] := (Rgas*tc)/(8*pc);
latentK = 769*10^2*39*167*10^-29*6022*10^20;
(* Note that you do not use latentK *)
tcK = 2223;(*Kelvin*)
pcK =
16*10^6;(*Pa*)

eqn = (p - Rgas*t)*v^3 - bcoef[tcK, pcK]*p*v^2 + acoef[tcK, pcK]*v -
acoef[tcK, pcK]*bcoef[tcK, pcK] == 0;

soln = Solve[eqn, v] // Simplify;


Verifying that soln satisfies eqn

And @@ (eqn /. soln // Simplify)

(*  True  *)

roots = v /. soln;
(*  long result  *)


Only the first root is real (for appropriate values of p and t)

fd = FunctionDomain[#, {p, t}] & /@ roots


Reduce[fd[[1]], {p, t}, Reals] // Simplify