# NDSolve third order PDE with boundary conditions

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!

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.

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]}


• thank you, but nothing changed for me. The Results don't make any sense. For example one would expect that a x^3 potential would reduce and then invert the average momentum over time. Also if you use ps ~ 10^-20, then there shouldn't be any oscillations. But there are a lot. – Luke Feb 27 '20 at 14:48
• @Luke You are mistaken in your assumptions. The solution to equation  -D[u[t, p], t] == I p^2/4 u[t, p] has the form Exp[-I p^2/4 t]. These are oscillations with a frequency of p^2/4. This is what we see in the figures. – Alex Trounev Feb 27 '20 at 16:17
• We plot the absolute value, which should not oscillate, right? – Luke Feb 27 '20 at 16:25
• @Luke There is initial data. This is not a constant, but a package, which is also visible in the figures. – Alex Trounev Feb 27 '20 at 16:56
• @ Alex I don't think initial data matters here as the time dependency is just a product with the initial wave packet. If you solve the Equation -D[u[t, p], t] == I p^2/4 u[t, p] , and plot Abs[u[t,p]], you will see that the wave packet just stays put with respect to time. That also makes physical sense since there is no potential present which could affect any momenta. – Luke Feb 27 '20 at 17:14