形如 X ′ = A X + b X'=AX+b X′=AX+b的矩阵微分方程目前在scipy中是没有接口直接处理的,需要自己手动转换一下格式,但对应这种格式MATLAB的ode求解器是支持直接求解的。
本文基于https://wenku.baidu.com/view/8fc76d6bfbd6195f312b3169a45177232f60e4f9.html?re=view
的例子(原例子是用MATLAB求解的,采用scipy求解一组矩阵微分方程:
python程序如下:
from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
k1 = k2 = k3 = 1
m1 = m2 = m3 = 1
# f1 = f2 = f3 = 5*np.sin(1.25*t)
A = np.array([[0,0,0,1,0,0],
[0,0,0,0,1,0],
[0,0,0,0,0,1],
[-(k1+k2)/m1, k2/m1, 0,0,0,0],
[k2/m2,-(k2+k3)/m2,k3/m2,0,0,0],
[0,k3/m3,-k3/m3,0,0,0]]
)
B = np.array([[0,0,0],
[0,0,0],
[0,0,0],
[1/m1,0,0],
[0,1/m2,0],
[0,0,1/m3]
])
# solve dX/dt = A*X+b
def func(t,y):
#dX/dt = A*X+b
# x1,x2,x3,x4,x5,x6 = y
f = np.array([5*np.sin(1.25*t),5*np.sin(1.25*t),5*np.sin(1.25*t)])
# 主要下面直接用到了列表推导
return [np.dot(A[i,:], y)+np.dot(B[i,:],f) for i in range(6)]
N = 30
t_span = (0,N)
t_eval = np.linspace(0,N,1000)
y0 = [0,0,0,0,0,0]
sol = solve_ivp(func, t_span, y0, t_eval=t_eval)
plt.plot(sol.t, sol.y.T)
plt.grid('on')
plt.show()
结果如下,和原文结果一致:
总结:scipy也是可以求解矩阵微分方程的,只不过需要转换一下格式