フィルターのクリア

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 Center および File ExchangeLoops and Conditional Statements についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by