Python SciPy integrate.quad用法及代碼示例

本文簡要介紹 python 語言中 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： {函數，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[0] == x 的 xx 數組的長度，其餘項是包含在 quad 的 args 參數中的數字。

此外，為了向後兼容，支持某些 ctypes 調用簽名，但不應在新代碼中使用這些調用簽名。

a： 浮點數

積分的下限(使用 -numpy.inf 表示 -infinity)。

b： 浮點數

積分上限(使用 numpy.inf 表示 +infinity)。

args： 元組，可選

傳遞給 func 的額外參數。

full_output： 整數，可選

非零返回集成信息字典。如果非零，則警告消息也會被抑製，並且消息會附加到輸出元組。

complex_func： 布爾型，可選

指示函數的 (函數) 返回類型為實數 (complex_func=False：默認)或複雜(complex_func=True)。在這兩種情況下，函數的參數都是實數。如果 full_output 也不為零，則信息詞典,信息， 和解釋實數和複數分量在帶有鍵 “real output” 和 “imag output” 的字典中返回。

返回 ：：

y： 浮點數

func 從 a 到 b 的積分。

abserr： 浮點數

結果中絕對誤差的估計。

infodict： dict

包含附加信息的字典。

信息

收斂消息。

解釋

僅附加 ‘cos’ 或 ‘sin’ 權重和無限積分限製，它包含 infodict[‘ierlst’] 中代碼的解釋

其他參數 ：：

epsabs： float 或 int，可選

絕對容錯。默認值為 1.49e-8。quad試圖獲得準確度abs(i-result) <= max(epsabs, epsrel*abs(i))其中i= 積分函數從a到b， 和result是數值近似。看埃普雷爾以下。

epsrel： float 或 int，可選

相對誤差容限。默認值為 1.49e-8。如果epsabs <= 0,埃普雷爾必須大於 5e-29 和50 * (machine epsilon).看易勝寶更多。

limit： float 或 int，可選

自適應算法中使用的子區間數的上限。

points： (浮點數，整數序列)，可選

有界積分區間中的一係列斷點，其中可能出現被積函數的局部困難(例如，奇異點、不連續點)。序列不必排序。請注意，此選項不能與 weight 一起使用。

weight： float 或 int，可選

表示權重函數的字符串。對此和其餘論點的完整解釋可以在下麵找到。

wvar： 可選的

與加權函數一起使用的變量。

wopts： 可選的

用於重用切比雪夫矩的可選輸入。

maxp1： float 或 int，可選

切比雪夫矩數的上限。

limlst： 整數，可選

用於正弦加權和無限 end-point 的周期數上限 (>=3)。

注意：

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

quad() 輸入和輸出的額外信息

如果 full_output 不為零，則第三個輸出參數 (infodict) 是一個字典，其條目如下表所示。對於無限限製，範圍轉換為 (0,1)，並針對此轉換範圍給出可選輸出。令 M 為輸入參數限製，令 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 整數數組，如果為 K<=M/2+2 或 L=M+1-K，則其中的前 L 個元素是指向帶有 L=K 的子間隔上的錯誤估計的指針。讓 I 為序列 infodict['iord'] ，讓 E 為序列 infodict['elist'] 。然後E[I[1]], ..., E[I[L]] 形成一個遞減序列。

如果提供了輸入參數點(即，它不是無)，則以下附加輸出將放置在輸出字典中。假設點序列的長度為 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) 如果 |bb-aa| = |pts[2]-pts[1]| * 2**(-l) 具有級別 l。

‘ndin’：

長度為 P+2 的 rank-1 整數數組。在對區間 (pts[1], pts[2]) 進行第一次積分後，一些區間的誤差估計可能已經人為地增加了，以推進它們的細分。該數組在與發生這種情況的子間隔相對應的插槽中包含一個。

加權被積函數

輸入變量 weight 和 wvar 用於通過選擇的函數列表對被積函數進行加權。不同的積分方法用於計算這些加權函數的積分，這些方法不支持指定斷點。權重的可能值和相應的權重函數是。

weight	使用的權重函數	wvar
‘cos’	cos(w*x)	wvar = w
‘sin’	罪(w*x)	wvar = w
‘alg’	g(x) = ((x-a)**alpha)*((b-x)**beta)	wvar =(阿爾法，貝塔)
‘alg-loga’	g(x)*log(x-a)	wvar =(阿爾法，貝塔)
‘alg-logb’	g(x)*log(b-x)	wvar =(阿爾法，貝塔)
‘alg-log’	g(x)*log(x-a)*log(b-x)	wvar =(阿爾法，貝塔)
‘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’：

包含計算出的 Chebyshev 矩的 rank-2 形狀數組 (25, maxp1)。通過將此數組作為序列 wopts 的第二個元素傳遞並將 infodict[‘momcom’] 作為第一個元素傳遞，這些可以在相同的時間間隔內傳遞給積分。

如果積分極限之一是無窮大，則計算傅裏葉積分(假設 w neq 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 中的例程。本節提供有關調用每個例程的條件的詳細信息以及每個例程的簡短說明。調用的例程取決於重量,積分和積分限製a和b.

四包例程	重量	積分	無限界限
qagse	None	No	No
qagie	None	No	Yes
qagpe	None	Yes	No
qawoe	‘sin’, ‘cos’	No	No
qawfe	‘sin’, ‘cos’	No	A或B
qawse	‘阿爾格*’	No	No
qawce	‘cauchy’	No	No

下麵提供了 [1] 中每個例程的簡短說明。

察格色：

是一個基於全局自適應區間細分與外推法相結合的積分器，它將消除幾種類型的被積函數奇異性的影響。

卡吉：

處理無限間隔內的積分。無限範圍映射到有限區間，隨後應用與 QAGS 中相同的策略。

察格佩：

與 QAGS 的用途相同，但還允許用戶提供有關 trouble-spots 的位置和類型的明確信息，即內部奇點、不連續性和被積函數的其他困難的橫坐標。

卡沃伊：

是用於在有限區間 [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)\) 中的奇點。

誇夫：

計算用戶提供的 \(\omega\) 和 \(f\) 的傅裏葉變換 \(\int^\infty_a \cos(\omega x)f(x)dx\) 或 \(\int^\infty_a \sin(\omega x)f(x)dx\) 。 QAWO 的過程應用於連續的有限間隔，並且通過\(\varepsilon\) 算法的收斂加速應用於一係列積分近似。

誇瑟：

近似 \(\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\)。應用全局自適應細分策略，並對包含以下內容的子區間進行修改後的 Clenshaw-Curtis 積分：a或者b.

誇斯：

計算 \(\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通過分別集成實部和虛部來集成complex-valued函數。

參考：

[1]

羅伯特·皮森斯； de Doncker-Kapenga，伊莉斯；克裏斯托夫·W·尤伯胡貝爾；大衛·卡哈納 (1983)。 QUADPACK：用於自動集成的子程序包。 Springer-Verlag。 ISBN 978-3-540-12553-2。

[2]

麥卡洛，托馬斯；菲利普斯·基思 (1973)。複平麵分析基礎。霍爾特·萊因哈特·溫斯頓。國際標準書號 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)

計算 \(a = 1, 3\) 的 \(\int^1_0 a x \,dx\)

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