# Getting a stable solution for a simple first-order PDE

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:

• What do you want to describe? – Alex Trounev Nov 20 at 19:30
• @AlexTrounev I’m not sure what you mean; I’m trying to numerically solve the given PDE but I think the solution output by NDSolve is incorrect. The physical setup is a gaussian ensemble of classical particles colliding with a potential barrier $V(x) \sim 1/\cosh(x)$. – Diffycue Nov 20 at 19:34
• How do you know what the solution should be? – Alex Trounev Nov 20 at 19:42
• @AlexTrounev I don't know exactly what it should look like but I know from Liouville's theorem some qualitative properties. For instance, it can never be negative. (This is directly contradicted by the NDSolve solution). It should also look like an incompressible flow, which in my opinion it doesn't. Also in the early-time dynamics the initial gaussian clump evolves into three clumps before many of the particles have even impacted the barrier, which is strange as well. – Diffycue Nov 20 at 19:50
• In an animation the solution should look something like this: The blob approaches $x=0$ with the top of the blob (higher $p$) getting there earlier than the bottom of the blob. Some of the blob (the parts with $p^2/2m > v_0$) will make it over the barrier and the rest will be reflected. And the whole affair should be pretty uncomplicated; no rapid oscillations in the $(x,p)$ plane or any regions of negative $P(x,p)$ at any time. – Diffycue Nov 20 at 19:59

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]

• Wow that looks great! (Also your computer must be much more powerful than mine...) – Diffycue Nov 20 at 20:42
• We have to pay for quality :) – Alex Trounev Nov 20 at 20:52