Is this integral integrable?

  • #1
1,367
61
After following the logical steps to derive something, I reached to the following integral:

[tex]\int_0^{\infty}\frac{1}{x^2}\exp\left(\frac{A}{x}-x\right)E_1\left(B+\frac{A}{x}\right)\,dx[/tex]

where ##E_1(.)## is the exponential integral function, and ##A## and ##B## are non-zero positive constants. I tried to find this integral in the table of integrals, but couldn't find it. Next best thing to do is to evaluate this integral numerically, I guess. So, I resorted to Mathematica for that purpose, and used the Nintegrate command, but I got a message that the integrand has evaluated to overflow, indeterminate, or infinity. It suggested to increase the number of recursive refinements, but it didn't work. Is there something fundamentally wrong in the integral, or it's just a technical issue, and how to overcome it? I have to mention that when I replaced the lower bound of the integral by a very small value that is greater than 0, the integral was calculated perfectly. Is this the problem maybe?

Thanks
 

Answers and Replies

  • #2
347
134
Many symbolic algebra systems have this issue with their numerical functions. When there are [itex]\frac{1}{x}[/itex] terms and such inside the expression, even if the limit of the entire expression is finite, the numerical algorithm fails early when evaluating the expression from inside out. There are slightly kludgy looking solutions that usually work. I personally like to add a while-type conditional to the expression which forces it to evaluate to the proper limit at the points in question.
 
  • #3
1,367
61
Many symbolic algebra systems have this issue with their numerical functions. When there are [itex]\frac{1}{x}[/itex] terms and such inside the expression, even if the limit of the entire expression is finite, the numerical algorithm fails early when evaluating the expression from inside out. There are slightly kludgy looking solutions that usually work. I personally like to add a while-type conditional to the expression which forces it to evaluate to the proper limit at the points in question.

How do you add "a while-type conditional to the expression"?
 
  • #4
347
134
It depends on the particular CAS you're using. It's the same as an if-then-else structure if that helps.

For example, in some CAS languages, while(x=0,0,exp(-1/x^2)) will evaluate to 0 when x is 0 and to exp(-1/x^2) everywhere else. Look up how to do piecewise functions in your language, it's the same thing. In this case, one of the pieces of the domain is a single point instead of an interval, is all.
 
  • #5
George Jones
Staff Emeritus
Science Advisor
Gold Member
7,544
1,385
In order to evaluate the integral numerically, what values of A and B do you use?
 
  • #6
1,367
61
In order to evaluate the integral numerically, what values of A and B do you use?

I need to evaluate the integral at different values of A and B, but let's say A=50 and B=5 as one of the cases.
 
  • #8
177
61
The integral diverges at ##0##, the integrand behaves near ##0## as ##1/(Ax)##, which is not integrable.
 
  • #9
1,367
61
Can we approximate it such that it becomes integrable?
 
  • #10
177
61
What do you mean? Your function is positive, non integrable, so integral is ##+\infty##, there is no other reasonable interpretation. For non-integrable functions than change sign, there are some possibilities to define a finite value of the integral, for example, the principal value integral, but for positive functions there is no alternative interpretation, the integral is ##+\infty##.

As for the "approximation" you can change your integrand to get a finite number, but this will not have any relation to the original integral.
 
  • #11
1,367
61
OK. Now if I followed the mathematical derivations for something correctly, such that this something is between 0.5 and 1 for sure (I can get this using Monte-Carlo simulations), is it possible to reach to such an integral that is not integrable?
 
  • #12
177
61
You probably made mistake in your calculations.
 
  • #13
1,367
61
You probably made mistake in your calculations.

OK. Could I walk you through the derivation. Would you patient with me?
 
  • #15
1,367
61
OK, I'll try to follow.

OK. perfect. Now I have something I need to finish. Tomorrow, I will start the derivation from the beginning to see where I erred. Thanks
 
  • #16
1,367
61
Hi. Sorry, I was busy yesterday. Basically, I have this expession

[tex]\varepsilon(\alpha_1,\,\alpha_2,\,\alpha_3)=\exp\left(-\frac{\frac{\alpha_1}{\alpha_2}A}{\alpha_3 B+1}\right)[/tex]

