Phase Method basic eigenvalue search

Grant Bunker May 2024 
example code for book “The Phase Method”

This code, while still very basic, is a little more complex than the python example since it does extended precision, parallel, 
more graphics, and tests, and much of it is comments like these. 
This example code does not do wave functions, just energies.  A different notebook does everything (see web site).

initial setup

ClearAll["Global`*"]

starttime=AbsoluteTime[ ];

Choose desired digits of working precision $MachinePrecision ≤ workingprecision ≤ 40 or so;
$MachinePrecision is one option (fastest). More digits makes the computation slower but more accurate

workingprecision=31;

estimated ultimate accuracy, based on experience;  your mileage may vary

ngooddigits=workingprecision/2;

set nominal value and slope initial conditions at xmin

Phi0=1/10;Phi0p=1/10;

define the potential here however you want to.  Here we just define a finite square well using an If statement

a=3;V0=-10;
U=If[Abs[x/a]<1,V0,0];

set etest to the upper end of your energy search range

etest=-1/10;

set the range -- should be well outside turning points at etest.  Nearly-unbound states will require wider x-range for accuracy

{xmin,xmax}={-15,15};

verify that xmin, xmax are OK by checking the phaseplot at etest

{umin,uminxsoln}=NMinimize[U,x,Method->"RandomSearch"]//Chop;
uminx=x/.uminxsoln;
Print["here is the potential plot and etest - the turning points are where they intersect:"]
Plot[{U,etest},{x,xmin,xmax},PlotRange->All,Frame->True]

set energy search range for eigenvalue search.  
Set this to whatever makes sense for your problem.

{eminsearch,emaxsearch}={umin,etest};
searchinterval=SetPrecision[{eminsearch,emaxsearch},workingprecision];

the "nbs" (number of bound states aka "Ω") function estimates the number of bound states below specified energy;it also calculates the phase plot in the process. Clear out previously memo-ized values

ClearAll[nbs];

Sampling grid for counting transitions, just has to be small enough. This could be estimated from max momentum at potential minimum, but generally there is no need.

dx=.001;

Make the variables extended-precision so as to not degrade precision in calculation
A nicer way to count transitions is to use "WhenEvents" in NDSolve, but the stone-axe approach shown here is portable to other implementations that may not offer that capability.

