NDSolve third order PDE with boundary conditions

Question

I am trying to solve following PDE involving a third derivative with respect to p:

a = 40; 
ps = 2*10^-5;
ic = u[0, p] == 0.3 E^(-0.06 (-(7/2) + p/2)^2);
bc = {u[t, -a] == 0, u[t, a] == 0, Derivative[0, 1][u][t, -a] == 0};
pde = - D[u[t, p], t] == ps D[u[t, p], {p, 3}] + I p^2/4 u[t, p];

solp = NDSolve[{pde, ic, bc}, u, {t, 0, 60}, {p, -a, a}];

Then I animate the solution:

Animate[Plot[Evaluate[Abs[u[t, p] /. First[solp]]^2], {p, - a/2, a},  
PlotRange -> {0, 0.2}], {t, 0, 60, 0.05}]

The solution does not look like it should. For example, for very small ps, the solution should just be the initial function. Instead it is oscillating all over the place. How can I improve this code? My guess is, that I should choose different boundary conditions. I played around with them a lot, without success.

The physics behind this: I want to solve the time-dependent Schrödinger equation in momentum space. Situation: a particle wave packet collides with a weak x^3 potential. Any help very appreciated!

Answer 1

1. You need to make the grid size small enough.
2. Based on the statement, I believe what you really need is periodic b.c..

The following code generates the same result as that in Henrik's answer:

a = 10; 
ps = 2 10^-5;
ic = u[0, p] ==E^(-(p^2/4)) (2/π)^(1/4);
(*
bc = {u[t, -a] == 0, u[t, a] == 0, 
      Derivative[0, 1][u][t, -a] == 0, Derivative[0, 1][u][t, a] == 0};
 *)
pde = - D[u[t, p], t] == ps D[u[t, p], {p, 3}] + I p^2/4 u[t, p];

solp = NDSolveValue[{pde, ic, u[t,-a]==u[t,a](*bc*)}, u, {t, 0, 2Pi 10}, {p, -a, a}, 
         Method->{MethodOfLines, SpatialDiscretization -> 
           {TensorProductGrid, MinPoints->1200, MaxPoints->1200, DifferenceOrder->2}}];
Plot[(solp[2 Pi 10, p]//Abs)^2, {p,-3,2}, PlotRange->All]

BTW, it's interesting that NDSolveValue accepts the bc above i.e. one can set 4 b.c.s in p direction, the output is almost the same.

Answer 2

We change the boundary conditions a bit, then we get

a = 40;
ps = 2*10^-5; f = 0.3 E^(-0.06 (-(7/2) + p/2)^2); f1 = D[f, p];
ic = u[0, p] == f;
bc = {u[t, -a] == f /. p -> -a, 
   Derivative[0, 1][u][t, -a] == f1 /. p -> -a, 
   Derivative[0, 1][u][t, a] == f1 /. p -> a};
pde = -D[u[t, p], t] == ps D[u[t, p], {p, 3}] + I p^2/4 u[t, p];

solp = NDSolve[{pde, ic, bc}, u, {t, 0, 60}, {p, -a, a}];

Large and small scale visualization

{Plot3D[Evaluate[Abs[u[t, p] /. First[solp]]^2], {p, -a, a}, {t, 0, 
   60}, ColorFunction -> Hue, Mesh -> None, PlotRange -> All, 
  PlotPoints -> 50, AxesLabel -> Automatic, Boxed -> False], 
 Plot3D[Evaluate[Abs[u[t, p] /. First[solp]]^2], {p, -a, a}, {t, 0, 
   1}, ColorFunction -> Hue, Mesh -> None, PlotRange -> All, 
  PlotPoints -> 50, AxesLabel -> Automatic, Boxed -> False]}