评估mpmath gammainc函数的困难之处

Question

评估mpmath gammainc函数的困难之处

pythonnumerical-integrationmpmathgamma-function

4

我正在使用Python库mpmath，特别是在评估不完整的Gamma函数。这是根查找例程的一部分，但对于某些复值参数的组合，其评估非常缓慢。

import mpmath as mp
from mpmath import gammainc
def Gamma(a,z0,z1):
    return gammainc(a,a=z0,b=z1,regularized=False)

这里 mpmath.gammainc 函数的评估被卡住了：

>> Gamma(mp.mpc(12.5+17.5j), mp.mpf(0.0), mp.mpf(-12.5))

另一方面，Mathematica 几乎瞬间返回了结果。

In[1]:= Gamma[12.5 + 17.5 I, 0, -12.5]
Out[1]:= 2.38012*10^-7 + 5.54827*10^-7 I

在其他情况下，对于不同的参数，mpmath 和 Mathematica 返回相同的输出结果:

Mathematica

In[2]: Gamma[3.5 I, 0, 10]
Out[2]:= 0.0054741 + 0.000409846 I

Python mpmath

>> Gamma(3.5j,0,10)
mpc(real='0.0054741038497352953', imag='0.00040984640431357779')

您是否有关于这种行为的原因的一些想法？这可以被视为mpmath的问题还是数学上的积分问题？不幸的是，scipy没有提供复杂参数的gamma函数实现，所以这不是一个选择。

- linello

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- user6655984 · Accepted Answer

显然，mpmath中的一个错误在某些情况下会导致其在评估gammainc时进入无限循环。值得在mpmath跟踪器上报告mpmath tracker。但是至少对于您提到的情况，一种解决方法是将三个参数的不完全伽玛函数写为两个上不完全伽玛函数的差reference。上不完全伽玛函数通过向gammainc传递两个参数来计算（即，z1被隐式地视为正无穷大）。

def Gamma(a, z0, z1):
    return gammainc(a, z0) - gammainc(a, z1)

print(Gamma(12.5+17.5j, 0.0, -12.5))

打印出(2.3801203496987e-7 + 5.54827238374855e-7j)，与WolframAlpha一致。