对具有无限上限的值数组进行积分

Posted 2023-03-24

【中文标题】对具有无限上限的值数组进行积分

【英文标题】：Integrating an array of values with infinite upper limit

【发布时间】：2020-10-02 21:03:23

【问题描述】：

我正在研究一个化学动力学程序，该程序需要将函数从 0 积分到无穷大。有问题的变量是从数据库中获取的数组。变量的简化版本是 x、F(x) 和 G(x)，其中 F(x) 和 G(x) 是从所述数据库中插值的。

所以我有我的三个数组可以传入被积函数（我不包括用于读取和插值 F 和 G 的文件解析函数）-

def integrand(x,Y,Z):
   return x*Y*Z 
#integrate (x*F*G)dx 0 to infinity
integrated=scipy.integrate.quad(integrand,0,np.inf,args=(F,G)) 

我知道这会返回错误，因为 F 和 G 是数组。对于具有 3 个变量但相对于单个参数（在本例中为 dx）进行积分的被积函数，是否有更好的积分方法？ 注意：在某些情况下，G 收敛速度足够快，可以估计为浮点数。

尝试的解决方案

使用 np.trapz 与 x 数组的任意大上限来代替无穷大

x_n=np.linspace(0,10000,500)
#now, F(x_n) and G(x_n) from np.interp
vars=integrand(x_n,F,G)
answer=np.trapz(vars,x_n,dx=.01)

这种方法可能不适用于其他数据集，也不完全准确。

【参考方案1】：

您可能可以只对向量的每个元素使用此函数。或者，如果对 Y 和 Z 的依赖如此简单，只需将乘数从积分中提取出来，这样您只需积分一次，然后将此结果乘以 Y * Z。

from scipy import integrate
import numpy as np

def integrand(x,Y,Z):
    return (1 / x ** 2) * Y * Z


Y, Z = [1, 2], [3, 4]
integrated = [integrate.quad(integrand,1,np.inf,args=(y,z)) for y, z in zip(Y, Z)]
print(integrated)
# [(3.0, 3.3306690738754696e-14), (8.0, 8.881784197001252e-14)]

【讨论】：

返回列表中的第一个和第二个元组有什么区别？如果我有 len(Y)=len(Z)=100，则该函数返回一个填充了元组的等长列表。积分的结果会是 Integrated[:,0] 的总和吗？

此外，积分的下限 a 在很大程度上影响总和（一个数量级）。我会继续研究您提出的方法。