Getting a stable solution for a simple first-order PDE

Question

I have what is in my estimation a pretty simple PDE. It's the Liouville equation for the density of points in phase space with a hyperbolic secant potential. But when I try to solve it with NDSolve, I get nonphysical results. In particular the solution tends to ripple and diffract; this shouldn't be happening for a classical problem with non-negative initial data. Can anyone help me figure out how to fix this?

Here is some code that first defines the pieces that appear in the problem then tries to solve it using NDSolve. The initial data is a gaussian distribution P0[x,p]

poisson[f_, g_] := Inactivate[D[f, x] D[g, p] - D[f, p] D[g, x], D];
liouvilleEqn = 
  Inactivate[D[P[x, p, t], t] + poisson[P[x, p, t], H[x, p]] == 0, D];

P0[x_, p_] := 
 Module[{σp}, σp = ℏ/(2 σ); 
  1/(2 π Sqrt[σ σp])
    Exp[-(x - x0)^2/(2 σ)^2 - (p - p0)^2/(2 σp)^2]]

v[x_] := v0/Cosh[x/ξ0]; h[x_, p_] := p^2/(2 m) + v[x];

sol1 = With[{vars := {m = 2, ℏ = 1, v0 = 1, ξ0 = 1, L = 20, 
     Lp = 50, x0 = -10, p0 = 1, σ = 1}},
    Block[vars,
           classicalProblem = {Activate[liouvilleEqn /. H[x, p] :> h[x, p]]};
           boundaryConditions = { P[-L, p, t] == P[L, p, t], 
                 P[x, -Lp, t] == P[x, Lp, t] == 0};
           initialCondition = {P[x, p, 0] == P0[x, p]};
    NDSolveValue[{classicalProblem, boundaryConditions, 
        initialCondition}, P, {x, -L, L}, {p, -Lp, Lp}, {t, 0, 40}, 
           MaxStepSize -> .8]
   ]]

But as you can see from the plot of the solution at the initial and final times it cannot be correct; it starts off as a gaussian:

but isn't even positive everywhere at $t=40$:

This is how it evolves:

And here's a plot of the trajectories of 100 particles drawn from the initial distribution, to give a hint of what the correct solution should suggest in the end:

Answer 1

It is necessary to reduce the size of the integration region and to establish homogeneous boundary conditions. Then one can observe the process of scattering of particles by the potential

pars = {m = 2, \[HBar] = 1, v0 = 1, \[Xi]0 = 1, L = 10, Lp = 5, 
  x0 = -2, p0 = 1, \[Sigma] = 1, \[Sigma]p = \[HBar]/(2 \[Sigma])}; 
poisson[f_, g_] := Inactivate[D[f, x] D[g, p] - D[f, p] D[g, x], D];
liouvilleEqn = 
  Inactivate[D[P[x, p, t], t] + poisson[P[x, p, t], H[x, p]] == 0, D];

P0[x_, p_] := 
 1/(2 \[Pi] Sqrt[\[Sigma] \[Sigma]p]) Exp[-(x - 
         x0)^2/(2 \[Sigma])^2 - (p - p0)^2/(2 \[Sigma]p)^2]

v[x_] := v0/Cosh[x/\[Xi]0]; h[x_, p_] := p^2/(2 m) + v[x];

classicalProblem = {Activate[liouvilleEqn /. H[x, p] :> h[x, p]]};
boundaryConditions = {P[-L, p, t] == 0, P[L, p, t] == 0, 
   P[x, -Lp, t] == P[x, Lp, t] == 0};
initialCondition = {P[x, p, 0] == P0[x, p]};
sol = NDSolveValue[{classicalProblem, boundaryConditions, 
   initialCondition}, P, {x, -L, L}, {p, -Lp, Lp}, {t, 0, 4}, 
  Method -> {"MethodOfLines", 
    "SpatialDiscretization" -> {"TensorProductGrid", 
      "MinPoints" -> 900, "MaxPoints" -> 1000, 
      "DifferenceOrder" -> 2}}, MaxSteps -> 10^6]

frame = Table[
   ContourPlot[sol[x, p, t], {x, -5, 5}, {p, -5, 5}, 
    ColorFunction -> Hue, Frame -> False, PlotRange -> All, 
    Contours -> 20], {t, 0, 4, .05}];

ListAnimate[frame]