I am trying to simulate the movement of a coherent state in a quantum harmonic oscilator, but for some reason the answer diverges and there is a warning about not enought boundary conditions.

Also, since I plan on moving to more complex situations in the future, I would like to avoid b.c.s of the kind:

p[t, L] == 0

The reason is that for higher energy situations it will be more taxing to keep L much greater than my area of interest an having the wave get too close to it would cause undesired reflexions.

My code, for a 1D case, is as follows:

w = 1;
L = 3;
c = 1;
usol = NDSolveValue[{I D[p[t, x], t] + 1/2 D[p[t, x], x, x] == 
    1/2 w^2 x^2 p[t, x], p[0, x] == E^(-w (x - c)^2/2)*(w/Pi)^(1/4)}
  , p, {t, 0, 10}, {x, -L, L}]
  • 1
    $\begingroup$ Yes, the b.c. isn't enough, for a initial value problem defined on $-\infty<x<\infty$ you need to add artificial b.c. approximating "infinity". There're many posts about this type of issue, see e.g. mathematica.stackexchange.com/q/128516/1871 $\endgroup$
    – xzczd
    May 8 '17 at 5:54
  • $\begingroup$ @xzczd, thank you for the hint! I adapted the b.c. and it worked, but I must admit to not undertanding why it did... abc = D[p[t,x], x] + direction p[t,x] == 0 /. {{x -> lb, direction -> -1}, {x -> rb, direction -> 1}}; $\endgroup$
    – ivbc
    May 8 '17 at 16:00
  • $\begingroup$ I'm a bit surprised :D . I posted the example just to show that artificial boundary is needed, I didn't know ABC is also applicable for 1D Schrodinger's equation. $\endgroup$
    – xzczd
    May 8 '17 at 16:48
  • $\begingroup$ PLenty of references about it, like: sciencedirect.com/science/article/pii/S0021999108004804. But the subject is complex and I dont know if what I tried is a good answer... $\endgroup$
    – ivbc
    May 8 '17 at 17:48
  • 2
    $\begingroup$ Maybe you can ask a question e.g. "What's the state-of-art/most popular artificial boundary condition for 1D Schrodinger's equation" in scicomp.stackexchange.com ? $\endgroup$
    – xzczd
    May 9 '17 at 5:50

Since OP has found this interesting post, let me try to implement the exterior complex scaling method mentioned there.

