无需近似即可求解不正确积分

Question

我在用 python 求解这个积分时遇到了麻烦。被集成的功能没有定义在集成的边界上。 我发现了更多与此类似的问题，但所有问题都是针对该问题的非常具体的答复。 我不想过于近似积分，如果可能的话根本不想，因为我首先进行积分的原因是为了避免近似。 有什么办法可以解这个积分吗？

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)

第一块变量和定义很好。问题在于“functionz(z,zvar)”及其集成。 非常感谢任何帮助！

Answer 1

除非您可以解析地求解积分，否则无法在不超过其边界的近似值的情况下求解它。这不是一个 Python 问题，而是一个一般的微积分问题，这就是为什么数学课会费尽心思向你展示数字近似值。

如果你不想让它相差太大，可以选择一个小的 epsilon，并采用快速收敛的方法。

编辑-最后一句话的清晰度：

Epsilon - ɛ - 指的是通过积分边界的步长 - delta x - 请记住，数值近似方法都将积分切成条子并将它们加起来，将其视为每个条子的宽度，越小条子越接近，近似值就越好。您可以在数字包中指定这些。

快速收敛的方法意味着该方法能够快速接近积分的真实值，并且对于每个条子的近似误差很小。例如，黎曼和是一种简单的方法，它假设每个条子都是矩形，而梯形则用一条线连接条子的开头和结尾以形成梯形。在这两者中，梯形通常收敛得更快，因为它试图解释形状内的变化。 （这两种方法通常都不使用，因为对大多数函数都有更好的猜测）

这两个变量都会改变计算的计算费用。通常，改变 epsilon 的成本最高，因此选择一种好的近似方法很重要（对于相同的 epsilon，有些方法可能会存在一个数量级的误差）。

所有这些都取决于您的计算可以容忍多少误差。

Answer 2

通过重新缩放变量通常有助于消除可能的数值不稳定性。在您的情况下，

zvar

从

1e6

开始，这可能由于

quad()

中的一些实现细节而导致问题。如果将其缩放为

y = zvar / z

，那么积分从

1

开始，对于

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

（我设置了

m_Earth_air = 28.8e-3

——你的代码中缺少这个常数，我假设它是空气的摩尔质量，单位为（edit）kg/mole）。

至于

z = 5e6

，

scale_height0(z)

为负数，在指数下给出了巨大的正值，使得积分在无穷大上发散。

Answer 3

我遇到了类似的问题，发现 SciPy 四边形需要您指定另一个参数，

epsabs=1e-1000

，

limit=1000

（步长限制），

epsrel=1e1

适用于我尝试过的所有内容。 IE。在这种情况下：

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

似乎具有很高的绝对误差容限，但对于不能快速收敛的积分，它似乎可以解决问题。只是为其他有类似问题的人发帖，因为这篇文章已经过时了。其他包中的算法收敛速度更快，但我在 SciPy 中没有找到。结果基于发布的代码（而不是所选的答案）。