Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Python Integrating to Infinity Numerically

  1. Sep 14, 2017 #1

    joshmccraney

    User Avatar
    Gold Member

    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)

    # 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
     
     
  2. jcsd
  3. Sep 14, 2017 #2

    jim mcnamara

    User Avatar

    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
     
  4. Sep 14, 2017 #3
    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})$$
     
  5. Sep 14, 2017 #4

    Dr Transport

    User Avatar
    Science Advisor
    Gold Member

  6. Sep 14, 2017 #5
    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.
     
  7. Sep 15, 2017 at 3:24 PM #6

    George Jones

    User Avatar
    Staff Emeritus
    Science Advisor
    Gold Member

    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.
     
  8. Sep 21, 2017 at 12:17 PM #7

    joshmccraney

    User Avatar
    Gold Member

    Thanks a lot everyone! I always find your advise very helpful!
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: Integrating to Infinity Numerically
Loading...