scipy.integrate.

quad#

scipy.integrate.quad(func, a, b, args=(), full_output=0, epsabs=1.49e-08, epsrel=1.49e-08, limit=50, points=None, weight=None, wvar=None, wopts=None, maxp1=50, limlst=50, complex_func=False)[原始碼]#

計算定積分。

使用 Fortran 程式庫 QUADPACK 的技術，對 func 從 a 到 b 進行積分 (可能是無限區間)。

參數:

func{function, scipy.LowLevelCallable}

要積分的 Python 函數或方法。如果 func 接受多個參數，則會沿著對應於第一個參數的軸進行積分。

如果使用者希望提高積分效能，則 f 可以是具有以下簽章之一的 scipy.LowLevelCallable

double func(double x)
double func(double x, void *user_data)
double func(int n, double *xx)
double func(int n, double *xx, void *user_data)

user_data 是 scipy.LowLevelCallable 中包含的資料。在帶有 xx 的呼叫形式中，n 是 xx 陣列的長度，其中包含 xx[0] == x，其餘項目是 quad 的 args 參數中包含的數字。

此外，為了向後相容性，還支援某些 ctypes 呼叫簽章，但不應在新程式碼中使用這些簽章。

afloat

積分下限 (無限小使用 -numpy.inf)。

bfloat

積分上限 (無限大使用 numpy.inf)。

argstuple，選用

要傳遞給 func 的額外參數。

full_outputint，選用

非零值會傳回包含積分資訊的字典。如果是非零值，也會抑制警告訊息，並且訊息會附加到輸出元組。

complex_funcbool，選用

指出函數 (func) 的傳回類型是實數 (complex_func=False：預設值) 還是複數 (complex_func=True)。在這兩種情況下，函數的引數都是實數。如果 full_output 也為非零值，則實數和複數成分的 infodict、message 和 explain 會在字典中傳回，索引鍵分別為 “real output” 和 “imag output”。

傳回值:

yfloat: func 從 a 到 b 的積分。
abserrfloat: 結果中絕對誤差的估計值。
infodictdict: 包含其他資訊的字典。
message: 收斂訊息。
explain: 僅在使用 ‘cos’ 或 ‘sin’ 加權和無限積分限制時附加，其中包含 infodict[‘ierlst’] 中程式碼的說明

其他參數:

epsabsfloat 或 int，選用: 絕對誤差容忍度。預設值為 1.49e-8。quad 嘗試取得 abs(i-result) <= max(epsabs, epsrel*abs(i)) 的準確度，其中 i = func 從 a 到 b 的積分，而 result 是數值近似值。請參閱下方的 epsrel。
epsrelfloat 或 int，選用: 相對誤差容忍度。預設值為 1.49e-8。如果 epsabs <= 0，則 epsrel 必須大於 5e-29 和 50 * (機器 epsilon)。請參閱上方的 epsabs。
limitfloat 或 int，選用: 自適應演算法中使用的子區間數量的上限。
points(float 或 int 序列，選用): 有界積分區間中的斷點序列，被積分函數的局部困難可能發生在這些斷點 (例如，奇異點、不連續點)。序列不必排序。請注意，此選項不能與 weight 一起使用。
weightfloat 或 int，選用: 字串，表示加權函數。下方可以找到此參數和其餘參數的完整說明。
wvar選用: 用於加權函數的變數。
wopts選用: 用於重複使用切比雪夫矩的選用輸入。
maxp1float 或 int，選用: 切比雪夫矩數量的上限。
limlstint，選用: 用於正弦加權和無限端點的週期數上限 (>=3)。

另請參閱

dblquad: 二重積分
tplquad: 三重積分
nquad: n 維積分 (遞迴使用 quad)
fixed_quad: 固定階高斯求積法
simpson: 取樣資料的積分器
romb: 取樣資料的積分器
scipy.special: 用於正交多項式的係數和根

註解

