Python Integrating to Infinity Numerically

1. Sep 14, 2017

joshmccraney

Hi PF!

I am trying to integrate functions over an infinite domain. One example is $$\int_0^\infty \frac{e^{-x}}{\sqrt{x}}\,dx$$ I know the substitution $u = \sqrt{x}$ reduces this problem to integrating $\exp(-x^2)$, but if I want to integrate the function as is, how would I do this?

I've already tried Gauss-Legendre quadrature and Romberg integration. GL reports NaN and Romberg is evidently unable to handle the infinite limits.

Code (Python):

import numpy as np
import scipy
import scipy.linalg# SciPy Linear Algebra Library
from matplotlib import pyplot as plt# plotting
from scipy import integrate

f = lambda x: np.exp(-x)/np.sqrt(x)# function to integrate
a = 0# lower bound
b = np.inf# upper bound

toler = 10e-3# tolerance

exact = 1.772453850# exact value of integral

# Romberg Integration
I = integrate.romberg(f, a, b, rtol=toler, show=True, divmax=25)

deg = 1# degree of Legendre poly
gauss = 0# initial guess
while abs(exact-gauss) > toler:
x, w = np.polynomial.legendre.leggauss(deg)
# Translate x values from the interval [-1, 1] to [a, b]
t = 0.5*(x + 1)*(b - a) + a
gauss = sum(w * f(t)) * 0.5*(b - a)
deg = deg + 1

print gauss
print deg

2. Sep 14, 2017

Staff: Mentor

NaN or INF are what you get when you go beyond the range of floating point. As you wrote it, you will not get an answer using python. You can always use wolfram alpha, it will give you an answer. There are other approximations and workarounds.

This will help you to bypass some FP issues in the future:
https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html

3. Sep 14, 2017

NFuller

Right, because actually integrating from 0 to infinity would take an infinite amount of time. You should instead integrate to a very large number. The integral should converge very fast since this is a Gaussian integral; you wont have to go that far out to get a very good approximation. If you are conserned with just how accurate the numerical solution is, you should look up the error estimates for the numerical integration methods you are using. I think for Gaussian Quadrature, it goes something like
$$\Delta f(\eta_{n})=\frac{2^{2n+1}(n!)^{4}}{(2n+1)[(2n)!]^{3}}f^{(2n)}(\eta_{n})$$

4. Sep 14, 2017

Dr Transport

5. Sep 14, 2017

NFuller

Looking more closely at your code. It is not clear to me what is happening with the variables named t and gauss. They both have a term $(b-a)$ in them but since $b=\infty$ these values will also be infinity.

6. Sep 15, 2017

George Jones

Staff Emeritus
The book "Numerical Recipes" explains how to handle this type of improper integral. There are two "problems" for numerical integration: 1) the integrand blows up at $x=0$; the region of integration is infinite. Separate the problems, i.e., write
$$I = \int_0^\infty \frac{e^{-x}}{\sqrt{x}}dx = \int_0^1 \frac{e^{-x}}{\sqrt{x}}dx + \int_1^\infty \frac{e^{-x}}{\sqrt{x}}dx = I_1 + I_2.$$

Since $I_1$ blows up like $x^{-1/2}$ as $x$ goes to zero, "Numerical Recipes" says to make the substitution $x=u^2$ in $I_1$. Because the region of integration is infinite in $I_2$, "Numerical Recipes" says to make the substitution $u = e^{-x}$ in $I_2$. These substitutions easily give
$$I = 2 \int_0^1 e^{-u^2}du + \int_0^{1/e} \frac{du}{\sqrt{- \ln u}}.$$

It is easy to integrate numerically each of the integrals. As you say, the first substitution turns the whole question into the integral of a Gaussian, but the idea here is to illustrate techniques that can be used on improper integrals.

7. Sep 21, 2017