Skip to main content
Mathematica

Return to Answer

Added potential energy
Source Link
Jens
Jens
  • 97.9k
  • 7
  • 215
  • 510

Edit 2: Adding a potential energy term.

The above numerical solutions are basically for a free particle, except that the spatial grid is forcing us to choose some boundary conditions on the sides of the square. Periodic boundary conditions are a common choice. But the whole effort is overkill for a free particle because the solutions can be obtained analytically. It gets more interesting if we add an arbitrary potential energy to see how the wave packet is deflected over time.

The periodic boundary conditions in this calculation allow you to add a potential energy to the Hamiltonian, as long as it doesn't conflict with the periodicity of box. Here is an example where I added the potential

$$V(x, y) = - 20 \cos(\frac{\pi x}{10}) \cos(\frac{\pi y}{10})$$

with a box of side length $10$. This potential vanishes on the box boundaries, and has an attractive center at the origin.

Also, I started the Gaussian slightly offset from the center, with a momentum tangential to the equipotential lines, so we expect it to go around with some angular momentum:

ψ = u /. 
   First@NDSolve[{I D[u[t, x, y], t] == -D[u[t, x, y], {x, 2}] - 
        D[u[t, x, y], {y, 2}] - 
        20 Cos[Pi x/10] Cos[Pi y/10] u[t, x, y], 
      u[0., x, y] == Exp[-((x - 1)^2. + y^2.)] Exp[I y], 
      u[t, -5., y] == u[t, 5., y], u[t, x, -5.] == u[t, x, 5.]}, 
     u, {t, 0., 3.}, {x, -5., 5.}, {y, -5., 5.}, 
     Method -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"TensorProductGrid", 
         "DifferenceOrder" -> "Pseudospectral"}}];

pl = Table[
   Plot3D[Abs[ψ[t, x, y]], {x, -5, 5}, {y, -5, 5}, 
    PlotRange -> {0, 1}], {t, 0, 3, .1}];

Export["revolve.gif", pl, AnimationRepetitions -> Infinity, 
 "DisplayDurations" -> .1]

revolve

The packet still disperses but is clearly trapped in the potential minimum, as expected.

Edit 2: Adding a potential energy term.

The above numerical solutions are basically for a free particle, except that the spatial grid is forcing us to choose some boundary conditions on the sides of the square. Periodic boundary conditions are a common choice. But the whole effort is overkill for a free particle because the solutions can be obtained analytically. It gets more interesting if we add an arbitrary potential energy to see how the wave packet is deflected over time.

The periodic boundary conditions in this calculation allow you to add a potential energy to the Hamiltonian, as long as it doesn't conflict with the periodicity of box. Here is an example where I added the potential

$$V(x, y) = - 20 \cos(\frac{\pi x}{10}) \cos(\frac{\pi y}{10})$$

with a box of side length $10$. This potential vanishes on the box boundaries, and has an attractive center at the origin.

Also, I started the Gaussian slightly offset from the center, with a momentum tangential to the equipotential lines, so we expect it to go around with some angular momentum:

ψ = u /. 
   First@NDSolve[{I D[u[t, x, y], t] == -D[u[t, x, y], {x, 2}] - 
        D[u[t, x, y], {y, 2}] - 
        20 Cos[Pi x/10] Cos[Pi y/10] u[t, x, y], 
      u[0., x, y] == Exp[-((x - 1)^2. + y^2.)] Exp[I y], 
      u[t, -5., y] == u[t, 5., y], u[t, x, -5.] == u[t, x, 5.]}, 
     u, {t, 0., 3.}, {x, -5., 5.}, {y, -5., 5.}, 
     Method -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"TensorProductGrid", 
         "DifferenceOrder" -> "Pseudospectral"}}];

pl = Table[
   Plot3D[Abs[ψ[t, x, y]], {x, -5, 5}, {y, -5, 5}, 
    PlotRange -> {0, 1}], {t, 0, 3, .1}];

Export["revolve.gif", pl, AnimationRepetitions -> Infinity, 
 "DisplayDurations" -> .1]

revolve

The packet still disperses but is clearly trapped in the potential minimum, as expected.

Source Link
Jens
Jens
  • 97.9k
  • 7
  • 215
  • 510

I think it's worth pointing out that the problem can be solved "straightforwardly" (i.e., really using only NDSolve) once you know the options that Stefan used in ProcessEquations (which I upvoted because those options are the main ingredient):

Below I show the original problem of a Gaussian wave packet with no initial momentum, and then a modified case where an initial momentum has been imparted, making the initial condition complex as well. I call the complex wave function $\psi$ and plot its absolute value:

ψ = u /. 
   First@NDSolve[{I D[u[t, x, y], t] == -D[u[t, x, y], {x, 2}] - 
        D[u[t, x, y], {y, 2}], u[0., x, y] == Exp[-(x^2. + y^2.)], 
      u[t, -5., y] == u[t, 5., y], u[t, x, -5.] == u[t, x, 5.]}, 
     u, {t, 0., 2.}, {x, -5., 5.}, {y, -5., 5.}, 
     Method -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"TensorProductGrid", 
         "DifferenceOrder" -> "Pseudospectral"}}];

pl =
  Table[Plot3D[Abs[ψ[t, x, y]], {x, -5, 5}, {y, -5, 5}, 
    PlotRange -> {0, 1}], {t, 0, 2, .1}];

Export["spreading.gif", pl, AnimationRepetitions -> Infinity, 
 "DisplayDurations" -> .4]

spreading

ψ = 
 u /. First@
   NDSolve[{I D[u[t, x, y], t] == -D[u[t, x, y], {x, 2}] - 
       D[u[t, x, y], {y, 2}], 
     u[0., x, y] == Exp[-(x^2. + y^2.)] Exp[3 I x], 
     u[t, -10., y] == u[t, 10., y], u[t, x, -10.] == u[t, x, 10.]}, 
    u, {t, 0., 1.}, {x, -10., 10.}, {y, -10., 10.}, 
    Method -> {"MethodOfLines", 
      "SpatialDiscretization" -> {"TensorProductGrid", 
        "DifferenceOrder" -> "Pseudospectral"}}]

pl =
  Table[Plot3D[Abs[ψ[t, x, y]], {x, -10, 10}, {y, -10, 10}, 
    PlotRange -> {0, 1}], {t, 0, 1, .1}];

Export["moving.gif", pl, AnimationRepetitions -> Infinity, 
 "DisplayDurations" -> .4]

moving