為了獲得有效結果，積分必須收斂；不保證發散積分的行為。

quad() 輸入和輸出的額外資訊

如果 full_output 為非零值，則第三個輸出引數 (infodict) 是一個字典，其中包含如下表所示的條目。對於無限限制，範圍會轉換為 (0,1)，並且相對於此轉換範圍提供選用輸出。設 M 為輸入引數 limit，設 K 為 infodict[‘last’]。這些條目是

‘neval’: 函數評估的次數。
‘last’: 在細分過程中產生的子區間數 K。
‘alist’: 長度為 M 的 rank-1 陣列，其中前 K 個元素是積分範圍分割中子區間的左端點。
‘blist’: 長度為 M 的 rank-1 陣列，其中前 K 個元素是子區間的右端點。
‘rlist’: 長度為 M 的 rank-1 陣列，其中前 K 個元素是子區間上的積分近似值。
‘elist’: 長度為 M 的 rank-1 陣列，其中前 K 個元素是子區間上絕對誤差估計值的模數。
‘iord’: 長度為 M 的 rank-1 整數陣列，其中前 L 個元素是指向子區間上誤差估計值的指標，如果 K<=M/2+2 則 L=K，否則 L=M+1-K。設 I 為序列 infodict['iord']，設 E 為序列 infodict['elist']。則 E[I[1]], ..., E[I[L]] 形成遞減序列。

如果提供了輸入引數 points (即，它不是 None)，則以下額外輸出會放置在輸出字典中。假設 points 序列的長度為 P。

‘pts’: 長度為 P+2 的 rank-1 陣列，其中包含積分限制和區間的斷點 (依遞增順序排列)。這是一個陣列，提供將在其中進行積分的子區間。
‘level’: 長度為 M (=limit) 的 rank-1 整數陣列，其中包含子區間的細分層級，即，如果 (aa,bb) 是 (pts[1], pts[2]) 的子區間，其中 pts[0] 和 pts[2] 是 infodict['pts'] 的相鄰元素，則 (aa,bb) 的層級為 l，如果 |bb-aa| = |pts[2]-pts[1]| * 2**(-l)。
‘ndin’: 長度為 P+2 的 rank-1 整數陣列。在對區間 (pts[1], pts[2]) 進行第一次積分後，某些區間的誤差估計值可能會被人為地增加，以便提前進行細分。此陣列在對應於發生這種情況的子區間的插槽中包含 1。

加權被積分函數

輸入變數 weight 和 wvar 用於透過選定的函數清單對被積分函數進行加權。不同的積分方法用於計算具有這些加權函數的積分，並且這些方法不支援指定斷點。weight 的可能值和對應的加權函數為。

`weight`	使用的加權函數	`wvar`
‘cos’	cos(w*x)	wvar = w
‘sin’	sin(w*x)	wvar = w
‘alg’	g(x) = ((x-a)*alpha)((b-x)**beta)	wvar = (alpha, beta)
‘alg-loga’	g(x)*log(x-a)	wvar = (alpha, beta)
‘alg-logb’	g(x)*log(b-x)	wvar = (alpha, beta)
‘alg-log’	g(x)log(x-a)log(b-x)	wvar = (alpha, beta)
‘cauchy’	1/(x-c)	wvar = c

wvar 保留參數 w、(alpha, beta) 或 c，具體取決於選取的權重。在這些表達式中，a 和 b 是積分限制。

對於 ‘cos’ 和 ‘sin’ 加權，可以使用其他輸入和輸出。

對於有限積分限制，積分是使用 Clenshaw-Curtis 方法執行的，該方法使用切比雪夫矩。對於重複計算，這些矩會儲存在輸出字典中

