积分函数误差

Question

我想通过 MATLAB 求解以下积分。我写了以下代码

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

    %% System parmaeters 
    lambda_m= 10^(-2);
    lambda_f= kk.*lambda_m; %can be change based on senario
    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);% chose proper value

    Pdbm_req=0; % the optimum one should be calculate
    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));
    %%
    %% Interfernce 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;
    %% Expectiation over h in lemmma 3 :  Interfernce 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;

当

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

我可以为不同的 kk 运行它，但对于某些 kk，它给了我错误。例如当 kk=10 时，我可以通过

运行这个程序

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

但是对于 kk=100，它给了我一些错误（错误使用 .* 矩阵尺寸必须一致。）

你能帮我解决这个错误吗？谢谢

Answer 1

问题在于行

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扩展为水平向量，内积分将与h做同样的事情，这意味着只有当两个大小恰好相同时才能进行乘法运算，并给出错误的结果。我会对h进行转置，让它在垂直维度上扩展，并将乘法与bsxfun而不是.*匹配，如下所示：

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

对两个内积分做同样的事情。

注意：我没有integral 函数，但这种方法对quadv 有效（您需要一个能够处理向量函数的集成命令）。