Shooting method

Quantum Mechanics
 K. Konishi, G. Paffuti

Solve the eigenvalue problem

H ψ = E ψ ;      H = ((p^2)/(2m)) + V[x];

by the shooting method.

Packages

ClearAll["Global`*"];
Off[General::spell, General::spell1,InterpolatingFunction::"dmval"];

General remarks

A typical eigenvalue problem in Quantum Mechanics is solving an equation of the form

y''[x] + F[x] y[x] = λ y[x]

with homogeneous boundary conditions:

y[xL]=0;           y[xR]=0;

xL, xR are the boundary points. Boundary conditions can be satisfied only for special values of the parameter λ, called eigenvalues of the equation.
A simple method to find eigenvalues is the shooting method, in which the boundary values problem is traded for a initial value problem, easily solved with standard algorithms. The basic idea is shown in the figure below:

Plot[{-0.4x^2 + x,-0.4x^2+0.9x +0.01x^3,-0.4x^2+1.1x},{x,0,2},
PlotRange-> {{0,3.5},{-0.2,1}},
Axes-> False,Ticks->None ,
Epilog->{ Line[{{0,-0.2},{2,-0.2},{2,1},{0,1},{0,-0.2}}] ,
Text["desired",{3,0.5}],
Text["boundary",{3,0.4}],
Text["condition",{3,0.3}],
Arrow[{{2.6,0.4},{2.1,0.4}}]
}]

Trial integrations that satisfy boundary conditions at one endpoint are performed. The discrepancies at the other endpoint are used to adjust the starting conditions or other parameters until both conditions are ultimately satisfied.
Equations () and () are homogeneous, i.e. if y[x] is a solution the same is true for  c y[x]. We can choose the overall scale putting, for example, y'[xL]=1. In this case the value y[xR] is a function of λ,  δ[λ] = y[xR], and the condition to be satisfied is δ[λ] = 0.

The function Sin[x] satisfies the equation y''[x] =  λ y[x]  with y[0] = y[π] = 0 for λ = -1. Let us solve the equation for an arbitrary value of λ with initial conditions y[0] = 0, y'[0] = 1. The function δ[λ] is defined as

