solve_ivp#
- scipy.integrate.solve_ivp(fun, t_span, y0, method='RK45', t_eval=None, dense_output=False, events=None, vectorized=False, args=None, **options)[source]#
求解常微分方程式系統的初值問題。
此函數以數值方式積分給定初值的常微分方程式系統
dy / dt = f(t, y) y(t0) = y0
此處 t 是一維獨立變數（時間），y(t) 是 N 維向量值函數（狀態），而 N 維向量值函數 f(t, y) 決定微分方程式。目標是找到近似滿足微分方程式的 y(t)，給定初值 y(t0)=y0。
一些求解器支援在複數域中積分，但請注意，對於剛性 ODE 求解器，右側必須是複數可微分的（滿足柯西-黎曼方程式 [11]）。若要在複數域中求解問題，請傳遞具有複數資料類型的 y0。另一個始終可用的選項是分別重寫實部和虛部的問題。
- 參數:
- funcallable
系統的右側：時間 t 時狀態
y的時間導數。呼叫簽名為fun(t, y)，其中t是純量，而y是具有len(y) = len(y0)的 ndarray。如果使用args，則需要傳遞額外引數（請參閱args引數的文件）。fun必須傳回與y相同形狀的陣列。請參閱 vectorized 以取得更多資訊。- t_span2 成員序列
積分區間 (t0, tf)。求解器從 t=t0 開始積分，直到達到 t=tf。t0 和 tf 都必須是浮點數或可由浮點數轉換函數解釋的值。
- y0array_like，形狀 (n,)
初始狀態。對於複數域中的問題，請傳遞具有複數資料類型的 y0（即使初始值是純實數）。
- method字串或
OdeSolver，選用 要使用的積分方法
‘RK45’ (預設)：5(4) 階顯式龍格-庫塔法 [1]。誤差控制假設四階方法的準確性，但步驟是使用五階精確公式（完成局部外推）。四次內插多項式用於密集輸出 [2]。可以應用於複數域。
‘RK23’：3(2) 階顯式龍格-庫塔法 [3]。誤差控制假設二階方法的準確性，但步驟是使用三階精確公式（完成局部外推）。三次埃爾米特多項式用於密集輸出。可以應用於複數域。
‘DOP853’：8 階顯式龍格-庫塔法 [13]。最初以 Fortran [14] 撰寫的 “DOP853” 演算法的 Python 實作。7 階內插多項式精確到 7 階，用於密集輸出。可以應用於複數域。
‘Radau’：Radau IIA 系列的 5 階隱式龍格-庫塔法 [4]。誤差以三階精確嵌入式公式控制。滿足配置條件的三次多項式用於密集輸出。
‘BDF’：基於導數近似的後向微分公式的隱式多步變階（1 到 5）方法 [5]。[6] 中描述了實作。使用準恆定步長方案，並使用 NDF 修改來增強準確性。可以應用於複數域。
‘LSODA’：具有自動剛性偵測和切換的 Adams/BDF 方法 [7], [8]。這是 ODEPACK 中 Fortran 求解器的包裝器。
顯式龍格-庫塔法（‘RK23’、‘RK45’、‘DOP853’）應用於非剛性問題，而隱式方法（‘Radau’、‘BDF’）應用於剛性問題 [9]。在龍格-庫塔法中，建議使用 ‘DOP853’ 進行高精度求解（rtol 和 atol 的值較低）。
如果不確定，請先嘗試執行 ‘RK45’。如果迭代次數異常多、發散或失敗，則您的問題很可能是剛性的，您應該使用 ‘Radau’ 或 ‘BDF’。‘LSODA’ 也可以是一個不錯的通用選擇，但使用起來可能不太方便，因為它包裝了舊的 Fortran 程式碼。
您也可以傳遞從
OdeSolver衍生的任意類別，該類別實作求解器。- t_evalarray_like 或 None，選用
要儲存計算解的時間點，必須排序且位於 t_span 內。如果為 None（預設），則使用求解器選取的點。
- dense_outputbool，選用
是否計算連續解。預設為 False。
- eventscallable，或 callable 列表，選用
要追蹤的事件。如果為 None（預設），則不會追蹤任何事件。每個事件都發生在時間和狀態的連續函數的零點。每個函數都必須具有簽名
event(t, y)，如果使用args，則必須傳遞額外引數（請參閱args引數的文件）。每個函數都必須傳回一個浮點數。求解器將使用尋根演算法找到 t 的準確值，在此值處event(t, y(t)) = 0。預設情況下，將找到所有零點。求解器在每個步驟中尋找符號變更，因此如果一個步驟中發生多個過零點，則可能會遺漏事件。此外，每個 event 函數可能具有以下屬性- terminal: bool 或 int，選用
當為布林值時，是否在此事件發生時終止積分。當為整數時，在指定次數的此事件發生後終止。如果未指派，則隱含為 False。
- direction: float，選用
過零點的方向。如果 direction 為正數，則 event 僅在從負數變為正數時觸發，反之亦然，如果 direction 為負數。如果為 0，則任一方向都會觸發事件。如果未指派，則隱含為 0。
您可以將屬性（例如
event.terminal = True）指派給 Python 中的任何函數。- vectorizedbool，選用
fun 是否可以向量化方式呼叫。預設為 False。
如果
vectorized為 False，則始終使用形狀為(n,)的y呼叫 fun，其中n = len(y0)。如果
vectorized為 True，則可以使用形狀為(n, k)的y呼叫 fun，其中k是整數。在這種情況下，fun 的行為必須如此，使得fun(t, y)[:, i] == fun(t, y[:, i])（即傳回陣列的每一列都是與y的一列對應的狀態的時間導數）。設定
vectorized=True允許使用方法 ‘Radau’ 和 ‘BDF’ 更快地進行 Jacobian 的有限差分近似，但對於其他方法以及在某些情況下（例如，小的len(y0)）的 ‘Radau’ 和 ‘BDF’，將導致執行速度較慢。- argstuple，選用
要傳遞給使用者定義函數的其他引數。如果給定，則其他引數會傳遞給所有使用者定義函數。因此，舉例來說，如果 fun 具有簽名
fun(t, y, a, b, c)，則 jac（如果給定）和任何事件函數都必須具有相同的簽名，並且 args 必須是長度為 3 的元組。- **options
傳遞給所選求解器的選項。以下列出了已實作求解器的所有可用選項。
- first_stepfloat 或 None，選用
初始步長。預設值為 None，表示應由演算法選擇。
- max_stepfloat，選用
允許的最大步長。預設值為 np.inf，即步長不受限制，僅由求解器決定。
- rtol, atolfloat 或 array_like，選用
相對和絕對容錯。求解器使局部誤差估計值小於
atol + rtol * abs(y)。此處 rtol 控制相對準確度（正確位數），而 atol 控制絕對準確度（正確小數位數）。若要達到所需的 rtol，請將 atol 設定為小於可以從rtol * abs(y)預期的最小值，以便 rtol 主導允許的誤差。如果 atol 大於rtol * abs(y)，則不保證正確位數。相反地，若要達到所需的 atol，請設定 rtol，使得rtol * abs(y)始終小於 atol。如果 y 的成分具有不同的尺度，則可能有利於為不同的成分設定不同的 atol 值，方法是為 atol 傳遞形狀為 (n,) 的 array_like。預設值為 rtol 的 1e-3 和 atol 的 1e-6。- jacarray_like、sparse_matrix、callable 或 None，選用
系統右側關於 y 的 Jacobian 矩陣，‘Radau’、‘BDF’ 和 ‘LSODA’ 方法需要。Jacobian 矩陣的形狀為 (n, n)，其元素 (i, j) 等於
d f_i / d y_j。定義 Jacobian 有三種方法如果為 array_like 或 sparse_matrix，則 Jacobian 假定為常數。‘LSODA’ 不支援。
如果為 callable，則 Jacobian 假定為取決於 t 和 y；將在必要時呼叫為
jac(t, y)。如果使用args，則必須傳遞額外引數（請參閱args引數的文件）。對於 ‘Radau’ 和 ‘BDF’ 方法，傳回值可能是稀疏矩陣。如果為 None（預設），則 Jacobian 將通過有限差分來近似。
通常建議提供 Jacobian，而不是依賴有限差分近似。
- jac_sparsityarray_like、sparse matrix 或 None，選用
定義 Jacobian 矩陣的稀疏結構，用於有限差分近似。其形狀必須為 (n, n)。如果 jac 不是 None，則忽略此引數。如果 Jacobian 在每一列中只有少數非零元素，則提供稀疏結構將大大加快計算速度 [10]。零條目表示 Jacobian 中的對應元素始終為零。如果為 None（預設），則 Jacobian 假定為密集矩陣。‘LSODA’ 不支援，請參閱 lband 和 uband。
- lband, ubandint 或 None，選用
參數定義 ‘LSODA’ 方法的 Jacobian 的頻寬，即
jac[i, j] != 0 only for i - lband <= j <= i + uband。預設值為 None。設定這些參數需要您的 jac 常式以壓縮格式傳回 Jacobian：傳回的陣列必須具有n列和uband + lband + 1行，其中寫入 Jacobian 對角線。具體而言，jac_packed[uband + i - j , j] = jac[i, j]。相同的格式用於scipy.linalg.solve_banded（檢查範例）。這些參數也可以與jac=None一起使用，以減少有限差分估計的 Jacobian 元素數量。- min_stepfloat，選用
‘LSODA’ 方法允許的最小步長。預設情況下，min_step 為零。
- 傳回值:
- 具有以下已定義欄位的 Bunch 物件
- tndarray，形狀 (n_points,)
時間點。
- yndarray，形狀 (n, n_points)
t 處的解的值。
- sol
OdeSolution或 None 找到的解作為
OdeSolution實例；如果 dense_output 設定為 False，則為 None。- t_eventsndarray 列表或 None
針對每種事件類型，包含偵測到該類型事件的陣列列表。如果 events 為 None，則為 None。
- y_eventsndarray 列表或 None
對於每個 t_events 值，解的對應值。如果 events 為 None，則為 None。
- nfevint
右側的評估次數。
- njevint
Jacobian 的評估次數。
- nluint
LU 分解的次數。
- statusint
演算法終止的原因
-1：積分步驟失敗。
0：求解器成功達到 tspan 的末尾。
1：發生終止事件。
- messagestring
終止原因的人工可讀描述。
- successbool
如果求解器達到區間末尾或發生終止事件，則為 True (
status >= 0)。
參考文獻
[1]J. R. Dormand, P. J. Prince, “A family of embedded Runge-Kutta formulae”, Journal of Computational and Applied Mathematics, Vol. 6, No. 1, pp. 19-26, 1980.
[2]L. W. Shampine, “Some Practical Runge-Kutta Formulas”, Mathematics of Computation,, Vol. 46, No. 173, pp. 135-150, 1986.
[3]P. Bogacki, L.F. Shampine, “A 3(2) Pair of Runge-Kutta Formulas”, Appl. Math. Lett. Vol. 2, No. 4. pp. 321-325, 1989.
[4]E. Hairer, G. Wanner, “Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems”, Sec. IV.8.
[5][6]L. F. Shampine, M. W. Reichelt, “THE MATLAB ODE SUITE”, SIAM J. SCI. COMPUTE., Vol. 18, No. 1, pp. 1-22, January 1997.
[7]A. C. Hindmarsh, “ODEPACK, A Systematized Collection of ODE Solvers,” IMACS Transactions on Scientific Computation, Vol 1., pp. 55-64, 1983.
[8]L. Petzold, “Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations”, SIAM Journal on Scientific and Statistical Computing, Vol. 4, No. 1, pp. 136-148, 1983.
[9][10]A. Curtis, M. J. D. Powell, and J. Reid, “On the estimation of sparse Jacobian matrices”, Journal of the Institute of Mathematics and its Applications, 13, pp. 117-120, 1974.
[11][13]E. Hairer, S. P. Norsett G. Wanner, “Solving Ordinary Differential Equations I: Nonstiff Problems”, Sec. II.
範例
基本指數衰減，顯示自動選擇的時間點。
>>> import numpy as np >>> from scipy.integrate import solve_ivp >>> def exponential_decay(t, y): return -0.5 * y >>> sol = solve_ivp(exponential_decay, [0, 10], [2, 4, 8]) >>> print(sol.t) [ 0. 0.11487653 1.26364188 3.06061781 4.81611105 6.57445806 8.33328988 10. ] >>> print(sol.y) [[2. 1.88836035 1.06327177 0.43319312 0.18017253 0.07483045 0.03107158 0.01350781] [4. 3.7767207 2.12654355 0.86638624 0.36034507 0.14966091 0.06214316 0.02701561] [8. 7.5534414 4.25308709 1.73277247 0.72069014 0.29932181 0.12428631 0.05403123]]
指定所需解的點。
>>> sol = solve_ivp(exponential_decay, [0, 10], [2, 4, 8], ... t_eval=[0, 1, 2, 4, 10]) >>> print(sol.t) [ 0 1 2 4 10] >>> print(sol.y) [[2. 1.21305369 0.73534021 0.27066736 0.01350938] [4. 2.42610739 1.47068043 0.54133472 0.02701876] [8. 4.85221478 2.94136085 1.08266944 0.05403753]]
向上發射的砲彈，在撞擊時終止事件。
terminal和direction事件的欄位通過猴子補丁函數應用。此處y[0]是位置，y[1]是速度。拋射物從位置 0 開始，速度為 +10。請注意，積分永遠不會達到 t=100，因為事件是終止事件。>>> def upward_cannon(t, y): return [y[1], -0.5] >>> def hit_ground(t, y): return y[0] >>> hit_ground.terminal = True >>> hit_ground.direction = -1 >>> sol = solve_ivp(upward_cannon, [0, 100], [0, 10], events=hit_ground) >>> print(sol.t_events) [array([40.])] >>> print(sol.t) [0.00000000e+00 9.99900010e-05 1.09989001e-03 1.10988901e-02 1.11088891e-01 1.11098890e+00 1.11099890e+01 4.00000000e+01]
使用 dense_output 和 events 找到砲彈軌跡頂點的位置 100。頂點未定義為終止，因此會找到頂點和著地。t=20 沒有資訊，因此使用 sol 屬性來評估解。sol 屬性通過設定
dense_output=True傳回。或者，可以使用 y_events 屬性來存取事件發生時的解。>>> def apex(t, y): return y[1] >>> sol = solve_ivp(upward_cannon, [0, 100], [0, 10], ... events=(hit_ground, apex), dense_output=True) >>> print(sol.t_events) [array([40.]), array([20.])] >>> print(sol.t) [0.00000000e+00 9.99900010e-05 1.09989001e-03 1.10988901e-02 1.11088891e-01 1.11098890e+00 1.11099890e+01 4.00000000e+01] >>> print(sol.sol(sol.t_events[1][0])) [100. 0.] >>> print(sol.y_events) [array([[-5.68434189e-14, -1.00000000e+01]]), array([[1.00000000e+02, 1.77635684e-15]])]
作為具有其他參數的系統範例，我們將實作 Lotka-Volterra 方程式 [12]。
>>> def lotkavolterra(t, z, a, b, c, d): ... x, y = z ... return [a*x - b*x*y, -c*y + d*x*y] ...
我們使用 args 引數傳入參數值 a=1.5、b=1、c=3 和 d=1。
>>> sol = solve_ivp(lotkavolterra, [0, 15], [10, 5], args=(1.5, 1, 3, 1), ... dense_output=True)
計算密集解並繪製它。
>>> t = np.linspace(0, 15, 300) >>> z = sol.sol(t) >>> import matplotlib.pyplot as plt >>> plt.plot(t, z.T) >>> plt.xlabel('t') >>> plt.legend(['x', 'y'], shadow=True) >>> plt.title('Lotka-Volterra System') >>> plt.show()
使用 solve_ivp 求解具有複數矩陣
A的微分方程式y' = Ay的幾個範例。>>> A = np.array([[-0.25 + 0.14j, 0, 0.33 + 0.44j], ... [0.25 + 0.58j, -0.2 + 0.14j, 0], ... [0, 0.2 + 0.4j, -0.1 + 0.97j]])
求解具有上述
A和作為 3x1 向量的y的 IVP>>> def deriv_vec(t, y): ... return A @ y >>> result = solve_ivp(deriv_vec, [0, 25], ... np.array([10 + 0j, 20 + 0j, 30 + 0j]), ... t_eval=np.linspace(0, 25, 101)) >>> print(result.y[:, 0]) [10.+0.j 20.+0.j 30.+0.j] >>> print(result.y[:, -1]) [18.46291039+45.25653651j 10.01569306+36.23293216j -4.98662741+80.07360388j]
求解具有上述
A和作為 3x3 矩陣的y的 IVP>>> def deriv_mat(t, y): ... return (A @ y.reshape(3, 3)).flatten() >>> y0 = np.array([[2 + 0j, 3 + 0j, 4 + 0j], ... [5 + 0j, 6 + 0j, 7 + 0j], ... [9 + 0j, 34 + 0j, 78 + 0j]])
>>> result = solve_ivp(deriv_mat, [0, 25], y0.flatten(), ... t_eval=np.linspace(0, 25, 101)) >>> print(result.y[:, 0].reshape(3, 3)) [[ 2.+0.j 3.+0.j 4.+0.j] [ 5.+0.j 6.+0.j 7.+0.j] [ 9.+0.j 34.+0.j 78.+0.j]] >>> print(result.y[:, -1].reshape(3, 3)) [[ 5.67451179 +12.07938445j 17.2888073 +31.03278837j 37.83405768 +63.25138759j] [ 3.39949503 +11.82123994j 21.32530996 +44.88668871j 53.17531184+103.80400411j] [ -2.26105874 +22.19277664j -15.1255713 +70.19616341j -38.34616845+153.29039931j]]