实现这个二维数值积分的计算更快的方法是什么？

Question

我想做一个二维的数值积分。目前我在用 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(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

一般来说，通过矩阵运算来进行求和要比使用scipy.integrate.quad（或dblquad）快得多。你可以把你的函数f(q, z, t)改写成接收q、z和t的向量，并返回一个三维数组的f值。然后用np.tensordot来计算这些值，再把你的面积元素（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的内存）。所以如果你在进行非常大的高精度积分时，最好先快速检查一下所需的内存。

Answer 2

你可以使用Numba或者低级可调用函数

几乎是你的例子

我直接把函数传给了scipy.integrate.dblquad，而不是像你那样用lambda表达式来生成函数。

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编译函数

要让这个运行，你还需要额外安装Numba和numba-scipy。numba-scipy的目的是提供来自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

使用低级可调用函数

scipy.integrate的函数还可以接受C语言的回调函数，而不是Python函数。这些函数可以用C、Cython或者Numba来编写，就像我在这个例子中使用的那样。主要的好处是，在调用函数时不需要和Python解释器进行交互。

一个很棒的回答来自@Jacques Gaudin，展示了如何轻松做到这一点，包括额外的参数。

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