δ[λ_?NumberQ]:= ψ[π]/.(NDSolve[{ψ''[x]== λ ψ[x],ψ[0]== 0,ψ'[0]==1},ψ,{x,0,π}][[1]])

and the eigenvalue can be found using FindRoot

λ/.FindRoot[δ[λ]==0,{λ,0.}]

We can plot δ[λ] and look for higher eigenvalues

Plot[δ[λ],{λ,-10,0}]

The eigenvalues are -1,-4,... - (n^2) as is easily seen.They can be obtained choosing a different starting point for the function FindRoot:

λ/.FindRoot[δ[λ]==0,{λ,-5.}]

The syntax δ[λ_?NumberQ] is due to the non-standard evaluation of FindRoot (see Mathematica Help). With this definition the function is evaluated only when his argument, λ, is numeric.
The output of NDSolve is in general a vector whose first argument is the solution given as an interpolating function. The syntax NDSolve[....][[1]]  select this first element.

Quantum Mechanics

In usual quantum mechanics problems the above procedure needs sometimes a little modification: the end points are at infinity and/or the equation the coefficients in the equation are singular at the endpoints. Let us consider the first case and let us assume  xL = - ∞, xR = + ∞. The simplest way  is to put the system in a box with finite boundaries xL, xR, ex. xl = - L, xR = +L , and then to evaluate the solution for large enough L. 
In simple cases the potential is growing as |x| gets large and one can easily choose a value of L deep inside the classically forbidden region, where the wave function goes rapidly to 0. Let us note some points:

Example

Let us apply the procedure sketched above to a  potential defined in the range -∞, + ∞. We use harmonic potential to check the results but the procedure is completely general.

We choose the potential

V[x_]:= x^2/2;

and write condition, allowing a box of size L:

δ[λ_?NumberQ,L_?NumberQ]:= ψ[L]/.(NDSolve[{-1/2ψ''[x]+V[x] ψ[x]== λ ψ[x],ψ[-L]== 0,ψ'[-L]==1},ψ,{x,-L,L}][[1]])

It is easily shown that the asymptotic behavior of the eigenfunction is exp[-x^3/3] so, for L = 4 we have already a small factor Exp[-64/3] = 5 10^-10.

λ/.FindRoot[δ[λ,4.]==0,{λ,0.}]

We get the correct result. For large value of L we have a numerical instability, as shown by the output of FindRoot

λ/.FindRoot[δ[λ,10.]==0,{λ,0.}]

Let us plot the result as a function of L (we put off the Warning of FindRoot)

Off[FindRoot::"lstol"]

values=λ/.FindRoot[δ[λ,#]==0,{λ,0.}]&/@Range[2,10,0.2];

ListPlot[Transpose[{Range[2,10,0.2`],values}],PlotRange->{{1,10},{0.49`,0.51`}},AxesOrigin->{1,0.49`}]

For L > 4 the result is stable.  We give some details on Mathematica commands:

The expression for values given above uses two typical Mathematica constructions: abstract functions and Map. 
body& is the abstract form of a function, #, #1, #2 etc. are the parameters. Eg.   Sin[#]& 2.5 is equivalent to Sin[2.5],  (1 + #^2)& x is equivalent to  1+x^2  etc. In the expression values we want to compute the root for a generic L, this is achieved leaving in abstract form the second argument of δ[ℰ,#].
The second typical construction is the function Map, or /@.  f /@ List applies f  to each element of List. Eg. Sin /@ {1.1, 1.2} gives {Sin[1.1], Sin[1.2]}.
The reader can consult the online Help for more detailed explanations.

Excited states

We can start with a plot of the eigenvalue at a fixed value of L, say 8, changing the starting value for FindRoot:

values=λ/.FindRoot[δ[λ,8]==0,{λ,#}]&/@Range[0,3,0.2]

The reader recognizes the Hamiltonian eigenvalues for excited states. He/she can do a stability study.

Animation

We can follow the step by step the results of  FindRoot using the Option EvaluationMonitor (see Help in Mathematica ). The option is given in the form :> instead of the usual -> to prevent the expression  ( here ℰ ) to be immediately evaluated.

rList = {};rr=FindRoot[δ[λ,4.]==0,{λ,0.},EvaluationMonitor:>AppendTo[rList,λ]]

The partial results are:

rList

The number of steps is:

Length[rList]

We can create an animation showing the approach to the solution using Animate. The reader have to consult Mathematica Help for this function.

Animate[Module[{L=4,V},
V[x_]:= x^2/2;
solution = NDSolve[{-1/2ψ''[x]+V[x] ψ[x]== λ  ψ[x],ψ[-L]== 0,ψ'[-L]==1},ψ,{x,-L,L}];
Plot[(ψ[x]/.solution)/(ψ[0]/.solution),{x,-L,L},PlotRange-> {{-L,L},{-1,2}}]],
{λ, rList},AnimationRunning->False]

The simplest syntax for Animation is
Animate[expr, {u, (Subscript[u, min]), (Subscript[u, max])}]for a continuous range of the parameter u
Animate[expr, {u, { (Subscript[u, 1]), (Subscript[u, 2]) , ... (Subscript[u, k])}}] for u varying in a list.

Clear[L];

Example 2

A more efficient, and usually more stable, procedure to compute the eigenvalues is to solve the equation starting from both endpoints and imposing a matching condition at some  intermediate value. Using the homogeneity  of the problem we can impose the continuity of the logarithmic derivative ψ'/ψ.

Let us take x=0 as the matching point. We define the discrepancy function as

δ[λ_?NumberQ,L_?NumberQ]:=  (solLeft =NDSolve[{-1/2ψ''[x]+V[x] ψ[x]== λ ψ[x],ψ[-L]== 0,ψ'[-L]==1},ψ,{x,-L,0}][[1]];
solRight =NDSolve[{-1/2ψ''[x]+V[x] ψ[x]== λ ψ[x],ψ[L]== 0,ψ'[L]==-1},ψ,{x,L,0}][[1]];
((ψ'[0]/ψ[0])/.solLeft) - ((ψ'[0]/ψ[0])/.solRight))

Operating as in the previous case:

rList = {};rr=FindRoot[δ[λ,4.]==0,{λ,0.},EvaluationMonitor:>AppendTo[rList,λ]];

{Length[rList],rList}

The solution is found in less steps, but we have solved 2 equations per step.

Example 3

The procedure can be refined in particular cases. For example for even potentials we can limit ourselves to the half line [0,∞] taking as left boundary condition ψ[0] = 1, ψ'[0] = 0 for even states and ψ[0] = 0, ψ'[0] = 1 for odd states., i.e. defining:

δEven[λ_?NumberQ,L_?NumberQ]:= ψ[L]/.(NDSolve[{-1/2ψ''[x]+V[x] ψ[x]== λ ψ[x],ψ[0]== 1,ψ'[0]==0},ψ,{x,0,L}][[1]]);

δOdd[λ_?NumberQ,L_?NumberQ]:= ψ[L]/.(NDSolve[{-1/2ψ''[x]+V[x] ψ[x]== λ ψ[x],ψ[0]== 0,ψ'[0]==1},ψ,{x,0,L}][[1]])

For the harmonic oscillator we have now:

rList = {};rr=FindRoot[δEven[λ,4.]==0,{λ,0.},EvaluationMonitor:>AppendTo[rList,λ]];
{Length[rList],rList}

On[FindRoot::"lstol"];

Exercises

Compute the lowest eigenvalue for the Schrödinger equation in the potential V[x] = x^2/2 + g x^4/2 (anharmonic oscillator) with g = 10.. The approximate result is

ℰ = 1.224587036059193459134396953093865213110138999308673549845104365885027406655096801688754290649704430

Compute the lowest eigenvalue for the Schrödinger equation in the potential V[x] = - U for |x| < a. Choose values for a and U.

Compute the lowest eigenvalue for the Schrödinger equation in the potential V[x] = - U/Cosh[x]^2. The exact value is [1]

eneTh = -((((Sqrt[1+ 8 V0]) - 1))^2) /8

Boundary value problems in an infinite interval can be reduced to  finite interval problems by a change of coordinates. Some useful transformations are

x = ArcTanh[ξ/L];  x = ξ/(Sqrt[(L^2)- (ξ^2)]);  ...

The reader can explore the shooting method using this technique and study the dependence on the scale parameter L.

References

1) L.Landau, E. Lifchitz, Quantum Mechanics: Non Relativistic Theory, Butterworth-Heinemann; 3 edition (1981).