Python中的快速二维数组数值积分答案

【问题标题】：Fast 2d array numerical integration in PythonPython中的快速二维数组数值积分
【发布时间】：2020-11-01 17:03:49
【问题描述】：

（已编辑）

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

其中 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，计算需要大约 100 分钟。 10 秒，我必须针对不同的 n-s 执行一系列计算。并且做很多很多次。所以 10 秒太长了。

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

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

谢谢

【问题讨论】：

一般来说，您应该始终尝试对您的代码进行矢量化处理，即确保它适用于任意 shpae 的输入参数。 quadpy 做到了，这也许是一个很好的参考。（免责声明：我写的。）
jn/jv 评估支配着整个运行时。你的被积函数的并行化 cython 版本看起来是要走的路。如果您提供一个可运行的示例（这还包括一些正确形状和数据类型的随机输入），这将非常有用。

标签： python math numerical-integration

【解决方案1】：

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

【讨论】：