How to can I solve a cubic equation with wide ranges of T and P?

I want to calculate the volume using VdW. It works properly with assigned values, but I want to calculate using ranges of P and T values.

R = 8.314;
Tc = 304.1282;
Pc = 7.3773*1000;
a = (27/64)*(R^2*Tc^2/Pc);
b = (R*Tc)/(8*Pc);
T = 350; %instead of 350, I want to calculate from 250 to 350
P = 5; %Same here, I want P values from 5 to 70

A = 1 ;
B = -(b + ((R*T)/P)) ;
C = a/P ;
D = -((a*b)/P) ;

% A*V^3 - B*V^2 + C*V - D = 0
vdw = [A B C D];
V = real(roots(vdw));
Vg = max(V);

Answer 1

[John D'Errico]

MATLAB does not offer an interval arithmetic solution for roots. (Nor will it ever do so, I predict.) This would get pretty complicated if you even tried to use interval arithmetic in some form, since the roots will be often complex for some parameter values.

Could you just use various values for your parameters, then call roots multiple times, and record the different sets of roots found? Yes. But be careful, as the roots can easily "swap" places if that matters to you. Regardless, this requires nothing more than a loop. (Could you get fancy, and use arrayfun, creating an implicit loop internally? Yes. But a simple for loop is trivial to write, use and debug, and it would be no slower in the end.)

Answer 2

[Walter Roberson]

Q = @(v) sym(v);

R = Q(8314)/Q(1000);
Tc = Q(3041282)/Q(10000);
Pc = Q(73773)/Q(10000)*Q(1000);
a = (Q(27)/Q(64))*(R^2*Tc^2/Pc);
b = (R*Tc)/(8*Pc);
syms T   %250 to 350
syms P   %5 to 70

A = Q(1) ;
B = -(b + ((R*T)/P)) ;
C = a/P ;
D = -((a*b)/P) ;

% A*V^3 - B*V^2 + C*V - D = 0
syms V
vdw = [A B C D];
VDW = poly2sym(vdw, V)
V = solve(VDW, 'maxdegree', 3)

You can now study diff(V,T) looking for critical points for T, and diff(V,P) looking for critical points for P. When I looked at diff(V,T) I find there are not critical points in the range 250 to 350. That means that you can study the range of values by substituting the endpoints of the range.

I did not examine diff(V,P)