时间戳

首页 > 解决方案 > Python中的快速二维数组数值积分

python - Python中的快速二维数组数值积分

问题描述

（已编辑）

我需要执行以下数值积分：

公式

其中 Pn - Legendre 多项式，λ - 波长，ρ 和 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 秒，我必须针对不同的 ns 执行一系列计算。并且做很多很多次。所以10秒太长了。

这个数学表达式是更大表达式的一部分。当我运行这个问题中上面列出的代码时 - 它的运行速度是我在整个程序中运行它时的两倍。

你能告诉我 - 如何在 Python 中与二维数组进行快速数值积分？一些软件包，使用 cython，任何提示。

谢谢你

标签： pythonmathnumerical-integration

解决方案

我能找到解决问题的最好方法是 quadpy 库 - https://pypi.org/project/quadpy/。将您的数据分组到一个 N 维数组中，例如（batch_size、Image_width、Image_height、Legendre 顺序）并使用 np.multiply.outer。在我尝试的所有选项中，这提供了最快的计算。

推荐阅读