实现这种2D数值积分的计算速度更快的方法是什么？

Question

我对2D数值积分感兴趣。现在我正在使用，但它非常慢。请参阅下面的代码。我需要用完全不同的参数来评估这个积分的100次。因此，我想使处理尽可能快速和高效。代码为：scipy.integrate.dblquadimport numpy as npfrom scipy import integratefrom scipy.special import erffrom scipy.special import j0import timeq = np.linspace(0.03, 1.0, 1000)start = time.time()def f(q, z, t):    return t * 0.5 * (erf((t - z) / 3) - 1) * j0(q * t) * (1 / (np.sqrt(2 * np.pi) * 2)) * np.exp(        -0.5 * ((z - 40) / 2) ** 2)y = np.empty([len(q)])for n in range(len(q)):    y[n] = integrate.dblquad(lambda t, z: f(q[n], z, t), 0, 50, lambda z: 10, lambda z: 60)[0]end = time.time()print(end - start)所需时间是212.96751403808594这太过分了。请建议一个更好的方法来实现我想做的事情。在来这里之前，我试图做一些搜索，但没有找到任何解决方案。我已经阅读过可以更好，更快地完成这项工作，但我不知道如何实现相同的工作。请帮忙。quadpy

Answer 1

您可以使用努巴或低级可调用

几乎是你的例子

我只是直接将函数传递给，而不是使用lambdas生成函数的方法。scipy.integrate.dblquad

import numpy as np

from scipy import integrate

from scipy.special import erf

from scipy.special import j0

import time

q = np.linspace(0.03, 1.0, 1000)

start = time.time()

def f(t, z, q):

    return t * 0.5 * (erf((t - z) / 3) - 1) * j0(q * t) * (1 / (np.sqrt(2 * np.pi) * 2)) * np.exp(

        -0.5 * ((z - 40) / 2) ** 2)

def lower_inner(z):

    return 10.

def upper_inner(z):

    return 60.

y = np.empty(len(q))

for n in range(len(q)):

    y[n] = integrate.dblquad(f, 0, 50, lower_inner, upper_inner,args=(q[n],))[0]

end = time.time()

print(end - start)

#143.73969149589539

这已经快了一点点（143与151），但唯一的用途是有一个简单的例子来优化。

使用Numba简单地编译函数

要让它运行，你还需要努姆巴和努姆巴西皮。努巴西比的目的是提供 来自 的包装函数。scipy.special

import numpy as np

from scipy import integrate

from scipy.special import erf

from scipy.special import j0

import time

import numba as nb

q = np.linspace(0.03, 1.0, 1000)

start = time.time()

#error_model="numpy" -> Don't check for division by zero

@nb.njit(error_model="numpy",fastmath=True)

def f(t, z, q):

    return t * 0.5 * (erf((t - z) / 3) - 1) * j0(q * t) * (1 / (np.sqrt(2 * np.pi) * 2)) * np.exp(

        -0.5 * ((z - 40) / 2) ** 2)

def lower_inner(z):

    return 10.

def upper_inner(z):

    return 60.

y = np.empty(len(q))

for n in range(len(q)):

    y[n] = integrate.dblquad(f, 0, 50, lower_inner, upper_inner,args=(q[n],))[0]

end = time.time()

print(end - start)

#8.636585235595703

使用低级可调用

这些函数还提供了传递 C 回调函数而不是 Python 函数的可能性。例如，这些函数可以编写在C，Cython或Numba中，我在此示例中使用。主要优点是，在函数调用时不需要Python解释器交互。scipy.integrate

@Jacques高丁的出色答案显示了一种简单的方法来做到这一点，包括其他参数。

import numpy as np

from scipy import integrate

from scipy.special import erf

from scipy.special import j0

import time

import numba as nb

from numba import cfunc

from numba.types import intc, CPointer, float64

from scipy import LowLevelCallable

q = np.linspace(0.03, 1.0, 1000)

start = time.time()

def jit_integrand_function(integrand_function):

    jitted_function = nb.njit(integrand_function, nopython=True)

    #error_model="numpy" -> Don't check for division by zero

    @cfunc(float64(intc, CPointer(float64)),error_model="numpy",fastmath=True)

    def wrapped(n, xx):

        ar = nb.carray(xx, n)

        return jitted_function(ar[0], ar[1], ar[2])

    return LowLevelCallable(wrapped.ctypes)

@jit_integrand_function

def f(t, z, q):

    return t * 0.5 * (erf((t - z) / 3) - 1) * j0(q * t) * (1 / (np.sqrt(2 * np.pi) * 2)) * np.exp(

        -0.5 * ((z - 40) / 2) ** 2)

def lower_inner(z):

    return 10.

def upper_inner(z):

    return 60.

y = np.empty(len(q))

for n in range(len(q)):

    y[n] = integrate.dblquad(f, 0, 50, lower_inner, upper_inner,args=(q[n],))[0]

end = time.time()

print(end - start)

#3.2645838260650635

Answer 2

一般来说，通过矩阵运算求和比使用scipy.integrate.quad（或dblquad）要快得多。您可以重写 f（q， z， t） 以取 q、z 和 t 向量，并使用 np.tensordot 返回 f 值的 3D 数组，然后将面积元素 （dtdz） 与函数值相乘，并使用 np.sum 对它们求和。如果您的面积元素不是常量，则必须创建一个面积元素数组并使用 np.einsum 要考虑您的集成限制，您可以在汇总之前使用屏蔽的数组来屏蔽集成限制之外的函数值。请注意，np.einsum 忽略了掩码，因此，如果您使用 einsum，则可以使用 np.where 将积分限制之外的函数值设置为零。示例（使用常量面积元素和简单的积分限制）：

import numpy as np

import scipy.special as ss

import time

def f(q, t, z):

    # Making 3D arrays before computation for readability. You can save some time by

    # Using tensordot directly when computing the output

    Mq = np.tensordot(q, np.ones((len(t), len(z))), axes=0)

    Mt = np.tensordot(np.ones(len(q)), np.tensordot(t, np.ones(len(z)), axes = 0), axes = 0)

    Mz = np.tensordot(np.ones((len(q), len(t))), z, axes = 0)

    return Mt * 0.5 * (ss.erf((Mt - Mz) / 3) - 1) * (Mq * Mt) * (1 / (np.sqrt(2 * np.pi) * 2)) * np.exp(

     -0.5 * ((Mz - 40) / 2) ** 2)

q = np.linspace(0.03, 1, 1000)

t = np.linspace(0, 50, 250)

z = np.linspace(10, 60, 250)

#if you have constand dA you can shave some time by computing dA without using np.diff

#if dA is variable, you have to make an array of dA values and np.einsum instead of np.sum

t0 = time.process_time()

dA = np.diff(t)[0] * np.diff(z)[0]

func_vals = f(q, t, z)

I = np.sum(func_vals * dA, axis=(1, 2))

t1 = time.process_time()

这在我的2012年macbook pro（2.5GHz i5）上花费了18.5秒，dA = 0.04。以这种方式做事还允许您在精度和效率之间轻松进行选择，并将 dA 设置为在了解函数行为方式时有意义的值。

但是，值得注意的是，如果您想要更大的点数，则必须拆分积分，否则您将冒着最大化内存（1000 x 1000 x 1000）的风险，需要8GB的RAM。因此，如果您正在执行具有高优先级的非常大的集成，那么在运行之前对所需的内存进行快速检查是值得的。