# how to avoid dividing by zeros.

114 views (last 30 days)
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
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;
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
end
f=f(:);
function fun=sumin(Friction,i)
% global Protrusion;
% Protrusion(i)=evalin(Protrusion(i));
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
else
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);
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
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));
Sum(m+1)=p(m)*Pr(i)+Sum(m);
end;
fun=Sum(m+1);

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
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
Error in ==> quad at 71
y = f(x, varargin{:});
Error in ==> double_int at 11
"
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;
try
MeanBedVelocity(i)
catch
end
lasterr;
else
end

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
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 on 7 Oct 2011

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