# Simulating Phase Space Evolution

I am interested in modeling the time evolution of phase-space $\rho(\vec{q},\vec{p},t)$. I have attempted to use Liouville's theorem $\partial_t\rho=-\sum_{i=1}^{3}(\partial_{q_i}\rho)\dot q_i+(\partial_{p_i}\rho)\dot p_i$, but even in the simple 1D case with no potential my numerical simulation fails.

I would like help with the simplest case, where there is one spatial dimension ($q=x$), and for the sake of simplicity we can assume m=1, so $p=v$, and Liouille's equation becomes $\partial_t \rho = v \partial_x \rho$, with $\rho=\rho(x,v,t)$. What I want is to numerically simulate $\rho(x,v,t)$ given $\rho(x,v,0)$, but I'm not sure how to make this happen. My attempt to do this was by using Mathematica, calling the NDSolve command and passing the simplified liouville equation and an initial distribution, but no sensible results are obtained. There is a warning thrown that I don't have enough boundary conditions on $x$, but I don't know what additional boundary conditions should look like.

• You may want to look at finite difference methods of solving PDEs. However, as it stands this question is really a programming question and probably should be asked at Computational Science. Commented Dec 15, 2014 at 3:44
• I would agree with Kyle that for purposes of simulation of Hamiltonian systems solving the Liouville equation directly may not be a good choice. Basically every non-trivial Hamiltonian system has a very complex (usually chaotic) dynamic. This means that the global solution to the Liouville equation should be a fractal. I am not aware of any general purpose PDE solver (like the one in Mathematica) which would be able to produce sensible results in such a case. Did you look for tools that are specifically created to deal with this kind of problem? Commented Dec 15, 2014 at 4:44

## 1 Answer

You might be best off creating a bunch of trial systems and evolving each one independently using the (much easier) Hamilton's equations, which are good old ordinary differential equations. After all this is how the Liouville equation was constructed in the first place by people like Gibbs---in the limit of an infinite number of trial systems, you get the continuous version for a point density $\rho$.

That's the trick I used in creating the following animation, where 20000 dots represent trial systems [ https://commons.wikimedia.org/wiki/File:Hamiltonian_flow_classical.gif ]:

It won't work for all purposes but it is likely good enough to give you a feeling for the evolution.