Scipy求解常微分方程組有scipy.integrate.solve_ivp和scipy.integrate.odeint，后者是較老的版本主要是采用 FORTRAN 的odepack庫里面的lsoda 方法，而前者是后面更新的函數，支持的方法也更多，按照官方的文檔介紹大致有如下的方法。

import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
def fun(t, y):
    dydt = [1]
    return dydt# 初始條件
y0 = [2]
yy = solve_ivp(fun, (0, 5000), y0, method='Radau')
t = yy.t
data = yy.y
plt.plot(t, data[0, :])
plt.xlabel("時間s")
plt.legend(["求解變量"])
plt.show()

然后逐漸復雜一點，假設方程為

發現得到的結果是這樣的，很明顯不合理

這是因為求解的時候輸出的插值的問題，我們改一下代碼，結果如下

是不是就合理多了完整的代碼如下

import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import numpy as np
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']

def fun(t, y):
    print(t)
    dydt = [t]
    return dydt

# 初始條件
y0 = [0]

yy = solve_ivp(fun, (0, 5000), y0, method='Radau',t_eval = np.arange(1,5000,1) )
xx = solve_ivp(fun, (0, 5000), y0, method='Radau')
t = yy.t
data = yy.y
t2 = xx.t
data2 = xx.y
plt.plot(t, data[0, :])
plt.plot(t2, data2[0, :])
plt.xlabel("時間s")
plt.legend(["求解變量"])
plt.show()

現在看著還不高端，可以多定義兩個函數讓他看起來更高端一些，方程如下

import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import numpy as np
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']

def fun(t, y):
    y1 = y[0]
    y2 = y[1]
    dydt = [y2, t*t-y2+t]
    return dydt

# 初始條件
y0 = [0,1]

yy = solve_ivp(fun, (0, 500), y0, method='Radau',t_eval = np.arange(1,500,1) )
t = yy.t
data = yy.y
plt.plot(t, data[0, :])
plt.plot(t, data[1, :])
plt.xlabel("時間s")
plt.legend(["求解變量"])
plt.show()

喜歡的朋友可以給個關注或者聯系我