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.

  • 2
    $\begingroup$ I've given you a starting point. If you are looking for a complete solution, people are unlikely to put in the effort to write one, unless you demonstrate that you put in the effort to learn a bit of Mathematica and at least try to solve the problem. $\endgroup$
    – Szabolcs
    Jul 4, 2014 at 17:20
  • $\begingroup$ Thank you for the help, I try to use NDsolve for the same equations, but many initial values. $\endgroup$
    – user16320
    Jul 4, 2014 at 17:35
  • $\begingroup$ Be careful not to use E as a symbol for energy in your code. Mathematica reserves this symbol for the base of the natural logarithms. $\endgroup$
    – m_goldberg
    Jul 4, 2014 at 18:46

1 Answer 1


I suggest going with the second approach: choose a set of points and follow their trajectories.

To solve differential equations numerically, use NDSolve (as a documentation search quickly reveals).

To plot a set of points in 2D, use ListPlot.

This information, together the the documentation of the mentioned functions and some of the basic tutorials is sufficient to solve the problem.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.