计算科学 - 单积分的数值积分 - 吾爱随笔录

单积分的数值积分

计算科学 Python

2021-12-17 02:12:00

我正在尝试使用 Mathematica 和 Python 中的数值积分函数直接评估这个积分。

\int_{0}^{\infty} y^{- (n + 1)} \prod_{i = 1}^{n} γ (a_{i} + 1, b_{i} y) d y

$\int_0^\infty {y^{-(n+1)}\prod_{i=1}^n \gamma(a_i+1,b_i y) }dy$

在哪里 $\gamma(a,x)=\int_0^x dt\ t^{a-1}e^{-t}$ 是下不完全 Gamma 函数。在我的实际问题中 $n$ 可以取 10 到 100 之间的值，其中 $a_i$ 可以取 100 到 1000 之间的值。

这是Python代码

import numpy as np
import scipy.integrate as integrate
import scipy.special as special


def integrand(a,b,y):
    l=len(a)
    x=pow(y,-(l+1))
    for ai,bi in zip(a,b):
        # x*=special.gammainc(ai+1,bi*y) 
        x*=1/bi*special.gammainc(ai+1,bi*y)
    return x

ns=[2,4,5,6]

print integrate.quad(lambda x: integrand([n for n in ns],[0.1,0.8,0.3,0.4],x),0,np.inf)
print integrate.quad(lambda x: integrand([10*n for n in ns],[0.1,0.8,0.3,0.4],x),0,np.inf)
print integrate.quad(lambda x: integrand([100*n for n in ns],[0.1,0.8,0.3,0.4],x),0,np.inf)

相应的输出是

$(8.231566916434924e-05, 2.7894916135933814e-10)$
$(5.884546215247849e-09, 1.1634128468162877e-08)$
$(0.0, 0.0)$

考虑到谦虚，这令人失望 $n$ =4。在第二个输出中 $a_i$ ~ $60$ 并且误差比得到的积分大一个数量级。不需要评论第三个审判案例 $a_i$ ~ $600$ ，零错误！

当然这不是魔术，我使用的是完整的积分区间 $[0,\infty]$ ，并要求数字技巧，如减少积分间隔和任何其他处理以获得合理的结果。

我通过乘以一个常数（ $\prod_r 1/b_i$ ) 并且数字在我的应用程序的所有范围内变得更加稳定。

例如，通过绘制被积函数，我看到了一条跨越范围的单峰曲线 $[5,60]$ , $[50,600]$ 和 $[500,6000]$ 分别消失在外面。如果我分别截断积分区间，结果是完全可以接受的，并且与 Wolfram Mathematica 的结果一致。我尝试了其他更大的集合，并通过试验/错误找到间隔。所以如果我能找到合适的间隔，我会很高兴。

我怎样才能找到任何组的间隔 $a_i$ 和 $b_i$ 和一个使数字稳定的乘法常数？

2个回答

像这样的不正确积分通常无法与双指数求积相匹配。我可以用mpmath计算你的积分。对于较大的 n，似乎需要使用更高的精度。

from mpmath import mp, gammainc, quad, power, inf, nstr

def integrand(a,b,y):
    l=len(a)
    x=power(y,-(l+1))
    for ai,bi in zip(a,b):
        x*=gammainc(ai+1,0,bi*y,regularized=True) / bi
    return x

ns=[2,4,5,6]

mp.pretty = True
mp.dps = 30
print quad(lambda x: integrand([n for n in ns],[0.1,0.8,0.3,0.4],x),[0,inf],error=True)
print quad(lambda x: integrand([10*n for n in ns],[0.1,0.8,0.3,0.4],x),[0,inf],error=True)
mp.dps = 100
print nstr(quad(lambda x: integrand([100*n for n in ns],[0.1,0.8,0.3,0.4],x),[0,inf],error=True),20)

输出是：

(0.0000823156691641082138376342855035, 1.0e-41)
(0.0000000156728839596991659411718340463, 1.0e-28)
(1.6721149546712801741e-12, 1.0e-32)

在某一点拆分间隔也可以解决问题：

>>> mp.dps = 30
>>> print quad(lambda x: integrand([100*n for n in ns],[0.1,0.8,0.3,0.4],x),[0,5000,inf],error=True)
(1.6721149546712801740899041739e-12, 1.0e-32)

如果您不知道要使用哪些参数，则可以使其自动化：检查 quad() 中的误差估计，如果它太大，则将精度加倍或将间隔一分为二并递归。

不正确积分的求积总是有问题的。如果您发现大错误，我会建议以下两件事之一：

正如 choward 建议的那样，使用变量变化技巧将积分更改为具有有限界限，然后像往常一样使用求积。
使用正交方法，该方法使用适合您的问题的基函数。在这种情况下，拉盖尔多项式（或半无限区间的有理切比雪夫）工作（参考博伊德：切比雪夫和傅里叶谱方法。）拉盖尔正交应该更好地工作。甚至还有一个 numpy 模块。

其它你可能感兴趣的问题

上一篇什么决定了通常的化学教科书原子轨道图？下一篇选择不连续 Galerkin 的惩罚