So I'm doing some quantum mechanics and I have these datapoints for a wave function:

wfData = {{15.3,0.},{15.6,0.},{15.9,0.},{16.2,0.},{16.5,0.},{16.8,0.},{17.1,0.},{17.4,0.},{17.7,0.},{18.,0.},{18.3,0.},{18.6,0.},{18.9,0.},{19.2,0.},{19.5,0.},{19.8,0.},{20.1,0.},{20.4,0.},{20.7,0.},{21.,0.},{21.3,0.},{21.6,0.},{21.9,0.},{22.2,0.},{22.5,0.},{22.8,0.},{23.1,0.},{23.4,0.},{23.7,0.},{24.,0.},{24.3,0.},{24.6,0.},{24.9,0.},{25.2,0.},{25.5,0.},{25.8,0.},{26.1,0.},{26.4,0.},{26.7,0.},{27.,-0.0001},{27.3,-0.0001},{27.6,-0.0002},{27.9,-0.0005},{28.2,-0.001},{28.5,-0.002},{28.8,-0.0041},{29.1,-0.0084},{29.4,-0.0173},{29.7,-0.0355},{30.,-0.073},{30.3,-0.1019},{30.6,-0.1186},{30.9,-0.1214},{31.2,-0.1098},{31.5,-0.0854},{31.8,-0.051},{32.1,-0.0108},{32.4,0.0308},{32.7,0.0688},{33.,0.0991},{33.3,0.1184},{33.6,0.1249},{33.9,0.1179},{34.2,0.0985},{34.5,0.0689},{34.8,0.0323},{35.1,-0.0076},{35.4,-0.0466},{35.7,-0.0813},{36.,-0.1084},{36.3,-0.1258},{36.6,-0.1322},{36.9,-0.1273},{37.2,-0.112},{37.5,-0.0877},{37.8,-0.0566},{38.1,-0.0213},{38.4,0.0155},{38.7,0.0512},{39.,0.0836},{39.3,0.1105},{39.6,0.1306},{39.9,0.1429},{40.2,0.147},{40.5,0.143},{40.8,0.1314},{41.1,0.1132},{41.4,0.0894},{41.7,0.0614},{42.,0.0306},{42.3,-0.0015},{42.6,-0.0335},{42.9,-0.0642},{43.2,-0.0924},{43.5,-0.1169},{43.8,-0.137},{44.1,-0.152},{44.4,-0.1615},{44.7,-0.1651},{45.,-0.1627},{45.3,-0.1545},{45.6,-0.1408},{45.9,-0.1219},{46.2,-0.0986},{46.5,-0.0717},{46.8,-0.042},{47.1,-0.0106},{47.4,0.0211},{47.7,0.052},{48.,0.0807},{48.3,0.1058},{48.6,0.1262},{48.9,0.1406},{49.2,0.1481},{49.5,0.1481},{49.8,0.1402},{50.1,0.1247},{50.4,0.1019},{50.7,0.0731},{51.,0.0398},{51.3,0.0039},{51.6,-0.0323},{51.9,-0.0662},{52.2,-0.0953},{52.5,-0.1174},{52.8,-0.1303},{53.1,-0.1329},{53.4,-0.1245},{53.7,-0.1056},{54.,-0.0774},{54.3,-0.0422},{54.6,-0.0031},{54.9,0.0362},{55.2,0.0721},{55.5,0.1007},{55.8,0.1192},{56.1,0.1253},{56.4,0.1182},{56.7,0.0984},{57.,0.0679},{57.3,0.0298},{57.6,-0.0116},{57.9,-0.0516},{58.2,-0.0858},{58.5,-0.11},{58.8,-0.1214},{59.1,-0.1186},{59.4,-0.1018},{59.7,-0.073},{60.,-0.0355},{60.3,-0.0173},{60.6,-0.0084},{60.9,-0.0041},{61.2,-0.002},{61.5,-0.001},{61.8,-0.0005},{62.1,-0.0002},{62.4,-0.0001},{62.7,-0.0001},{63.,0.},{63.3,0.},{63.6,0.},{63.9,0.},{64.2,0.},{64.5,0.},{64.8,0.},{65.1,0.},{65.4,0.},{65.7,0.},{66.,0.},{66.3,0.},{66.6,0.},{66.9,0.},{67.2,0.},{67.5,0.},{67.8,0.},{68.1,0.},{68.4,0.},{68.7,0.},{69.,0.},{69.3,0.},{69.6,0.},{69.9,0.},{70.2,0.},{70.5,0.},{70.8,0.},{71.1,0.},{71.4,0.},{71.7,0.},{72.,0.},{72.3,0.},{72.6,0.},{72.9,0.},{73.2,0.},{73.5,0.},{73.8,0.},{74.1,0.},{74.4,0.},{74.7,0.},{75.,0.}};

