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)[原始碼]#
解常微分方程式系統的初始值問題。
此函數給定初始值，用數值方法積分常微分方程式系統。
dy / dt = f(t, y) y(t0) = y0
這裡 t 是一個一維的自變數（時間），y(t) 是一個 N 維的向量值函數（狀態），而一個 N 維的向量值函數 f(t, y) 決定了微分方程式。目標是在給定初始值 y(t0)=y0 的情況下，找到近似滿足微分方程式的 y(t)。
某些求解器支援在複數域中積分，但請注意，對於剛性常微分方程式求解器，右側必須是複可微的（滿足柯西-黎曼方程式 [11]）。要在複數域中解決問題，請傳遞具有複數資料類型的 y0。另一個始終可用的選項是將您的問題的實部和虛部分開重寫。
- 參數：
- fun可呼叫物件
系統的右側：狀態
y在時間t的時間導數。呼叫簽章為fun(t, y)，其中t是一個純量，而y是一個 ndarray，其len(y) = len(y0)。如果使用了args，則需要傳遞額外的參數（請參閱args參數的說明文件）。fun必須返回一個與y形狀相同的陣列。有關更多資訊，請參閱向量化。- t_span二元序列
積分區間 (t0, tf)。求解器從 t=t0 開始積分，直到達到 t=tf。t0 和 tf 都必須是浮點數或可由浮點數轉換函式解讀的值。
- y0類似陣列，形狀 (n,)
初始狀態。對於複數域中的問題，即使初始值是純實數，也要以複數資料類型傳遞 y0。
- method字串或
OdeSolver，可選 要使用的積分方法
‘RK45’（預設）：5(4) 階的顯式 Runge-Kutta 方法 [1]。在假設四階方法準確的情況下控制誤差，但使用五階準確公式進行步驟計算（進行局部外推）。使用四次插值多項式進行密集輸出 [2]。可應用於複數域。
‘RK23’：3(2) 階的顯式 Runge-Kutta 方法 [3]。在假設二階方法準確的情況下控制誤差，但使用三階準確公式進行步驟計算（進行局部外推）。使用三次 Hermite 多項式進行密集輸出。可應用於複數域。
‘DOP853’：8 階的顯式 Runge-Kutta 方法 [13]。最初以 Fortran 編寫的 “DOP853” 演算法的 Python 實作 [14]。使用精度達 7 階的 7 階插值多項式進行密集輸出。可應用於複數域。
‘Radau’：5 階 Radau IIA 家族的隱式 Runge-Kutta 方法 [4]。使用三階準確的嵌入式公式控制誤差。使用滿足配置條件的三次多項式進行密集輸出。
‘BDF’：基於導數近似的反向微分公式的多步變階（1 至 5）隱式方法 [5]。實作遵循 [6] 中描述的方法。使用準恆定步長方案，並使用 NDF 修正來提高精度。可應用於複數域。
「LSODA」：具有自動剛性檢測和切換的 Adams/BDF 方法 [7]、[8]。這是 ODEPACK 中 Fortran 解算器的包裝器。
顯式 Runge-Kutta 方法（「RK23」、「RK45」、「DOP853」）應用於非剛性問題，而隱式方法（「Radau」、「BDF」）應用於剛性問題 [9]。在 Runge-Kutta 方法中，建議使用「DOP853」進行高精度求解（rtol 和 atol 值較低）。
如果不確定，請先嘗試運行「RK45」。如果它進行了異常多次的迭代、發散或失敗，則您的問題很可能是剛性的，您應該使用「Radau」或「BDF」。「LSODA」也可以作為一個通用的好選擇，但由於它包裝了舊的 Fortran 代碼，使用起來可能不太方便。
您也可以傳遞任何繼承自
OdeSolver並實作解算器的類別。- t_evalarray_like 或 None，可選
儲存計算結果的時間點，必須經過排序且位於 t_span 範圍內。如果為 None（預設值），則使用解算器選擇的點。
- dense_output布林值，可選
是否計算連續解。預設值為 False。
- events可呼叫物件，或可呼叫物件的列表，可選
要追蹤的事件。如果為 None（預設值），則不會追蹤任何事件。每個事件發生在時間和狀態的連續函數的零點處。每個函數必須具有簽章
event(t, y)，如果使用了args，則必須傳遞額外的參數（請參閱args參數的說明文件）。每個函數必須返回一個浮點數。解算器將使用尋根演算法找到event(t, y(t)) = 0的精確 t 值。預設情況下，將找到所有零點。解算器會在每個步驟中尋找符號變化，因此如果在一個步驟中發生多個零交叉，則可能會遺漏事件。此外，每個 event 函數可能具有以下屬性：- terminal：布林值或整數，可選
當為布林值時，表示如果發生此事件是否終止積分。當為整數時，在該事件發生指定次數後終止。如果未指定，則隱式為 False。
- direction：浮點數，可選
零交叉的方向。如果 direction 為正，則 event 只會在從負到正時觸發，反之亦然，如果 direction 為負。如果為 0，則任一方向都會觸發事件。如果未指定，則隱式為 0。
您可以在 Python 中將
event.terminal = True等屬性賦值給任何函數。- vectorized布林值，可選
是否可以向量化方式呼叫 fun。預設值為 False。
如果
vectorized為 False，則 fun 將始終以形狀為(n,)的y作為參數被呼叫，其中n = len(y0)。如果
vectorized為 True，則 fun 可能會以形狀為(n, k)的y作為參數被呼叫，其中k為整數。在此情況下，fun 的行為必須滿足fun(t, y)[:, i] == fun(t, y[:, i])（即，返回陣列的每一欄都是與y的一欄對應的狀態的時間導數）。設定
vectorized=True可讓 ‘Radau’ 和 ‘BDF’ 方法更快地進行雅可比矩陣的有限差分近似，但在某些情況下（例如，len(y0)很小），會導致其他方法以及 ‘Radau’ 和 ‘BDF’ 方法執行速度變慢。- argstuple，選用
要傳遞給使用者定義函數的額外參數。如果提供，額外參數會傳遞給所有使用者定義的函數。因此，例如，如果 fun 的簽章為
fun(t, y, a, b, c)，則 jac（如果提供）和任何事件函數都必須具有相同的簽章，且 args 必須是長度為 3 的 tuple。- **options
傳遞給所選求解器的選項。所有已實作求解器可用的選項如下所示。
- first_step浮點數或 None，選用
初始步長。預設值為 None，表示演算法應自行選擇。
- max_step浮點數，選用
允許的最大步長。預設值為 np.inf，即步長不受限制，僅由求解器決定。
- rtol, atol浮點數或類似陣列，選用
相對和絕對容差。求解器會將局部誤差估計值保持在
atol + rtol * abs(y)以下。其中，rtol 控制相對精度（正確數字的位數），而 atol 控制絕對精度（正確小數位數）。要達到所需的 rtol，請將 atol 設定為小於rtol * abs(y)可能產生的最小值，以便 rtol 主導允許的誤差。如果 atol 大於rtol * abs(y)，則無法保證正確數字的位數。反之，要達到所需的 atol，請設定 rtol，使rtol * abs(y)永遠小於 atol。如果 y 的組成元具有不同的尺度，則透過傳遞形狀為 (n,) 的類陣列物件作為 atol，為不同的組成元設定不同的 atol 值可能會有益。rtol 的預設值為 1e-3，atol 的預設值為 1e-6。- jac類陣列物件、稀疏矩陣、可呼叫物件或 None，可選
系統右側相對於 y 的雅可比矩陣，這是「Radau」、「BDF」和「LSODA」方法所需的。雅可比矩陣的形狀為 (n, n)，其元素 (i, j) 等於
d f_i / d y_j。有三種方法可以定義雅可比矩陣：如果是類陣列物件或稀疏矩陣，則假設雅可比矩陣為常數。「LSODA」不支援此種方式。
如果是可呼叫物件，則假設雅可比矩陣取決於 t 和 y；它將根據需要以
jac(t, y)的形式被呼叫。如果使用了args，則必須傳遞額外的參數（請參閱args參數的文件）。對於「Radau」和「BDF」方法，返回值可能是一個稀疏矩陣。如果是 None（預設值），則雅可比矩陣將透過有限差分來近似。
通常建議提供雅可比矩陣，而不是依賴有限差分近似。
- jac_sparsity類陣列物件、稀疏矩陣或 None，可選
為有限差分近似定義雅可比矩陣的稀疏結構。其形狀必須為 (n, n)。如果 jac 不是 None，則會忽略此參數。如果雅可比矩陣在*每一*列中只有少數非零元素，則提供稀疏結構將大大加快計算速度 [10]。零項表示雅可比矩陣中對應的元素始終為零。如果是 None（預設值），則假設雅可比矩陣是稠密的。「LSODA」不支援此參數，請改用 lband 和 uband。
- lband, uband整數或 None，可選
定義「LSODA」方法之雅可比矩陣頻寬的參數，即
jac[i, j] != 0 僅適用於 i - lband <= j <= i + uband。預設值為 None。設定這些參數需要您的 jac 常式以打包格式返回雅可比矩陣：返回的陣列必須具有n個欄和uband + lband + 1個列，其中寫入了雅可比矩陣的對角線。具體來說，jac_packed[uband + i - j , j] = jac[i, j]。scipy.linalg.solve_banded中使用了相同的格式（查看圖示）。這些參數也可以與jac=None一起使用，以減少透過有限差分估計的雅可比矩陣元素的數量。- min_step浮點數，可選
「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。
- nfev整數
右側的評估次數。
- njev整數
雅可比矩陣的評估次數。
- nlu整數
LU 分解的次數。
- status整數
演算法終止的原因
-1：積分步驟失敗。
0：求解器成功到達 tspan 的末尾。
1：發生終止事件。
- message字串
終止原因的人類可讀描述。
- success布林值
如果求解器到達區間末尾或發生終止事件 (
status >= 0)，則為 True。
參考文獻
[1]J. R. Dormand, P. J. Prince，“A family of embedded Runge-Kutta formulae”，Journal of Computational and Applied Mathematics，第 6 卷，第 1 期，第 19-26 頁，1980 年。
[2]L. W. Shampine，“Some Practical Runge-Kutta Formulas”，Mathematics of Computation，第 46 卷，第 173 期，第 135-150 頁，1986 年。
[3]P. Bogacki, L.F. Shampine，“A 3(2) Pair of Runge-Kutta Formulas”，Appl. Math. Lett. 第 2 卷，第 4 期，第 321-325 頁，1989 年。
[4]E. Hairer, G. Wanner，“Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems”，第四章第八節。
[6]L. F. Shampine, M. W. Reichelt，“THE MATLAB ODE SUITE”，SIAM J. SCI. COMPUTE.，第 18 卷，第 1 期，第 1-22 頁，1997 年 1 月。
[7]A. C. Hindmarsh，“ODEPACK, A Systematized Collection of ODE Solvers”，IMACS Transactions on Scientific Computation，第 1 卷，第 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，第 4 卷，第 1 期，第 136-148 頁，1983 年。
[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，第 117-120 頁，1974 年。
[12]維基百科上的洛特卡－沃爾泰拉方程式。
[13]E. Hairer, S. P. Norsett, G. Wanner，“Solving Ordinary Differential Equations I: Nonstiff Problems”，第二章。
範例
基本指數衰減，顯示自動選擇的時間點。
>>> 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 屬性來評估解。透過設定
dense_output=True可以返回 sol 屬性。或者，可以使用 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]])]
作為具有額外參數的系統範例，我們將實作洛特卡－沃爾泰拉方程式 [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 解微分方程式
y' = Ay的例子，其中A為複數矩陣。>>> 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解初始值問題>>> 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解初始值問題>>> 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]]