how to avoid dividing by zeros.

114 views (last 30 days)
Lin LI
Lin LI on 7 Oct 2011
Answered: cui,xingxing on 24 Oct 2023
i run a code but it gave me warnings that Warning that dividing by zero and index exceeding matrix dimensions. I looked into it and found probaby it is caused by this sentence in function sumin
if Y2(i)<=Thickness
MeanBedVelocity(i)=quad(ff1,Y1,Y2(i))./(quad(ff2,Y1,Y2(i)));
how to avoid dividing by zeros, how can i fix this problem, thank you very much for your help.
here is the code. it includes one main code and 3 function codes.
Di=4e-3;
PR=double_int(0,Di,0,Di)
function SS=double_int(innlow,innhi,outlow,outhi)
Protrusion1=outlow;Protrusion2=outhi;Friction1=innlow;Friction2=innhi;
SS=quad(@G_yi,Protrusion1,Protrusion2,[],[],Friction1,Friction2);
function f=G_yi(Protrusion,Friction1,Friction2)
Protrusion=Protrusion(:);n=length(Protrusion);
% if isnumeric(Friction1)==1;FrictionF1=Friction1*ones(size(Protrusion));
% else FrictionF1=feval(Friction1,Protrusion);
% end
% if isnumeric(Friction2)==1;FrictionF2=Friction2*ones(size(Protrusion));
% else FrictionF2=feval(Friction2,Protrusion);
% end
save G_yi.mat;
for i=1:n
f(i)=quad(@sumin,Friction1,Friction2,[],[],i);
end
f=f(:);
function fun=sumin(Friction,i)
% global Protrusion;
% Protrusion(i)=evalin(Protrusion(i));
load G_yi.mat;
MeanShearStress=5 ; %unit Pa
Density=1e3;
Ustar=sqrt(MeanShearStress/Density);
N=15;
D=zeros(N,1);
p=zeros(N,1);
for k=1:N-1;
D(1)=0.5e-3;
p(1)=1/N;
D(k+1)=D(k)+0.5e-3;
p(k+1)=p(k);
end;
%p1=1/3, p2=1/3,p3=1/3
%D1=40e-6,D2=60e-6, D3=80e-6
KinematicViscosity=1.004e-6;
D50=4e-3;
Roughness=2*D50;
RoughnessRenolds=Ustar*Roughness/KinematicViscosity;
if RoughnessRenolds>100
Intensity=Ustar*2.14;
Skewness=0.43;
Flatness=2.88;
else
Intensity=Ustar*(-0.187*log(RoughnessRenolds)+2.93);
Skewness=0.102*log(RoughnessRenolds);
Flatness=0.136*log(RoughnessRenolds)+2.30;
end
CoefficientC=-0.993*log(RoughnessRenolds)+12.36;
SpecificWeightofSand=1.8836e4;
SpecificWeightofWater=9.789e3 ;
Di=4e-3;
% Protrusion=1e-3;
Thickness=1.5*D50;
Y1=0.25*Thickness;
Y2(i)=0.25*Thickness+Protrusion(i);
Sum=zeros(N,1);
for m=1:N-1
% syms y;
Kapa=0.4;
ff1=@(y) (Ustar*CoefficientC.*y/Thickness).*sqrt((0.5*Di)^2-(y-Protrusion(i)-Y1+0.5*Di).^2);
ff2=@(y) sqrt((0.5*Di)^2-(y-Protrusion(i)-Y1+0.5*Di).^2);
if Y2(i)<=Thickness
MeanBedVelocity(i)=quad(ff1,Y1,Y2(i))./(quad(ff2,Y1,Y2(i)));
else
MeanBedVelocity(i)=(quad(ff1,Y1,Thickness)+quad(ff1,Thickness,Y2(i)))./(quad(ff2,Y1,Y2(i)));
end if MeanBedVelocity(i)<=Ustar*CoefficientC
Yb(i)=(MeanBedVelocity(i)*Thickness)/(Ustar*CoefficientC);
else
Yb(i)=Thickness*exp(Kapa*(MeanBedVelocity(i)/Ustar-CoefficientC));
end; ParticleRenolds(i)=MeanBedVelocity(i).*Protrusion(i)/KinematicViscosity;
if ParticleRenolds(i)<=1754
Cd(i)=(24./ParticleRenolds(i)).*(1+0.15*ParticleRenolds(i).^0.687);
else
Cd(i)=0.36;
end;
(Protrusion,Friction,Dk,MeanBedVelocity,Yb,Cd, Cl,SpecificWeightofSand,SpecificWeightofWater,MeanShearStress,Ustar,RoughnessRenolds,N,Y1,Di)
h1(i)=Yb(i)-Y1-Protrusion(i)+0.5*Di;
h2=Di*(Friction+0.5*D(m)-0.5*Di)/(Di+D(m));
Ld=h1(i)+h2;
Ll=sqrt((0.5*Di)^2-h2.^2);
Lw=Ll;
Pei=zeros(N,1);
Phi=zeros(N,1);
for j=2:N;
Pei(1)=p(1)*Di/(Di+D(1));
Pei(j)=p(j)*Di/(Di+D(j))+ Pei(j-1);
Phi(j)=1-Pei(j);
end;
HidingFactor=(Pei(N)/Phi(N))^0.6;
EffectiveShearStress=HidingFactor*MeanShearStress;
% SpecificWeightofSand=1.8836e4;
% SpecificWeightofWater=9.789e3 ; %at 20 degree centigrade
DimensionlessEffectiveShearStress= EffectiveShearStress/((SpecificWeightofSand-SpecificWeightofWater)*Di);
ff3=@(y) sqrt((0.5*Di)^2-(y-Protrusion(i)-Y1+0.5*Di).^2);
A(i)=quad(ff3,Y1,Y2(i));
Br=Ustar*sqrt((2*Lw*pi*Di^2)./((Cd(i).*Ld+Cl(i).*Ll)*6.*A(i)*DimensionlessEffectiveShearStress)); % Rolling Threshold
Bl=Ustar*sqrt((2*pi*Di^2)/(Cl(i)*6.*A(i)*DimensionlessEffectiveShearStress)); %Lifting Threshold %syms Ub ; % Instantaneous velocity
% U=((Ub-MeanBedVelocity)/Intensity);% PDF=@(Ub) (exp(-((Ub-MeanBedVelocity)/Intensity).^2/2)/(Intensity*sqrt(2*pi))).*(1+(Skewness/factorial(3))*(((Ub-MeanBedVelocity)/Intensity).^3-3*((Ub-MeanBedVelocity)/Intensity))+(Flatness-3)*(((Ub-MeanBedVelocity)/Intensity).^4-6*((Ub-MeanBedVelocity)/Intensity).^2+3)/factorial(4))/Intensity; %PDF of Velocity Fluctuation
% Pr(i)=quad(PDF,-Bl,-Br)+quad(PDF,Br,Bl);
% Pr(i)=sum(arrayfun(@(BL,BR) quad(PDF, -BL, -BR) + quad(PDF, BR, BL), Bl, Br));
syms Ub ; % Instantaneous velocity
PD(i)=int((exp(-((Ub-MeanBedVelocity)/Intensity).^2/2)/(Intensity*sqrt(2*pi))).*(1+(Skewness/factorial(3))*(((Ub-MeanBedVelocity)/Intensity).^3-3*((Ub-MeanBedVelocity)/Intensity))+(Flatness-3)*(((Ub-MeanBedVelocity)/Intensity).^4-6*((Ub-MeanBedVelocity)/Intensity).^2+3)/factorial(4))/Intensity,-Bl,-Br)...
+int((exp(-((Ub-MeanBedVelocity)/Intensity).^2/2)/(Intensity*sqrt(2*pi))).*(1+(Skewness/factorial(3))*(((Ub-MeanBedVelocity)/Intensity).^3-3*((Ub-MeanBedVelocity)/Intensity))+(Flatness-3)*(((Ub-MeanBedVelocity)/Intensity).^4-6*((Ub-MeanBedVelocity)/Intensity).^2+3)/factorial(4))/Intensity,Br,Bl);
Pr(i)= double(PD(i));
% Pl=1-quad(PDF,-Bl,Bl); %end Sum(1)=p(1)*Pr(i);
Sum(m+1)=p(m)*Pr(i)+Sum(m);
end;
fun=Sum(m+1);

