在 Python 中模仿 ode45() 函式

Aditya Raj 2022年7月18日
在 Python 中模仿 ode45() 函式

MatLab 中使用常微分方程來解決許多科學問題。ode45() 在 MatLab 中用於求解微分方程。

本文將展示我們如何在 Python 中模仿 ode45() 函式。

在 Python 中模仿 ode45() 函式

為了模仿 python 中的 ode45() 函式,我們可以使用 scipy 模組中定義的 solve_ivp() 方法。solve_ivp() 方法整合了一個常微分方程 (ODE) 系統。

  • solve_ivp() 方法將函式作為其第一個輸入引數。輸入引數中給出的函式必須返回一個包含微分方程係數的陣列。
  • 在第二個輸入引數中,solve_ivp() 方法採用包含兩個數值的元組或列表。這些值表示積分割槽間,其中元組中的第一個值表示區間的開始,元組的第二個值表示區間中的最大值。
  • 在第三個輸入引數中,solve_ivp() 方法採用一個表示初始值的陣列。
  • 執行後,solve_ivp() 方法返回具有各種屬性的一堆物件。
    1. t 屬性包含一個包含時間點的 numpy 陣列。
    2. y 屬性包含一個 numpy 陣列,其值和時間點在 t 中。
    3. sol 屬性包含一個 Odesolution 物件,其中包含微分方程的解。如果在 solve_ivp() 方法中 dense_output 引數設定為 false,則 sol 屬性包含 None

為了更好地理解這一點,請參閱以下示例。

from scipy.integrate import solve_ivp

def function(t, y):
    return 2 * y

interval = [0, 10]
initial_values = [10, 15, 25]
solution = solve_ivp(function, interval, initial_values)
print("Time:", solution.t)
print("Y:", solution.y)

輸出:

Time: [ 0.          0.07578687  0.56581063  1.18741382  1.85887096  2.55035821
  3.25007544  3.95320486  4.65775424  5.36289544  6.06828346  6.77377445
  7.47930839  8.18486026  8.89041961  9.59598208 10.        ]
Y: [[1.00000000e+01 1.16366412e+01 3.10073783e+01 1.07492109e+02
  4.11689241e+02 1.64114780e+03 6.65071446e+03 2.71362627e+04
  1.11036049e+05 4.54874443e+05 1.86437495e+06 7.64300835e+06
  3.13352156e+07 1.28474398e+08 5.26752964e+08 2.15973314e+09
  4.84541488e+09]
 [1.50000000e+01 1.74549617e+01 4.65110674e+01 1.61238163e+02
  6.17533861e+02 2.46172171e+03 9.97607169e+03 4.07043941e+04
  1.66554074e+05 6.82311665e+05 2.79656243e+06 1.14645125e+07
  4.70028233e+07 1.92711598e+08 7.90129446e+08 3.23959970e+09
  7.26812231e+09]
 [2.50000000e+01 2.90916029e+01 7.75184457e+01 2.68730272e+02
  1.02922310e+03 4.10286951e+03 1.66267862e+04 6.78406569e+04
  2.77590123e+05 1.13718611e+06 4.66093739e+06 1.91075209e+07
  7.83380389e+07 3.21185996e+08 1.31688241e+09 5.39933284e+09
  1.21135372e+10]]

在上面的例子中,我們首先定義了一個名為 function 的函式,它以 ty 作為其輸入引數,並根據 y 返回一個值。

然後,我們分別使用變數 intervalinitial_values 為 ODE 定義了一個區間和初始值。我們將 functionintervalinitial_values 作為輸入引數傳遞給 solve_ivp() 函式,最後,我們在變數解決方案中獲得輸出。

在輸出中,你可以觀察到時間值分佈在 0 到 10 的區間內。同樣,輸出包含與每個時間值對應的 y 值。

我們還可以在解決方案的屬性 t 中明確指定時間點。為此,我們需要將包含所需時間值的陣列傳遞給 solve_ivp() 方法的 t_eval 引數,如下所示。

from scipy.integrate import solve_ivp

def function(t, y):
    return 2 * y

interval = [0, 10]
initial_values = [10, 15, 25]
time_values = [1, 2, 3, 6, 7, 8]
solution = solve_ivp(function, interval, initial_values,t_eval=time_values)
print("Time:", solution.t)
print("Y:", solution.y)

輸出:

Time: [1 2 3 6 7 8]
Y: [[7.38683416e+01 5.46053271e+02 4.03089733e+03 1.62618365e+06
  1.20160156e+07 8.87210156e+07]
 [1.10802512e+02 8.19079906e+02 6.04634600e+03 2.43927547e+06
  1.80240234e+07 1.33081523e+08]
 [1.84670854e+02 1.36513318e+03 1.00772433e+04 4.06545912e+06
  3.00400390e+07 2.21802539e+08]]

你可以看到時間值僅包含那些作為輸入引數傳遞給 t_eval 引數的值。類似地,屬性 y 僅包含指定 t 值的值。

這種方法可以幫助你獲取區間中某些點的值。

相關文章 - Python Function