積分 (scipy.integrate)#
scipy.integrate 子套件提供了數種積分技術，包含常微分方程式積分器。模組的概觀可透過 help 指令取得。
>>> help(integrate)
Methods for Integrating Functions given function object.
quad -- General purpose integration.
dblquad -- General purpose double integration.
tplquad -- General purpose triple integration.
fixed_quad -- Integrate func(x) using Gaussian quadrature of order n.
quadrature -- Integrate with given tolerance using Gaussian quadrature.
romberg -- Integrate func using Romberg integration.
Methods for Integrating Functions given fixed samples.
trapezoid -- Use trapezoidal rule to compute integral.
cumulative_trapezoid -- Use trapezoidal rule to cumulatively compute integral.
simpson -- Use Simpson's rule to compute integral from samples.
romb -- Use Romberg Integration to compute integral from
-- (2**k + 1) evenly-spaced samples.
See the special module's orthogonal polynomials (special) for Gaussian
quadrature roots and weights for other weighting factors and regions.
Interface to numerical integrators of ODE systems.
odeint -- General integration of ordinary differential equations.
ode -- Integrate ODE using VODE and ZVODE routines.
通用積分 (quad)#
函式 quad 用於積分單變數函數在兩點之間的值。這些點可以是 \(\pm\infty\) (以 \(\pm\) inf 表示) 以表示無限的極限。例如，假設您想要沿著區間 \([0, 4.5].\) 積分貝索函數 jv(2.5, x)。
這可以使用 quad 計算。
>>> import scipy.integrate as integrate
>>> import scipy.special as special
>>> result = integrate.quad(lambda x: special.jv(2.5,x), 0, 4.5)
>>> result
(1.1178179380783249, 7.8663172481899801e-09)
>>> from numpy import sqrt, sin, cos, pi
>>> I = sqrt(2/pi)*(18.0/27*sqrt(2)*cos(4.5) - 4.0/27*sqrt(2)*sin(4.5) +
... sqrt(2*pi) * special.fresnel(3/sqrt(pi))[0])
>>> I
1.117817938088701
>>> print(abs(result[0]-I))
1.03761443881e-11
quad 的第一個參數是「可呼叫的」Python 物件（即，函式、方法或類別實例）。請注意在此案例中使用 lambda 函式作為參數。接下來兩個參數是積分的極限。回傳值是一個 tuple，第一個元素包含積分的估計值，第二個元素包含絕對積分誤差的估計值。請注意，在此案例中，此積分的真值為
其中
是菲涅耳正弦積分。請注意，數值計算的積分值與精確結果的誤差在 \(1.04\times10^{-11}\) 之內 — 遠低於報告的誤差估計。
如果要積分的函式需要額外的參數，可以在 args 參數中提供。假設需要計算以下積分
此積分可以使用以下程式碼評估
>>> from scipy.integrate import quad
>>> def integrand(x, a, b):
... return a*x**2 + b
...
>>> a = 2
>>> b = 1
>>> I = quad(integrand, 0, 1, args=(a,b))
>>> I
(1.6666666666666667, 1.8503717077085944e-14)
quad 也允許無限輸入，方法是使用 \(\pm\) inf 作為其中一個參數。例如，假設需要指數積分的數值
（並且忘記了此積分可以計算為 special.expn(n,x)）。函式 special.expn 的功能可以透過根據例程 quad 定義一個新函式 vec_expint 來複製。
>>> from scipy.integrate import quad
>>> import numpy as np
>>> def integrand(t, n, x):
... return np.exp(-x*t) / t**n
...
>>> def expint(n, x):
... return quad(integrand, 1, np.inf, args=(n, x))[0]
...
>>> vec_expint = np.vectorize(expint)
>>> vec_expint(3, np.arange(1.0, 4.0, 0.5))
array([ 0.1097, 0.0567, 0.0301, 0.0163, 0.0089, 0.0049])
>>> import scipy.special as special
>>> special.expn(3, np.arange(1.0,4.0,0.5))
array([ 0.1097, 0.0567, 0.0301, 0.0163, 0.0089, 0.0049])
被積分的函式甚至可以使用 quad 參數（儘管由於被積分函式中使用 quad 可能會產生數值誤差，誤差界限可能會低估誤差）。在此案例中的積分為
>>> result = quad(lambda x: expint(3, x), 0, np.inf)
>>> print(result)
(0.33333333324560266, 2.8548934485373678e-09)
>>> I3 = 1.0/3.0
>>> print(I3)
0.333333333333
>>> print(I3 - result[0])
8.77306560731e-11
最後這個範例顯示，可以使用重複呼叫 quad 來處理多重積分。
警告
數值積分演算法在有限數量的點上對被積分函式取樣。因此，它們無法保證任意被積分函式和積分極限的精確結果（或精確度估計）。例如，考慮高斯積分
>>> def gaussian(x):
... return np.exp(-x**2)
>>> res = integrate.quad(gaussian, -np.inf, np.inf)
>>> res
(1.7724538509055159, 1.4202636756659625e-08)
>>> np.allclose(res[0], np.sqrt(np.pi)) # compare against theoretical result
True
由於被積分函式除了在原點附近幾乎為零之外，我們預期大的但有限的積分極限會產生相同的結果。然而
>>> integrate.quad(gaussian, -10000, 10000)
(1.975190562208035e-203, 0.0)
發生這種情況是因為 quad 中實作的自適應正交積分例程，雖然設計上運作正常，但沒有注意到在如此大的有限區間內函數的微小但重要的部分。為了獲得最佳結果，請考慮使用緊密圍繞被積分函式重要部分的積分極限。
>>> integrate.quad(gaussian, -15, 15)
(1.772453850905516, 8.476526631214648e-11)
具有多個重要區域的被積分函式可以根據需要分成多個部分。
通用多重積分 (dblquad、tplquad、nquad)#
雙重和三重積分的機制已封裝到函式 dblquad 和 tplquad 中。這些函式分別接受要積分的函式以及四個或六個參數。所有內部積分的極限都需要定義為函式。
以下顯示使用雙重積分計算 \(I_{n}\) 數個值的範例
>>> from scipy.integrate import quad, dblquad
>>> def I(n):
... return dblquad(lambda t, x: np.exp(-x*t)/t**n, 0, np.inf, lambda x: 1, lambda x: np.inf)
...
>>> print(I(4))
(0.2500000000043577, 1.29830334693681e-08)
>>> print(I(3))
(0.33333333325010883, 1.3888461883425516e-08)
>>> print(I(2))
(0.4999999999985751, 1.3894083651858995e-08)
作為非常數極限的範例，請考慮積分
此積分可以使用以下表達式評估（請注意，非常數 lambda 函式用於內部積分的上限）
>>> from scipy.integrate import dblquad
>>> area = dblquad(lambda x, y: x*y, 0, 0.5, lambda x: 0, lambda x: 1-2*x)
>>> area
(0.010416666666666668, 1.1564823173178715e-16)
對於 n 重積分，scipy 提供了函式 nquad。積分界限是可迭代物件：可以是常數界限的列表，也可以是非常數積分界限的函式列表。積分的順序（以及因此界限）是從最內層積分到最外層積分。
上面的積分
可以計算為
>>> from scipy import integrate
>>> N = 5
>>> def f(t, x):
... return np.exp(-x*t) / t**N
...
>>> integrate.nquad(f, [[1, np.inf],[0, np.inf]])
(0.20000000000002294, 1.2239614263187945e-08)
請注意，f 的引數順序必須與積分界限的順序一致；也就是說，關於 \(t\) 的內部積分在區間 \([1, \infty]\) 上，而關於 \(x\) 的外部積分在區間 \([0, \infty]\) 上。
非常數積分界限可以用類似的方式處理；上面的範例
可以使用以下方式評估
>>> from scipy import integrate
>>> def f(x, y):
... return x*y
...
>>> def bounds_y():
... return [0, 0.5]
...
>>> def bounds_x(y):
... return [0, 1-2*y]
...
>>> integrate.nquad(f, [bounds_x, bounds_y])
(0.010416666666666668, 4.101620128472366e-16)
這與之前的結果相同。
高斯求積法#
fixed_quad 在固定區間上執行固定階數的高斯求積法。此函式使用 scipy.special 提供的正交多項式集合，它可以計算各種正交多項式的根和求積權重（多項式本身可以作為特殊函數使用，傳回多項式類別的實例 — 例如，special.legendre）。
使用樣本進行積分#
如果樣本是等間隔的，並且可用樣本數對於某些整數 \(k\) 為 \(2^{k}+1\)，則可以使用 Romberg romb 積分，以使用可用樣本獲得積分的高精度估計值。Romberg 積分在步長大小以 2 的冪相關時使用梯形法則，然後對這些估計值執行 Richardson 外插法，以更高精確度逼近積分。
在任意間隔樣本的情況下，可以使用兩個函式 trapezoid 和 simpson。它們分別使用 1 階和 2 階的 Newton-Coates 公式來執行積分。梯形法則將函數近似為相鄰點之間的直線，而辛普森法則將函數在三個相鄰點之間近似為拋物線。
對於奇數個等間隔的樣本，如果函數是 3 階或更低階的多項式，則辛普森法則是精確的。如果樣本不是等間隔的，則只有當函數是 2 階或更低階的多項式時，結果才是精確的。
>>> import numpy as np
>>> def f1(x):
... return x**2
...
>>> def f2(x):
... return x**3
...
>>> x = np.array([1,3,4])
>>> y1 = f1(x)
>>> from scipy import integrate
>>> I1 = integrate.simpson(y1, x=x)
>>> print(I1)
21.0
這完全對應於
然而，積分第二個函數
>>> y2 = f2(x)
>>> I2 = integrate.simpson(y2, x=x)
>>> print(I2)
61.5
不對應於
因為 f2 中多項式的階數大於 2。
使用低階回呼函式進行更快速的積分#
希望縮短積分時間的使用者可以透過 scipy.LowLevelCallable 將 C 函式指標傳遞給 quad、dblquad、tplquad 或 nquad，它將被積分並在 Python 中傳回結果。此處效能的提升來自兩個因素。主要的改進是更快的函式評估，這是透過編譯函式本身提供的。此外，我們還透過移除 quad 中 C 和 Python 之間的函式呼叫來加速。對於諸如正弦波之類的簡單函式，此方法可以提供約 2 倍的速度提升，但對於更複雜的函式，可以產生更顯著的改進（10 倍以上）。因此，此功能適用於需要數值密集積分且願意編寫少量 C 程式碼以大幅縮短計算時間的使用者。
例如，此方法可以透過 ctypes 以幾個簡單的步驟使用
1.) 用 C 語言編寫一個被積分函式，其函式簽名為 double f(int n, double *x, void *user_data)，其中 x 是一個陣列，包含評估函式 f 的點，而 user_data 則是要提供的任意額外資料。
/* testlib.c */
double f(int n, double *x, void *user_data) {
double c = *(double *)user_data;
return c + x[0] - x[1] * x[2]; /* corresponds to c + x - y * z */
}
2.) 現在將此檔案編譯為共享/動態函式庫（快速搜尋將有助於此，因為它與作業系統相關）。使用者必須連結任何使用的數學函式庫等。在 linux 上，它看起來像這樣
$ gcc -shared -fPIC -o testlib.so testlib.c
輸出函式庫將被稱為 testlib.so，但它可能具有不同的檔案副檔名。現在已建立一個函式庫，可以使用 ctypes 載入到 Python 中。
3.) 使用 ctypes 將共享函式庫載入到 Python 中，並設定 restypes 和 argtypes - 這允許 SciPy 正確解譯函式
import os, ctypes
from scipy import integrate, LowLevelCallable
lib = ctypes.CDLL(os.path.abspath('testlib.so'))
lib.f.restype = ctypes.c_double
lib.f.argtypes = (ctypes.c_int, ctypes.POINTER(ctypes.c_double), ctypes.c_void_p)
c = ctypes.c_double(1.0)
user_data = ctypes.cast(ctypes.pointer(c), ctypes.c_void_p)
func = LowLevelCallable(lib.f, user_data)
函式中的最後一個 void *user_data 是選用的，如果不需要，可以省略（在 C 函式和 ctypes argtypes 中都可以）。請注意，座標是以雙精度浮點數陣列而不是單獨的引數傳遞的。
4.) 現在像平常一樣積分函式庫函式，此處使用 nquad
>>> integrate.nquad(func, [[0, 10], [-10, 0], [-1, 1]])
(1200.0, 1.1102230246251565e-11)
Python tuple 會在更短的時間內如預期傳回。所有選用參數都可以與此方法一起使用，包括指定奇異點、無限界限等。
常微分方程式 (solve_ivp)#
積分一組給定初始條件的常微分方程式 (ODE) 是另一個有用的範例。SciPy 中提供了函式 solve_ivp，用於積分一階向量微分方程式
給定初始條件 \(\mathbf{y}\left(0\right)=y_{0}\)，其中 \(\mathbf{y}\) 是長度為 \(N\) 的向量，而 \(\mathbf{f}\) 是從 \(\mathcal{R}^{N}\) 到 \(\mathcal{R}^{N}.\) 的映射。透過將中間導數引入 \(\mathbf{y}\) 向量，高階常微分方程式始終可以簡化為這種型別的微分方程式。
例如，假設需要找到以下二階微分方程式的解
具有初始條件 \(w\left(0\right)=\frac{1}{\sqrt[3]{3^{2}}\Gamma\left(\frac{2}{3}\right)}\) 和 \(\left.\frac{dw}{dz}\right|_{z=0}=-\frac{1}{\sqrt[3]{3}\Gamma\left(\frac{1}{3}\right)}.\) 已知具有這些邊界條件的此微分方程式的解是 Airy 函數
這提供了一種使用 special.airy 檢查積分器的手段。
首先，透過設定 \(\mathbf{y}=\left[\frac{dw}{dz},w\right]\) 和 \(t=z\)，將此 ODE 轉換為標準形式。因此，微分方程式變為
換句話說，
作為一個有趣的提醒，如果 \(\mathbf{A}\left(t\right)\) 在矩陣乘法下與 \(\int_{0}^{t}\mathbf{A}\left(\tau\right)\, d\tau\) 可交換，則此線性微分方程式具有使用矩陣指數的精確解
然而，在這種情況下，\(\mathbf{A}\left(t\right)\) 及其積分不可交換。
此微分方程式可以使用函式 solve_ivp 求解。它需要導數 fprime、時間跨度 [t_start, t_end] 和初始條件向量 y0 作為輸入引數，並傳回一個物件，其 y 欄位是一個陣列，其中包含連續的解值作為列。因此，初始條件在第一個輸出欄中給出。
>>> from scipy.integrate import solve_ivp
>>> from scipy.special import gamma, airy
>>> y1_0 = +1 / 3**(2/3) / gamma(2/3)
>>> y0_0 = -1 / 3**(1/3) / gamma(1/3)
>>> y0 = [y0_0, y1_0]
>>> def func(t, y):
... return [t*y[1],y[0]]
...
>>> t_span = [0, 4]
>>> sol1 = solve_ivp(func, t_span, y0)
>>> print("sol1.t: {}".format(sol1.t))
sol1.t: [0. 0.10097672 1.04643602 1.91060117 2.49872472 3.08684827
3.62692846 4. ]
如所見，solve_ivp 會自動決定其時間步長，除非另有指定。為了將 solve_ivp 的解與 airy 函數進行比較，由 solve_ivp 建立的時間向量會傳遞給 airy 函數。
>>> print("sol1.y[1]: {}".format(sol1.y[1]))
sol1.y[1]: [0.35502805 0.328952 0.12801343 0.04008508 0.01601291 0.00623879
0.00356316 0.00405982]
>>> print("airy(sol.t)[0]: {}".format(airy(sol1.t)[0]))
airy(sol.t)[0]: [0.35502805 0.328952 0.12804768 0.03995804 0.01575943 0.00562799
0.00201689 0.00095156]
solve_ivp 使用其標準參數的解顯示與 airy 函數有很大的偏差。為了盡量減少這種偏差，可以使用相對和絕對容差。
>>> rtol, atol = (1e-8, 1e-8)
>>> sol2 = solve_ivp(func, t_span, y0, rtol=rtol, atol=atol)
>>> print("sol2.y[1][::6]: {}".format(sol2.y[1][0::6]))
sol2.y[1][::6]: [0.35502805 0.19145234 0.06368989 0.0205917 0.00554734 0.00106409]
>>> print("airy(sol2.t)[0][::6]: {}".format(airy(sol2.t)[0][::6]))
airy(sol2.t)[0][::6]: [0.35502805 0.19145234 0.06368989 0.0205917 0.00554733 0.00106406]
為了為 solve_ivp 的解指定使用者定義的時間點，solve_ivp 提供了兩種可能性，也可以互補使用。透過將 t_eval 選項傳遞給函式呼叫，solve_ivp 會在其輸出中傳回 t_eval 這些時間點的解。
>>> import numpy as np
>>> t = np.linspace(0, 4, 100)
>>> sol3 = solve_ivp(func, t_span, y0, t_eval=t)
如果已知函式的 Jacobian 矩陣，則可以將其傳遞給 solve_ivp 以獲得更好的結果。但是請注意，預設積分方法 RK45 不支援 Jacobian 矩陣，因此必須選擇另一種積分方法。支援 Jacobian 矩陣的積分方法之一是以下範例的 Radau 方法。
>>> def gradient(t, y):
... return [[0,t], [1,0]]
>>> sol4 = solve_ivp(func, t_span, y0, method='Radau', jac=gradient)
求解具有帶狀 Jacobian 矩陣的系統#
可以告知 odeint Jacobian 是帶狀的。對於已知是剛性的微分方程式大型系統，這可以顯著提高效能。
作為範例，我們將使用線法 [MOL] 求解一維 Gray-Scott 偏微分方程式。區間 \(x \in [0, L]\) 上函數 \(u(x, t)\) 和 \(v(x, t)\) 的 Gray-Scott 方程式為
其中 \(D_u\) 和 \(D_v\) 分別是組件 \(u\) 和 \(v\) 的擴散係數，而 \(f\) 和 \(k\) 是常數。（有關該系統的更多資訊，請參閱 http://groups.csail.mit.edu/mac/projects/amorphous/GrayScott/）
我們將假設 Neumann（即「無通量」）邊界條件
為了應用線法，我們透過定義 \(N\) 個點 \(\left\{x_0, x_1, \ldots, x_{N-1}\right\}\) 的均勻間隔網格來離散化 \(x\) 變數，其中 \(x_0 = 0\) 且 \(x_{N-1} = L\)。我們定義 \(u_j(t) \equiv u(x_k, t)\) 和 \(v_j(t) \equiv v(x_k, t)\)，並用有限差分取代 \(x\) 導數。也就是說，
然後我們有一個 \(2N\) 個常微分方程式的系統
為了方便起見，已省略 \((t)\) 引數。
為了強制執行邊界條件，我們引入「虛擬」點 \(x_{-1}\) 和 \(x_N\)，並定義 \(u_{-1}(t) \equiv u_1(t)\)、\(u_N(t) \equiv u_{N-2}(t)\)；\(v_{-1}(t)\) 和 \(v_N(t)\) 以類似方式定義。
然後
以及
我們完整的 \(2N\) 個常微分方程式系統是 (1)，適用於 \(k = 1, 2, \ldots, N-2\)，以及 (2) 和 (3)。
我們現在可以開始在程式碼中實作此系統。我們必須將 \(\{u_k\}\) 和 \(\{v_k\}\) 組合成長度為 \(2N\) 的單個向量。兩個明顯的選擇是 \(\{u_0, u_1, \ldots, u_{N-1}, v_0, v_1, \ldots, v_{N-1}\}\) 和 \(\{u_0, v_0, u_1, v_1, \ldots, u_{N-1}, v_{N-1}\}\)。從數學上講，這無關緊要，但選擇會影響 odeint 求解系統的效率。原因在於順序如何影響 Jacobian 矩陣的非零元素的模式。
當變數排序為 \(\{u_0, u_1, \ldots, u_{N-1}, v_0, v_1, \ldots, v_{N-1}\}\) 時，Jacobian 矩陣的非零元素的模式為
變數交錯排序為 \(\{u_0, v_0, u_1, v_1, \ldots, u_{N-1}, v_{N-1}\}\) 時，Jacobian 模式為
在這兩種情況下，都只有五個非平凡對角線，但是當變數交錯排列時，帶寬會小得多。也就是說，主對角線以及主對角線正上方和正下方緊鄰的兩條對角線是非零對角線。這點非常重要，因為 odeint 的輸入 mu 和 ml 是 Jacobian 矩陣的上下帶寬。當變數交錯排列時，mu 和 ml 都是 2。當變數堆疊，\(\{v_k\}\) 接著 \(\{u_k\}\) 時，上下帶寬為 \(N\)。
做出這個決定後，我們就可以編寫實現微分方程組的函數。
首先，我們定義系統的源項和反應項函數
def G(u, v, f, k):
return f * (1 - u) - u*v**2
def H(u, v, f, k):
return -(f + k) * v + u*v**2
接下來，我們定義計算微分方程組右側項的函數
def grayscott1d(y, t, f, k, Du, Dv, dx):
"""
Differential equations for the 1-D Gray-Scott equations.
The ODEs are derived using the method of lines.
"""
# The vectors u and v are interleaved in y. We define
# views of u and v by slicing y.
u = y[::2]
v = y[1::2]
# dydt is the return value of this function.
dydt = np.empty_like(y)
# Just like u and v are views of the interleaved vectors
# in y, dudt and dvdt are views of the interleaved output
# vectors in dydt.
dudt = dydt[::2]
dvdt = dydt[1::2]
# Compute du/dt and dv/dt. The end points and the interior points
# are handled separately.
dudt[0] = G(u[0], v[0], f, k) + Du * (-2.0*u[0] + 2.0*u[1]) / dx**2
dudt[1:-1] = G(u[1:-1], v[1:-1], f, k) + Du * np.diff(u,2) / dx**2
dudt[-1] = G(u[-1], v[-1], f, k) + Du * (- 2.0*u[-1] + 2.0*u[-2]) / dx**2
dvdt[0] = H(u[0], v[0], f, k) + Dv * (-2.0*v[0] + 2.0*v[1]) / dx**2
dvdt[1:-1] = H(u[1:-1], v[1:-1], f, k) + Dv * np.diff(v,2) / dx**2
dvdt[-1] = H(u[-1], v[-1], f, k) + Dv * (-2.0*v[-1] + 2.0*v[-2]) / dx**2
return dydt
我們不會實作計算 Jacobian 矩陣的函數，但我們會告訴 odeint Jacobian 矩陣是帶狀矩陣。這允許底層求解器 (LSODA) 避免計算已知為零的值。對於大型系統，這可以顯著提高效能，如下面的 ipython 會議所示。
首先，我們定義所需的輸入
In [30]: rng = np.random.default_rng()
In [31]: y0 = rng.standard_normal(5000)
In [32]: t = np.linspace(0, 50, 11)
In [33]: f = 0.024
In [34]: k = 0.055
In [35]: Du = 0.01
In [36]: Dv = 0.005
In [37]: dx = 0.025
計時在不利用 Jacobian 矩陣的帶狀結構的情況下的計算
In [38]: %timeit sola = odeint(grayscott1d, y0, t, args=(f, k, Du, Dv, dx))
1 loop, best of 3: 25.2 s per loop
現在設定 ml=2 和 mu=2，讓 odeint 知道 Jacobian 矩陣是帶狀矩陣
In [39]: %timeit solb = odeint(grayscott1d, y0, t, args=(f, k, Du, Dv, dx), ml=2, mu=2)
10 loops, best of 3: 191 ms per loop
速度快了不少！
讓我們確保它們計算出相同的結果
In [41]: np.allclose(sola, solb)
Out[41]: True