在不近似的情况下求解不正确的积分答案

【问题标题】：Solving improper integral without approximating在不近似的情况下求解不正确的积分
【发布时间】：2015-08-26 12:09:24
【问题描述】：

我无法在 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)”及其集成。非常感谢任何帮助！

【问题讨论】：

您希望 SciPy 以分析的方式解决这个问题吗？它不会那样做； scipy.integrate 例程都计算数值近似值。
m_Earth_air 的值是多少？

标签： python scipy integration

【解决方案1】：

除非您可以通过解析求解积分，否则如果不对其边界进行近似，则无法求解。这不是 Python 问题，而是一般的微积分问题，这就是为什么数学课如此费力地向您展示数值近似值的原因。

如果您不希望它相差太大，请选择一个小的 epsilon，并采用快速收敛的方法。

编辑-最后陈述的清晰性：

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

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

这两个变量都会改变计算的计算开销。通常 epsilon 是最昂贵的更改，因此选择一种好的近似方法很重要（对于相同的 epsilon，有些方法的误差可能相差一个数量级）。

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

【讨论】：

当你说“选择一个收敛速度快的小ε”时，你在想什么方法？

【解决方案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)是负数，在指数下给出一个巨大的正值，使得积分在无穷大上发散。

【讨论】：

尽管它现在确实给出了一个值，但问题仍然是当 y=1.0 时 function(z,y) 是无限的，而当 y=inf 时是 nan。因此，积分的值远小于应有的值。附言。你是对的 scale_height0(z)，实际上对于不同的 z 范围，我有不同的比例高度函数，而 scale_height0(z) 仅适用于长达 11 公里（以我的单位为 11e5）。 PS2。 m_Earth_air 的好猜测！它是 28.8 kg/kmol，空气的平均分子量 :)
quad() 在积分收敛时没有无穷大问题。你可以尝试将exp(-x)/sqrt(x)从0整合到inf，你会得到sqrt(pi)。你的积分也没有什么不妥之处。您怎么知道该值小于应有的值？ Wolframalpha 毫无问题地集成它并给出了相同的答案。

【解决方案3】：

我遇到了类似的问题，发现 SciPy quad 需要您指定另一个参数，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 中没有找到。结果基于发布的代码（不是所选答案）。

【讨论】：