where ##\alpha_i## for i=1, 2, 3 are independent and identically distributed exponential random variables with mean 1, and A and B are non-zero positive constants. I need to find the average value of ##\varepsilon(\alpha_1,\,\alpha_2,\,\alpha_3)##. To do so, I let ##X=\frac{\frac{\alpha_1}{\alpha_2}A}{\alpha_3 B+1}## and then found the cumulative distribution function (CDF) of X denoted by ##F_X(x)##, and then found its derivative with respect to ##x## to find the probability density function (PDF), denoted by ##f_X(x)##. After finding the PDF, the average value of ##\varepsilon(\alpha_1,\,\alpha_2,\,\alpha_3)## can be found as

[tex]\varepsilon=\int_0^{\infty}e^{-x}\,f_X(x)\,dx[/tex]

I will stop here to see everything I did so far is correct, and to make sure you're following. I'll provide the details of ##F_X(x)## and ##f_X(x)## in the next post.
 
  • #17
177
61
The expectation should be $$\int_{-\infty}^\infty x f_X(x) dx,$$ see Wikipedia. Your quantity is positive, so ##f_X(x)=0## for ##x<0##, so you can integrate only from ##0## to ##\infty##, but ##e^{-x}## is wrong.
 
  • #18
1,367
61
But I need to find the expected value of ##\varepsilon## not ##X##!! Right?
 
  • #19
177
61
Yes, you are right. Sorry, my bad here.
 
  • #20
1,367
61
OK, perfect. So first I need to find the CDF as

[tex]F_X(x)=\text{Pr}\left[X\leq x\right]=\text{Pr}\left[\frac{\frac{\alpha_1}{\alpha_2}A}{\alpha_3 B+1}\leq x\right]=\int_{\alpha_2=0}^{\infty}\int_{\alpha_3=0}^{\infty}\text{Pr}\left[\alpha_1\leq\frac{x\alpha_2}{A}\left(\alpha_3 B+1\right)\right]f_{\alpha_2}(\alpha_2)f_{\alpha_3}(\alpha_3)\,d\alpha_2d\alpha_3=\int_{\alpha_2=0}^{\infty}\int_{\alpha_3=0}^{\infty}\left[1-\exp\left(-\frac{x\alpha_2}{A}\left(\alpha_3 B+1\right)\right)\right]e^{-\alpha_2}e^{-\alpha_3}d\alpha_2d\alpha_3=1-\frac{A}{xB}\int_{\alpha_3=0}^{\infty}\frac{e^{-\alpha_3}}{\alpha_3+\frac{x+1}{B}}\,d\alpha_3[/tex]

From the table of integral we have

[tex]\int_0^{\infty}\frac{e^{-\mu x}}{x+\beta}\,dx=-e^{\beta\mu}\text{Ei}(-\beta\mu)[/tex]

where ##\text{Ei}## is the exponential integral function. But then I used the property that ##E_1(z)=-\text{Ei}(-z)##, which results in

[tex]F_X(x)=1-\frac{A}{xB}\exp((x+1)/B)E_1((x+1)/B)[/tex].

Is everything OK so far?
 
Last edited:
  • #21
177
61
I am not following the last identity in the chain. Didn't you forget term ##e^{-\alpha_2}## when you integrated?
 
  • #22
1,367
61
I am not sure which one you mean, but I have double integrals; one over ##\alpha_2## and one over ##\alpha_3##. The last integral is after integrating the integrad over ##\alpha_2## which is a straightforward integral. This will result in one final integral over ##\alpha_3##. Did I answer your question?
 
  • #23
177
61
I mean the last identity in this chain:
$$F_X(x)=\text{Pr}\left[X\leq x\right]=\text{Pr}\left[\frac{\frac{\alpha_1}{\alpha_2}A}{\alpha_3 B+1}\leq x\right]=\int_{\alpha_2=0}^{\infty}\int_{\alpha_3=0}^{\infty}\text{Pr}\left[\alpha_1\leq\frac{x\alpha_2}{A}\left(\alpha_3 B+1\right)\right]f_{\alpha_2}(\alpha_2)f_{\alpha_3}(\alpha_3)\,d\alpha_2d\alpha_3=\int_{\alpha_2=0}^{\infty}\int_{\alpha_3=0}^{\infty}\left[1-\exp\left(-\frac{x\alpha_2}{A}\left(\alpha_3 B+1\right)\right)\right]e^{-\alpha_2}e^{-\alpha_3}d\alpha_2d\alpha_3=1-\frac{A}{xB}\int_{\alpha_3=0}^{\infty}\frac{e^{-\alpha_3}}{\alpha_3+\frac{x+1}{B}}\,d\alpha_3$$
 
  • #24
