How would I solve for following, else any implemented algorithms or solvers in matlab (even ways to solve it) for Lyapunov differential equation of form:
$$P'(t) + A(t)^TP(t) + P(t)A(t) + Q(t) = 0,$$
where $Q(t)> 0$, with $P(t)$ being symmetric periodic positive definite, $A(t)$ is linearized dynamics of the system.
On
Recently, (1) proposed novel formulations of solution algorithms for differential matrix equations based on an $LDL^{T}$ decomposition that keep the computations in real arithmetic.
According to the authors, the proposed methods are a more practical alternative for large-scale problems arising in applications, and the numerical results from the paper show better performance of the proposed methods compared to earlier formulations.
Albeit the methods in the paper focus more on the differential Riccati equation, they naturally restrict to the differential Lyapunov equation:
However, all our methods and techniques naturally restrict to the differential Lyapunov equation (DLE) (see comments in Section 2.2). A more detailed explanation and extensive numerical experiments for the DLE will be presented elsewhere to keep the presentation within usual page limits.
(1): Lang, Norman; Mena, Hermann; Saak, Jens (September 2015). "On the benefits of the $LDL^{T}$ factorization for large-scale differential matrix equation solvers". Linear Algebra and its Applications 480: 44–71. doi:10.1016/j.laa.2015.04.006.
On
The analytical solution to this problem is
\begin{equation} P(t) = \int_0^t \Phi^T(\tau,t)Q(\tau)\Phi(\tau,t)d\tau \end{equation} where $\Phi(\tau,t)$ is the state transition matrix from time $t$ to time $\tau$ associated with the linear system $\dot{x}(t) = A(t) x(t)$. (In your case, $A(t)$ is coming from a linearization.) You may or may not be able to obtain a closed-form representation of the state transition matrix $\Phi(\tau,t)$. Sometimes, an analytical expression for $\Phi(\tau,t)$ can be easily obtained (e.g. in the linear time-invariant case, $\Phi(\tau,t)=e^{A(\tau-t)}$). Even if you do not have an exact expression for $\Phi(\tau,t)$, you can still verify that this form is correct by differentiation. Since $t$ appears both in the upper limit and within the integral, you need to use Leibniz's rule. You also need to use the fact that $\frac{d}{dt} \Phi(\tau,t) = -\Phi(\tau,t)A$. Substitution completes the proof.
Suppose that $A$ is constant. Let us take $$P(t)=\exp(-tA^T)R(t)\exp(-tA),$$ where $R(t)$ is some matrix. We can write the derivative of $P$: $$P'(t)=-A^TP(t)-P(t)A+\exp(-tA^T)R'(t)\exp(-tA) = -A^TP(t)-P(t)A-Q(t).$$ Now we have $$R'(t) = - \exp(tA^T)Q(t)\exp(tA),$$ which is easy to solve.
This method can be generalized to $A=A(t)$, but it will require strong hypothesis on $A(t)$: we will need that $A(t_1)A(t_2)=A(t_2)A(t_1)$ for all $t_1$ and $t_2$.