Solving an implicit equation in matlab

Dear all,

I need to solve the following equation:

and the equation that coorelates by the following equation( the term 5 was ignored):

and the variables b, A, B, C, k, Tc can be founded in this table:

As the variavle V can not be isolated, we should solve it implicit. I am trying witht he following code:

Tc = 367.85;
Pc = 3374.84;
R = 7.290839*10.^-2;
b = 3.973395*10.^-4;
k = 5.033;
A2 = -0.1398755;
A3 = 4.031561*10^-4;
A4 = -1.508838*10.^-7;
B2 = 1.677524*10.^-4;
B3 = -6.739375*10.^-7;
B4 = 0;
C2 = -1.518125;
C3 = 5.521381*10^-5;
C4 = 0;
Vapor_pressure_ref = 100
syms F P V T
F = -P + ((R*T)/(V-b)) + (((A2) + (B2*T) + (C2*exp(-k*T/Tc)))/(V-b)^2)  +  (((A3) + (B3*T) + (C3*exp(-k*T/Tc)))/(V-b).^3) + (((A4) + (B4*T) + (C4*exp(-k*T/Tc)))/(V-b).^4);

dV_dT = - diff(F,T)/diff(F,V)
Eq = (V - (dV_dT*T))
int(Eq,P,Vapor_pressure_ref,1)

However,the output is:

-(17211847986098922057056496272843303410021431145154622943954950247403281191336156397430653844910382510721336111875*exp(3956657/1051000) + 171221942040754850037313460159985281792124894590638846549063088630653073472919905693201072792645916265515881332736)/(1065840624671895280778428055080140800000*(64450689896085148577946341147025193190720307559496600271150138269013626089*exp(3956657/1051000) + 236272089264744816453435291457063368284234553216075863541003706881470889984))

How to simplify it?

[Torsten]

How should "int" know that V is the solution of your quartic equation in V

P - RT/(V-b) - sum(...) = 0

and that this V has to be inserted in the expression for dV_dT ?

And the integration must be for constant T, I guess - so T has to be specified.

[Star Strider]

First, plot the function —

    syms F P V T
    Tc = 367.85;
    Pc = 3374.84;
    R = 7.290839*10.^-2;
    b = 3.973395*10.^-4;
    k = 5.033;
    A2 = -0.1398755;
    A3 = 4.031561*10^-4;
    A4 = -1.508838*10.^-7;
    B2 = 1.677524*10.^-4;
    B3 = -6.739375*10.^-7;
    B4 = 0;
    C2 = -1.518125;
    C3 = 5.521381*10^-5;
    C4 = 0;
    Vapor_pressure_ref = 100;
%     syms F P V T
    F = -P + ((R*T)/(V-b)) + (((A2) + (B2*T) + (C2*exp(-k*T/Tc)))/(V-b)^2)  +  (((A3) + (B3*T) + (C3*exp(-k*T/Tc)))/(V-b).^3) + (((A4) + (B4*T) + (C4*exp(-k*T/Tc)))/(V-b).^4);
    
    dV_dT = - diff(F,T)/diff(F,V);
    Eq = (V - (dV_dT*T));
    IntP(T,V) = int(Eq,P,Vapor_pressure_ref,1);
    
    figure
    fsurf(IntP, [-5000 5000 -5000 5000], 'MeshDensity',150)
    hold on
    patch([-5000 5000 5000 -5000], [-5000 -5000 5000 5000], zeros(1,4), 'r')
    hold off
    grid on
    zlim([-1 1]*5E-1)
    xlabel('T')
    ylabel('V')
    view(35,30)

What do you want to do?

.

[Torsten]

How to simplify it?

Use double(expression).

But this is not the output of the above code.

[Yogan]

Lets imagem that our P-V-T equation is PV = RT.

In this case we can isolate the variavle V and we will have: V = RT/P

So we can Derivate V regarding T (DV/DT) and we will have: V = R/P or V/T (because I can write P = RT/V)

But in my case, my PVT equation in something that I can not isolate the variable V (as I have done here). For this reason, I thought to use implicit deviratives.... I don't wanna extend this comment, but basically if I have a F(x,y) and I am interested on DX/DY, I can use the following equation:

DX/DY = - (DF/DY)/ (DF/DX)

This link you can have more details about Implicit Differentiation

Once I have the DV/DT done, theoreticaly I need just to multiply by -T, add V and integrate from Pressure 1 to pressure 2.

[Torsten]

Once I have the DV/DT done, theoreticaly I need just to multiply by -T, add V and integrate from Pressure 1 to pressure 2.

Yes, but what are these T and V in the integrand ? Both dF/dT and dF/dV depend on P,T and V. And T and V appear directly in the expression of the integrand. So you must know which T and V values are to be inserted in the integrand for the pressures between 1 and 100 in order to evaluate the integral.

Answer:

As I wrote, T must be prescribed and V must be derived as the solution of your equation-of-state for given P within the integration process and the specified T.

[Bruno Luong]

"T must be prescribed"

Correct. OP forgets to tell us.

Answer 1

[Bruno Luong]

The numerical answer is -0.284179329144572

format long
-(17211847986098922057056496272843303410021431145154622943954950247403281191336156397430653844910382510721336111875*exp(3956657/1051000) + 171221942040754850037313460159985281792124894590638846549063088630653073472919905693201072792645916265515881332736)/(1065840624671895280778428055080140800000*(64450689896085148577946341147025193190720307559496600271150138269013626089*exp(3956657/1051000) + 236272089264744816453435291457063368284234553216075863541003706881470889984))

ans = 
  -0.284179329144572

Answer 2

[Bruno Luong]

D=1/(V-b) is a root of fifth order polynomial, where the coefficients depend on P, T, etc...

It looks like it is small non-linear correction of the linear equation.

P = RT*D

I assume the equation has single real root in the range you consider, so actually it is not really an implicit equation, it is explicit equation if you write it as root of the polynomial.

It seems int is smart enough to figure this fact out, and able to gie back the expression you look for.