# About numerical Integration for a non-solvable improper integral

I have an improper integral which does not have a closed-form expression.

$$p1 = 2\pi \lambda_M \int_{r=0}^{\infty} \frac{e^{- \mathcal{Z} r^2}}{1 + \mathcal{Y}r^\alpha} r \, dr$$

If I used numerical-integration using simpson's method, how can I find the ideal value of the upper limit of $r$. I am having different results for different values of $r$

• How about a substitution to make it a definite integral? For instance $r=\tan\frac{\pi x}{2}$? – Klaas van Aarsen Feb 11 '18 at 4:13
• Sorry, it is not indefinite integral, but improper integral instead. My bad ! – Kashan Feb 11 '18 at 4:18
• But $\tan \pi/2$ is undefined? – Kashan Feb 11 '18 at 4:28
• I think this can actually be done in terms of a sum by expanding through the geometric series and then changing variables to relate it to the Gamma function or maybe the incomplete Gamma functions. But anyway one can turn $(0,\infty)$ into a bounded interval using a substitution like $u=1/r$. – Ian Feb 11 '18 at 4:42
• @I like Serena Sorry, I was sloppy. We can map $(1,\infty)$ to $(0,1)$ using $1/x$. So $x=1/(1+r)$ works fine here. – Ian Feb 11 '18 at 5:01

## 2 Answers

One way is to not use an upper truncation limit at all and instead map to an integral over $(0,1)$ using a substitution such as $x=1/(1+r)$, or $r=\tan(\pi x/2)$.

Another way is to note that the integrand is bounded by $e^{-\mathcal{Z} r^2}$ (assuming $\mathcal{Y} \geq 0$ so that everything makes sense). This combined with some familiar estimates for tails of the normal distribution allows you to show that the tail of your integral (without the coefficient in front) is bounded by $e^{-\mathcal{Z} R^2}$.

Still another way would be to use Gauss-Hermite quadrature (which would converge quite rapidly if $\mathcal{Y}$ isn't too large).

• Gauss-Hermite might not work so well because of the $r^{\alpha}$ factor in the integrand. – Ron Gordon Feb 11 '18 at 18:57
• @RonGordon I doubt it will matter, because unless $\mathcal{Y}$ is rather large compared to $\mathcal{Z}^{-1}$, $\frac{1}{1+\mathcal{Y} r^\alpha}$ will be pretty flat where $e^{-\mathcal{Z} r^2}$ is appreciable. Of course you could do direct analysis of this by adding-and-subtracting $\int_0^\infty e^{-\mathcal{Z} r^2} dr$. – Ian Feb 11 '18 at 19:13
• If I use $x = 1/(1+r)$, the upper limit of $r \= \infty$ changes to $x\ = 0$ and $r = 0$ changes to $x = 1$, i.e. upper and lower limits are changed, do I have to change the sign of the expression...pardon me if that is a silly ask :( – Kashan Feb 15 '18 at 11:48
• @Sjaffry You probably should, but if you actually integrate from right to left then your $\Delta x$ will be negative, which will get the sign to match. – Ian Feb 15 '18 at 12:48

Sub $r=\sqrt{x/Z}$ and the integral becomes

$$\frac{\pi \lambda_M}{\sqrt{Z}} \int_0^{\infty} dx \, \frac{e^{-x}}{1+Y Z^{-\alpha/2} x^{\alpha/2}}$$

Now you can use something called Gauss-Laguerre quadrature. This is where one uses the orthogonality properties of Laguerre polynomials $L_n$ to form the approximation

$$\int_0^{\infty} dx \, f(x) e^{-x} \approx \sum_{i=1}^n w_i f(x_i)$$

where

$$L_n(x_i) = 0$$

and

$$w_i = \frac{x_i}{(n+1)^2 L_{n+1}(x_i)^2}$$

Note that the Laguerre polynomial is defined as

$$L_n(x) = e^x \frac{d^n}{dx^n} (e^{-x} x^n)$$

A table of zeros may be found here.

• What shoud be the value of $n$ in this case? Is it any value of choice and the greater the $n$, the lesser the error between approximated and actual result? – Kashan Feb 11 '18 at 22:07
• @Sjaffry: yes, the higher the $n$ the smaller the error in theory. See the table of zeroes at the link - for large $n$, the weights get smaller for larger $k$. There is a trade off, not just in terms of computation but in terms of potential roundoff error. – Ron Gordon Feb 12 '18 at 0:56