我對2D數值積分感興趣?，F在我正在使用，但它非常慢。請參閱下面的代碼。我需要用完全不同的參數來評估這個積分的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
2 回答
海綿寶寶撒
TA貢獻1809條經驗 獲得超8個贊
您可以使用努巴或低級可調用
幾乎是你的例子
我只是直接將函數傳遞給，而不是使用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
絕地無雙
TA貢獻1946條經驗 獲得超4個贊
一般來說，通過矩陣運算求和比使用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。因此，如果您正在執行具有高優先級的非常大的集成，那么在運行之前對所需的內存進行快速檢查是值得的。
添加回答
舉報
0/150
提交
取消