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?
rho1
. Also complete error messages including line numbers are much more helpful. Have you tried debugging? Finally, StackOverflow would probably be a more appropriate place for this question. $\endgroup$