-1
$\begingroup$

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$

$\endgroup$
  • $\begingroup$ Show us what u have tried. Post the actual code so that we can copy it. See documentation of NIntegrate $\endgroup$ – Hubble07 May 17 '16 at 6:10
  • $\begingroup$ @Hubble07, Please have a look at the Edited Question with the MMA. But It does not work. $\endgroup$ – Srestha Narayanan May 17 '16 at 6:21
  • $\begingroup$ 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. $\endgroup$ – Hubble07 May 17 '16 at 6:38
  • 1
    $\begingroup$ Related / Duplicate: mathematica.stackexchange.com/q/114943/6588 $\endgroup$ – Quantum_Oli May 17 '16 at 6:58
  • $\begingroup$ @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... $\endgroup$ – Srestha Narayanan May 17 '16 at 7:07
2
$\begingroup$

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}
]
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Srestha Narayanan May 17 '16 at 7:58
  • $\begingroup$ 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). $\endgroup$ – Quantum_Oli May 17 '16 at 8:30
  • $\begingroup$ Thank you very much. However, when I run your above MMA code, I still get some error. Would you please check it again? $\endgroup$ – Srestha Narayanan May 17 '16 at 12:45
  • $\begingroup$ 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. $\endgroup$ – Quantum_Oli May 17 '16 at 13:25
  • $\begingroup$ I really appreciate your effort. Thank you very much for all your support. I hope you will be able to fix it. $\endgroup$ – Srestha Narayanan May 17 '16 at 14:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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