python - 评估 mpmath gammainc 函数的困难

Question

我正在使用 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 返回相同的输出:

数学

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 函数的实现，所以它不是一个选项。

Answer 1

显然，在某些情况下，mpmath 中的错误导致它在评估 gammainc 时进入无限循环。值得举报mpmath tracker .但至少对于您提到的情况，解决方法是将三参数不完全 Gamma 函数编写为两个 upper 不完全 Gamma 函数 (reference) 的差值。上不完全 Gamma 函数是通过将 两个 参数传递给 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 一致。