while loop convergence criteria
4 次查看(过去 30 天)
显示 更早的评论
im writing a program to compute the definite integral of a function over a given interval using the midpoint rule. i have the program doing what i want but instead of just looping 1000 times to get a fairly accurate answer i want it to loop until the answer is increasing by less than .0001 im having trouble implementing that criteria.
heres my code so far
clc
clear
f = input('Gimme the function: ', 's');
f = str2func (['@(x)' f]);
n = 0;
a = input('lower bound ');
b = input('upper bound ');
while n<=1000
n=n+1;
delta = (b-a)/n;
x = a + delta*(1:2:2*n-1)/2;
m = delta*sum(f(x));
end
disp(['The integral is ' num2str(m
0 个评论
采纳的回答
Walter Roberson
2012-11-6
Before the loop,
old_fx = infinity;
new_fx = -infinity;
Change the loop condition to
while abs(old_fx - new_fx) > 0.0001
And just inside the loop finishing:
old_fx = new_fx;
new_fx = f(x);
更多回答(1 个)
Rajeev Upadhyay
2020-4-11
I need to resolve the convergence issue with While loop. Please suggest what should I do to run the program.
T_Res=582.0;
Molecular_Weight=16.64035;
Molecular_Wt_Air=28.97;
Gas_Gravity=Molecular_Weight/Molecular_Wt_Air;
K1=(0.00094+0.000002*Molecular_Weight)*(T_Res^1.5)/(200+19*Molecular_Weight+T_Res);
X=3.5+(986/T_Res)+(0.01*Molecular_Weight);
Y=2.4-0.2*X;
Sp_gr_oil=0.85;
Sp_gr_gas=0.6;
API_gr=(141.5/Sp_gr_oil)-131.5;
a=0.0125*API_gr-0.00091*(T_Res-460);
% for Z Calculation
T_crit=341.8116;
t=T_crit/T_Res;
P_crit=666.541;
E=0.06125*t*exp(-1.2*(1-t)^2)/P_crit;
% Dead Oil Viscosity
Z=3.0324-0.02023*Sp_gr_oil;
y=10^Z;
X=y*T_Res^(-1.163);
Dead_Oil_Viscosity=10^X-1;
% For Rel_Perm
Swcon=0.3;
Sorg=0.2;
Sgcon=0.05;
Krogcg=0.8;
Krgcl=0.9;
nog=2;
ng=2;
N=100;
P=3550;
GOR=0;
Np=0;
Gp=0;
Delta_P=50;
nstep=70;
Rsi=Sp_gr_gas*(((P/18.2)+1.4)*10^a);
Boi=0.9759+0.00012*((Rsi*(Sp_gr_gas/Sp_gr_oil)^0.5)+(1.25*(T_Res-460)))^1.2;
results=zeros(nstep,12);
f=inline('-E*P+((x+x^2+x^3-x^4)/(1-x)^3)-((14.76*t-9.76*t^2+4.58*t^3)*x^2)+((90.7*t-242.2*t^2+42.4*t^3)*x^(2.18+2.82*t))');
syms x
g=inline(diff(f(E,P,t,x)));
for i=1:nstep
P=P-Delta_P;
Rs=Sp_gr_gas*(((P/18.2)+1.4)*10^a);
Bo=0.9759+0.00012*((Rs*(Sp_gr_gas/Sp_gr_oil)^0.5)+(1.25*(T_Res-460)))^1.2;
A=10.715*(Rs+100)^(-0.515);
B=5.44*(Rs+150)^(-0.338);
Oil_Viscosity=A*Dead_Oil_Viscosity^B;
Gas_Density=0.00149406*P*Molecular_Weight/Z/T_Res;
Gas_Viscosity=K1*exp(X*Gas_Density^Y);
x=0.1;
while abs(f(E,P,t,x))>1e-05
x1=x-f(E,P,t,x)/g(x);
x=x1;
end
z=E*P/x;
Bg=0.00503676*z*T_Res/P;
m=0.1;
Np1=(N*m-Np);
Np=Np+Np1;
Gp1=(N*((Bo-Boi)+(Rsi-Rs)*Bg)-Np*(Bo-Rs*Bg))/Bg;
So=(1-Swcon)*(1-m)*Bo/Boi;
Sl=So+Swcon;
if Sl<(Sorg+Swcon)
Krog=0;
Krg=Krgcl;
elseif Sl>=(1-Sgcon);
Krog=Krogcg;
Krg=0;
else Krg=Krgcl*((1-Sl-Sgcon)/(1-Sgcon-Sorg-Swcon))^ng;
Krog=Krogcg*((Sl-Sorg-Swcon)/(1-Sgcon-Sorg-Swcon))^nog;
end
Ratio_Rel_Perm=Krg/Krog;
GOR1=Rs+(Ratio_Rel_Perm*Oil_Viscosity*Bo/Gas_Viscosity/Bg);
GOR_Avg=(GOR1+GOR)/2;
GOR=GOR_Avg;
Gp2=Gp+(GOR*Np1);
u=inline('((N*((Bo-Boi)+(Rsi-Rs)*Bg)-m*N*(Bo-Rs*Bg))/Bg)-(Gp+(GOR*Np1))')
syms m
v=inline(diff(u(Bg,Bo,Boi,GOR,Gp,m,N,Np1,Rs,Rsi)))
m=0.1;
while abs(u(Bg,Bo,Boi,GOR,Gp,N,Np1,Rs,Rsi,m))>0.5
m1=m-u(Bg,Bo,Boi,GOR,Gp,N,Np1,Rs,Rsi,m)/v(m);
m=m1;
end
Np=N*m;
Gp=(N*((Bo-Boi)+(Rsi-Rs)*Bg)-Np*(Bo-Rs*Bg))/Bg;
results(i,1)=P;
results(i,2)=Rs;
results(i,3)=Bo;
results(i,4)=Oil_Viscosity;
results(i,5)=Gas_Viscosity;
results(i,6)=Bg;
results(i,7)=Sl;
results(i,8)=Ratio_Rel_Perm;
results(i,9)=GOR;
results(i,10)=Np;
results(i,11)=Gp;
results(i,12)=m;
end
plotyy(results(:,1),results(:,2),results(:,1),results(:,3))
plot(results(:,1),results(:,10))
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!