Accepted Answer

Walter Roberson
Walter Roberson on 7 Oct 2011
Looking at your ff2, I do not see any obvious reason why the integral cannot ever come out as zero. If Protrusion(i) can ever be 0, then the bounds of the integral will be identical, leading to an integral of 0.
If you want to avoid division by 0, you need to make sure the denominator cannot possibly be 0 for any useful data range (and you have to test the data to be sure it does not violate the constraint); OR, you need to calculate the denominator first and test whether it is non-zero enough before you go ahead with the division.
If you are looking for a bad-programming-but-it-works short-cut, you could use a try/catch block to substitute some other value if division by 0 is detected.
You really should put in a breakpoint and examine the variables when it happens. If you break the numerator and denominator in to two statements then you can put in a conditional breakpoint that looks something like
dbstop at 123 if abs(ff2_quad) < 1e-5
where ff2_quad had been assigned the denominator.
  1 Comment
Lin LI
Lin LI on 7 Oct 2011
thank you very much, walter, I added the try-catch, and changed the code,trying to replace the denominator to (0+1e-6) if denominator is 0 is detected. Does this work? it still gives me error message
"??? Index exceeds matrix dimensions.
Error in ==> quad at 79
if ~isfinite(y(7))
Error in ==> G_yi at 19
f(i)=quad(@sumin,Friction1,Friction2,[],[],i);
Error in ==> quad at 71
y = f(x, varargin{:});
Error in ==> double_int at 11
SS=quad(@G_yi,Protrusion1,Protrusion2,[],[],Friction1,Friction2);
"
what kind of codes will work?
ff1=@(y) (Ustar*CoefficientC.*y/Thickness).*sqrt((0.5*Di)^2-(y-Protrusion(i)-Y1+0.5*Di).^2);
ff2=@(y) sqrt((0.5*Di)^2-(y-Protrusion(i)-Y1+0.5*Di).^2);
if Y2(i)<=Thickness
% dbstop at 49 if abs(quad(ff2,Y1,Y2(i)))<1e-5;
% display(abs(quad(ff2,Y1,Y2(i))));
MeanBedVelocity(i)=quad(ff1,Y1,Y2(i))./(quad(ff2,Y1,Y2(i)));
try
MeanBedVelocity(i)
catch
MeanBedVelocity(i)=quad(ff1,Y1,Y2(i))./(quad(ff2,Y1,Y2(i))+1e-6);
end
lasterr;
else
MeanBedVelocity(i)=(quad(ff1,Y1,Thickness)+quad(ff1,Thickness,Y2(i)))./(quad(ff2,Y1,Y2(i)));
end

