Time-dependent Schrödinger equation and density plot

Question

I'm interested in the density plot of the solution of a 1D time-dependent Schrödinger equation with a given potential. So I have:

ψ = NDSolve[{I D[u[t, x], t] == -D[u[t, x], {x, 2}] + 
 I Sin[x] u[t, x], u[0, x] == Exp[-x^2] Exp[I x], u[t, -Pi] == u[t, Pi]}, u, {t, 0, 5}, {x, -Pi, Pi}]

DensityPlot[Evaluate[Abs[u[t, x]] /. ψ], {x, -Pi, Pi}, {t, 0, 5}, PlotPoints -> 200, Mesh -> False]

However, Mathematica spits out NDSolve::mxsst warning and does not return anything. Any help is appreciated!

Not the answer to your question but are you sure that there is I in front of the potential? — , Dec 17 '18 at 13:12
@Buddha_the_Scientist Thanks for your comment. Yes, it's a complex potential. — SomeOne, Dec 17 '18 at 13:13
There is this post: https://mathematica.stackexchange.com/questions/60386/ndsolve-wave-equation-triangular-wave-pulse-inital-condition which can help. It is talking about the same error. — , Dec 17 '18 at 13:15
Try this command: \[Psi] = NDSolve[{ D[u[t, x], t] == I*D[u[t, x], {x, 2}] + Sin[x] u[t, x], u[0, x] == Exp[-x^2] Exp[I x], u[t, -\[Pi]] == u[t, \[Pi]]}, u, {t, 0, 5}, {x, -\[Pi], \[Pi]}, Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> 101}}] Then it plots something. — , Dec 17 '18 at 13:18
@Buddha_the_Scientist This is a correct solution, why not turn it into an answer? — xzczd, Dec 18 '18 at 06:27
@xzczd It is already presented in the post that I have mentioned. — , Dec 18 '18 at 07:56
@Buddha_the_Scientist It's always better to post it as an answer: https://mathematica.meta.stackexchange.com/q/76/1871 — xzczd, Dec 18 '18 at 08:05
Your code worked for me as soon as I added Method -> "MethodOfLines" to NDSolve. It took a few minutes though. — Alexei Boulbitch, Dec 18 '18 at 09:00

score 1 · Answer 1 · edited Dec 18 '18 at 10:31

There is this post which is talking about the same error. Following their advice I believe this must work. Please check it and tell me if it is correct or not.

ψ = NDSolve[{ D[u[t, x], t] == I*D[u[t, x], {x, 2}] + Sin[x] u[t, x], 
              u[0, x] == Exp[-x^2] Exp[I x], u[t, -π] == u[t, π]}, 
            u, {t, 0, 5}, {x, -π, π}, Method -> {"MethodOfLines", 
            "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> 101}}]

and then your plot

DensityPlot[Evaluate[Abs[u[t, x]] /. ψ], {x, -Pi, Pi}, {t, 0, 5}, 
            PlotPoints -> 200, Mesh -> False]

Alex Trounev · Answer 2 · 2018-12-18T12:55:18.300

It is possible that the condition of periodicity 'u[t, -π] == u[t, π]' is wrong. A solution with such a potential cannot be symmetric about zero. On the other hand, we do not know which method to solve is correct. We can compare solutions

sol = NDSolveValue[{I D[u[t, x], t] == -D[u[t, x], {x, 2}] + 
     I Sin[x] u[t, x], u[0, x] == Exp[-(x - Pi)^2] Exp[I (x - Pi)], 
   PeriodicBoundaryCondition[u[t, x], x == 2*\[Pi], 
    Function[x, x - 2*\[Pi]]]}, u, {t, 0, 5}, {x, 0, 2*Pi}];
sol1 = NDSolveValue[{D[u[t, x], t] == 
    I*D[u[t, x], {x, 2}] + Sin[x] u[t, x], 
   u[0, x] == Exp[-(x - Pi)^2] Exp[I (x - Pi)], 
   u[t, 0] == u[t, 2*\[Pi]]}, u, {t, 0, 5}, {x, 0, 2*\[Pi]}, 
  Method -> {"MethodOfLines", 
    "SpatialDiscretization" -> {"TensorProductGrid", 
      "MaxPoints" -> 10000}}];
sol2 = NDSolveValue[{I D[u[t, x], t] == -D[u[t, x], {x, 2}] + 
     I Sin[x] u[t, x], u[0, x] == Exp[-(x - Pi)^2] Exp[I (x - Pi)], 
   u[t, 0] == u[t, 2*\[Pi]]}, u, {t, 0, 5}, {x, 0, 2*\[Pi]}, 
  Method -> {"MethodOfLines", 
    "SpatialDiscretization" -> {"TensorProductGrid", 
      "MaxPoints" -> 100}}]
    {DensityPlot[Abs[sol[t, x]], {x, 0, 2*Pi}, {t, 0, 5}, 
      PlotPoints -> 50, MaxRecursion -> 2, PlotLegends -> Automatic, 
      ColorFunction -> Hue, FrameLabel -> {"x", "t"}, PlotRange -> All], 
     DensityPlot[Abs[sol1[t, x]], {x, 0, 2*Pi}, {t, 0, 5}, 
      PlotPoints -> 50, MaxRecursion -> 2, PlotLegends -> Automatic, 
      ColorFunction -> Hue, FrameLabel -> {"x", "t"}, PlotRange -> All], 
     DensityPlot[Abs[sol1[t, x]], {x, 0, 2*Pi}, {t, 0, 5}, 
      PlotPoints -> 50, MaxRecursion -> 2, PlotLegends -> Automatic, 
      ColorFunction -> Hue, FrameLabel -> {"x", "t"}, PlotRange -> All]}

Obviously, one of the methods is erroneous. But only sol1 and sol2 have a message NDSolveValue::mxsst: Using maximum number of grid points 10000 allowed by the MaxPoints or MinStepSize options for independent variable x. and NDSolveValue::mxsst: Using maximum number of grid points 100 allowed by the MaxPoints or MinStepSize options for independent variable x.

The result is incorrect because: 1. Since you've moved the domain to $[0,2π]$, the PDE should be modified accordingly, you forgot to modify the I Sin[x] u[t, x] term. 2. There's a bug in PeriodicBoundaryCondition: https://mathematica.stackexchange.com/q/124737/1871 — xzczd, Dec 18 '18 at 10:45
@xzczd 1. It is possible that the condition of periodicity u[t, -π] == u[t, π] is wrong. A solution with such a potential cannot be symmetric about zero. 2. We can compare two solutions. — Alex Trounev, Dec 18 '18 at 12:13
The mxsst warning can be eliminated with e.g. Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> 200, "MinPoints" -> 200, "DifferenceOrder" -> 4}}, and the solution still looks the same. — xzczd, Dec 18 '18 at 13:21
@xzczd I agree with the use of the method "MethodOfLines", but the condition u[t, -π] == u[t, π] is obviously erroneous. It should also be understood why the first method does not agree with the second. — Alex Trounev, Dec 18 '18 at 16:48
I've started a new question for the inconsistency between FiniteElement and TensorProductGrid: https://mathematica.stackexchange.com/q/188109/1871 — xzczd, Dec 18 '18 at 17:44

Time-dependent Schrödinger equation and density plot

2 Answers2

Linked