根据Maxwell Construction，可以求出气液共存密度，参考博主Maxwell Construction解法 - 知乎 (zhihu.com)，根据等面积法则逼近气液共存状态，这对于远离临界压力的求解是相对准确的，然而如果要求解临界压力附近的气液密度，比如说我们就要求解临界密度，这样的求解方式存在局限性，存在数值误差，难以求得临界密度。
某一压力下的P-R状态方程P-V图
我们以P-R方程为例：
，整理得
当时，存在三个不相等的实根，当
时，存在三个完全相等的实根，此时
的一阶导只存在一个零点，即二次方程
存在两个完全相同的实根，根据二元一次方程求根公式，我们很容易求得
，这就是临界点的比容，
，其中A和B与系数a无关，及系数a不会影响临界密度。
Matlab求解如下：
clear, clc;
a = 2/49;
b = 2/21;
R = 1;
ome = 0.344;
Tc = 0.07779607 / 0.45723553 * a / (b * R);
Pc = 0.07779607 * R * Tc / b;
T = Tc;
phi = (1 + (0.37464 + 1.54226*ome - 0.26992*ome^2)*(1 - sqrt(T/Tc)))^2;
% p = R*T./(V - b) - a*phi./(V.^2 + 2*b.*V - b^2);
A=Pc*b^3+R*T*b^2-a*phi*b;
B=-3*Pc*b^2-2*b*R*T+a*phi;
C=Pc*b-R*T;
D=Pc;
%rho=Solve3Polynomial(A,B,C,D)
A1=Pc;
B1=b*Pc-R*T;
C1=-3*Pc*b^2-2*R*T*b+a*phi;
D1=Pc*b^3+R*T*b^2-a*phi*b;
%V=Solve3Polynomial(A1,B1,C1,D1)
A2=3*A1;
B2=2*B1;
C2=C1;
Delta=B2^2-4*A2*C2;
vc=-B2/(2*A2)
rhoc=1.0/vc
结果：
vc =
0.3763
rhoc =
2.6573
以下两篇文献的rhoc我觉得是错了：
参考代码：(12条消息) MATLAB实现一元三次方程求解/盛金公式_matlab解一元三次方程_TocI的博客-CSDN博客
