I am facing some problems solving a time-dependent matrix differential equation. The problem is that the time-dependent coefficient is not just following some time-dependent shape, rather it is the solution of another differential equation.
Up until now, I have considered the trivial case where my coefficient G(t) is just G(t)=pulse(t) where this pulse(t) is a function I define. Here is the code:
# Matrix differential equation
def Leq(t,v,pulse):
v=v.reshape(4,4) #covariance matrix
M=np.array([[-kappa,0,E_0*pulse(t),0],\. #coefficient matrix
[0,-kappa,0,-E_0*pulse(t)],\
[E_0*pulse(t),0,-kappa,0],\
[0,-E_0*pulse(t),0,-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,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,))
And I get a nice result. The problem is that is it not exactly what I need. Want I need is that dot(G(t))=-kappa*G(t)+pulse(t), i.e. the coefficient is the solution of a differential equation.
I have tried to implement this differential equation in a sort of vectorized way in Leq by returning another parameter G(t) that would be fed to M, but I was getting problems with the dimensions of the arrays.
Any idea of how should I proceed?
Thanks,
In principle you have the right idea, you just have to split and join the state and derivative vectors.
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()])
The initial vector has likewise to be such a composite.