Is my integral integrable using Mathematica or is there a fundamental error?

Click For Summary
SUMMARY

The integral under discussion is defined as \int_0^{\infty}\frac{1}{x^2}\exp\left(\frac{A}{x}-x\right)E_1\left(B+\frac{A}{x}\right)\,dx, where E_1(.) is the exponential integral function and A and B are positive constants. Attempts to evaluate this integral using Mathematica's NIntegrate resulted in overflow errors due to the \frac{1}{x} term causing numerical instability. The discussion concluded that the integral diverges at zero, making it non-integrable without modification. A workaround involves using piecewise functions to handle the singularity at zero.

PREREQUISITES
  • Understanding of integral calculus, specifically improper integrals.
  • Familiarity with the exponential integral function E_1(z).
  • Experience with Mathematica, particularly the NIntegrate function.
  • Knowledge of conditional expressions in computer algebra systems (CAS).
NEXT STEPS
  • Learn how to implement piecewise functions in Mathematica to handle singularities.
  • Explore the properties of the exponential integral function E_1(z) in detail.
  • Study numerical integration techniques for handling improper integrals.
  • Investigate alternative methods for evaluating divergent integrals, such as regularization techniques.
USEFUL FOR

Mathematicians, physicists, and engineers who are working with complex integrals, particularly those involving singularities and numerical methods for integration.

  • #31
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##.
 
Physics news on Phys.org
  • #32
Oh, really? So, I think you can see now where the integral in my first post came from. I will follow on this.
 
  • #33
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

\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.

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

\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

Is the above integral inetgrable now?
 
  • #35
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,

\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

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!
 
  • #36
Is it an error margin in the numerical evaluation of the integral?
 
  • #37
Can you show the figures? What error do you have?
 
  • #38
How to upload a figure from my PC?
 
  • #39
Upload image to a hosting site (dropbox, google photos, tumblr, flickr, etc) then click to "image" on the toolbar and insert the image url.
 
  • #40
It doesn't work. I tried both dropbox and google photos. Isn't it "get a link" that I need to insert here?
 
  • #41
Attached is the figure. Sim=Monte-Carlo simulations
 

Attachments

  • untitled.jpg
    untitled.jpg
    24.6 KB · Views: 473
  • #42
I think Monte-Carlo simulations are accurate, so, I suspect, given that everything else done properly, the accuracy of NIntegrate in Mathematica is the issue. Is there any possible reason and I cannot see it?
 
  • #43
Monte-Carlo convergence is quite slow, ##C/\sqrt N##, maybe the error is due to that. Also there could be some details that Monte-Carlo misses.
On the other hand how Mathematica does NIntegrate is hidden, it might also introduce some systematic error here.
 
  • #44
I remember using Mathematica to evaluate some numerical integral in my master thesis. There was no difference between the numerical integration and Monte-Carlo simulations then!
 

Similar threads

  • · Replies 21 ·
Replies
21
Views
4K
  • · Replies 3 ·
Replies
3
Views
2K
  • · Replies 8 ·
Replies
8
Views
3K
  • · Replies 6 ·
Replies
6
Views
3K
  • · Replies 1 ·
Replies
1
Views
2K
  • · Replies 9 ·
Replies
9
Views
2K
  • · Replies 3 ·
Replies
3
Views
3K
  • · Replies 6 ·
Replies
6
Views
5K
  • · Replies 6 ·
Replies
6
Views
2K
  • · Replies 5 ·
Replies
5
Views
2K