# Can Mathematica solve this integral numerically?

I need to evaluate numerically the following

$$\int_0^{\infty}\exp\left(-\pi\lambda\mathbb{E}_H\left(\int_r^\infty\left(1-\exp\left(-sP_{\text{M}}Hr^{-1/\delta}\right)\right){\rm{d}}r\right)\right)2\pi\lambda r\exp\left(-\pi\lambda r^2\right)dr$$

with $s=\mu r^4/P_{\text{M}}$

The PDF of $H$ is given by

$f_H(h)=\frac{K+1}{\tau}\exp\left(-K-(K+1)h/\tau\right)I_0\left(\sqrt{4K(K+1)h/\tau}\right)$

Here $I_0$ is the 0-th order Bessel function of first kind

Let, $P_{\text{M}}=20$, $\mu=10$, $\lambda=20$, $K=1$, $\tau=1$, $\delta=0.5$

• Show us what u have tried. Post the actual code so that we can copy it. See documentation of NIntegrate – Hubble07 May 17 '16 at 6:10
• @Hubble07, Please have a look at the Edited Question with the MMA. But It does not work. – Srestha Narayanan May 17 '16 at 6:21
• Have you defined your int correctly. Is the integration limit correct in int. NExpectation will work only when int is just a function of h but you have made it a function of both h and r. – Hubble07 May 17 '16 at 6:38
• Related / Duplicate: mathematica.stackexchange.com/q/114943/6588 – Quantum_Oli May 17 '16 at 6:58
• @Hubble07, I believe int is defined correctly. Note that I get expected result when $K$=0. For other values of $K$, I do not get results. Note that I have changed NExpectation to Expectation... – Srestha Narayanan May 17 '16 at 7:07

Most recent edit

This seems to work (tested exactly what is posted here in a fresh kernel), however it does complain a little.

α = 4;
μ = 10;
k = 1;
τ = 1;
δ = 1/2;

s[r_] := (μ r^α)/Pm

f[h_] = FullSimplify[
Exp[2 k] (k + 1)/τ Exp[-k - (k + 1) h/τ] BesselJ[0, Sqrt[4 k (k + 1) h/τ]]
];

int[r_, h_] = Integrate[(1 - Exp[-s[r] Pm h rd^(-(1/δ))]), {rd, r, Infinity}, Assumptions -> {r > 0, h >= 0}];

AverageProbSuccess[λ_] := NIntegrate[
Exp[-π λ NIntegrate[int[r, h] f[h], {h, 0, Infinity}]] 2 π λ r Exp[-π λ r^2],
{r, 0, Infinity}]

AverageProbSuccess[20]


0.979115

None of the below actually appears to run... Should I delete it?

Edit

It is unclear in your question (or at least was initially to me) which r is which at certain points, often dummy variables are used for integration to avoid ambiguities like this. My s[rd] should have been s[r] which gives AverageProbSuccess[λ = 20, k = 1] = 0.979.

Pm = 20;
μ = 10;
τ = 1;
δ = 1/2;

s[r_] := (μ r^4)/Pm

f[h_, k_] := Exp[2 k] (k + 1)/τ Exp[-k - (k + 1) h/τ] BesselJ[0, Sqrt[4 k (k + 1) h/τ]]

int[r_, h_] = Integrate[(1 - Exp[-s[r] Pm h rd^(-(1/δ))]), {rd, r, Infinity}, GenerateConditions -> False]

AverageProbSuccess[λ_, k_] := NIntegrate[
Exp[-π λ NIntegrate[int[r, h] f[h, k], {h, 0, Infinity}]] 2 π λ r Exp[-π λ r^2],
{r, 0, Infinity}
]


Old, Unchanged, Incorrect

Do you know what result you expect? For the values you give this spits out AverageProbSuccess = 4.13 * 10^33:

Pm = 20;
μ = 10;
λ = 20;
k = 1;
τ = 1;
δ = 1/2;

s[r_] := (μ r^4)/Pm

f[h_] := (k + 1)/τ Exp[-k - (k + 1) h/τ] BesselJ[0, Sqrt[4 k (k + 1) h/τ]]

A = Integrate[f[h], {h, 0, Infinity}]

int[r_, h_] = Integrate[(1 - Exp[-s[rd] Pm h rd^(-(1/δ))]), {rd, r, Infinity}, GenerateConditions -> False]

AverageProbSuccess = NIntegrate[Exp[-π λ NIntegrate[int[r, h] 1/A f[h], {h, 0, Infinity}]] 2 π λ r Exp[-π λ r^2],
{r, 0, Infinity}
]

• Thanks for your effort. I expect that "0<AverageProbSuccess<1" . Just for checking, putting $K=0$, you get exponential distribution and the MMA code I provided works perfectly. When $K$ is not equal to zero, its a different distribution for $h$, where $h$ is the square of a Rician distributed variable. – Srestha Narayanan May 17 '16 at 7:58
• I guess the clue is in the name that AverageProbSuccess should likely be between 0 and 1. But I don't think your PDF f[h] is normalised for general k, which is why I calculate and include a normalisation constant A. You need a additional prefactor of E^(2 k). – Quantum_Oli May 17 '16 at 8:30
• Thank you very much. However, when I run your above MMA code, I still get some error. Would you please check it again? – Srestha Narayanan May 17 '16 at 12:45
• You're right, not sure what's going on. I must have had something defined in the background, I thought I had done a check on a clean kernel. I'll try to fix it, otherwise I'll delete my (not) answer. – Quantum_Oli May 17 '16 at 13:25
• I really appreciate your effort. Thank you very much for all your support. I hope you will be able to fix it. – Srestha Narayanan May 17 '16 at 14:19