1
$\begingroup$

I have written a MATLAB code as following, in order to calculate an integration involving Modified Bessel Function with variable limits:

function num_exact_JD()
    rho = 0.5 * ones(1, 5); 
    N = 4; 
    sigma = ones(1, 5); 

    Omega = sigma .* sqrt((1 - rho.^2) / 2); 
    R_c = 0.5; 

    sys_SNR_dB = -5:2:25; 
    sys_SNR = 10.^(sys_SNR_dB / 10); 

    out_probs = zeros(1, size(sys_SNR, 2)); 

    for i = 1:size(sys_SNR, 2)
    out_probs(i) = outage_prob(rho, sys_SNR(i), Omega, sigma, N, R_c);
    end

end

function res = phi(r)
    res = log10(r + 1) / log(2);
end

function res = outage_prob(rho, snr, omega, sigma, N, R_c)

    fac = snr / N;

    f = @(om, gam, sig, R, lam) 1/(2*om^2*fac) * ...
        exp(-(gam/fac+(lam*sig*R).^2)./(2*om^2)) .* ...
        besseli(0, sqrt(gam/fac)*lam*sig*R./om^2); 

    f1 = @(g4, g3, g2, g1, R) f(omega(1), g1, sigma(1), R, rho(1)) .* ...
        f(omega(2), g2, sigma(2), R, rho(2)) .* ...
        f(omega(3), g3, sigma(3), R, rho(3)) .* ...
        f(omega(4), g4, sigma(4), R, rho(4));

    m1 = @(g3, g2, g1, R) integral(@(g4) f1(g4, g3, g2, g1, R),...
        0, 2.^(R_c-phi(g1)-phi(g2)-phi(g3))-1, 'ArrayValued',true);

    m2 = @(g2, g1, R) integral(@(g3) m1(g3, g2, g1, R),...
        0, 2.^(R_c-phi(g1)-phi(g2))-1, 'ArrayValued',true);

    m3 = @(g1, R) integral(@(g2) m2(g2, g1, R),...
        0, 2.^(R_c-phi(g1))-1, 'ArrayValued',true);

    m4 = @(R) integral(@(g1) m3(g1, R), 0, 2.^(R_c)-1, 'ArrayValued',true) 
        *2.*exp(-R.^2).*R;

    res = integral(@(R) m4(R), 0, 20);

end

It costs lots of time to execute, so I'd like to convert it to Mathematica and check whether it'll run faster. The following is the corresponding Mathematica code:

num = 4;
rho = 0.5*Table[1, num];
sigma = Table[1, num];
omega = sigma*Sqrt[0.5*(1 - rho^2)];
r = 0.5;
snr = 10^(Table[i, {i, -5, 25, 2}]/10);
probs = Table[0, 16];

phi[r_] := Log[r + 1]/Log[2];

fac = snr[1]/num;

f[om, gam, sig, y, lam] := 1/(2*om^2*fac)*Exp[-(gam/fac + 
   (lam*sig*y)^2)/(2*om^2)]*BesselI[0, Sqrt[gam/fac*lam*sig*y/om^2]];

f1[g4, g3, g2, g1, y] := f[omega[[1]], g1, sigma[[1]], y, rho[[1]]]*
   f[omega[[2]], g2, sigma[[2]], y, rho[[2]]]*
   f[omega[[3]], g3, sigma[[3]], y, rho[[3]]]*
   f[omega[[4]], g4, sigma[[4]], y, rho[[4]]];

m1[y] := Integrate[f1, {g4, 0, 2^(r - phi[g1] - phi[g2] - phi[g3]) - 1}, 
   {g3, 0, 2^(r - phi[g1] - phi[g2]) - 1}, {g2, 0, 2^(r - phi[g1]) - 1}, 
   {g1,0, 2^(r) - 1}];

m4 = Integrate[m1*2*Exp[-y^2]*y, {y, 0, Infinity}]

First, I don't know how to convert the for-loop part, so in the Mathematica code I just write down the first iteration. Moreover, I know the multiple dimensional integration is wrong, but I get stuck in this.

Could anyone check my Mathematica code and give me some suggests? Thanks very much indeed!

$\endgroup$
0

1 Answer 1

3
$\begingroup$

This is how the code looks after some adjustments. In some of the function declarations, you missed the underscore _. Moreover, I changed the argument pattern of m1 to y_?NumericQ in order to prevent the NIntegrate in the definition of m4 from trying to evaluate m1[y] symbolically. Note that I used NIntegrate, not Integrate. Moreover, Method -> "GaussKronrodRule" speeds up the m1 (m1[1.] needs about 2.7 seconds to evaluate). I haven't tried to evaluate m4. It will take very long.

num = 4;
ρ = Table[0.5, num];
σ = Table[1., num];
ω = σ Sqrt[0.5 (1 - ρ^2)];
r = 0.5;
snr = 10.^(Table[i, {i, -5, 25, 2}]/10);
probs = Table[0, 16];
ϕ[r_] := Log[r + 1]/Log[2];

fac = snr[[1]]/num;

f[ω_, γ_, σ_, y_, λ_] :=  1/(2 ω^2 fac) Exp[-(γ/fac + (λ σ y)^2)/(2 ω^2)] BesselI[0, Sqrt[γ/fac λ σ y/ω^2]];
f1[g_, y_] := Times @@ MapThread[f[#1, #2, #3, y, #4] &, {ω, g, σ, ρ}];

m1[y_?NumericQ] := NIntegrate[
   f1[{g1, g2, g3, g4}, y],
   {g1, 0, 2^(r) - 1},
   {g2, 0, 2^(r - ϕ[g1]) - 1},
   {g3, 0, 2^(r - ϕ[g1] - ϕ[g2]) - 1},
   {g4, 0, 2^(r - ϕ[g1] - ϕ[g2] - ϕ[g3]) - 1},
   Method -> "GaussKronrodRule"
   ];

m4 = NIntegrate[m1[y] 2 Exp[-y^2]*y, {y, 0, ∞}]
$\endgroup$
1
  • $\begingroup$ Thanks for your help. However, when I run your code, the output is NIntegrate[m1[y] 2 Exp[-y^2] y, {y, 0, 15}] rather than a certain value. And a Warning would come out that pointed out the calculation of integrand will get a NaN. Does the problem locate in that the multi-dimensional integral order was reversed? $\endgroup$
    – yuhou CHEN
    Jun 15, 2018 at 20:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.