2
$\begingroup$

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

enter image description here

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}}]
$\endgroup$
5
  • $\begingroup$ Where is the question? Is your approach not working? Is the solution wrong? $\endgroup$
    – Domen
    Apr 29, 2023 at 15:29
  • $\begingroup$ Why do you think the kinetic energy should be constant. You have a non constant potential. $\endgroup$ Apr 29, 2023 at 15:40
  • $\begingroup$ Is this problem connect to mathematica.stackexchange.com/questions/266736/… ? $\endgroup$ Apr 30, 2023 at 2:51
  • $\begingroup$ Yes @AlexTrounev But with a well present from -0.5 to 0.5 in x. $\endgroup$ Apr 30, 2023 at 4:29
  • $\begingroup$ Well, my psi(x,t)=exp(i k x-i w t)psi(x). My norm is 40 and if I put this expression into the energy expression it should come to (N*k^2)/2. But energy keeps on increasing with time. I am not sure about the code for calculating the expression.@DanielHuber @Domen $\endgroup$ Apr 30, 2023 at 4:41

1 Answer 1

3
$\begingroup$

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

Figure 1

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

Figure 2

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

Figure 3

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.

$\endgroup$
1
  • 1
    $\begingroup$ Thanks for the answers and also giving me an explanation.Also thanks @AlexTrounev for providing another way of integration process, it was taking me long time to integrate. $\endgroup$ Apr 30, 2023 at 9:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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