# More elegant/efficient numerical WKB implementation for 1D Schrödinger Equation

I'm trying to figure out a better way to implement numerical WKB quantization of a spectrum. WKB quantization is the condition that

$$\oint_{E}p \, \mathrm{d}q = \pi \hbar (n+1/2)$$

So the idea to implement it numerically is to do the integral over and over again for different values of $$E$$ and when it turns out to be a half-integer multiple of $$2 \pi \hbar$$ I take the corresponding value of $$E$$ to be in the spectrum.

Here's how it's implemented. The integral is rescaled so that it is always on the interval $$(-1,1)$$ and the energy can be related to the value a used to do this rescaling:

h = 1; m = 2; xmax = 10; nMax = 87;
Vv[x_] := x^2/2 + 0.375 x^4;

action = Table[{a,
Re[NIntegrate[a Sqrt[2 m (Vv[a] - Vv[a x])], {x, -1, 1}]]},
{a, 0.000001, xmax/2, 0.005}];

curve = Interpolation @ action;

allowedEnergies = a /. Table[
FindRoot[ curve[a] == π h (n - 1/2), {a, 1, 2}], {n, 1, nMax}];

ListPlot @ allowedEnergies


Is there a better (or even built-in) way to implement this? I think this way of going about it is pretty hacky and inelegent.

• Note, by the way, that your allowedEnergies are not actually the allowed energies. You've got to "unscale" them by working backwards to get the turning points, and then plug the turning points back into V to get the real energies. That will turn your square-root-ish looking energies into something else. – march Nov 15 '19 at 0:14
• Actually, the way you've parameterized things, a is the turning point, and so you just need to plug it back into the potential to get the actual energy. – march Nov 15 '19 at 0:49
• I actually like your implementation. What do you want to improve? – yarchik Nov 15 '19 at 9:09

I think we can do direct root-finding by defining a numerical function for your integration variable, like so:

ClearAll@f
f[a_?NumericQ] := NIntegrate[a Sqrt[2 m (Vv[a] - Vv[a x])], {x, -1, 1}]


Then we can root-find on this function using FindRoot, i.e.

Table[a /. FindRoot[f[a] == π (n + 1/2), {a, n}], {n, 0, 5}]
(* {0.75565, 1.18227, 1.43771, 1.62971, 1.78691, 1.92165} *)


If you compare these values to what you get, it's actually very similar, but this method is much faster. Note that you might want to tweak the initial guesses to get better convergence, and Mathematica does spit out some errors, so you want to be careful with the method to make sure that you can trust the results.