For a deterministic dynamical system $\dot{\mathbf{x}} = F(\mathbf{x})$ with $\mathbf{x} \in \mathbb{R}^n$ the Liouville equation governs the evolution of probability densities $\rho(\mathbf{x}, t)$:
$$ \frac{\partial \rho(\mathbf{x}, t)}{\partial t} = -\nabla \mathbf{\cdot}[F(\mathbf{x}) \rho(\mathbf{x}, t)] = \hat{L}\rho(\mathbf{x}, t) $$
where $\hat{L}$ is the Liouville operator. For an initial distribution $\rho_0 = \rho(\mathbf{x}, t=0)$ this equation has the solution
$$ \rho_t = \mathrm{exp} \left(t \hat{L} \right) \rho_0. $$
For a given vector field $F$ and initial distribution $\rho_0$ the first equation can be solved numerically for $\rho_t$ for example with finite difference methods - however I would like to solve for the operator $\hat{L}$ (approximated as a matrix). How can I go about doing this?
Best regards!
One possibility is found in this article: first estimate the Perron-Frobenius operator for small $t$, and then differentiate numerically (Equation 13).