Multiple Output for a function made for a single output.

Hi

I'm currently having issue with how to set up this problem and use the current matlab function i have. I currently have a function designed to calculate Vapour pressure. The problem is i'm currently trying to calculate several vapour pressures in a mixture. So when i run it a single vapour pressure is calculated..

My input currently looks like:

Pc=[4599000 4872000 4248000 3796000 3370000 3025000 2490000 1817000 1401000]; 
Tc=[190.6 305.3 369.8 425.1 469.7 507.6 568.7 658.1 722];
P=0.1E6;
T=320;
Pso=0.2e6;

Where the Function calulates Ps value, instead of calculating several i get a single value.

This is the function for vapour pressure:

function [Ps,VL,VV] = Ps_VDW(Pc, Tc, Pso, T)
%The funtion Ps_VDW calculates the vapor pressure and the specific molar
%volume V for a pure fluid at given temperature T using
%VDW-CEoS. This function requires Z_VDW and FUGACITY_VDW funtions to run.
%All input/output data are expressed in SI.
%P[Pa], T[K], V[m^3/mol]

%First value for Ps
Ps=Pso;

%calculation of VL and VV using Z_VDW function
[Z,V,ZL,ZV,VL,VV] = Z_VDW(Pc, Tc, Ps, T);

%comparing VL and VV values to performed calculations
if VL==VV
    disp('WRONG INITIAL Pso INPUT VALUE');
else
    [fiL,fiV,fL,fV] = FUGACITY_VDW(Pc, Tc, Ps, T);
    
    iter=0;
    while abs(1-fL/fV)>0.00001
        Ps=Ps*(fL/fV);
        [fiL,fiV,fL,fV] = FUGACITY_VDW(Pc, Tc, Ps, T);
        iter=iter+1;
    end
    
    %calculation of VL and VV using Z_VDW function to be sure
    [Z,V,ZL,ZV,VL,VV] = Z_VDW(Pc, Tc, Ps, T);
    
    if VL==VV
        disp('ALGORITHM FAILED')
    else
        Ps=Ps;
        VL=VL;
        VV=VV;
    end
        
end

The compressibility subprogram used:

function [Z,V,ZL,ZV,VL,VV] = Z_VDW(Pc, Tc, P, T)
%The funtion Z_VDW calculates the compressibility factor z and the specific 
%molar volume V for a pure fluid at given pressure P, temperature T using 
%VDW-CEoS. This function requires REALROOTS_CARDANO funtion to run.
%All input/output data are expressed in SI.
    %P[Pa], T[K], V[m^3/mol], Z[dimensionless]

%universal gas constant
R=8.314;

%VDW attractive constant and covolume at critical point
a=27*R.^2*Tc.^2/(64*Pc);
b=R*Tc/(8*Pc);

%calculation of A and B
A=a*P/(R*T)^2;
B=b*P/(R*T);

%calculation of polynomial coefficients a2, a1, a0 for VDW_CEoS
a2=-1-B;
a1=A;
a0=-A*B;

%calculation of Z and V using REALROOTS_CARDANO function
Z=REALROOTS_CARDANO(a2,a1,a0);
V=Z*R*T/P;

%selection of Z and V values for liquid and vapor
if min(V)<b
    ZL=max(Z);
    ZV=max(Z);
    VL=max(V);
    VV=max(V);
else
    ZL=min(Z);
    ZV=max(Z);
    VL=min(V);
    VV=max(V);
end

end

The fugacity subprogram used:

function [fiL,fiV,fL,fV] = FUGACITY_VDW(Pc, Tc, P, T)
%The funtion FUGACITY_VDW calculates the fugacity (f) and fugacity 
%coefficient (fi) for a pure fluid at given pressure P, temperature T using 
%VDW-CEoS. This function requires Z_VDW funtion to run.
%All input/output data are expressed in SI.
    %P[Pa], T[K], f[Pa], fi[dimensionless]

%universal gas constant
R=8.314;

%VDW attractive constant and covolume at critical point
a=27*R.^2*Tc.^2/(64*Pc);
b=R*Tc/(8*Pc);

%calculation of A and B
A=a*P/(R*T)^2;
B=b*P/(R*T);

%calculation of ZL and ZV using Z_VDW function
[Z,V,ZL,ZV,VL,VV] = Z_VDW(Pc, Tc, P, T);

%calculation of fugacity and fugacity coefficient for liquid and vapor 
%fL, fV, fiL and fiV
fiL=exp(ZL-1-log(ZL-B)-A/ZL);
fiV=exp(ZV-1-log(ZV-B)-A/ZV);
fL=fiL*P;
fV=fiV*P;

end

Real root subprogram:

function [RR]=REALROOTS_CARDANO(a2,a1,a0)
%This function calculates the real roots of a polynomial of degree 3 using 
%an algorithm based on the Cardano's formula. 
%The polynomial must be written as: x^3 + a2*x^2 + a1*x + a0. 
%Input data: a2, a1, a0.
%Output data: array containing the real roots. 
%The output array can be 1x1 or 1x3. 

Q=(a2^2-3*a1)/9;
R=(2*a2^3-9*a1*a2+27*a0)/54;

    if Q==0 && R==0 %Case of single real root with multiplicity 3
        RR=-a2/3;
    else
        if Q^3-R^2>=0 
            theta=acos(R/(Q^(3/2)));
            RR(1)=-2*sqrt(Q)*cos(theta/3)-a2/3;
            RR(2)=-2*sqrt(Q)*cos((theta+2*pi)/3)-a2/3;
            RR(3)=-2*sqrt(Q)*cos((theta+4*pi)/3)-a2/3;
        else
            RR=-sign(R)*((sqrt(R^2-Q^3)+abs(R))^(1/3)+Q/(sqrt(R^2-Q^3)+abs(R))^(1/3))-a2/3;
        end
    end
end

[Stephen23]

Most likely you need to vectorize your code:

https://www.mathworks.com/help/matlab/matlab_prog/vectorization.html

In particular:

https://www.mathworks.com/help/matlab/matlab_prog/array-vs-matrix-operations.html

Not that your code throws an error:

>> [Ps,VL,VV] = Ps_VDW(Pc(1), Tc(1), Pso, T)
Undefined function 'REALROOTS_CARDANO' for input arguments of type 'double'.
Error in Z_VDW (line 20)
Z=REALROOTS_CARDANO(a2,a1,a0);
Error in Ps_VDW (line 10)
[Z,V,ZL,ZV,VL,VV] = Z_VDW(Pc, Tc, Ps, T); 

[Abdullahi Geedi]

Thanks. I've added the final subprogram, which is the error you see.

Do i need to vectorize the function or my input when using the function.

[Stephen23]

"Do i need to vectorize the function or my input when using the function."

Either:

• vectorize the function code, or
• call the function in a loop. Use indexing so that the function is called on scalar values on each loop iteration. If you want to store the output values in arrays then read this: https://www.mathworks.com/help/matlab/matlab_prog/preallocating-arrays.html