Сложности в оценке mpmath gammainc function - PullRequest
Новая
       88

Сложности в оценке mpmath gammainc function

0 голосов
/ 24 августа 2018

Я играю с библиотекой python mpmath, в частности, оценивая неполную гамма-функцию.Это часть процедуры поиска корня, но ее оценка очень медленная для некоторых комбинаций комплексных аргументов. 

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 для сложных аргументов, поэтому это не вариант.

1 Ответ

0 голосов
/ 26 августа 2018

По-видимому, ошибка в mpmath приводит к тому, что при некоторых значениях gammainc он входит в бесконечный цикл. Стоит сообщить о mpmath tracker . Но по крайней мере для случая, который вы упомянули, обходной путь - написать неполную гамма-функцию с тремя аргументами как разность двух верхних неполных гамма-функций ( ссылка ). Верхняя неполная гамма-функция вычисляется путем передачи двух аргументов в 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.

...