Python scipy.integrate.quad со сложными границами интеграции

Я вычисляю верхнюю неполную гамма-функцию, используя Python mpmath.gammainc:

import numpy as np
from mpmath import gammainc    

z = 0 # define power of t as t^(-1)
a = 0.5+0.4j # integral lower limit
b = np.inf # integral upper limit

myfun = np.array([gammainc(z,a,b,regularized=False)], dtype=complex)

Это одномерный интеграл, определенный в документации mpmath. Я хочу сравнить этот результат myfun с помощью scipy quad функция:

myfun2 = scipy.integrate.quad(exp(-t)/t, a, inf)[0]

Однако я не думаю, что quad принимает сложные аргументы для верхних/нижних границ интегрирования. Я не знаю, возможно ли разделить проблему на реальную/воображаемую части. Есть идеи?


python scipy complex-numbers mpmath
person Medulla Oblongata    schedule 08.08.2018    source источник
comment
Не могли бы вы дать более подробную информацию о математическом определении интегрирования со сложными границами? Я не думаю, что формула exp(-t)/t на самом деле используется для сложного случая, я пытаюсь понять Incomplete_gamma_function. quad не будет работать со сложным аргументом: путь интегрирования должен быть определен (или область) на комплексной плоскости...   -  person xdze2    schedule 08.08.2018

Ответы (1)


arrow_upward
4
arrow_downward

Интеграл нужно брать по горизонтальной полупрямой, идущей от a вправо. (Это не удастся, если a является отрицательным действительным числом, но в любом случае это территория сокращения ветвей). Эта полулиния параметризуется числом a+t, где t является действительным и изменяется от 0 до бесконечности. Итак, интегрируем exp(-(a+t))/(a+t) от 0 до бесконечности.

Кроме того, quad из SciPy требует функций с действительными значениями, поэтому разделите их на действительную и мнимую части. И не забывайте, что функция должна передаваться как вызываемая, например lambda t: np.exp(-t), а не просто exp(-t).

from scipy import integrate
myfun2_re = integrate.quad(lambda t: np.real(np.exp(-(a+t))/(a+t)), 0, np.inf)[0]
myfun2_im = integrate.quad(lambda t: np.imag(np.exp(-(a+t))/(a+t)), 0, np.inf)[0]
myfun2 = myfun2_re + 1j*myfun2_im
print(myfun2)

Это печатает

(0.3411120086192922-0.36240971724285814j)

по сравнению с myfun,

array([ 0.34111201-0.36240972j])

Кстати, нет необходимости оборачивать myfun в массив, если вы просто хотите преобразовать одно число: подойдет complex(gammainc(...)).

person Community    schedule 08.08.2018
comment
Интеграл надо брать по горизонтальной полупрямой, идущей от а вправо. это как-то общее или специфичное для неполной гамма-функции? - person xdze2; 09.08.2018
comment
Специфическая для этой функции, которая определяется как интеграл от a до положительной бесконечности, когда a действительно. Это задает общее направление, по которому должен идти путь интеграции. Существует большая свобода в выборе его формы, потому что такой интеграл не зависит от пути, пока мы не обходим сингулярность (есть один в 0) или сталкиваемся с большими значениями на бесконечности (в правой половине их нет). самолет). Горизонтальная линия — самый простой доступный путь. Итог: для работы со специальными функциями нужно знать какой-то сложный анализ. - person ; 09.08.2018