1,367
61
OK, let me write it this way:

[tex]
\int_{0}^{\infty}\underbrace{\left[\int_{0}^{\infty}\left[1-\exp\left(-\frac{x\alpha_2}{A}\left(\alpha_3 B+1\right)\right)\right]e^{-\alpha_2}\,d\alpha_2\right]}_{\frac{1}{\frac{x}{A}(\alpha_3 B+1)+1}}e^{-\alpha_3}\,d\alpha_3
[/tex]

After some mathematical manipulations we reached to the final identity. Is it clear now?
 
  • #25
177
61
I see that $$ \int_{\alpha_2=0}^{\infty}\int_{\alpha_3=0}^{\infty} e^{-\alpha_2}e^{-\alpha_3} d\alpha_2d\alpha_3=1.$$
And I see that $$\int_{\alpha_2=0}^{\infty} \exp\left(-\frac{x\alpha_2}{A}\left(\alpha_3 B+1\right)\right)e^{-\alpha_2} d\alpha_2 = \left[ \frac{x}{A} (\alpha_3 B +1) + 1\right]^{-1}. $$ But I do not follow the algebra you use to get your result. And it looks like the function
$$ \left[ \frac{x}{A} (\alpha_3 B +1) + 1\right]^{-1}$$ and your integrand
$$\frac{e^{-\alpha_3}}{\alpha_3+\frac{x+1}{B}}$$ has poles at different values of ##\alpha_3##.
 
  • #26
177
61
From my calculations I get that $$F_X(x) = 1-\int_{\alpha_3=0}^\infty \left[ \frac{x}{A} (\alpha_3 B +1) + 1\right]^{-1} e^{-\alpha_3} d\alpha_3. $$ And I do not see how to get your integral from here.
 
  • #27
177
61
Yes, your formula
$$F_X(x)=1-\frac{A}{xB}\int_{\alpha_3=0}^{\infty}\frac{e^{-\alpha_3}}{\alpha_3+\frac{x+1}{B}}\,d\alpha_3$$ is wrong, because for small ##x## it will give you negative values. Namely, if you fix ##A,B>0##, then when ##x## is decreasing your integral $$\int_{\alpha_3=0}^{\infty}\frac{e^{-\alpha_3}}{\alpha_3+\frac{x+1}{B}}\,d\alpha_3>0$$ increases. But you also have a factor ##\frac{A}{xB}## in front of the integral, so $$\lim_{x\to 0+} \frac{A}{xB}\int_{\alpha_3=0}^{\infty}\frac{e^{-\alpha_3}}{\alpha_3+\frac{x+1}{B}}\,d\alpha_3 =+\infty,$$ so for small ##x>0## you will get ##P_X(x)<0##.
 
  • #28
1,367
61
I got it like this

[tex]\frac{x}{A}(\alpha_3B+1)+1=\frac{x}{A}\alpha_3B+\frac{x}{A}+1=\frac{xB}{A}(\alpha_3+\underbrace{\frac{1}{B}+\frac{A}{xB}}_{\frac{x+A}{xB}})[/tex]

So, yes instead of 1 it should be A. Why it's wrong again? By the way I did the Monte-Carlo Simulations and compared it with the above numerical expression, and both gave the same curves. No infinities!!
 
  • #29
177
61
Yes, it should be ##A## instead of ##1##, and also there should be ##xB## instead of ##B## in the denominator.
 
  • #30
1,367
61
Right. My bad. Are we good again? If yes, I will continue tomorrow the PDF ##f_X(x)##, where, if there is any error, it will be there. Thanks for following the derivations.

PS.: I used different notations here for simplicity, that's why I sometimes get something wrong as a typo. But in my original derivations they are correct.
 
  • #31
177
61
This might simplify your calculations a bit: integrating by parts you get $$\int_0^\infty e^{-x} f_X(x) dx = \int_0^\infty F_X(x) e^{-x} dx, $$ so you do not need to compute derivative of ##F_X##.
 
  • #32
