how to avoid dividing by zeros.

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);

Answer 1

[Walter Roberson]

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.

[Lin LI]

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

Answer 2

[William]

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]

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]

http://www.mathworks.com/matlabcentral/answers/13205-tutorial-how-to-format-your-question-with-markup

Answer 3

[cui,xingxing]

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.