scipy.integrate.

odeint

scipy.integrate.odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None, rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0, tfirst=False)

積分常微分方程式系統。

注意

對於新程式碼，請使用 scipy.integrate.solve_ivp 來求解微分方程式。

使用 FORTRAN 程式庫 odepack 中的 lsoda 求解常微分方程式系統。

求解剛性或非剛性一階 ode-s 系統的初始值問題

dy/dt = func(y, t, ...)  [or func(t, y, ...)]

其中 y 可以是向量。

注意

預設情況下，func 的前兩個引數的必要順序與 scipy.integrate.ode 類別和函式 scipy.integrate.solve_ivp 使用的系統定義函式中的引數順序相反。若要使用簽名為 func(t, y, ...) 的函式，則必須將引數 tfirst 設定為 True。

參數:

func可呼叫物件(y, t, …) 或 可呼叫物件(t, y, …)

計算 y 在 t 的導數。如果簽名為 callable(t, y, ...)，則必須將引數 tfirst 設定為 True。func 不得修改 y 中的資料，因為它是 ODE 求解器在內部使用的資料的檢視。

y0陣列

y 的初始條件（可以是向量）。

t陣列

要求解 y 的時間點序列。初始值點應為此序列的第一個元素。此序列必須單調遞增或單調遞減；允許重複值。

argstuple，選用

要傳遞給函式的額外引數。

Dfun可呼叫物件(y, t, …) 或 可呼叫物件(t, y, …)

func 的梯度（雅可比矩陣）。如果簽名為 callable(t, y, ...)，則必須將引數 tfirst 設定為 True。Dfun 不得修改 y 中的資料，因為它是 ODE 求解器在內部使用的資料的檢視。

col_derivbool，選用

如果 Dfun 定義了沿著欄的導數（較快），則為 True，否則 Dfun 應定義沿著列的導數。

full_outputbool，選用

若為 True，則傳回選用輸出的字典作為第二個輸出

printmessgbool，選用

是否列印收斂訊息

tfirstbool，選用

若為 True，則 func（和 Dfun，如果給定）的前兩個引數必須是 t, y 而不是預設的 y, t。

在版本 1.1.0 中新增。

傳回:

y陣列，形狀 (len(t), len(y0))

包含 t 中每個所需時間的 y 值的陣列，初始值 y0 位於第一列。

infodictdict，僅在 full_output == True 時傳回

包含額外輸出資訊的字典

key	意義
‘hu’	每個時間步長成功使用的步長向量
‘tcur’	每個時間步長達到的 t 值向量（始終至少與輸入時間一樣大）
‘tolsf’	公差比例因子向量，大於 1.0，在偵測到要求過高精確度時計算
‘tsw’	上次方法切換時的 t 值（針對每個時間步長給定）
‘nst’	時間步長的累積次數
‘nfe’	每個時間步長的函式評估累積次數
‘nje’	每個時間步長的雅可比矩陣評估累積次數
‘nqu’	每個成功步長的方法階數向量
‘imxer’	在錯誤傳回時，加權局部誤差向量 (e / ewt) 中最大量級的元件索引，否則為 -1
‘lenrw’	所需的雙精度工作陣列長度
‘leniw’	所需的整數工作陣列長度
‘mused’	每個成功時間步長的方法指示器向量：1：adams（非剛性），2：bdf（剛性）

其他參數:

ml, muint，選用

如果這些值都不是 None 或非負數，則假設雅可比矩陣為帶狀矩陣。這些值給出了此帶狀矩陣中下對角線和上對角線的非零對角線數。對於帶狀矩陣情況，Dfun 應傳回一個矩陣，其列包含非零帶（從最低對角線開始）。因此，當 ml >=0 或 mu >=0 時，從 Dfun 傳回的矩陣 jac 的形狀應為 (ml + mu + 1, len(y0))。jac 中的資料必須以 jac[i - j + mu, j] 保存第 i 個方程式對於第 j 個狀態變數的導數的方式儲存。如果 col_deriv 為 True，則必須傳回此 jac 的轉置。

rtol, atolfloat，選用

輸入參數 rtol 和 atol 決定求解器執行的錯誤控制。求解器將根據 max-norm of (e / ewt) <= 1 形式的不等式來控制 y 中估計局部錯誤的向量 e，其中 ewt 是計算為 ewt = rtol * abs(y) + atol 的正錯誤權重向量。rtol 和 atol 可以是與 y 長度相同的向量或純量。預設值為 1.49012e-8。

tcritndarray，選用

應注意積分的臨界點向量（例如，奇異點）。

h0float，（0：求解器決定），選用

要在第一步嘗試的步長。

hmaxfloat，（0：求解器決定），選用

允許的最大絕對步長。

hminfloat，（0：求解器決定），選用

允許的最小絕對步長。

ixprbool，選用

是否在方法切換時產生額外列印。

mxstepint，（0：求解器決定），選用

每個 t 中的積分點允許的最大（內部定義）步數。

mxhnilint，（0：求解器決定），選用

列印的最大訊息數。

mxordnint，（0：求解器決定），選用

非剛性（Adams）方法允許的最大階數。

mxordsint，（0：求解器決定），選用

剛性（BDF）方法允許的最大階數。

另請參閱

solve_ivp

求解 ODE 系統的初始值問題

ode

更物件導向的積分器，基於 VODE

quad

用於尋找曲線下的面積

範例

在瀏覽器中試試看！

受重力和摩擦力作用的擺錘角度 theta 的二階微分方程式可以寫成

theta''(t) + b*theta'(t) + c*sin(theta(t)) = 0

其中 b 和 c 是正常數，而撇號（’）表示導數。若要使用 odeint 求解此方程式，我們必須先將其轉換為一階方程式系統。透過定義角速度 omega(t) = theta'(t)，我們得到系統

theta'(t) = omega(t)
omega'(t) = -b*omega(t) - c*sin(theta(t))

令 y 為向量 [theta, omega]。我們在 Python 中將此系統實作為

>>> import numpy as np
>>> def pend(y, t, b, c):
...     theta, omega = y
...     dydt = [omega, -b*omega - c*np.sin(theta)]
...     return dydt
...

我們假設常數為 b = 0.25 和 c = 5.0

>>> b = 0.25
>>> c = 5.0

對於初始條件，我們假設擺錘幾乎垂直，theta(0) = pi - 0.1，並且最初處於靜止狀態，因此 omega(0) = 0。則初始條件向量為

>>> y0 = [np.pi - 0.1, 0.0]

我們將在區間 0 <= t <= 10 中生成 101 個均勻間隔樣本的解。因此，我們的時間陣列為

>>> t = np.linspace(0, 10, 101)

呼叫 odeint 以產生解。若要將參數 b 和 c 傳遞給 pend，我們使用 args 引數將它們提供給 odeint。

>>> from scipy.integrate import odeint
>>> sol = odeint(pend, y0, t, args=(b, c))

解是一個形狀為 (101, 2) 的陣列。第一欄是 theta(t)，第二欄是 omega(t)。以下程式碼繪製兩個元件。

>>> import matplotlib.pyplot as plt
>>> plt.plot(t, sol[:, 0], 'b', label='theta(t)')
>>> plt.plot(t, sol[:, 1], 'g', label='omega(t)')
>>> plt.legend(loc='best')
>>> plt.xlabel('t')
>>> plt.grid()
>>> plt.show()

返回在新分頁中開啟