1
$\begingroup$

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!

Improve this question
$\endgroup$
12
  • $\begingroup$ Not the answer to your question but are you sure that there is I in front of the potential? $\endgroup$
    – user59583
    Dec 17 '18 at 13:12
  • $\begingroup$ @Buddha_the_Scientist Thanks for your comment. Yes, it's a complex potential. $\endgroup$
    – SomeOne
    Dec 17 '18 at 13:13
  • 1
    $\begingroup$ You did not mention that your code would return a warning. $\endgroup$
    – Αλέξανδρος Ζεγγ
    Dec 17 '18 at 13:14
  • $\begingroup$ There is this post: mathematica.stackexchange.com/questions/60386/… which can help. It is talking about the same error. $\endgroup$
    – user59583
    Dec 17 '18 at 13:15
  • $\begingroup$ It tries to calculate, but nothing happens. $\endgroup$
    – SomeOne
    Dec 17 '18 at 13:16
1
$\begingroup$

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]

enter image description here

Improve this answer
$\endgroup$
0
$\begingroup$

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

fig1

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.

Improve this answer
$\endgroup$
5
  • $\begingroup$ 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: mathematica.stackexchange.com/q/124737/1871 $\endgroup$
    – xzczd
    Dec 18 '18 at 10:45
  • $\begingroup$ @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. $\endgroup$
    – Alex Trounev
    Dec 18 '18 at 12:13
  • $\begingroup$ 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. $\endgroup$
    – xzczd
    Dec 18 '18 at 13:21
  • $\begingroup$ @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. $\endgroup$
    – Alex Trounev
    Dec 18 '18 at 16:48
  • $\begingroup$ I've started a new question for the inconsistency between FiniteElement and TensorProductGrid: mathematica.stackexchange.com/q/188109/1871 $\endgroup$
    – xzczd
    Dec 18 '18 at 17:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.