Solving improper integral without approximating

Question

I'm having trouble solving this integral in python. The function being integrated is not defined on the boundaries of integration. I've found a few more questions similar to this, but all were very specific replies to the issue in particular. I don't want to approximate the integral too much, if possible not at all, as the reason I'm doing this integral in the first place is to avoid an approximation. Is there any way to solve this integral?

import numpy as np
from pylab import *
import scipy
from math import *
from scipy import integrate

m_Earth_air = (28.0134*0.78084)+(31.9988*0.209476)+(39.948*0.00934)+(44.00995*0.000314)+(20.183*0.00001818)+(4.0026*0.00000524)+(83.80*0.00000114)+(131.30*0.000000087)+(16.04303*0.000002)+(2.01594*0.0000005)
Tb0 = 288.15
Lb0 = -6.5
Hb0 = 0.0
def Tm_0(z):
    return Tb0+Lb0*(z-Hb0)
k = 1.38*10**-19 #cm^2.kg/s^2.K   #Boltzmann cst
mp = 1.67262177*10**-27 #kg
Rad= 637100000.0 #radius planet #cm
g0 = 980.665 #cm/s^2
def g(z):
    return (g0*((Rad/(Rad+z))**2.0))
def scale_height0(z):
    return k*Tm_0(z*10**-5)/(m_Earth_air*mp*g(z))



def functionz(z,zvar):
    return np.exp(-zvar/scale_height0(z))*((Rad+zvar)/(Rad+z))/((np.sqrt(((Rad+zvar)/(Rad+z))**2.0-1.0)))

def chapman0(z):
    return (1.0/(scale_height0(z)))*((integrate.quad(lambda zvar: functionz(z,zvar), z, np.inf))[0])

print chapman0(1000000)
print chapman0(5000000)

The first block of variables and definitions are fine. The issue is in the "functionz(z,zvar)" and its integration. Any help very much appreciated !

Answer 1

It often helps to eliminate possible numerical instabilities by rescaling variables. In your case zvar starts from 1e6, which is probably causing problems due to some implementation details in quad(). If you scale it as y = zvar / z, so that the integration starts from 1 it seems to converge pretty well for z = 1e6:

def functiony(z, y):
    return np.exp(-y*z/scale_height0(z))*(Rad+y*z)/(Rad+z) / np.sqrt(((Rad+y*z)/(Rad+z))**2.0-1.0)

def chapman0y(z):
    return (1.0/(scale_height0(z)))*((integrate.quad(lambda y: functiony(z,y), 1, np.inf))[0])

>>> print(chapman0y(1000000))

1.6217257661844094e-06

(I set m_Earth_air = 28.8e-3 — this constant is missing in your code, I assumed it is the molar mass of air in (edit) kg/mole).

As for z = 5e6, scale_height0(z) is negative, which gives a huge positive value under the exponent, making the integral divergent on the infinity.

Answer 2

Unless you can solve the integral analytically there is no way to solve it without an approximation over its bounds. This isn't a Python problem, but a calculus problem in general, thus why math classes take such great pains to show you the numeric approximations.

If you don't want it to differ too much, choose a small epsilon with a method that converges fast.

Edit- Clarity on last statement:

Epsilon - ɛ - refers to the step size through the bounds of integration- the delta x- remember that the numeric approximation methods all slice the integral into slivers and add them back up, think of it as the width of each sliver, the smaller the sliver the better the approximation. You can specify these in numerical packages.

A method that converges fast implies the method approaches the true value of the integral quickly and the error of approximation is small for each sliver. For example, the Riemann sum is a naive method which assumes each sliver is a rectangle, while a trapezoid connects the beginning and the end of the sliver with a line to make a trapezoid. Of these 2, trapezoid typically converges faster as it tries to account for the change within the shape. (Neither is typically used as there are better guesses for most functions)

Both of these variables change the computational expense of the calculation. Typically epsilon is the most expensive to change, thus why it is important you choose a good method of approximation (some can differ by an order of magnitude in error for the same epsilon).

All of this will depend on how much error your calculation can tolerate.

Answer 3

I had a similar issue and found that SciPy quad needs you to specify another parameter, epsabs=1e-1000, limit=1000 (stepsize limit), epsrel=1e1 works for everything I've tried. I.e. in this case:

def chapman0(z):    
    return (1.0/(scale_height0(z)))*((integrate.quad(lambda zvar: functionz(z,zvar), z, np.inf, limit=1000, epsabs=1e-1000, epsrel=1e1))[0])[0])
#results:
0.48529410529321887
-1.276589093231806e+21

Seems to be a high absolute error tolerance but for integrals that don't rapidly converge it seems to fix the issue. Just posting for others with similar problems as this post is quite dated. There are algorithms in other packages that converge faster but none that I've found in SciPy. The results are based on the posted code (not the selected answer).