Что плохого в том, чтобы просто разделить его на реальные и мнимые части?scipy.integrate.quad требует наличия встроенной функции возвратных чисел (иначе вещественных чисел) для алгоритма, который она использует.
import scipy
from scipy.integrate import quad
def complex_quadrature(func, a, b, **kwargs):
def real_func(x):
return scipy.real(func(x))
def imag_func(x):
return scipy.imag(func(x))
real_integral = quad(real_func, a, b, **kwargs)
imag_integral = quad(imag_func, a, b, **kwargs)
return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])
Например,
>>> complex_quadrature(lambda x: (scipy.exp(1j*x)), 0,scipy.pi/2)
((0.99999999999999989+0.99999999999999989j),
(1.1102230246251564e-14,),
(1.1102230246251564e-14,))
, что является ожидаемой ошибкой округления - интегралexp (ix) из 0, pi / 2 - это (1 / i) (e ^ i pi / 2 - e ^ 0) = -i (i - 1) = 1 + i ~ (0.99999999999999989 + 0.99999999999999989j).
И для записи, если не всем понятно на 100%, интеграция - это линейный функционал, означающий, что ∫ {f (x) + kg (x)} dx = ∫ f (x) dx +k ∫ g (x) dx (где k - постоянная по отношению к x).Или для нашего конкретного случая ∫ z (x) dx = ∫ Re z (x) dx + i ∫ Im z (x) dx при z (x) = Re z (x) + i Im z (x).
Если вы пытаетесь выполнить интеграцию по пути в комплексной плоскости (отличной от реальной оси) или области в комплексной плоскости, вам потребуется более сложный алгоритм.
Примечание:Scipy.integrate не будет напрямую обрабатывать сложную интеграцию.Зачем?Он выполняет тяжелую работу в библиотеке FORTRAN QUADPACK , особенно в qagse.f , которая явно требует, чтобы функции / переменные были реальными, прежде чем выполнять свою "глобальную адаптивную квадратуру на основе 21-точечнойКвадратура Гаусса-Кронрода в каждом подинтервале с ускорением по алгоритму Питера Уинна "Эпсилон".Поэтому, если вы не захотите изменить базовый FORTRAN, чтобы он мог обрабатывать сложные числа, скомпилировать его в новую библиотеку, вы не сможете заставить его работать.
Если вы действительно хотите использовать метод Гаусса-Кронрода с комплексными числами ровно в одной интеграции, посмотрите на страницу википедии и реализуйте непосредственно, как сделано ниже (используя 15-пт, 7-пт)правило).Обратите внимание, я запомнил функцию повторения общих вызовов общих переменных (при условии, что вызовы выполняются медленно, как если бы функция была очень сложной).Также применили только правило 7 и 15 пунктов, так как мне не хотелось самому вычислять узлы / веса, а те, которые перечислены в википедии, но получать разумные ошибки для тестовых случаев (~ 1e-14)
import scipy
from scipy import array
def quad_routine(func, a, b, x_list, w_list):
c_1 = (b-a)/2.0
c_2 = (b+a)/2.0
eval_points = map(lambda x: c_1*x+c_2, x_list)
func_evals = map(func, eval_points)
return c_1 * sum(array(func_evals) * array(w_list))
def quad_gauss_7(func, a, b):
x_gauss = [-0.949107912342759, -0.741531185599394, -0.405845151377397, 0, 0.405845151377397, 0.741531185599394, 0.949107912342759]
w_gauss = array([0.129484966168870, 0.279705391489277, 0.381830050505119, 0.417959183673469, 0.381830050505119, 0.279705391489277,0.129484966168870])
return quad_routine(func,a,b,x_gauss, w_gauss)
def quad_kronrod_15(func, a, b):
x_kr = [-0.991455371120813,-0.949107912342759, -0.864864423359769, -0.741531185599394, -0.586087235467691,-0.405845151377397, -0.207784955007898, 0.0, 0.207784955007898,0.405845151377397, 0.586087235467691, 0.741531185599394, 0.864864423359769, 0.949107912342759, 0.991455371120813]
w_kr = [0.022935322010529, 0.063092092629979, 0.104790010322250, 0.140653259715525, 0.169004726639267, 0.190350578064785, 0.204432940075298, 0.209482141084728, 0.204432940075298, 0.190350578064785, 0.169004726639267, 0.140653259715525, 0.104790010322250, 0.063092092629979, 0.022935322010529]
return quad_routine(func,a,b,x_kr, w_kr)
class Memoize(object):
def __init__(self, func):
self.func = func
self.eval_points = {}
def __call__(self, *args):
if args not in self.eval_points:
self.eval_points[args] = self.func(*args)
return self.eval_points[args]
def quad(func,a,b):
''' Output is the 15 point estimate; and the estimated error '''
func = Memoize(func) # Memoize function to skip repeated function calls.
g7 = quad_gauss_7(func,a,b)
k15 = quad_kronrod_15(func,a,b)
# I don't have much faith in this error estimate taken from wikipedia
# without incorporating how it should scale with changing limits
return [k15, (200*scipy.absolute(g7-k15))**1.5]
Контрольный пример:
>>> quad(lambda x: scipy.exp(1j*x), 0,scipy.pi/2.0)
[(0.99999999999999711+0.99999999999999689j), 9.6120083407040365e-19]
Я не доверяю оценке ошибок - я взял что-то из вики для рекомендуемой оценки ошибок при интегрировании от [-1 до 1] и значений donмне не кажется разумным.Например, ошибка выше по сравнению с истиной составляет ~ 5e-15, а не ~ 1e-19.Я уверен, что если кто-то сверился с num рецептами, вы могли бы получить более точную оценку.(Вероятно, придется умножить на (a-b)/2 до некоторой степени или что-то подобное).
Напомним, версия на python менее точна, чем просто дважды вызывать интеграцию на основе QUADPACK от scipy.(Вы можете улучшить его, если хотите).