what functions can be used to do a non-scalar limits of integration?

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;

Answer 1

[Walter Roberson]

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))

[Lin LI]

Thank you very much,Walter. In fact, I want to keep thoes BR,BL remain unchanged in the integration Pr(i)=quad(PDF,-Bl,-Br)+quad(PDF,Br,Bl)

, and then I will integrate them in the next integration in G_yi.m.f(i)=quad(@sumin,Friction1,Friction2,[],[],i) for the variable Friction. that's what I want to do.

I have tried your method, thank you very much, but it tells me that "??? Error using ==> arrayfun

All of the input arguments must be of the same size and shape.

Previous inputs had size 1 in dimension 2. Input #3 has size 2." . so I change the code to the int form. and then messages tells me "??? 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); "

I could find the reason why Index exceeds matrix dimensions.

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

% Sum=@(Protrusion,Friction) p(1)*Pr;

[Lin LI]

and how to avoid dividing by zeros, because it will create infinity.

[Lin LI]

I could not find the reason why Index exceeds matrix dimensions.