Python - 求解瞬态矩阵微分方程

编程问答 2022-05-29

问题描述

我在求解瞬态矩阵微分方程时遇到了一些问题。 问题在于瞬态系数不只是遵循某种瞬态形状，而是另一个微分方程的解。

到目前为止，我已经考虑了一个微不足道的情况，其中我的系数 G(t) 只是 G(t)=pulse(t)，其中这个 pulse(t) 是我定义的函数。代码如下：

# Matrix differential equation
def Leq(t,v,pulse): 
    v=v.reshape(4,4) #covariance matrix 
    M=np.array([[-kappa,E_0*pulse(t),0],\.  #coefficient matrix 
                [0,-kappa,-E_0*pulse(t)],\
                [E_0*pulse(t),\
                [0,-E_0*pulse(t),-kappa]])
    
    Driff=kappa*np.ones((4,4),float) #constant term
    
    dv=M.dot(v)+v.dot(M)+Driff #solve dot(v)=Mv+vM^(T)+D
    return dv.reshape(-1)  #return vectorized matrix

#pulse shape
def Gaussian(t):
    return np.exp(-(t - t0)**2.0/(tau ** 2.0))

#scipy solver
cov0=np.zeros((4,float) ##initial vector
cov0 = cov0.reshape(-1);   ## vectorize initial vector
Tmax=20 ##max value for time
Nmax=10000 ##number of steps
dt=Tmax/Nmax  ##increment of time
t=np.linspace(0.0,Tmax,Nmax+1)

Gaussian_sol=solve_ivp(Leq,[min(t),max(t)],cov0,t_eval= t,args=(Gaussian,))

我得到了一个不错的结果。问题是这不正是我所需要的。我需要的是 dot(G(t))=-kappa*G(t)+pulse(t)，即系数是微分方程的解。

我试图通过返回另一个将被馈送到 M 的参数 G(t) 在 Leq 中以一种矢量化的方式实现这个微分方程，但我遇到了数组维度的问题。

知道我应该如何进行吗？

谢谢，

解决方法

原则上您的想法是正确的，您只需拆分并连接状态向量和导数向量即可。

def Leq(t,u,pulse): 
    v=u[:16].reshape(4,4) #covariance matrix 
    G=u[16:].reshape(4,4) 
    ... # compute dG and dv
    return np.concatenate([dv.flatten(),dG.flatten()])

初始向量也必须是这样的复合。

matrixmatrixnumpyodepythonscipyscipy