在 Python 中模仿 ode45() 函式
Aditya Raj
2022年7月18日
MatLab 中使用常微分方程來解決許多科學問題。ode45() 在 MatLab 中用於求解微分方程。
本文將展示我們如何在 Python 中模仿 ode45() 函式。
在 Python 中模仿 ode45() 函式
為了模仿 python 中的 ode45() 函式，我們可以使用 scipy 模組中定義的 solve_ivp() 方法。solve_ivp() 方法整合了一個常微分方程 (ODE) 系統。
solve_ivp()方法將函式作為其第一個輸入引數。輸入引數中給出的函式必須返回一個包含微分方程係數的陣列。- 在第二個輸入引數中，
solve_ivp()方法採用包含兩個數值的元組或列表。這些值表示積分割槽間，其中元組中的第一個值表示區間的開始，元組的第二個值表示區間中的最大值。 - 在第三個輸入引數中，
solve_ivp()方法採用一個表示初始值的陣列。 - 執行後，
solve_ivp()方法返回具有各種屬性的一堆物件。t屬性包含一個包含時間點的 numpy 陣列。y屬性包含一個 numpy 陣列，其值和時間點在t中。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 的函式，它以 t 和 y 作為其輸入引數，並根據 y 返回一個值。
然後，我們分別使用變數 interval 和 initial_values 為 ODE 定義了一個區間和初始值。我們將 function、interval 和 initial_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 值的值。
這種方法可以幫助你獲取區間中某些點的值。