First, make the transform $x= \left\{\begin{array}{cc} & \begin{array}{cc} -R_0+e^{i \theta } (\xi +R_0) & -\xi >R_0 \\ R_0+e^{i \theta } (\xi -R_0) & \xi >R_0 \\ \xi & -R_0\leq \xi\leq R_0 \\ \end{array} \\ \end{array}\right. $, I'll use DChange for the task:

set = {I D[p[t, x], t] + 1/2 D[p[t, x], x, x] == 1/2 w^2 x^2 p[t, x], 
   p[0, x] == E^(-w (x - c)^2/2)*(w/Pi)^(1/4)};

right = Exp[I θ] (ξ - R0) + R0;
left = Exp[I θ] (ξ + R0) - R0;
middle = ξ;

coevalue = CoefficientList[
    Piecewise[{{right, ξ > R0}, {left, -ξ > R0}}, middle], ξ];

neweq = DChange[set, x == coe@0 + coe@1 ξ, x, ξ, p[t, x]] /. 
   Thread[(coe /@ {0, 1}) -> coevalue];

(* Alternative approach for deducing neweq: *)
help = DChange[#, x == #2, x, ξ, p[t, x]] &;
change = Piecewise[{{help[#, right], ξ > R0}, {help[#, left], -ξ > R0}}, 
    help[#, middle]] &;  

neweq = Simplify`PWToUnitStep@Map[change@# &, set, {2}];

Notice currently DChange doesn't directly support piecewise function so the coding is a little roundabout, but I think it's still simpler than transforming by hand.


Simplify`PWToUnitStep is an undocumented function that expands Piecewise into a combination of UnitStep. I use this function to "extract" ξ from the piecewise part, or CoefficientList in first approach and NDSolveValue in second approach will fail.

The next step, which is also the final step, is to solve the equation:

w = 1;
L = 3;
c = 1;
tend = 20;
boundarylayer = L/5; thvalue = 1/2;
mol[n_Integer, o_: "Pseudospectral"] := {"MethodOfLines", 
    "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, 
        "MinPoints" -> n, "DifferenceOrder" -> o}}
newsol = NDSolveValue[{neweq, 
      p[t, -L - boundarylayer] == p[t, L + boundarylayer] == 0} /. {R0 -> L, θ -> 
       thvalue}, p, {ξ, -L - boundarylayer, L + boundarylayer}, {t, 0, tend}, 
    Method -> mol[25, 4]]; // AbsoluteTiming
DensityPlot[Norm@newsol[t, x], {t, 0, tend}, {x, -L, L}, PlotRange -> All, 
 PlotPoints -> 100, ColorFunction -> "AvocadoColors", FrameLabel -> {"t", "x"}]

Mathematica graphics

Notice I've manually set the number of spatial grid points, or NDSolveValue will automatically choose a too large one because the initial condition is no longer smooth.

Besides, the choosing of R0 ($R_0$), boundarylayer (distance from $R_0$ to the boundary of computational domain) and thvalue ($\theta$) turns out to be a kind of art to make the reflection small enough. There might be some hint in the original paper about the method, but I simply find the proper value by trial and error.

  • $\begingroup$ (+1) You have transformed x and tto xi and theta, right? $\endgroup$
    – zhk
    May 10 '17 at 7:49
  • $\begingroup$ @zhk Nope, in this transformation, only x is transformed to ξ. θ is a parameter. You can have a look at the linked post for more information. $\endgroup$
    – xzczd
    May 10 '17 at 8:03
  • $\begingroup$ Does this simplify the problem? In what way this transformation helps? $\endgroup$
    – zhk
    May 10 '17 at 8:04
  • 3
    $\begingroup$ @zhk This transform makes the equation more complicated, of course :) , but by solving this equation instead, we obtain a better approximation for infinite range. (OP wants to solve the equation in $-\infty<x<\infty$. ) $\endgroup$
    – xzczd
    May 10 '17 at 8:10
  • 1
    $\begingroup$ @ivbc Actually this is a function I place in SystemOpen@"init.m" file, because I adjust these options frequently when using NDSolve. (See Method -> mol[25, 4]? It's the same as writing Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> 25, "MinPoints" -> 25, "DifferenceOrder" -> 4}}. ) I adjust these options because, as mentioned above, if I don't manually set "MinPoints" and "MaxPoints", NDSolveValue will choose a too large one because the initial condition is not smooth. (You can take away the option and see what will happen. ) $\endgroup$
    – xzczd
    May 11 '17 at 3:10

Let's remember Schrodinger's equation:

$i\hbar\frac{\partial}{\partial t} \Psi(\mathbf{r},t) = \left [ \frac{-\hbar^2}{2\mu}\nabla^2 + V(\mathbf{r},t)\right ] \Psi(\mathbf{r},t)$

For the harmonic oscilator $V = x^2$, so your equation is missing a p[x,t] on the RHS besides the boundary conditions. L also needs to be larger too. This seems to work.

w = 2;
L = 10;
c = 3;
usol = NDSolveValue[{I D[p[t, x], t] + 1/2 D[p[t, x], x, x] == 
    1/2 w^2 x^2 p[t, x], 
   p[0, x] == Exp[(-w (x - c)^2/2)*(w/Pi)^(1/4)], p[t, L] == 0, 
   p[t, -L] == 0}, p, {t, 0, 10}, {x, -L, L}]

DensityPlot[Norm@usol[t, x], {t, 0, 10}, {x, -L/2, L/2}, 
 PlotRange -> All, PlotPoints -> 100, ColorFunction -> "Rainbow"]

enter image description here


For completeness purposes only, this is the solution with the abc boundary condition from @xzczd.

w = 1;
L = 3;
c = 1;
usol = NDSolveValue[{I D[p[t, x], t] + 1/2 D[p[t, x], x, x] - 
     1/2 w^2 x^2 p[t, x] == 0, 
   p[0, x] == E^(-w (x - c)^2/2)*(w/Pi)^(1/4), 
   Derivative[0, 1][p][t, -L] - p[t, -L] == 0, 
   Derivative[0, 1][p][t, L] + p[t, L] == 0}, 
  p, {t, 0, 20}, {x, -L, L}]

DensityPlot[Norm@usol[t, x], {t, 0, 20}, {x, -L/2, L/2}, 
 PlotRange -> All, PlotPoints -> 100, ColorFunction -> "Rainbow", 
 FrameLabel -> {"t", "x"}] 

enter image description here

Note how in the first case the wave function starts to diffuse as time passes, so it's not really a coherent state.

  • $\begingroup$ Thank you for spotting the mistake and for the solution. I improved the question a bit to tray to avoid the artificial b.c. p[t, L] == 0, as it must always be far from the area being analysed. $\endgroup$
    – ivbc
    May 8 '17 at 15:56
  • $\begingroup$ @ivb if with the ABC boundary condition it works, what is your question now? I don't understand. $\endgroup$ May 9 '17 at 2:04
  • $\begingroup$ @tsuresuregusa I guess it's because there's no guarentee that the simple ABC mentioned in my answer is always applicable. (It's a artificial boundary for 1D wave equation. ) Though currently it seems to work, it may be just a coincidence. (In the comment above, OP has found a paper about ABC for 1D Schrodinger's equation, in which the proposed ABC has a very different form. ) $\endgroup$
    – xzczd
    May 9 '17 at 5:47
  • $\begingroup$ @xzczd your last comment nails it, I think we are now out of the scope of this site. $\endgroup$ May 9 '17 at 14:53
  • 2
    $\begingroup$ Answer accepted since it solves the problem. The question about why it works is important, and can be further considered here: scicomp.stackexchange.com/questions/15973/… $\endgroup$
    – ivbc
    May 9 '17 at 22:51

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.