python - 如何提高积分速度？

Question

我有一个我想整合的形式的功能:

def f(z, t, q):
    return t * 0.5 * (erf((t - z) / 3) - 1) * j0(q * t) * np.exp(-0.5 * ((z - 40) / 2) ** 2)

我以这个函数为例来理解和展示不同积分方法之间的区别。
基于此平台上的各种答案以及文档中的早期知识，我尝试了不同的方法来提高函数的集成速度。我将在这里简要解释和比较它们:
1. 仅使用 python 和 scipy.quad: (需要 202.76 秒)

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)  * 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)

#202.76 s

正如预期的那样，这不是很快。
2. 使用低级可调用函数并使用 Numba 编译函数: (需要 6.46 秒)

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)

#6.46 s

这比之前的方法更快。
3. 使用 scipy.quad_vec : (需要 2.87 秒)

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()

def f(z, t, q):
    return t * 0.5 * (erf((t - z) / 3) - 1) * j0(q * t) * np.exp(-0.5 * ((z - 40) / 2) ** 2)

ans2 = integrate.quad_vec(lambda t: integrate.quad_vec(lambda z: f(z,t,q),10, 60)[0],0,50)[0]

end = time.time()
print(end - start)
#2.87s

方法涉及 scipy.quad_vec是最快的。我想知道如何让它更快？ scipy.quad_vec不幸的是不接受低级可调用函数；无论如何要矢量化 scipy.quad或 scipy.dblquad功能与 scipy.quad_vec 一样有效?或者如果我可以使用任何其他方法？任何形式的帮助将不胜感激。

Answer 1

这是一个快 10 倍的版本

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 fa(z, t, q):
    return (erf((t - z) / 3) - 1) * np.exp(-0.5 * ((z - 40) / 2) ** 2)

ans3 = 0.5 * integrate.quad_vec(lambda t: t * j0(q * t) * 
         integrate.quad_vec(lambda z: fa(z,t,q),10, 60)[0],0,50)[0]

end = time.time()
print(end - start)

这里的不同之处在于采取 j0(q*t)出内积分。