Calculating Kinetic Energy Using The Solution Obtained From Nonlinear Partial Differential Equation

Question

I want to calculate the kinetic energy of the dynamical wavefunction colliding with a well using

$$E = \int \mathrm{d}x \, \frac12 \left| \frac{\mathrm{d}\psi(x,t)}{\mathrm{d}x} \right|^2$$

which should be constant throughout the dynamics and the value comes close to $(NK^2)/2$.

For this particular case, $N=40$ (number of particles) and $K=0.1$ (speed). Here is my code for dynamics:

L = 200;
mu = -0.2222222222222222222222222211689169810351288660408099720909 - 
   1.2557735735320599617305095419979235954708106033*10^-31 I;
ϕ[x_, t_] := -3 mu Exp[-I mu t]/(1 + 
      Sqrt[1 + 9 mu/2] Cosh[Sqrt[-2 mu x^2]]);
particleno = NIntegrate[ϕ[x, 0]*Conjugate[ϕ[x, 0]], {x, -L, L}];(*N*)

W[x_] := If[0.5 > x > -0.5, -0.3, 0];(*potential*)

k0 = 0.1;(*speed or K*)
x0 = 65;(*initial position*)
eq = I D[u[x, t], t] == -1/2 D[u[x, t], x, x] + 
   Abs[u[x, t]]^2 u[x, t] - Abs[u[x, t]] u[x, t] + W[x] u[x, t];
ic = u[x, 0] == 
  Exp[I k0 x] ϕ[x + x0, 0]; bc = {u[L, t] == ic[[2]] /. {x -> L},
   u[-L, t] == ic[[2]] /. {x -> -L}};

ψ = NDSolveValue[{eq, bc, ic}, u, {x, -L, L}, {t, 0, 2000}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MinPoints" -> 400, "MaxPoints" -> 2001, 
       "DifferenceOrder" -> 4}}, MaxSteps -> 10^6];

(*Plotting the dynamics*)
DensityPlot[Abs[ψ[x, t]]^2, {t, 0, 2000}, {x, -L, L}, 
 AspectRatio -> 1/2, Frame -> True, FrameTicks -> Automatic, 
 PlotPoints -> 200, ImageSize -> 500, 
 ColorFunction -> "AvocadoColors", 
 FrameLabel -> {{Style["x", FontFamily -> "Times New Roman", 
     FontSlant -> "Italic", FontWeight -> Bold, FontSize -> 30], 
    None}, {Style["t", FontFamily -> "Times New Roman", 
     FontSlant -> "Italic", FontWeight -> Bold, FontSize -> 30], 
    None}}]

Here is my attempt to calculate and plot the energy throughout the dynamics

(*calculating ∫|(dψ(x,t))/dx|^2\[DifferentialD]x *)
data = Flatten[Table[{x, t, ψ[x, t]}, {x, -L, L}, {t, 0, 2000}], 1];
if = Interpolation[data]
dψdx = Derivative[1, 0][if]
finalkin = 
 Table[{t, NIntegrate[0.5 Abs[dψdx[x, t]]^2, {x, -L, L}]}, {t, 0,
    2000, 500}]

(*plotting energy*)
ListLinePlot[finalkin, AspectRatio -> 1/2, Frame -> True, 
 FrameTicks -> Automatic, ImageSize -> 500, 
 PlotStyle -> {Orange, Thickness[0.005]}, 
 LabelStyle -> {24, Bold, Large, Black}, 
 FrameLabel -> {{Style["E", FontFamily -> "Times New Roman", 
     FontSlant -> "Italic", FontWeight -> Bold, FontSize -> 30], 
    None}, {Style["t", FontFamily -> "Times New Roman", 
     FontSlant -> "Italic", FontWeight -> Bold, FontSize -> 30], 
    None}}]

Answer 1

Here we have an example of the scattering of a wave packet by a potential well. To compute kinetic energy and total energy we use code

L = 200;
mu = -0.2222222222222222222222222211689169810351288660408099720909 - 
   1.2557735735320599617305095419979235954708106033*10^-31 I;
\[Phi][x_, 
   t_] := -3 mu Exp[-I mu t]/(1 + 
      Sqrt[1 + 9 mu/2] Cosh[Sqrt[-2 mu x^2]]);
particleno = 
  NIntegrate[\[Phi][x, 0]*Conjugate[\[Phi][x, 0]], {x, -L, L}];(*N*)

W[x_] := If[0.5 > x > -0.5, -0.3, 0];(*potential*)

k0 = 0.1;(*speed or K*)
x0 = 65;(*initial position*)
eq = I D[u[x, t], t] == -1/2 D[u[x, t], x, x] + 
    Abs[u[x, t]]^2 u[x, t] - Abs[u[x, t]] u[x, t] + W[x] u[x, t];
ic = u[x, 0] == 
  Exp[I k0 x] \[Phi][x + x0, 0]; bc = {u[L, t] == ic[[2]] /. {x -> L},
   u[-L, t] == ic[[2]] /. {x -> -L}};

\[Psi] = 
  NDSolveValue[{eq, bc, ic}, u, {x, -L, L}, {t, 0, 2000}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MinPoints" -> 400, "MaxPoints" -> 2001, 
       "DifferenceOrder" -> 4}}, MaxSteps -> 10^6];

Kinetic energy computed with NIntegrate is given by

ekin[t_?NumericQ] := 
 1/2 NIntegrate[Abs[Derivative[1, 0][\[Psi]][x, t]]^2, {x, -L, L}];

This computation is very slow, for example

ekin[2000] // AbsoluteTiming

(*Out[]= {29.8448, 0.31453}*)

So, it takes about 30s to compute one point. To avoid this inconvenience let define numerical integration with using Gauss quadrature

Needs["NumericalDifferentialEquationAnalysis`"];

gg = GaussianQuadratureWeights[1000, -200, 200]; 

point = gg[[All, 1]];
weight = gg[[All, 2]];

With this definition we compute kinetic energy in 200 points into 2.5s

Ekin = Table[{t, 
     weight . 
      Table[Evaluate[1/2 Abs[Derivative[1, 0][\[Psi]][x, t]]^2], {x, 
        point}]}, {t, 0, 2000, 10}]; // AbsoluteTiming

Visualization

ListLinePlot[Ekin, Frame -> True, FrameLabel -> {"t", "Ekin"}]

Therefore, the kinetic energy is not saved. Let check the number of particles

Ntotal = 
   Table[{t, 
     weight . Table[Evaluate[Abs[\[Psi][x, t]]^2], {x, point}]}, {t, 
     0, 2000, 10}]; // AbsoluteTiming

ListLinePlot[Ntotal, PlotRange -> {0, 40}, 
 AxesLabel -> {"t", "Ntotal"}]

Therefore, the number of particles is very close to 40. Let check the total energy

ETotal = 
   Table[{t, 
     weight . 
      Table[Evaluate[
        1/2 Abs[Derivative[1, 0][\[Psi]][x, t]]^2 + 
         Abs[\[Psi][x, t]]^4 - Abs[\[Psi][x, t]]^3 + 
         W[x] Abs[\[Psi][x, t]]^2], {x, point}]}, {t, 0, 2000, 
     10}];
ListLinePlot[ETotal, PlotRange -> {-10, 0}, Frame -> True, 
 FrameLabel -> {"t", "ETotal"}]

The total energy is not saved as well. The reason is that dynamic term $\frac{i}{2}\int(\psi \psi_t^*-\psi^*\psi_t)dx$ is not a constant in this system.