pythonnumpyscipynumerical-integrationinfinity

Integration with infinite limit in python


I am facing difficulty in integrating within (0,infinite) limit. I am trying to to do this Jupyter notebook (Don't know if this is a relevant info). The code is as following:

import numpy as np
import matplotlib.pyplot as plt
import scipy as sc
import math
from scipy.integrate import quad

def integrand3(x,z,m,n):
    return ((np.cosh(x))**m)*((np.sinh(x))**n)*(np.exp(-z*np.cosh(x)))

def IgK0(z,m,n):
    IntK1 = quad(integrand3, 0, 1, args=(z,m,n))
    return IntK1

IgK0(1,1,1)

This gives the result:

(0.19224739693489237, 2.1343748650136926e-15)

This is alright. But when I replace the upper limit with infinity I get 'nan' as output. See the following code:

def IgKI(z,m,n):
    IntKI = quad(integrand3, 0, np.inf, args=(z,m,n))
    return IntKI

when I use (0,0) values for (m,n), although there is some error I don't understand, at least I get some answer.

IgKI(1,0,0)

<ipython-input-53-9b9d0956fa85>:2: RuntimeWarning: overflow encountered in cosh
  return ((np.cosh(x))**m)*((np.sinh(x))**n)*(np.exp(-z*np.cosh(x)))
<ipython-input-53-9b9d0956fa85>:2: RuntimeWarning: overflow encountered in sinh
  return ((np.cosh(x))**m)*((np.sinh(x))**n)*(np.exp(-z*np.cosh(x)))
(0.42102443824067737, 1.3628527263018718e-08)

But when I use any other value for (m,n), I get the following:

IgKI(1,0,1)

<ipython-input-53-9b9d0956fa85>:2: RuntimeWarning: overflow encountered in cosh
  return ((np.cosh(x))**m)*((np.sinh(x))**n)*(np.exp(-z*np.cosh(x)))
<ipython-input-53-9b9d0956fa85>:2: RuntimeWarning: overflow encountered in sinh
  return ((np.cosh(x))**m)*((np.sinh(x))**n)*(np.exp(-z*np.cosh(x)))
<ipython-input-53-9b9d0956fa85>:2: RuntimeWarning: invalid value encountered in double_scalars
  return ((np.cosh(x))**m)*((np.sinh(x))**n)*(np.exp(-z*np.cosh(x)))
<ipython-input-76-9bee7a5d6456>:2: IntegrationWarning: The occurrence of roundoff error is detected, which prevents 
  the requested tolerance from being achieved.  The error may be 
  underestimated.
  IntKI = quad(integrand3, 0, np.inf, args=(z,m,n))

(nan, nan)

So, what am I doing wrong?


Solution

  • The errors are clear - you are hitting overflow, and everything breaks. This is due to your functions which quickly diverge:

    >>> [np.cosh(10**x) for x in range(5)]
    __main__:1: RuntimeWarning: overflow encountered in cosh
    [1.5430806348152437, 11013.232920103324, 1.3440585709080678e+43, inf, inf]
    

    At 1000 already Python can't calculate cosh (for example). Indeed, integrating (1,1,1) until 100 is fine. Since this is numerical integration, the value at the bound (and infinity will be translated to a bound as well, but you can just test with 1000) needs to be computed. As you can see from the warning, this means each part of the function is computed separately, not to mention you raise this to a power, and exponentiate an exponent.

    This library will not meet your needs. You can try symbolic integration using sympy, or something more dedicated like Mathematica.