필터 지우기
필터 지우기

while loop convergence criteria

조회 수: 3 (최근 30일)
GreGG
GreGG 2012년 11월 6일
답변: Rajeev Upadhyay 2020년 4월 11일
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

채택된 답변

Walter Roberson
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
GreGG
GreGG 2012년 11월 6일
thank you that worked perfect

댓글을 달려면 로그인하십시오.

추가 답변 (1개)

Rajeev Upadhyay
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))

카테고리

Help CenterFile Exchange에서 Loops and Conditional Statements에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by