I run this code and it gave me error message, and it seems that this problem comes from the function quadl, because the variable y never appear in my code, y is from the function quadl.I set some breakpoint,still couldnt find why.If anyone can help fix this problem , i appreciate it very much. anyone can help find what the problem is. I am really annoyed by this problem.thank you so much!
here is the error:
"??? Attempted to access y(13); index out of bounds because numel(y)=1.Error in ==> quadl at 72if ~isfinite(y(13))Error in ==> G_yi at 20f(i)=quadl(@sumin,Friction1,Friction2,[],[],i);Error in ==> quadl at 64y = feval(f,x,varargin{:}); y = y(:).';Error in ==> double_int at 11SS=quadl(@G_yi,Protrusion1,Protrusion2,[],[],Friction1,Friction2); Error in ==> main at 4PR=double_int(0,Di,0,Di) "
here is the code, one main code and 3 functions.
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=quadl(@G_yi,Protrusion1,Protrusion2,[],[],Friction1,Friction2); function f=G_yi(Protrusion,Friction1,Friction2) Protrusion=Protrusion(:);n=length(Protrusion); save G_yi.mat; f=zeros(n,1); for i=1:n f(i)=quadl(@sumin,Friction1,Friction2,[],[],i); end f=f(:); function fun=sumin(Friction,i) load G_yi.mat; MeanShearStress=5 ; 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; 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; Thickness=1.5*D50; Y1=0.25*Thickness; Y2(i)=0.25*Thickness+Protrusion(i); Sum=zeros(N,1); for m=1:N-1 Kapa=0.4; ff1=@(t) (Ustar*CoefficientC.*t/Thickness).*sqrt((0.5*Di)^2-(t-Protrusion(i)-Y1+0.5*Di).^2); ff2=@(t) sqrt((0.5*Di)^2-(t-Protrusion(i)-Y1+0.5*Di).^2); MeanBedVelocity=zeros(n,1); MM=zeros(n,1); MM(i)= quadl(ff2,Y1,Y2(i)); if Y2(i)<=Thickness MeanBedVelocity(i)=quadl(ff1,Y1,Y2(i))./(MM(i)); else MeanBedVelocity(i)=(quadl(ff1,Y1,Thickness)+quadl(ff1,Thickness,Y2(i)))./(quadl(ff2,Y1,Y2(i))); end Yb=zeros(n,1); if MeanBedVelocity(i)<=Ustar*CoefficientC Yb(i)=(MeanBedVelocity(i)*Thickness)/(Ustar*CoefficientC); else Yb(i)=Thickness*exp(Kapa*(MeanBedVelocity(i)/Ustar-CoefficientC)); end; ParticleRenolds=zeros(n,1); ParticleRenolds(i)=MeanBedVelocity(i).*Protrusion(i)/KinematicViscosity; Cd=zeros(n,1); if ParticleRenolds(i)<=1754 Cd(i)=(24./ParticleRenolds(i)).*(1+0.15*ParticleRenolds(i).^0.687); else Cd(i)=0.36; end; Cl=zeros(n,1); Cl(i)=Cd(i); h1=zeros(n,1); 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; DimensionlessEffectiveShearStress= EffectiveShearStress/((SpecificWeightofSand-SpecificWeightofWater)*Di); ff3=@(t) sqrt((0.5*Di)^2-(t-Protrusion(i)-Y1+0.5*Di).^2); A=zeros(n,1); A(i)=quadl(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 ; Pr=zeros(n,1); PD=zeros(n,1); PD(i)=int((exp(-((Ub-MeanBedVelocity(i))/Intensity).^2/2)/(Intensity*sqrt(2*pi))).*(1+(Skewness/factorial(3))*(((Ub-MeanBedVelocity(i))/Intensity).^3-3*((Ub-MeanBedVelocity(i))/Intensity))+(Flatness-3)*(((Ub-MeanBedVelocity(i))/Intensity).^4-6*((Ub-MeanBedVelocity(i))/Intensity).^2+3)/factorial(4))/Intensity,-Bl,-Br)... +int((exp(-((Ub-MeanBedVelocity(i))/Intensity).^2/2)/(Intensity*sqrt(2*pi))).*(1+(Skewness/factorial(3))*(((Ub-MeanBedVelocity(i))/Intensity).^3-3*((Ub-MeanBedVelocity(i))/Intensity))+(Flatness-3)*(((Ub-MeanBedVelocity(i))/Intensity).^4-6*((Ub-MeanBedVelocity(i))/Intensity).^2+3)/factorial(4))/Intensity,Br,Bl); Pr(i)= double(PD(i)); Sum(1)=p(1)*Pr(i); Sum(m+1)=p(m)*Pr(i)+Sum(m); end; fun=Sum(m+1);
Best Answer