 Solution

Solving $\dot{x}=p/m$, $\dot{p}=-m{\omega }^{2}x$ gives the linear phase flowIts Jacobian determinant isHence every patch keeps its area $dxdp$, directly verifying Liouville's theorem. The energy ellipses are also invariant, but area preservation applies to arbitrary patches, not just whole energy shells.

In scaled canonical coordinates $X=\sqrt{m\omega }x$, $P=p/\sqrt{m\omega }$ the flow is a rigid clockwise rotation, $(X,P)=(X_0\cos \omega t+P_0\sin \omega t,-X_0\sin \omega t+P_0\cos \omega t)$. The following original plot follows an off-center patch and shows equal areas at successive times.