I then used interpolation in order to be able to work with these (e.g. take the derivative):

wf = Interpolation[wfData];

Now I want to calculate the kinetic energy as a function of x. The equation is

where T is the kinetic energy and psi is the wave function. Let's ignore the constant on the right side. How do I go about finding T while giving conditions like T being continuous? Obviously, I can't find an equation for T since I don't even have an equation for psi, but it should be possible to numerically find an interpolating function for T. How can I do this?

  • $\begingroup$ I have no idea of quantuum mechanics, but if the function $\psi(x)$ is given or approximated with your data, why dont you just compute the function $T(x)$ as $T = c \ \psi''(x)\ / \ (\psi(x))$ with $c$ as the proportionality constant of yours? $\endgroup$ – Mauricio Fernández May 16 '17 at 10:36
  • $\begingroup$ Because psi(x) gets 0 for several x values and thus simply dividing gives infinities. And even when I use Piecewise to eliminate these scenarios, I get ugly spikes in my function because the step size between my data points limits the accuracy. But I now found a way around it by just using MedianFilter to remove the spikes, so I guess my problem is solved now. $\endgroup$ – Keno May 16 '17 at 11:41
  • $\begingroup$ As an aside, I'd be interested to know a bit more context about the physics of this situation. Usually in quantum mechanics the kinetic energy isn't a function of position, but rather of the entire state, and you can't meaningfully talk about "the kinetic energy as a function of x". There are certain approximations (such as the WKB approximation) that relax this, so perhaps that's what you're doing here? $\endgroup$ – Michael Seifert May 16 '17 at 12:51
  • 1
    $\begingroup$ @MichaelSeifert I think you might be confusing the kinetic energy with the total energy here, but I only had quantum mechanics for a few weeks now, so I might be wrong. However, the time independent Schrödinger equation is $\left(E-V(r)\right)\psi(r)=-\frac{\hbar^2}{2m}\nabla^2\psi$ and since $T=E-V$ and $V=V(r)$, $T$ is indeed a function of $r$. Only $E$ is a constant. As to what it means - I have no clue. I guess quantum mechanics is often more about doing than about understanding. $\endgroup$ – Keno May 16 '17 at 13:09
  • $\begingroup$ You're right, of course. My apologies. $\endgroup$ – Michael Seifert May 16 '17 at 13:43

If we simply want an interpolation of the KE operator on the wavefunction, with $\hbar$ and $m$ set to unity, we can evaluate

wf = Interpolation[wfData];
ke[x_] = -(D[wf[λ], λ, λ]/2) /. λ -> x;
Plot[ke[x], {x, 20, 70}]

Then, to get the expectation value of the KE in the given state, we evaluate

region = Interval[{wfData[[1, 1]], wfData[[-1, 1]]}];
xpectKE = NIntegrate[wf[x] ke[x], x ∈ region] /
  NIntegrate[wf[x] wf[x], x ∈ region]

(*  {0.384641}  *)

What about the units? To get units of energy, we must multiply by $\hbar^2/m$ and divide by the square of the units of $x$. Also, the above applies for a real wavefunction. If the wavefunction had complex values, we would modify the formulas accordingly.


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.