how to avoid dividing by zeros.

102 次查看（过去 30 天）
Lin LI
Lin LI 2011-10-7
编辑： cui,xingxing 2024-4-27
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);

请先登录，再回答此问题。

采纳的回答

Walter Roberson
Walter Roberson 2011-10-7
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 个评论
Lin LI
Lin LI 2011-10-7
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

请先登录，再进行评论。

更多回答（2 个）

William
William 2011-10-7
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 个评论
Lin LI
Lin LI 2011-10-7
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 2011-10-7
http://www.mathworks.com/matlabcentral/answers/13205-tutorial-how-to-format-your-question-with-markup

请先登录，再进行评论。

cui,xingxing
cui,xingxing 2023-10-24
编辑：cui,xingxing 2024-4-27
Ran in:
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.
-------------------------Off-topic interlude, 2024-------------------------------
I am currently looking for a job in the field of CV algorithm development, based in Shenzhen, Guangdong, China,or a remote support position. I would be very grateful if anyone is willing to offer me a job or make a recommendation. My preliminary resume can be found at: https://cuixing158.github.io/about/ . Thank you!
Email: cuixingxing150@gmail.com

类别

在 Help Center 和 File Exchange 中查找有关 Matrix Indexing 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by