如何提高积分速度？

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)

。