Integrating to Infinity Numerically

[member 428835]

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.

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)

# Gauss-Legendre Quadrature Integration
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

 

[jim mcnamara]

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

 

[NFuller]

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

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 won't 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})$$

 

[Dr Transport]

use Gauss–Laguerre quadrature

 

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

 

[George Jones]

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?

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.

 

[member 428835]

Thanks a lot everyone! I always find your advise very helpful!