Matlab: Error in integral function

Question

I want to compute an integral with the following code written in Matlab.

function f = LO_scheme2a(r,kk)
% Run the Prog by :  q = integral(@(r)LO_scheme2a(r,kk), 0, inf)

%% System parameters 
lambda_m= 10^(-2);
lambda_f= kk.*lambda_m; %can be changed based on scenario
mu= 10^-4; %user density (find the proper user density or scenario based) 
Pdbm_m=43;
P_m=10.^(Pdbm_m./10);
Pdbm_f=20;
P_f=10.^(Pdbm_f./10);
gammadBm_m=0;
gamma_m=10.^(gammadBm_m./10);
gammadBm_f=0;
gamma_f=10.^(gammadBm_f./10);
alpha=4;
k=(P_f/P_m)^(1/alpha);% choose proper value
%k=.4
%% Busy Burst parameters 

Pdbm_req=0; % the optimum one should be calculated
P_req=10.^(Pdbm_req./10);
gammadBm_0=-40;
gamma_0=10.^(gammadBm_0./10);
%%
%% scheme 2
P_t=P_f;
C0=1; %By changing C_0 you can investigate the performance of system 
Theta_0=C0.* P_f;
%%
D1=(P_req./gamma_0).^(1./alpha);
D2=(Theta_0./gamma_0./P_t).^(1./alpha).*r;
D=D2;
lambda_fnew=lambda_f.*exp(-pi.*mu.*D.^2);
%%
PM_u=lambda_m./(lambda_m+k.^2.*lambda_fnew);
PF_u=lambda_fnew.*k.^2./(lambda_m+k.^2.*lambda_fnew);
%%
fm_M=2.*pi.*r.*(lambda_m+k.^2.*lambda_fnew).*exp(-pi.*r.^2.*(lambda_m+k.^2.*lambda_fnew));
ff_F=2.*pi.*r.*(lambda_m./k.^2+lambda_fnew).*exp(-pi.*r.^2.*(lambda_m./k.^2+lambda_fnew));
%%
%% Interference from MBSs to MU an FU
A_i= 1;
S_i=gamma_m .*r.^alpha;
LIm_m=exp(-pi.*lambda_m.*(A_i.*r).^2.*rho1(gamma_m,alpha));
%LIm_m=exp(-pi.*lambda_m.*(A_i.*r).^2.*rho0(S_i./(A_i.*r).^alpha,alpha));
A_i= 1./k;
S_i=gamma_f .*r.^alpha *P_m/P_f;
LIf_m=exp(-pi.*lambda_m.*(A_i.*r).^2.*rho1(gamma_f.*P_m.*k.^alpha./P_f,alpha));
%LIf_m=exp(-pi.*lambda_m.*(A_i.*r).^2.*rho0(S_i./(A_i.*r).^alpha,alpha));

%%
max_x= 10^5;
%% Expectation over h in lemmma 3 :  Interference from FAPs to MU an FU
S_i=gamma_f .*r.^alpha .*P_f./P_m;
z=(P_req/gamma_0).^(1/alpha);
%MIm_f=(exp(-h).*exp(-mu.*pi.*z^2.*h.^(2./alpha))).*((-((z.^2).*h^(2/alpha)).*(1-exp(-S_i.*z.^(-alpha))))+(S_i.*h).^(2./alpha).*real(gammainc(1-(2./alpha),S_i.*z.^(-alpha))));
for idx = 1:numel(r)
    MIm_f(idx)=  integral(@(h)exp(-h).*((exp(-mu.*pi.*z^2.*h.^(2./alpha))).*((-((z.^2).*h.^(2./alpha)).*(1-exp(-S_i.*z.^(-alpha))))+(S_i.*h).^(2./alpha).*real(gammainc(1-(2./alpha),S_i.*z.^(-alpha))))),0,max_x);
end
% % 
S_i=gamma_f .*r.^alpha ;
z=(P_req./gamma_0).^(1/alpha);

for idx = 1:numel(r)
    MIf_f(idx)=  integral(@(h)(exp(-h).*exp(-mu.*pi.*z^2.*h.^(2./alpha))).*((-((z.^2).*h.^(2./alpha)).*(1-exp(-S_i.*z.^(-alpha))))+(S_i.*h).^(2./alpha).*real(gammainc(1-(2./alpha),S_i.*z.^(-alpha)))),0,max_x);
end
%%

f=PM_u.*LIm_m.*exp(-2.*pi.*lambda_f.*MIm_f).*fm_M+PF_u.*LIf_m.*exp(-2.*pi.*lambda_f.*MIf_f).*ff_F;

where

function [ f_rho ] = rho1( x,alpha)
f_rho=x.^(2./alpha).*(pi./2-atan(x.^(-2./alpha)));
end

I can run it for various values of kk, but for some kk it gives me an error. For example, when kk=10, I can run this program with the command

 q = integral(@(r)LO_scheme2a(r,kk), 0, inf)

but for kk=100, it gives me an error.

(Error using .* Matrix dimensions must agree.)

Would you please help me to find the source of this error?

Answer 1

When you call q = integral(@(r)LO_scheme2a(r,kk), 0, inf), the integral function choosesr`, which is a vector that discretizes the integration bounds. The documentation states:

For scalar-valued problems, the function y = fun(x) must accept a vector argument, x, and return a vector result, y. This generally means that fun must use array operators instead of matrix operators. For example, use .* (times) rather than * (mimes).

Your error(s) occur where you call integral again from within your code. At this point r is a vector (a vector that can change size) and you don't have a scalar-valued problem. You have tried to deal with this with your two for loops. However, the right hand sides of both loops don't depend on your index variable idx, i.e., r(idx).

However, you don't even need for loops as the documentation above continues:

If you set the 'ArrayValued' option to true, fun must accept a scalar and return an array of fixed size.

Thus, using the 'ArrayValued' option, your for loops can be replaced with:

...
fun = @(h,S_i,z)exp(-h).*exp(-mu*pi*z^2*h.^(2/alpha)) ...
                       .*(-z^2*h.^(2/alpha).*(1-exp(-S_i*z^-alpha)) ...
                       +(S_i.*h).^(2/alpha).*real(gammainc(1-2/alpha,S_i*z^-alpha)));
MIm_f = integral(@(h)fun(h,S_i,z),0,max_x,'ArrayValued',true);

S_i = gamma_f.*r.^alpha ;
z = (P_req./gamma_0).^(1/alpha);

MIf_f = integral(@(h)fun(h,S_i,z),0,max_x,'ArrayValued',true);
...

I've also cleaned up your anonymous function and added two extra parameters so you can re-use it (see here).