Быстрое численное интегрирование двухмерного массива в Python - PullRequest
Новая
       65

Быстрое численное интегрирование двухмерного массива в Python

1 голос
/ 12 июля 2020

(EDITED)

Мне нужно выполнить численное интегрирование:

формула

где Pn - полином Лежандра, λ - длина волны, ρ и z - размер матрицы NxN, J - функция Бесселя нулевого вида.

Я написал такой код: 

import scipy.integrate as integrate
import datetime
from scipy.special import legendre, jn

image_size = 100
image_half_size = 50
scale = 10
cos_beta1 = 0.9999999
cos_beta2 = 0.01745
l_maxOrder = 10
wavelength = 660
k = 2 * np.pi / wavelength

x_center = 8 / scale
y_center = 14 / scale
x, y = np.meshgrid(np.arange(-image_half_size, image_half_size + 1) / scale,
                         np.arange(-image_half_size, image_half_size + 1) / scale,
                         sparse=False,
                         indexing='ij')

ro_p = np.sqrt((x-x_center)**2 + (y-y_center)**2 + 1e-4**2)

z = random.randint(-5, 5)
z_p = np.full((image_size + 1, image_size + 1), z)


def pi_plus_tau(n):
    def frst_deriv(arg):
        return  (n + 1) * (legendre(n+1)(arg) - arg*legendre(n)(arg)) / (arg**2-1)

    def sec_deriv(arg):
        return (n + 1) * (legendre(n)(arg)*((n-2)*(arg**4) + (3-n)*(arg**2)-1) / (arg**2-1) - legendre(n+1)(arg)*(5+2*n)*arg + legendre(n+2)(arg)*(n+2)) / (arg**2-1) ** 2

    def integrand(arg):
        coef = (-1 * jn(0, k * ro_p * (1-arg** 2) ** 0.5) * np.exp(1j * k * z_p * (1-arg)) * (arg**0.5))
        return (frst_deriv(arg) * (1 - arg) + (1 - arg ** 2) * sec_deriv(arg)) * coef

    return integrand


def intergration_plus(n, angle1, angle2):
    return (integrate.quad_vec(pi_plus_tau(int(n)), angle1, angle2, workers=1))[0]


for n in range(1, l_maxOrder):
    print(datetime.datetime.now())
    intergration_plus(n, cos_beta1, cos_beta2)
    print(datetime.datetime.now())

Это хорошо работает, но для N = 100 расчет занимает прибл. 10 секунд, и мне нужно выполнить серию вычислений для разных нс. И сделайте это много-много раз. Итак, 10 секунд - это слишком долго.

Это математическое выражение является частью более крупного выражения. Когда я запускаю код, указанный выше в этом вопросе, он выполняется вдвое быстрее, чем когда я запускаю его внутри всей программы. Некоторые пакеты для использования cython, любые подсказки.

Спасибо

...