I want to do an integration with non scalar limits of integration, quad and int seem not working,because Bl and Br are vectors. can anyone help me fix this problem?I will appreciate your help.
below are the codes, one main code and 3 function code . it gave message
"??? Error using ==> quadThe limits of integration must be scalars.Error in ==> sumin at 142 Pr(i)=quad(PDF,-Bl,-Br)+quad(PDF,Br,Bl);Error in ==> quad at 71y = f(x, varargin{:});" 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);
% endsave G_yi.mat;for i=1:nf(i)=quad(@sumin,Friction1,Friction2,[],[],i);endf=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;elseIntensity=Ustar*(-0.187*log(RoughnessRenolds)+2.93);Skewness=0.102*log(RoughnessRenolds);Flatness=0.136*log(RoughnessRenolds)+2.30;endCoefficientC=-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); % 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); % Sum=@(Protrusion,Friction) p(1)*Pr;
Best Answer