while loop convergence criteria

9 views (last 30 days)
GreGG
GreGG on 6 Nov 2012
Answered: Rajeev Upadhyay on 11 Apr 2020
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

Accepted Answer

Walter Roberson
Walter Roberson on 6 Nov 2012
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);

More Answers (1)

Rajeev Upadhyay
Rajeev Upadhyay on 11 Apr 2020
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))

Categories

Find more on Loops and Conditional Statements in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!