Issue
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 !
Solution
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.
Answered By - fjarri
0 comments:
Post a Comment
Note: Only a member of this blog may post a comment.