海绵宝宝撒
您可以使用努巴或低级可调用几乎是你的例子我只是直接将函数传递给，而不是使用lambdas生成函数的方法。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(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.specialimport numpy as npfrom scipy import integratefrom scipy.special import erffrom scipy.special import j0import timeimport numba as nbq = 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 npfrom scipy import integratefrom scipy.special import erffrom scipy.special import j0import timeimport numba as nbfrom numba import cfuncfrom numba.types import intc, CPointer, float64from scipy import LowLevelCallableq = 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_functiondef 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
绝地无双
一般来说，通过矩阵运算求和比使用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 npimport scipy.special as ssimport timedef 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.sumt0 = 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。因此，如果您正在执行具有高优先级的非常大的集成，那么在运行之前对所需的内存进行快速检查是值得的。