# Solving a linear system of ODEs in matrix form

I would like to solve the linear system of ordinary differential equations \begin{align*} \dot\mu&=R^\mathrm T\mu \\ \dot\Sigma&=R^\mathrm T\Sigma+\Sigma^\mathrm T R+\mathrm{diag}(R^\mathrm T\mu) \end{align*} where $\mu$ is a column vector, $\Sigma$ is a (symmetric) covariance matrix, and $R$ is a square rate matrix. Clearly $\mu(t)=e^{R^\mathrm Tt}\mu(0)$, but is there any way to solve for $\Sigma(t)$?

## 2 Answers

The matrix equation to solve reads $$\dot{\Sigma} = 2\,\mathrm{sym} (R^\mathrm{T} \Sigma) + \Lambda(t) \, ,$$ where $\Sigma$ is a $n\times n$ real symmetric matrix, $\Lambda(t) = \mathrm{diag}\!\left(R^\mathrm{T}\! \exp\left(t R^\mathrm{T}\right)\mu(0)\right)$ is a diagonal matrix, and $\mathrm{sym}(M) = \frac{1}{2}(M + M^\mathrm{T})$ denotes the symmetric part. One can note that $\dot\Sigma$ is symmetric, independently from the fact that $\Sigma$ is symmetric or not. In the case where $R$ and $\Sigma$ commute, e.g. $R$ is the identity matrix, then the homogeneous solutions are proportional to $\exp\left(2t\,\text{sym}(R)\right)$. In the general case, it seems hopeless to work on the matrix system itself. A strategy consists in rewriting the differential equation in vector form $$\dot{V} = A V + b(t) \, ,$$ where $V$ is a size-$n (n+1)/2$ vector of coordinates of $\Sigma$, $A$ is a square matrix, and $b(t)$ is a vector. To illustrate, let us consider the case of $2\times 2$ matrices $$\Sigma = \left(\begin{array}{cc} \sigma_1 & \sigma_2\\ \sigma_2 & \sigma_3\end{array}\right),\quad R = \left(\begin{array}{cc} r_1 & r_2\\ r_3 & r_4\end{array}\right),\quad \Lambda(t) = \left(\begin{array}{cc} \lambda_1(t) & 0\\ 0 & \lambda_2(t)\end{array}\right).$$ The differential equation in matrix form gives $$\frac{d}{dt} \left(\begin{array}{c} \sigma_1\\ \sigma_2\\ \sigma_3 \end{array}\right) = \left(\begin{array}{ccc} 2 r_1 & 2 r_3 & 0\\ r_2 & r_1 + r_4 & r_3\\ 0 & 2 r_2 & 2 r_4 \end{array}\right)\! \left(\begin{array}{c} \sigma_1\\ \sigma_2\\ \sigma_3 \end{array}\right) + \left(\begin{array}{c} \lambda_1(t)\\ 0\\ \lambda_2(t) \end{array}\right) ,$$ which can be solved with classical methods.

I have realized that the matrix equation is a continuous differential Lyapunov equation, which is studied in linear-quadratic optimal control. The system above has the following solution:

\begin{align*} \mu(t)&=e^{R^\text T t}\mu_0 \\ \Sigma(t)&=\int_0^t e^{R^\text T(t-\tau)} \,\text{diag}(R^\text T e^{R^\text T \tau}\mu_0)\,e^{R(t-\tau)}\,d\tau \end{align*}

assuming the usual initial conditions $\mu(0)=\mu_0$ and $\Sigma(0)=0$.

• There is a "proof" on page 6 of the book "Matrix Riccati Equations in Control and Systems Theory," but the authors simply state that it follows by differentiation. More useful discussions are given here and here. By the way, thank you for your previous answer -- it was very helpful! – Alex Jun 16 '17 at 16:05