I feel difficulty to find the roots of peng robinson eqution given below ..can any body tell me how to find the roots with respect to volume.

5 views (last 30 days)
P *V^3 + (b*P - R*T) *V ^2 - (3*b^2P + 2*b*R*T - a)*V + (b^3P + R*T*b^2 - a*b) = 0

Answers (1)

John D'Errico
John D'Errico on 5 Sep 2016
If the constants are known, then just use roots! WTP?
If the constants are not known, then it will get messy, since you have a non-constant coefficient cubic polynomial. Using solve will do it though, since the degree is low enough. It will look a bit messy though.
syms P V b R T a
vpa(solve(P *V^3 + (b*P - R*T) *V ^2 - (3*b^2*P + 2*b*R*T - a)*V + (b^3*P + R*T*b^2 - a*b) == 0,V),16)
ans =
((0.1111111111111111*(R*T - 1.0*P*b)^2)/P^2 + (0.3333333333333333*(3.0*P*b^2 + 2.0*R*T*b - 1.0*a))/P)/((0.03703703703703704*(R*T - 1.0*P*b)^3)/P^3 + (((0.03703703703703704*(R*T - 1.0*P*b)^3)/P^3 - (0.5*P*b^3 + 0.5*R*T*b^2 - 0.5*a*b)/P + (0.1666666666666667*(R*T - 1.0*P*b)*(3.0*P*b^2 + 2.0*R*T*b - 1.0*a))/P^2)^2 - 1.0*((0.1111111111111111*(R*T - 1.0*P*b)^2)/P^2 + (1.0*P*b^2 + 0.6666666666666667*R*T*b - 0.3333333333333333*a)/P)^3)^(1/2) - (0.5*P*b^3 + 0.5*R*T*b^2 - 0.5*a*b)/P + (0.1666666666666667*(R*T - 1.0*P*b)*(3.0*P*b^2 + 2.0*R*T*b - 1.0*a))/P^2)^(1/3) + ((0.03703703703703704*(R*T - 1.0*P*b)^3)/P^3 - (0.5*(P*b^3 + R*T*b^2 - 1.0*a*b))/P + (((0.03703703703703704*(R*T - 1.0*P*b)^3)/P^3 - (0.5*(P*b^3 + R*T*b^2 - 1.0*a*b))/P + (0.1666666666666667*(R*T - 1.0*P*b)*(3.0*P*b^2 + 2.0*R*T*b - 1.0*a))/P^2)^2 - 1.0*((0.1111111111111111*(R*T - 1.0*P*b)^2)/P^2 + (0.3333333333333333*(3.0*P*b^2 + 2.0*R*T*b - 1.0*a))/P)^3)^(1/2) + (0.1666666666666667*(R*T - 1.0*P*b)*(3.0*P*b^2 + 2.0*R*T*b - 1.0*a))/P^2)^(1/3) + (0.3333333333333333*(R*T - 1.0*P*b))/P
(((0.1111111111111111*(R*T - 1.0*P*b)^2)/P^2 + (0.3333333333333333*(3.0*P*b^2 + 2.0*R*T*b - 1.0*a))/P)*(- 0.5 - 0.8660254037844386i))/((0.03703703703703704*(R*T - 1.0*P*b)^3)/P^3 + (((0.03703703703703704*(R*T - 1.0*P*b)^3)/P^3 - (0.5*P*b^3 + 0.5*R*T*b^2 - 0.5*a*b)/P + (0.1666666666666667*(R*T - 1.0*P*b)*(3.0*P*b^2 + 2.0*R*T*b - 1.0*a))/P^2)^2 - 1.0*((0.1111111111111111*(R*T - 1.0*P*b)^2)/P^2 + (1.0*P*b^2 + 0.6666666666666667*R*T*b - 0.3333333333333333*a)/P)^3)^(1/2) - (0.5*P*b^3 + 0.5*R*T*b^2 - 0.5*a*b)/P + (0.1666666666666667*(R*T - 1.0*P*b)*(3.0*P*b^2 + 2.0*R*T*b - 1.0*a))/P^2)^(1/3) - ((0.03703703703703704*(R*T - 1.0*P*b)^3)/P^3 - (0.5*(P*b^3 + R*T*b^2 - 1.0*a*b))/P + (((0.03703703703703704*(R*T - 1.0*P*b)^3)/P^3 - (0.5*(P*b^3 + R*T*b^2 - 1.0*a*b))/P + (0.1666666666666667*(R*T - 1.0*P*b)*(3.0*P*b^2 + 2.0*R*T*b - 1.0*a))/P^2)^2 - 1.0*((0.1111111111111111*(R*T - 1.0*P*b)^2)/P^2 + (0.3333333333333333*(3.0*P*b^2 + 2.0*R*T*b - 1.0*a))/P)^3)^(1/2) + (0.1666666666666667*(R*T - 1.0*P*b)*(3.0*P*b^2 + 2.0*R*T*b - 1.0*a))/P^2)^(1/3)*(0.5 - 0.8660254037844386i) + (0.3333333333333333*(R*T - 1.0*P*b))/P
(((0.1111111111111111*(R*T - 1.0*P*b)^2)/P^2 + (0.3333333333333333*(3.0*P*b^2 + 2.0*R*T*b - 1.0*a))/P)*(- 0.5 + 0.8660254037844386i))/((0.03703703703703704*(R*T - 1.0*P*b)^3)/P^3 + (((0.03703703703703704*(R*T - 1.0*P*b)^3)/P^3 - (0.5*P*b^3 + 0.5*R*T*b^2 - 0.5*a*b)/P + (0.1666666666666667*(R*T - 1.0*P*b)*(3.0*P*b^2 + 2.0*R*T*b - 1.0*a))/P^2)^2 - 1.0*((0.1111111111111111*(R*T - 1.0*P*b)^2)/P^2 + (1.0*P*b^2 + 0.6666666666666667*R*T*b - 0.3333333333333333*a)/P)^3)^(1/2) - (0.5*P*b^3 + 0.5*R*T*b^2 - 0.5*a*b)/P + (0.1666666666666667*(R*T - 1.0*P*b)*(3.0*P*b^2 + 2.0*R*T*b - 1.0*a))/P^2)^(1/3) - ((0.03703703703703704*(R*T - 1.0*P*b)^3)/P^3 - (0.5*(P*b^3 + R*T*b^2 - 1.0*a*b))/P + (((0.03703703703703704*(R*T - 1.0*P*b)^3)/P^3 - (0.5*(P*b^3 + R*T*b^2 - 1.0*a*b))/P + (0.1666666666666667*(R*T - 1.0*P*b)*(3.0*P*b^2 + 2.0*R*T*b - 1.0*a))/P^2)^2 - 1.0*((0.1111111111111111*(R*T - 1.0*P*b)^2)/P^2 + (0.3333333333333333*(3.0*P*b^2 + 2.0*R*T*b - 1.0*a))/P)^3)^(1/2) + (0.1666666666666667*(R*T - 1.0*P*b)*(3.0*P*b^2 + 2.0*R*T*b - 1.0*a))/P^2)^(1/3)*(0.5 + 0.8660254037844386i) + (0.3333333333333333*(R*T - 1.0*P*b))/P

Categories

Find more on Mathematics and Optimization in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!