Sign in to comment.

More Answers (2)

William
William on 7 Oct 2011
First off, Use the code tool on the toolbar above when asking questions. It'll make things much easier to answer.
To fix this you first need to figure out where the divide by zero is occurring and where you are exceeding matrix dimensions.
The easiest way to do this (In my opinion) is to use the try and catch statements. This will attempt to run a bit of code and return the error. Here is an example of some code I wrote a while ago
try
rank(v)
catch
fprintf('Oops Won''t work "False"')
err = lasterror;
msg = err.message
end
fprintf('\nPart B')
It will try a bit of code and display the error if it occurs. You can also display warnings on it be it necessary. Hope this helps
  2 Comments
Lin LI
Lin LI on 7 Oct 2011
thank you for your suggestions, but I dont know how to use the code toolbar. should I put the code in {}? or use the numbered list or bulleted list or monospaced, how does it help? Sorry I dont know how to use it.
Walter Roberson
Walter Roberson on 7 Oct 2011
http://www.mathworks.com/matlabcentral/answers/13205-tutorial-how-to-format-your-question-with-markup

Sign in to comment.


cui,xingxing
cui,xingxing on 24 Oct 2023
suppose your denominator is "u" variable, this may be zero,you can do like this,use u+(u==0)*eps instead of u :
somevalues = [-1,0,1,2,3];
for i = 1:length(somevalues)
u = somevalues(i);
y = 1/(u+(u==0)*eps)
end
y = -1
y = 4.5036e+15
y = 1
y = 0.5000
y = 0.3333
If you don't write it that way, like the following:
for i = 1:length(somevalues)
u = somevalues(i);
y = 1/u
end
y = -1
y = Inf
y = 1
y = 0.5000
y = 0.3333
this can lead to unfinite/unbounded.

Community Treasure Hunt

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

Start Hunting!