1,367
61
Oh, really? So, I think you can see now where the integral in my first post came from. I will follow on this.
 
  • #33
1,367
61
Integration by parts goes like this. Let ##u=e^{-x}## and thus ##du=-e^{-x}\,dx##, and ##dv=f_X(x)\,dx##, and thus ##v=\int f_X(x)\,dx##. I know that ##F_X(x)=\int_0^xf_X(x)\,dx##, but does the above ##v=F_X(x)## and why? I will continue. Assuming the above is correct, we have

[tex]\int_0^{\infty}e^{-x}f_X(x)\,dx=\underbrace{\left. e^{-x}F_X(x)\right|_0^{\infty}}_{=0}+\int_0^{\infty}e^{-x}F_X(x)\,dx=\int_0^{\infty}e^{-x}F_X(x)\,dx[/tex].

Interesting!!
 
  • #34
1,367
61
Based on the above derivations, the average value of ##\varepsilon(\alpha_1,\alpha_2,\alpha_3)## is given by

[tex]\varepsilon=1-\frac{A}{B}\int_0^{\infty}x^{-1}\exp\left[\frac{x+A}{xB}-x\right]E_1\left[\frac{x+A}{xB}\right]\,dx[/tex]

Is the above integral inetgrable now?
 
  • #35
1,367
61
OK, now the integral is integrable, and Mathematica gives no complains. But when I compared the numerical results with Monte-Carlo simulations, where I generated ##10^5## samples for each point, I got very close but exactly the same curves. See attached Figure. My original equation is ##\varepsilon(\alpha_1,\alpha_2\alpha_3)=0.5\exp\left(-\frac{\frac{\alpha_1}{\alpha_2}G\gamma_Q}{\frac{1}{G}\alpha_3\gamma_p+1}\right)##. So,

[tex]\varepsilon=0.5-0.5\frac{G^2\gamma_Q}{\gamma_p}\int_0^{\infty}x^{-1}\exp\left[\frac{G\left(x+G\gamma_Q\right)}{x\gamma_p}-x\right]E_1\left[(\frac{G\left(x+G\gamma_Q\right)}{x\gamma_p}\right]\,dx[/tex]

The Mathematica code for the above equation is

Code:
yp = 10^(0/10);
GSS = 50;
For[yQdB = -10, yQdB <= 15, yQdB++;
yQ = 10^(yQdB/10);
A1 = 0.5 -
   0.5*((GSS^2)*yQ )/yp*
    NIntegrate[
     1/x*Exp[(GSS*(x + GSS*yQ))/(x*yp) - x]*
      ExpIntegralE[1, (GSS*(x + GSS*yQ))/(x*yp)], {x, 0, Infinity},
     PrecisionGoal -> 5, MaxRecursion -> 20];
Print[A1];]

I did Monte-Carlo simulations as follows

  1. For each ##\gamma_Q## generate three exponential random variables ##\alpha_i## for i=1,2,3.
  2. Find the value of ##\varepsilon(\alpha_1,\alpha_2\alpha_3,\gamma_Q)=\varepsilon(\alpha_1,\alpha_2\alpha_3,\gamma_Q)+0.5\exp\left(-\frac{\frac{\alpha_1}{\alpha_2}G\gamma_Q}{\frac{1}{G}\alpha_3\gamma_p+1}\right)##, where the initial value of ##\varepsilon(\alpha_1,\alpha_2\alpha_3,\gamma_Q)## is zero.
  3. Repeat for ##N=10^5## iterations.
  4. Find the average value as ##\varepsilon(\gamma_Q)=\varepsilon(\alpha_1,\alpha_2\alpha_3,\gamma_Q)/N##
Is there something wrong that gives me the difference between the two curves?

I didn't know how to attach the figure!!
 

Related Threads on Is this integral integrable?

  • Last Post
Replies
3
Views
2K
  • Last Post
Replies
2
Views
811
  • Last Post
Replies
4
Views
1K
  • Last Post
Replies
9
Views
597
  • Last Post
Replies
9
Views
3K
  • Last Post
Replies
1
Views
1K
Replies
2
Views
4K
Replies
1
Views
805
Replies
1
Views
770
  • Last Post
Replies
22
Views
2K
Top