‘momcom’: 已計算的切比雪夫矩的最大層級，即，如果 M_c 是 infodict['momcom']，則矩已針對長度為 |b-a| * 2**(-l) 的區間計算，l=0,1,...,M_c。
‘nnlog’: 長度為 M(=limit) 的 rank-1 整數陣列，其中包含子區間的細分層級，即，如果對應的子區間是 |b-a|* 2**(-l)，則此陣列的元素等於 l。
‘chebmo’: 形狀為 (25, maxp1) 的 rank-2 陣列，其中包含已計算的切比雪夫矩。這些矩可以傳遞到同一區間上的積分，方法是將此陣列作為序列 wopts 的第二個元素傳遞，並將 infodict[‘momcom’] 作為第一個元素傳遞。

如果其中一個積分限制是無限的，則會計算傅立葉積分 (假設 w 不等於 0)。如果 full_output 為 1 且遇到數值錯誤，除了附加到輸出元組的錯誤訊息之外，還會在輸出元組中附加一個字典，該字典將陣列 info['ierlst'] 中的錯誤代碼翻譯成英文訊息。輸出資訊字典包含以下條目，而不是 ‘last’、‘alist’、‘blist’、‘rlist’ 和 ‘elist’

‘lst’: 積分所需的子區間數 (稱為 K_f)。
‘rslst’: 長度為 M_f=limlst 的 rank-1 陣列，其前 K_f 個元素包含區間 (a+(k-1)c, a+kc) 上的積分貢獻，其中 c = (2*floor(|w|) + 1) * pi / |w| 且 k=1,2,...,K_f。
‘erlst’: 長度為 M_f 的 rank-1 陣列，其中包含與 infodict['rslist'] 中相同位置的區間對應的誤差估計值。
‘ierlst’: 長度為 M_f 的 rank-1 整數陣列，其中包含與 infodict['rslist'] 中相同位置的區間對應的錯誤旗標。請參閱說明字典 (輸出元組中的最後一個條目)，以瞭解程式碼的含義。

QUADPACK 層級常式的詳細資訊

quad 呼叫來自 FORTRAN 程式庫 QUADPACK 的常式。本節提供關於每個常式被呼叫的條件以及每個常式的簡短描述。呼叫的常式取決於 weight、points 和積分限制 a 和 b。

QUADPACK 常式	weight	points	無限邊界
qagse	無	否	否
qagie	無	否	是
qagpe	無	是	否
qawoe	‘sin’、‘cos’	否	否
qawfe	‘sin’、‘cos’	否	a 或 b
qawse	‘alg*’	否	否
qawce	‘cauchy’	否	否

以下提供來自 [1] 的每個常式的簡短描述。

qagse

是一種基於全域自適應區間細分並結合外推法的積分器，它將消除多種類型被積分函數奇異點的影響。

qagie

處理無限區間上的積分。無限範圍會映射到有限區間，然後套用與 QAGS 中相同的策略。

qagpe

與 QAGS 的用途相同，但也允許使用者提供關於問題點 (即，內部奇異點、不連續點和被積分函數其他困難的橫坐標) 的明確資訊。

qawoe

是用於評估有限區間 [a,b] 上 \(\int^b_a \cos(\omega x)f(x)dx\) 或 \(\int^b_a \sin(\omega x)f(x)dx\) 的積分器，其中 \(\omega\) 和 \(f\) 由使用者指定。規則評估元件基於修改後的 Clenshaw-Curtis 技術

自適應細分方案與外推程序結合使用，外推程序是 QAGS 中的修改版本，並允許演算法處理 \(f(x)\) 中的奇異點。

qawfe

計算使用者提供的 \(\omega\) 和 \(f\) 的傅立葉變換 \(\int^\infty_a \cos(\omega x)f(x)dx\) 或 \(\int^\infty_a \sin(\omega x)f(x)dx\)。QAWO 的程序應用於連續的有限區間，並且透過 \(\varepsilon\) 演算法將收斂加速應用於積分近似值序列。

qawse

