Shooting Methods

Shooting treats the energy in a stationary Schrödinger equation as an unknown parameter in an initial-value problem. One integrates from one boundary and adjusts the energy until the other boundary condition is satisfied. The method is economical for a few one-dimensional states, but forbidden-region growth can make one-sided shooting numerically unstable.

Required background. Boundary conditions in computation specifies the endpoint problem; numerical error and stability separates a root residual from solution accuracy.

The energy-dependent initial-value problem

For $-\alpha y''+Vy=Ey$ , with $\alpha=\hbar^2/(2m)>0$ , integrate

y'=z,\qquad z'=q(x,E)y, \qquad q(x,E)=\frac{V(x)-E}{\alpha}.

On a regular finite interval $[a,b]$ with Dirichlet endpoints, choose $y(a)=0$ , $z(a)=1$ . The derivative fixes an arbitrary nonzero scale; the final eigenfunction is normalized later. Define the shooting residual $F(E)=y(b;E)$ . Eigenvalues are zeros of $F$ .

Bracket a desired root and use a safeguarded scalar solver. The inner ODE accuracy must be tight enough that the sign and magnitude of $F$ are resolved. Reducing the root tolerance below the integration error creates a misleading precision claim.

An exact shooting residual

For $V=0$ , $a=0$ , and $b=L$ , positive energies give

y(x;E)=\frac{\sin(kx)}k, \qquad k=\sqrt{E/\alpha}, \qquad F(E)=\frac{\sin(kL)}k.

The nonzero roots are $kL=n\pi$ with $n\ge1$ . At $E=0$ , the continuous limit is $F(0)=L$ , not zero. Replacing $F$ by $\sin(kL)$ without treating $k=0$ separately introduces a spurious “ground state.” This example tests the initialization as well as the energy search.

At a regular separated Sturm–Liouville problem’s $n$ th eigenvalue, counting from $n=0$ , the eigenfunction has $n$ interior nodes. Node counting helps identify which root has been found and exposes skipped levels. Count resolved sign changes away from endpoint roundoff, not fluctuations in an underresolved tail.

Sensitivity equations and safeguarded Newton steps

When the initial data are independent of energy, differentiate the ODE with respect to $E$ . With $s=\partial_Ey$ ,

s''=q(x,E)s-\frac{y}{\alpha}, \qquad s(a)=s'(a)=0.

Then $F'(E)=s(b)$ can be computed alongside $F$ , yielding a Newton proposal $E-F/F'$ . Retain a bracket and reject proposals that leave it or rely on an unresolved derivative. If the boundary initial data themselves depend on energy, their derivatives must also be included; the zero sensitivity initial conditions no longer apply.

A small $F(E)$ depends on the arbitrary solution scale. Use a scale-aware residual or combine it with bracket width and refinement. Periodically rescaling the integrated solution prevents overflow, but the sensitivity bookkeeping must respect that rescaling if Newton derivatives are used.

What the residual says about energy error

For the real regular Dirichlet problem above, an exact eigenvalue $E_n$ has a nonzero residual derivative. To see why, form $W_s=ys'-y's$ from the solution and its energy sensitivity. Their equations give

W_s'=-\frac{y^2}{\alpha}, \qquad W_s(a)=0.

At an eigenvalue, $y(b)=0$ , so integration yields

F'(E_n)=s(b;E_n) =\frac{\displaystyle\int_a^b y(x;E_n)^2\,dx} {\alpha\,y'(b;E_n)}.

The endpoint derivative cannot vanish together with $y(b)$ for a nontrivial regular solution. Thus the root is simple. Close to it, $F/F'$ estimates the energy error in this fixed boundary-value problem; it is unchanged by an energy-independent rescaling of the shooting solution. An uncertainty $\delta F$ in the integrated residual correspondingly produces an energy uncertainty of order $|\delta F/F'|$ . Neither quantity measures the error caused by replacing an infinite domain with a finite interval.

Forbidden-region growth and matching

In a region with approximately constant $V-E>0$ , the solutions are $e^{\pm\kappa x}$ , $\kappa=\sqrt{(V-E)/\alpha}$ . A tiny contamination of the growing branch eventually overwhelms the desired decaying branch. Smaller time or space steps reduce local error but do not remove the exponential conditioning of a long outward integration.

For a bound state with two tails, integrate stable solutions inward from each side and match near the classically allowed region. Matching methods develops Wronskian and logarithmic-derivative residuals. A logarithmic derivative avoids amplitude overflow but has poles at wavefunction nodes; it is not a universally smooth scalar residual.

