Сципи: Ускорение расчета 2D комплексного интеграла
Я хочу многократно вычислить двумерный комплексный интеграл, используя dblquad из scipy.integrate. Поскольку количество оценок будет довольно большим, я хотел бы увеличить скорость оценки моего кода.
Dblquad, похоже, не в состоянии обрабатывать сложные интеграды. Таким образом, я разделил комплексное интегрирование на действительную и мнимую части:
def integrand_real(x, y):
R1=sqrt(x**2 + (y-y0)**2 + z**2)
R2=sqrt(x**2 + y**2 + zxp**2)
return real(exp(1j*k*(R1-R2)) * (-1j*z/lam/R2/R1**2) * (1+1j/k/R1))
def integrand_imag(x,y):
R1=sqrt(x**2 + (y-y0)**2 + z**2)
R2=sqrt(x**2 + y**2 + zxp**2)
return imag(exp(1j*k*(R1-R2)) * (-1j*z/lam/R2/R1**2) * (1+1j/k/R1))
y0, z, zxp, k и lam - переменные, определенные заранее. Чтобы оценить интеграл по площади круга с радиусом ra, я использую следующие команды:
from __future__ import division
from scipy.integrate import dblquad
from pylab import *
def ymax(x):
return sqrt(ra**2-x**2)
lam = 0.000532
zxp = 5.
z = 4.94
k = 2*pi/lam
ra = 1.0
res_real = dblquad(integrand_real, -ra, ra, lambda x: -ymax(x), lambda x: ymax(x))
res_imag = dblquad(integrand_imag, -ra, ra, lambda x: -ymax(x), lambda x: ymax(x))
res = res_real[0]+ 1j*res_imag[0]
Согласно профилировщику, два интеграта оцениваются около 35000 раз. Общий расчет занимает около одной секунды, что слишком долго для приложения, которое я имею в виду.
Я новичок в научных вычислениях с Python и Scipy и буду рад комментариям, которые указывают на способы повышения скорости оценки. Существуют ли способы переписать команды в функциях integrand_real и integrand_complex, которые могут привести к значительному улучшению скорости?
Имеет ли смысл компилировать эти функции, используя такие инструменты, как Cython? Если да: какой инструмент лучше всего подходит для этого приложения?