Solving the Liouville equation for classic harmonic oscillator via method of characteristics I'm interested in solving Liouville's equation 
$$\frac{\partial W}{\partial t} + \{ H,W\}=0$$
with $$H=\frac{p^2}{2m} + \frac{1}{2}m\omega^2x^2$$
using the method of characteristics. However I can't seem to find any resources on how to do it. Would appreciate any pointers.
 A: You can start by reading the wikipedia article on the method of characteristics. You will see that in our case the tangent of the characteristic is $(1,-\partial H/\partial q,\partial H/ \partial p)$ where the components are in order $t,p,q$. When you formulate the equation of the characteristic, you will actually find out you get equations of motion of a single particle with this Hamiltonian. 
Once you solve these equations for $t(s),p(s),q(s)$ with a parameter $s$, you will know that a solution will satisfy $W(t(s),p(s),q(s))=const.$. This is actually very intuitive - the probability distribution $W$ flows with the particles evolving according to their equations of motion.
The solutions will be determined by two constants corresponding to e.g. phase $\phi$ and amplitude $A$ of the oscillation, you will have $t^{\phi, A}(s),p^{\phi,A}(s),q^{\phi,A}$. Depending on the kind of initial data or problem you are solving, you should be able to get the solution by integrating over the family of characteristics parametrized by $A,\phi$ with a certain weight function $f(A,\phi)$ at a certain $s$ and then get the solution on the rest of the $t,p,q$ space by going through all values of $s$. 
Say you are given an initial $W_{0}(p,q)$ and want to evolve it through time. Then you have to find the map $p,q \to \phi,A$ substitute it into $W_0$ and this will be your weight function over the $\phi, A$.