近似 \(\int^b_a w(x)f(x)dx\)，其中 \(a < b\) 且 \(w(x) = (x-a)^{\alpha}(b-x)^{\beta}v(x)\)，其中 \(\alpha,\beta > -1\)，其中 \(v(x)\) 可能是以下函數之一：\(1\)、\(\log(x-a)\)、\(\log(b-x)\)、\(\log(x-a)\log(b-x)\)。

使用者指定 \(\alpha\)、\(\beta\) 和函數 \(v\) 的類型。全域自適應細分策略適用，並在包含 a 或 b 的那些子區間上使用修改後的 Clenshaw-Curtis 積分。

qawce

計算 \(\int^b_a f(x) / (x-c)dx\)，其中積分必須解釋為柯西主值積分，適用於使用者指定的 \(c\) 和 \(f\)。策略是全域自適應的。修改後的 Clenshaw-Curtis 積分用於包含點 \(x = c\) 的那些區間。

實數變數的複數函數積分

實數變數的複數值函數 \(f\) 可以寫成 \(f = g + ih\)。類似地，\(f\) 的積分可以寫成

\[\int_a^b f(x) dx = \int_a^b g(x) dx + i\int_a^b h(x) dx\]

假設 \(g\) 和 \(h\) 的積分在區間 \([a,b]\) [2] 上存在。因此，quad 透過分別積分實數和虛數成分來積分複數值函數。

參考文獻

[1]

Piessens, Robert; de Doncker-Kapenga, Elise; Überhuber, Christoph W.; Kahaner, David (1983). QUADPACK: A subroutine package for automatic integration. Springer-Verlag. ISBN 978-3-540-12553-2.

[2]

McCullough, Thomas; Phillips, Keith (1973). Foundations of Analysis in the Complex Plane. Holt Rinehart Winston. ISBN 0-03-086370-8

範例

計算 \(\int^4_0 x^2 dx\) 並與解析結果進行比較

>>> from scipy import integrate
>>> import numpy as np
>>> x2 = lambda x: x**2
>>> integrate.quad(x2, 0, 4)
(21.333333333333332, 2.3684757858670003e-13)
>>> print(4**3 / 3.)  # analytical result
21.3333333333

計算 \(\int^\infty_0 e^{-x} dx\)

>>> invexp = lambda x: np.exp(-x)
>>> integrate.quad(invexp, 0, np.inf)
(1.0, 5.842605999138044e-11)

計算 \(\int^1_0 a x \,dx\)，其中 \(a = 1, 3\)

>>> f = lambda x, a: a*x
>>> y, err = integrate.quad(f, 0, 1, args=(1,))
>>> y
0.5
>>> y, err = integrate.quad(f, 0, 1, args=(3,))
>>> y
1.5

使用 ctypes 計算 \(\int^1_0 x^2 + y^2 dx\)，將 y 參數設為 1

testlib.c =>
    double func(int n, double args[n]){
        return args[0]*args[0] + args[1]*args[1];}
compile to library testlib.*

from scipy import integrate
import ctypes
lib = ctypes.CDLL('/home/.../testlib.*') #use absolute path
lib.func.restype = ctypes.c_double
lib.func.argtypes = (ctypes.c_int,ctypes.c_double)
integrate.quad(lib.func,0,1,(1))
#(1.3333333333333333, 1.4802973661668752e-14)
print((1.0**3/3.0 + 1.0) - (0.0**3/3.0 + 0.0)) #Analytic result
# 1.3333333333333333

請注意，與積分區間大小相比，脈衝形狀和其他尖銳特徵可能無法使用此方法正確積分。此限制的一個簡化範例是積分 y 軸反射階梯函數，該函數在積分界限內具有許多零值。

>>> y = lambda x: 1 if x<=0 else 0
>>> integrate.quad(y, -1, 1)
(1.0, 1.1102230246251565e-14)
>>> integrate.quad(y, -1, 100)
(1.0000000002199108, 1.0189464580163188e-08)
>>> integrate.quad(y, -1, 10000)
(0.0, 0.0)