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