{Phi0,Phi0p,xmin,xmax,e}=SetPrecision[{Phi0,Phi0p,xmin,xmax,e},workingprecision];
nbs[e_]:=nbs[e]=(soln=Quiet@NDSolve[{-1/2Phi''[x]+(U-e)Phi[x]==0,Phi'[xmin]==Phi0p,Phi[xmin]==Phi0},Phi[x],{x,xmin,xmax},MaxSteps->10^7,Method->"LSODA",WorkingPrecision->workingprecision][[1]];
phaseplot[x_]=Tanh[Phi[x]/.soln];
transitions=DeleteCases[0]@Differences@Quiet@Table[Sign[Tanh[Phi[x]/.soln]],{x,xmin,xmax,dx}]/2;
Return[{e,nboundstates=Length[transitions]}]);

Examine the phaseplot before running the rest of calculation to verify the x-range is suitable for the chosen maximum energy. 
Transitions should be sharp, not rounded and there should regions of “saturation” at the near xmin and xmax.

nbs[etest];
dophaseplot:=Plot[phaseplot[x],{x,xmin,xmax},Frame->True,PlotLabel->"# bound states="<>ToString[nboundstates],PlotRange->All,PlotPoints->2000,ImageSize->8*72];

sanity check

Let's do a test run to make sure the x-range is set sensibly. If transitions are rounded instead of sharp, extend xmin and xmax further past the turning points; 
If overflow is encountered, reduce the range, or better, choose xmin and xmax by scaling the turning point values as described in the book, “The Phase Method”.

Print[Framed["Here is the phase plot at the upper energy limit etest="<>ToString[TraditionalForm@etest]]];
dophaseplot

start search

define "sift" and "refine" functions — accuracy seems to be limited to about half the working precision

maxiterations=Min[100,Ceiling[(workingprecision/2*Log[2,10]+5)]];

"sift" maps out where eigenvalues are located (serial process)

sift=bsearch/@(Partition[Flatten[#],2])&;

"refine" gets more accurate estimates (parallel process)

refine=FixedPoint[sift,#,maxiterations]&;

define tree - like binary search function

ns[e_]:=nbs[e][[2]];intervals={};
bsearch[intrvl_]:=Module[{vlow,vhigh,vmid,mid,newintervals},({elow,ehigh}=intrvl;emid=(elow+ehigh)/2;
{vlow,vmid,vhigh}={ns[elow],ns[emid],ns[ehigh]};newintervals={};
If[vlow<vmid,AppendTo[newintervals,{elow,emid}]];
If[vmid<vhigh,AppendTo[newintervals,{emid,ehigh}]];
Return[newintervals]);]

sifting (serial)

First, sift for eigenvalues (serially), enough to resolve the eigenvalues . 
This is mostly for visualization purposes. Don't worry if sifting is incomplete -- the refine process (in parallel) will still take care of the rest

Print[Framed@"for some potentials, sifting can take a few minutes, please be patient..."];

Rough estimate of number of bisections needed to separate the eigenvalue intervals

snit=Ceiling[Log[2,nboundstates]]+5;
{sifttime,intervals}=AbsoluteTiming[Flatten[Nest[sift,searchinterval,snit],1]];
Print["sifting complete: sifting time=",sifttime," sec,  # of intervals=",Length[intervals]];
ListPlot[Interval/@intervals,PlotRange->All,Frame->True]

This optional plot shows the interval search, using values already "memoized"

Print["the following plot shows the evolution of the binary search process"]
siftplot=ListPlot[Interval/@Partition[Flatten[NestList[sift,searchinterval,10]],2],PlotRange->All,ImageSize->Large,Frame->True,FrameLabel->{{Subscript[ℰ,n],None},{None,"U="<>ToString@TraditionalForm[U]}}]

refine (parallel)

Then "refine" the found eigenvalues

Print[Framed@"for some potentials, refining also can take a few minutes\n ... maybe enough time to get a cup of coffee or tea"];
{refinetime,intervals}=AbsoluteTiming[Flatten[ParallelMap[refine,intervals]//Quiet,2]];
eigs=Mean/@intervals;neigs=Length[eigs];roundedeigs=N[eigs,ngooddigits];
Print["sifting time=",sifttime,"; refining time=",refinetime,"; total=",sifttime+refinetime," (seconds)"];
Print["neigs=",neigs,";  ",(sifttime+refinetime)/neigs," seconds/eigenvalue"];

Here are the eigenvalue estimates, which are typically accurate to about workingprecision/2 digits;   your mileage may vary
If you want to be absolutely sure it’s right, drop a least-significant digit or two.

TableForm[Transpose[{Range[neigs],roundedeigs}]]

As a verification test, the following code checks accuracy of the calculated eigenvalues 
for a symmetric finite square well with the potential parameters above.
It can be ignored or deleted for other potentials.

opts={Frame->True,FrameLabel->{"n","ℰ"},GridLines->Full,ImageSize->6*72};

setpars = {m -> 1, hbar -> 1};
setvars = {k -> Sqrt[2 m (e - V0)]/hbar, Gamma -> Sqrt[2 m (-e)]/hbar};
zeta = hbar^2/(2 m a^2);
nbound[V0_] := 1 + Floor[(2/Pi) Sqrt[-V0/zeta]] /. setpars;
Print["nbound states=", nbound[V0]];

applying the following function to an eigenvalue should give zero for even parity states

eventest[e_]=Cot[k a]-k/Gamma/.setvars/.setpars;

applying the following function to an eigenvalue should give zero for odd parity states

oddtest[e_]=Tan[k a]+k/Gamma/.setvars/.setpars;

solve the exact equations for eigenvalue by numerical rootfinding using PM eigenvalues with random perturbations as guesses

rand:=.01Random[]

eroots=e/.FindRoot[eventest[e],{e,#},WorkingPrecision->workingprecision]&/@(eigs[[1;;;;2]]+rand);
oroots=e/.FindRoot[oddtest[e],{e,#},WorkingPrecision->workingprecision]&/@(eigs[[2;;;;2]]+rand);
exact=Riffle[eroots,oroots];

Precision[exact]

these are the maximum and minimum relative errors of the PM-computed eigenvalues relative to the “exact” rootfinding solutions to the exact equations

MinMax[Abs[eigs/exact-1 ]]

ListLogPlot[Abs[eigs/exact-1 ],PlotRange->All,opts]

for square well the levels are roughly quadratic in n

GraphicsRow[{ListPlot[Abs[eigs-V0],opts],ListLogLogPlot[Abs[eigs-V0],opts]},ImageSize->12*72]

elapsedtime=(AbsoluteTime[]-starttime)/60. min;
Print["total elapsed time= ",elapsedtime]