Integral Error message: first input argument must be a function handle

4 visualizzazioni (ultimi 30 giorni)
h= @(p)(p./DAK(0.6,p,150))
pseudop = integral(h,14.7,5000)
DAK is a function that's largely dependent on "p"(pressure), the other inputs are constant. The first line of code works fine but I keep getting an error for the second line of code. I tried fitting an equation with p as the variable for function DAK but the fit is bad and will negatively affect my solution..
This is the error message: Error using integral (line 82) First input argument must be a function handle.
This doesn't work either because the program keeps running...
integral(@(p)(p./DAK(0.6,p,150)),14.7,5000)
  2 Commenti
Temitope Amao
Temitope Amao il 20 Mar 2018
Here it is:
function Z = DAK(SGg,p,T) % gas compressibility factor
% Calculates the compressibility factor of natural gases for a range of pressures (Psia) 'minP' to
% 'maxP' in steps 'Pstep' (for a given temperature T(degree Farenheit) and specific gravity sg) using the Dranchuk-Abbou Kassem
% Equation of state. The range of validity is 0.2<Ppr<30 and 1.0<Tpr<3.0 (Craft
% &Hawkins, 1991)
%To calculate for a single pressure point P, set minP = maxP = Pstep
format short g
i = 1;
%for P = minP:Pstep:maxP %
for P = p %
% Estimation of Psuedocritical pressure and temperature based on
% suttons correlations (Craft & Hawkins, 1991)
Ppc = 677+15*SGg-37.5*SGg^2;
Tpc = 168+325*SGg-12.5*SGg^2;
%Calculation of Pseudo reduced pressure and temperature
Ppr = P / Ppc;
Tpr = (T+459.67)/ Tpc;
% Correlation Constants, c4 has to be updated within the root finding
% algorithm
A1 = 0.3265;
A2 = -1.0700;
A3 = -0.5339;
A4 = 0.01569;
A5 = -0.05165;
A6 = 0.5475;
A7 = -0.7361;
A8 = 0.1844;
A9 = 0.1056;
A10 = 0.6134;
A11 = 0.7210;
c1 = A1 + (A2/Tpr) + (A3/(Tpr^3))+ (A4/(Tpr^4))+ (A5/(Tpr^5));
c2 = A6 + (A7/Tpr) + (A8/(Tpr^2));
c3 = A9*((A7/Tpr) + (A8/(Tpr^2)));
% Implementation of bisection root finding method (The Newton-Raphson method becomes unstable and slower
% at low temperature and higher specific gravities, the bisection method
% appears more robust for all cases in this specific application)
brack1 = 1e-2;
brack2 = 4;
itr= 0;
while 2<3
zguess = (brack1 + brack2) / 2;
Rrguess = (0.27*Ppr) /(zguess*Tpr);
c4 = (A10)*(1+(A11*(Rrguess^2)))*((Rrguess^2)/(Tpr^3))*(exp(-A11*(Rrguess^2)));
syms z
Rr = (0.27*Ppr) /(z*Tpr);
f = z+ (c3*(Rr^5)) - (c2*(Rr^2)) - (c1*(Rr^1)) - c4 - 1;
fzguess = subs(f,z,zguess);
convergence = abs(fzguess-0);
conv(itr+1) = convergence;
disp_convergence = conv'; % remove ";" to display convergence progression for each pressure point
if convergence <=10^-4
itr = itr + 1;
break
end
if fzguess < 0
brack1 = zguess;
else
brack2 = zguess;
end
itr = itr + 1;
end
Z(i) = zguess;
Pressure(i) = P;
I__________P__________Z_____________I = [Pressure' Z'];
i=i+1;
%Uncomment following comments to view convergence progression for each
%pressure point
% plot(disp_convergence)
% xlabel('Convergence','fontsize',12)
% ylabel('Iteration','fontsize',12)
% hold on
% pause
clear conv
end
% hold off
% figure;plot(Pressure,Z, 'r*','MarkerSize',6)
% xlabel('Pressure','fontsize',12);
% ylabel('Z','fontsize',12)
end

Accedi per commentare.

Risposte (0)

Categorie

Scopri di più su Chemistry in Help Center e File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by