# Numerical solution of PDE - Liouville's equation

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.

• The syntax is wrong for your partial derivatives. 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; – Bob Hanlon May 28 '18 at 14:25

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.


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.

• Thanks! I definitely felt the need for boundary conditions but stubbornly refused to include any since I couldn't think of physically sensible ones. For other systems I expect the Gaussian to disperse, so there will definitely be cases where rho should be non-zero at the boundary. The problem, I suppose, is that physically "parts" of the solution leave the domain and then re-enter it later, however, when solving, those parts are lost as soon as they leave. Is there a clever way to choose boundary conditions to get around this? – Dogtard May 28 '18 at 15:09
• @Dogtard: Most of my own computational work uses periodic boundary conditions, but those wouldn't be suitable here. Imposing either Dirichlet or Neumann boundary conditions would lead to spurious reflections from the boundaries. I know that "absorbing boundary conditions" are sometimes used for wave-like equations to avoid these reflections, but I don't know how easy they are to implement in Mathematica. ... – Michael Seifert May 28 '18 at 15:26
• I'd recommend asking that question at Physics.SE or Computational Science.SE, since it's not really a Mathematica question per se. If you get a good answer at the other SE sites, feel free to open another question over here about how to implement them. – Michael Seifert May 28 '18 at 15:27
• That's true. Having thought about it a bit more, I think your approach of using Dirichlet boundaries and moving them a good distance away should work. If I make sure the solution doesn't "run into" the boundaries, the actual solution would be approximately to 0 there anyway, and so the result should be very close to correct. Thanks again. – Dogtard May 28 '18 at 15:45