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

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- 请记住，数值近似方法都将积分切成条状并将它们加起来，将其视为每个的宽度条子，条子越小，近似值越好。您可以在数字包中指定这些。

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

这两个变量都会改变计算的计算开销。通常 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 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

似乎是一个高绝对误差容忍度，但对于不能快速收敛的积分，它似乎解决了这个问题。只是 posting 给其他有类似问题的人，因为这个 post 已经过时了。其他包中的算法收敛速度更快，但我在 SciPy 中发现了 none。结果基于 posted 代码（不是选择的答案）。

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

Solving improper integral without approximating

python

scipy

integration