Solve for a non-differential equation where the solution is a function

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?

• 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? – Mauricio Fernández May 16 '17 at 10:36
• 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. – Keno May 16 '17 at 11:41
• 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? – Michael Seifert May 16 '17 at 12:51
• @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. – Keno May 16 '17 at 13:09
• You're right, of course. My apologies. – 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.