For a symmetric potential, parity reduces the problem to a half interval. Even states use $y'(0)=0$ with nonzero $y(0)$ ; odd states use $y(0)=0$ with nonzero derivative. Search the parity sectors separately, but do not impose this reduction on an asymmetric potential.

An oscillator with two independent numerical limits

In oscillator units $\hbar=m=\omega=1$ , the even ground state satisfies

y''=(x^2-2E)y, \qquad y(0)=1, \qquad y'(0)=0.

The whole-line energy is $E_0=1/2$ . Replacing its decaying tail by $y(L)=0$ instead defines an oscillator confined by walls at $\pm L$ . Wagner’s shooting examples discuss this parity reduction and the difficulty of the growing forbidden-region solution. The following calculation keeps the finite boundary and integration step separate.

Integrate the two first-order equations with classical fourth-order Runge–Kutta, using constant step $h=L/N$ . Bisect the ground-state bracket $[0.45,0.65]$ for 45 iterations, evaluating the integrated endpoint $F_h(E)=y_h(L;E)$ at each trial energy. Repeating for three boundaries gives:

Step $h$	Energy at $L=2$	Energy at $L=3$	Energy at $L=4$
$0.100$	$0.5374622816$	$0.5003916581$	$0.5000010083$
$0.050$	$0.5374612767$	$0.5003911190$	$0.5000005234$
$0.025$	$0.5374612135$	$0.5003910852$	$0.5000004929$

For each fixed $L$ , successive energy differences decrease by approximately sixteen, as expected for fourth-order integration in this resolved range. Yet the three columns approach different limits. The finite-boundary energies are approximately $0.5374612093$ , $0.5003910829$ , and $0.5000004909$ . In particular, $L=2$ retains an error of about $0.03746$ relative to the whole-line energy even though all the unrounded numerical roots have $|F_h|<10^{-11}$ with the stated initial scale. The small residual means that the numerical endpoint equation was solved accurately; rounding the displayed energies can increase that residual.

One independent check of these finite-boundary limits uses the exact even solution

y(x;E)=e^{-x^2/2}\, {}_1F_1\!\left(\frac{1-2E}{4};\frac12;x^2\right).

Here ${}_1F_1$ is the solution of Kummer’s equation analytic at zero, normalized to one there. With $t=x^2$ and $y=e^{-t/2}w(t)$ , the oscillator equation becomes $t w''+(1/2-t)w'-(1-2E)w/4=0$ , fixing the parameters above; the NIST definition also supplies a convergent series for evaluating it. Its zeros at $x=L$ give the quoted limits. An independently refined finite-difference Hamiltonian on the same interval approaches the same limits. The Harmonic Oscillator Spectrum lab provides that matrix calculation and its separate mesh and boundary studies.

Enlarging $L$ removes the artificial-wall error but worsens the conditioning of outward shooting. For longer forbidden regions, use inward propagation or two-sided matching and repeat both refinements. A smaller root tolerance cannot replace either one.

For a regular three-dimensional radial problem, the reduced function behaves as $u(r)\sim r^{\ell+1}$ near the origin. Begin at a small positive $r_0$ using a Frobenius expansion and its derivative, rather than evaluating a singular coefficient at zero. Coulomb and other singular potentials require the appropriate next terms and domain assumptions.

Vary $r_0$ , outer boundary, ODE tolerance, and root tolerance independently. A stable energy under one of these changes does not establish convergence in the others. Normalize with the correct radial measure: $u=rR$ has norm $\int|u|^2dr$ .

For each reported state, retain energy, node count, boundary residual, normalization, and an independent comparison with a matrix discretization or analytic limit. Shooting excels at low-dimensional boundary-value problems; it is not the natural method for a generic high-dimensional many-body Hamiltonian.

Exercises

Why does setting both $y(a)=0$ and $y'(a)=0$ make a linear shooting calculation useless?

Solution

Uniqueness of the regular linear initial-value problem gives the identically zero solution for every energy. Its right boundary residual vanishes trivially and contains no spectral information. One unconstrained amplitude or derivative must be nonzero.

For the free Dirichlet problem, use the sensitivity identity to find $F'(E_n)$ at $E_n=\alpha(n\pi/L)^2$ . How would multiplying the initial derivative by $10^{-8}$ affect $F$ and the local estimate $F/F'$ ?