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

4 次查看(过去 30 天)
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 个评论
Temitope Amao
Temitope Amao 2018-3-20
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

请先登录,再进行评论。

回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Chemistry 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by