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

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}}]
• Where is the question? Is your approach not working? Is the solution wrong? Commented Apr 29, 2023 at 15:29
• Why do you think the kinetic energy should be constant. You have a non constant potential. Commented Apr 29, 2023 at 15:40
• Is this problem connect to mathematica.stackexchange.com/questions/266736/… ? Commented Apr 30, 2023 at 2:51
• Yes @AlexTrounev But with a well present from -0.5 to 0.5 in x. Commented Apr 30, 2023 at 4:29
• 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 Commented Apr 30, 2023 at 4:41

## 1 Answer

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.

• 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. Commented Apr 30, 2023 at 9:31