Numerical solution of PDE - Liouville's equation

Question

I'm trying to solve Liouville's equation:

$\frac{\partial\rho}{\partial t}-f(q)\frac{\partial\rho}{\partial p}+p\frac{\partial\rho}{\partial q}=0$

I believe this should be solvable given an initial distribution, $\rho(q,p,0)$

For a harmonic oscillator $f(q)=q$ and an initial Gaussian distribution I have tried:

eqn = Derivative[0, 0, 1][ρ][q, p, t] - q*Derivative[0, 1, 0][ρ][q, p, t] + p*Derivative[1, 0, 0][ρ][q, p, t] == 0;


ρ0 = 10/Pi E^(-10*((q - 0.5)^2 + p^2)); 

sol = NDSolve[{eqn, ρ[q, p, 0] == ρ0}, ρ, {q, -1, 1}, {p, -1, 1}, {t, 0, 2*Pi}]

The result should have the whole Gaussian simply rotating around the origin over time.

I get two errors:

NDSolve::bcart: Warning: an insufficient number of boundary conditions have been specified for the direction of independent variable p. Artificial boundary effects may be present in the solution.

NDSolve::eerr: Warning: scaled local spatial error estimate of 51097.533173032294` at t = 6.283185307179586` in the direction of independent variable q is much greater than the prescribed error tolerance. Grid spacing with 35 points may be too large to achieve the desired accuracy or precision. A singularity may have formed or a smaller grid spacing can be specified using the MaxStepSize or MinPoints method options.

and clearly broken results:

I've tried using the finite difference method, which gives the error:

NDSolve::femcscd: The PDE is convection dominated and the result may not be stable. Adding artificial diffusion may help.

with less ugly but still incorrect plots.

Although I think this particular harmonic oscillator case can be solved with the method of characteristics, I'd like to be able to solve for more complex cases with a different $f(p)$ which would require numerical solving.

I'm not too fussed about the initial distribution being a Gaussian, if for some reason that's causing issues. In fact, I'd ideally like some kind of step function - although I feel that's only going to cause more problems.

Answer 1

Any numerical PDE solution technique for an initial-value problem (either method of lines or finite-element) is going to require spatial boundary conditions, not just initial conditions. This is what the error

NDSolve::bcart: Warning: an insufficient number of boundary conditions have been specified for the direction of independent variable p. Artificial boundary effects may be present in the solution.

is complaining about.

Mathematically, the simplest way to do this is to impose Dirichlet boundary conditions $\rho = 0$ at $q = \pm1$ and $p = \pm 1$, for all values of $t$:

eqn = Derivative[0, 0, 1][\[Rho]][q, p, t] - 
       q*Derivative[0, 1, 0][\[Rho]][q, p, t] + 
       p*Derivative[1, 0, 0][\[Rho]][q, p, t] == 0;

\[Rho]0 = 10/Pi E^(-10*((q - 0.5)^2 + p^2));

bcs = {\[Rho][-1, p, t] == 0, \[Rho][1, p, t] == 0, 
       \[Rho][q, -1, t] == 0, \[Rho][q, 1, t] == 0} ;
ics = {\[Rho][q, p, 0] == \[Rho]0};

sol = NDSolve[Join[{eqn}, bcs, ics], \[Rho], {q, -1, 1}, {p, -1, 1}, {t, 0, 2*Pi}]

Table[Plot3D[\[Rho][q, p, t] /. First[sol], {q, -1, 1}, {p, -1, 1}, PlotRange -> All], {t, 0, 2 \[Pi], \[Pi]/2}]

Whether or not the boundary condition I chose above is reasonable for your system is a question of the physics underlying the problem. These boundary conditions are probably not the greatest choice for this problem, if only because the Gaussian is not quite zero at the boundary; in fact, when you run this code, Mathematica returns the following warning because of this:

NDSolve::ibcinc: Warning: boundary and initial conditions are inconsistent. >>

If nothing else, I would move the boundary outwards (say to $q = \pm 3$ & $p = \pm 3$) for this Gaussian, just to minimize the level of disagreement between the initial conditions and the boundary conditions.