Proof

Liouville's theorem says that if you've got a set of possible states of a physical system, a "volume" in "phase space" which you might take to represent uncertainty, then that volume doesn't change as the system evolves. This applies to systems that obey conservation of energy. 

Liouville's theorem can be easily proven using the Hamiltonian formulation of classical mechanics, which we can phrase in terms of phase spaces. (Aside: Hamiltonian mechanics generalizes better to quantum mechanics.)

Phase space is just the space of co-ordinates $q$ and momenta $p$ of a physical system. If you add on some extra structure, you get dynamics in this space, where each point represents a state of a physical system and a parameterized curve represents how that state evolves over time. 

The core bit of structure you need to add here is the Hamiltonian $H(p,q)$, which is the energy as a function of position $q\left(t\right)$ and momentum $p\left(t\right)$. $H$ tells us everything about the state of our physical system, and in fact generates our equations of motion:

$\partial H/\partial p=\dot{q},\partial H/\partial q=-\dot{p}\left(1\right)$

This isn't anything mysterious: for a single particle in a potential well, we have  

$H(p,q)=m/2p^2+V\left(q\right)$

So 

$\partial H/\partial p=p/m=\dot{q}\partial H/\partial q=-\partial V/\partial q=\dot{p}$

But we won't need anything that specific for Liouville's theorem. The general equations of motion (1) are all we'll need for Liouville's theorem. 

As is usual for physics, start with an infinitesimal volume element, dp dq (fig 1). This can represent the possible states of a physical system. We wish to show this volume is conserved as the states evolve. 

Note that $dp\left(t\right)$ and $dq\left(t\right)$ are the differences in position/momenta along corners of the rectangle. To calculate how they change, we need to see how positions and momenta change. In an infinitesimal time step, by (1),

$p(t+dt)=p\left(t\right)+\dot{p}dt=p\left(t\right)-\partial H/\partial qdt,\left(2a\right)$

$q(t+dt)=q\left(t\right)+\dot{q}dt=q\left(t\right)+\partial H/\partial pdt.\left(2b\right)$

Using 2a and 2b, we can show that

$dp(t+dt)=p(t+dt)-p\left(t\right)=dp\left(t\right)-{\partial }^{2}H/\partial q\partial pdp\left(t\right)dt,$

$dq(t+dt)=q(t+dt)-q\left(t\right)=dq\left(t\right)+{\partial }^{2}H/\partial p\partial qdq\left(t\right)dt.$

This implies that for all t, 

$dq(t+dt)dp(t+dt)=dq\left(t\right)dp\left(t\right)(1+({\partial }^{2}H/\partial p\partial q-{\partial }^{2}H/\partial p\partial q\left)dt\right)$

So for all t

$dq(t+dt)dp(t+dt)=dp\left(t\right)dq\left(t\right).$

In words, infinitesimal rectangles in phase space have constant volume. This analysis generalizes to higher dimensions and different shapes. QED.

The Second Law

You might wonder how Liouville's theorem squares with the second law of thermodynamics, if volumes in phase space are meant to represent uncertainty. One theorem says uncertainty remains constant, one law says it must increase. The answer is that: 

1. We're using the full description of our physical system in Liouville's theorem, and in the context of thermodynamics we've thrown away almost all the info we have about our system. This grows the size of the phase space we're uncertain over.
2. In phase space, our system might evolve into very complex states that require fine-grained knowledge to keep track of. Instead of neat squares, we wind up with fractal, space filling shapes. But the thermodynamic model involves far coarser evolution, resulting in all the little gaps between thin "strands" of our phase space volume getting filled in because we can't keep track of the fine details of where our volume is/isn't.

So the thermodynamic uncertainty is larger and gets larger.