scipy求解矩阵微分方程

形如$$X^{\prime }=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也是可以求解矩阵微分方程的，只不过需要转换一下格式