Probability density in Hamiltonian Mechanics I am currently studying Liouville's theorem compare wikipedia
and there this mysterious probability density $\rho$ appears and I was wondering how one can determine this quantity analytically for a given problem. 
Let's think of the harmonic oscillator as a primitve example, how can I get an analytical representation of $\rho$?
 A: The probability density is the state of the system in this setting. While the hamiltonian contains all the information about how the system is set up, what forces are acting, and so on, the probability density contains all the physically measurable information about whereabouts in phase space your system actually is.
This is a complete shift in paradigm to other approaches to mechanics. In quantum mechanics, for example, you have a fundamental state - the wavefunction, $\psi$ - and you find $\rho=|\psi|^2$ in terms of it. In hamiltonian mechanics you might play the game of forgetting when you started the clock, and seeing where you might find the particle most often (see, e.g. "Quantum and classical probability distributions for position and momentum," RW Robinett, Am. J. Phys. 63 no. 9, p. 823 (1995), available here). This is different. Since you don't know where you've put your system, or where it might be now, the probability density $\rho$ is your fundamental handle on the state of the system, and you can't trade it for other objects.
Thus, given a specific initial state (and a specific hamiltonian), you have to solve Liouville's equation to obtain $\rho$ for later times. Depending on the system and the initial state this may or may not be possible analytically, or you may even find yourself in numerical difficulties. In general, though, this is relatively hard, since your equation of motion has changed from a system of ODEs (as in hamiltonian mechanics) to a relatively complex partial differential equation.

That said, the harmonic oscillator is in fact a very, very special case for liouvillian mechanics. This is because of its harmonicity: more precisely, because all the orbits have the same period. This implies that phase space moves rigidly, rotating around the equilibrium position at the oscillation frequency, and it carries the probability density with it.
Thus, for a harmonic hamiltonian $H=\frac12(p^2+q^2)$, the solution to Liouville's equation
$$\left[\frac\partial{\partial t}+p\frac\partial{\partial q}-q\frac\partial{\partial p} \right]\rho=0$$
under an arbitrary initial condition $\rho_0(q,p)$ is easily found to be
$$\rho(q,p,t)=\rho_0(q\cos(t)-p\sin(t),q\sin(t)+p\cos(t)).$$
(You should, of course, do the corresponding calculation! For a good exercise, add units and check that everything works out.) It should be pretty clear that what this solution does is trace the classical trajectory back from the point $(q,p)$ at time $t$ to where it would have started in at time 0, and retrieve the original particle density there.
I imagine there is an equivalent method to obtain (at least formal) solutions to Liouville's equations in terms of the corresponding hamiltonian trajectories. (A brief foray into the web didn't turn it up, but I'm pretty sure it exists.) Thus you'd reduce the hard PDE problem to the ODEs it was originally constructed from. Do note, though, that in general things will be at least slightly more complicated than the above, because if the phase-space motion isn't rigid - if trajectories bunch up or spread out - the probability density will correspondingly increase or decrease.
A: Liouville's theorem describes the time-evolution of the density, in canonical ensemble where the energy thereof is conservative, of system points in phase space. A common example is a closed system with Hamiltonian $\mathcal{H}=\frac{1}{2}p^2+V(q)$. According to Liouville's equation, the probability density distribution $\rho(q,p;t)$ of phase space remains invariant along trajectory. That is
$$\frac{d\rho}{dt}=\frac{\partial\rho}{\partial t}+\frac{\partial\rho}{\partial q}\frac{\partial\mathcal{H}}{\partial p}-\frac{\partial\rho}{\partial p}\frac{\partial\mathcal{H}}{\partial q}=\frac{\partial\rho}{\partial t}+\frac{\partial\rho}{\partial q}p-\frac{\partial\rho}{\partial p}\nabla V=0$$
It's a linear first-order pde and we can solve it by method of characteristics. The characteristic equation is
$$dt=\frac{dq}{p}=-\frac{dp}{\nabla V}$$
which gives two first integrals.
$$\mathcal{H}=\frac{1}{2}p^2+V(q)$$
$$q_0=q-\int pdt$$
Thus the solution is
$$\rho(q,p;t)=F(\mathcal{H},q-\int pdt)$$
for any function bivariate $F(x,y)$.
Assume the initial distribution is multivariate normal distribution $\mathcal{N}({\bf0,\mathcal{I}})$, namely
$$\rho(q,p;0)=F(\mathcal{H},q_0)=(2\pi)^{-\frac{n}{2}}\exp\{-\frac{1}{2}q_0^Tq_0\}$$
Then the final solution
$$\rho(q,p;t)\propto\exp\{-\mathcal{H(q(t),p(t))}\}$$
A: Consider some classical mechanical system with canonical coordinates and momenta $(p,q)$ on its phase space $\Gamma$.  We can imagine building a large number of identical copies of this system.  We call such a set of copies an ensemble.  
Next, we can prepare the initial conditions of each element of the ensemble however we like.  If we imagine that the state of each element of the ensemble is represented as a point in $\Gamma$, then the whole ensemble will look like a cloud of points in phase space.  If we imagine that the number of points is very large, then we can imagine representing the cloud in terms of a density on phase space.  If, for example, all of the members of the ensemble are initially prepared in nearly the same state, then the cloud will look like a clump around the corresponding point in phase space.
Once we specify the initial behavior of the ensemble however, the Hamiltonian evolution will determine how the cloud evolves with time.  In fact, the equation that governs the evolution of the phase density is
\begin{align}
  \frac{\partial\rho}{\partial t} = \{H,\rho\}
\end{align}

Constructing an example is easy.  Simply choose your favorite probability density on phase space at the initial time, and use the evolution equations to obtain it for later times!

