what functions can be used to do a non-scalar limits of integration?
Show older comments
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 ==> quad
The 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 71
y = 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);
% 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);
% 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;
Accepted Answer
Walter Roberson
on 6 Oct 2011
None of the quad* routines accept vectors for their integration limits.
I notice you assign the result of the quad() to a single scalar variable, which suggests that you are expecting a scalar result even though you would like to specify a vector of bounds. How did you want to aggregate those various results? Just add them together? Do quad integration across them?
Perhaps what you want is something like
sum(arrayfun(@(BL,BR) quad(PDF, -BL, -BR) + quad(PDF, BR, BL), Bl, Br))
More Answers (0)
Categories
Find more on Matrix Indexing in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
