My goal is to create a demonstration of the Liouville theorem in 2D phase space.
I made up an interesting potential energy function $U(x) = (x-4)x^3 + 27$, so that the minimal energy of system is exactly zero.
The Hamiltonian equations are (in dimensionless units):
$$dp/dt = -4(x-3)x^2,\ dx/dt = p$$
The equation for Contour plot is $(1/2)p^2 + U(x) = E$. Now I want to demonstrate that for given initial set of $(x, p)$ in phase plane, forming an area $D(0)$, which is time-evolving, and this area is conserved, i.e, $D(0) = D(t)$. The problem is that this non-linear system of ODEs is not solvable, so I cannot tell for a given initial $(x,p)$, how it will move in phase space after one second, it has to be solved numerically.
The goal is to draw several images -- first with chosen volume $D(0)$ formed by a large number of initial values of (x,p). Next image would be a system evolved in time, say in one second. It'd be nice to see the volume $D(0)$ being distorted due to the Hamiltonian equations and the non-linear force. One approach is 1) to parametrically describe the volume $D(0)$, 2) numerically solve the Hamilton equations for all those points in $D(0)$ to get their position for latter time 3) draw a new image with volume $D(t)$. It should be sufficient to solve the Hamiltonian equations only on the boundary of $D(0)$, which should preserve continuity in time. Another approach is to choose the initial volume $D(0)$, randomly choose, say, a hundred points in this volume, then solve Hamilton equations for all these points until time $t$, and draw these hundred points at their new position.
The problem is, I don't have the slightest idea how to do this practically in Mathematica 8. I'd appreciate any help.