PHYS 410/555 Computational Physics Solving Two-Point Boundary Value Problems Using “Shooting”

A Toy Model for the Deuteron

Recall that, by definition, two-point boundary value problems (BVPs), are ODE’s for which boundary conditions are supplied at two distinct points—typically the end points of the solution domain—rather than at some single point, as in the case of initial value problems (IVPs). In addition, we observe that two-point BVPs are often (but not always) eigenvalue problems, that is, BVPs are often characterized by one or more parameters such that only for specific parameter values (eigenvalues of the problem) will solutions satisfying the boundary conditions exist.

In such a case, the solution of the BVP also becomes a problem in search—in general we will not be able to construct an algorithm which in “one go” results in the eigenvalue and the associated eigenfunction. Rather, we will have to provide some initial estimate (guess) of the eigenvalue, and then successively refine it according to some criterion, until we have computed the eigensolution to some acceptable accuracy. Also, solutions of eigenvalue problems are generally not unique; typically an eigenvalue problem admits a countable infinity of eigenfunctions, each with an associated eigenvalue (the eigenvalues are often, but not always, distinct).

The technique we will briefly consider here is known as shooting. Shooting is based on the observation that any two-point BVP can also be solved as an IVP. Consider, for example, a BVP of the form:

u′′(x) = f(u, u′, x)  0 ≤ x ≤ 1  (1)

subject to the boundary conditions

u(0) = u_{0}  (2)
u(1) = u_{1}  (3)

where u_{0} and u_{1} are specified constants. Assume that we have somehow determined a function U(x) which satisfies (1), (2) and (3). We then note that if we now consider the initial value problem defined by (1) subject to the initial conditions:

u(0) = U(0) = u_{0}  (4)
u′(0) = U′(0)  (5)

then its solution must also be U(x). That is, we can solve the BVP as an IVP, by “guessing” what value of U′ we must specify at x = 0 so that when we are done integrating from x = 0 to x = 1, we have u(1) = u_{1}.

The term “shooting” comes from analogy with the problem of setting the elevation, θ, of an artillery gun (i.e. θ is the angle the gun’s barrel makes with the horizontal), so that the shell hits a target at some given range. Assuming that θ < π/4, if the elevation if too low/high, the shell will fall short/long of the target respectively. The gunner can use information about where his/her current shot lands to adjust θ so that the next shot comes closer to the target.

Similarly, when integrating a BVP such as (1) via shooting, we will typically find that if we specify u′(0) > U′(0) then we will have u(1) > u_{1}, while if we choose u′(0) < U′(0), then we will find u(1) < u_{1} (equally as likely is that u′(0) > U′(0) ⟶ u(1) < u_{1} and u′(0) < U′(0) ⟶ u(1) > u_{1}). As long as we can find an initial pair of “bracketing” values, [u′_{−}(0), u′_{+}(0)], such that separate integrations of (1) with these two values leads to values of u(1) which similarly bracket the desired boundary value, u_{1}, then we can narrow the bracket using, for example, the technique of bisection search described below, to determine U′(0), and hence the solution of the BVP, to whatever precision is desired.

We note that it is not always an initial value per se which we tune in a shooting problem. For some second order BVPs (such as the one considered below), we can deduce a second initial value at one of the boundary points (in addition to the one given by the boundary conditions) from mathematical or physical considerations, but, as discussed above, there is a parameter, λ, in the specification of the BVP which must have certain values in order that the boundary condition at the other boundary is satisfied. The idea is still the same; we look for an initial bracket [λ_{−}, λ_{+}] such that separate integration with the two parameter values gives end-point solution values which are too small and too large respectively. We can then refine the bracket until we determine the eigenvalue and eigenfunction to the required precision.

We now illustrate this technique by using it to solve the toy model for a deuteron which is discussed in Chapter 9 of Mathematical Methods for Physicists by Arfken. A deuteron, as most of you probably know, is a bound state of a proton and a neutron (the nucleus of deuterium, an isotope of hydrogen). The model is highly idealized; we assume that the deuteron wave function, ψ(r⃗), is a spherically symmetric solution of the time-independent Schrödinger equation, with the proton-neutron interaction described by a square wave potential.

Thus we wish to solve

−\frac{h̄^{2}}{2M} ∇^{2}ψ(r⃗) + V (r⃗)ψ(r⃗) = Eψ(r⃗)  (6)

where M is the deuteron “mass”, and E is the energy eigenvalue which will shortly become the “shooting parameter” in our BVP solution of the model. From the assumption of spherical symmetry, we have

ψ(r⃗) → ψ(r)  (7)

and

∇^{2}ψ(r) = \frac{1}{r^{2}} \frac{d}{dr} (r^{2} \frac{dψ}{dr})  (8)

Further, defining

u(r) ≡ rψ(r)  (9)

we have (as you should verify)

∇^{2}ψ(r) → \frac{1}{r} \frac{d^{2}u(r)}{dr^{2}}  (10)

Thus, (6) can be rewritten:

\frac{d^{2}u(r)}{dr^{2}} + \frac{2M}{h̄^{2}} (E − V (r)) u(r) = 0  (11)

As noted above, we will model the proton-neutron interaction as a finite-range square potential. Thus, we take

V (r) = V_{0}  0 ≤ r ≤ a  (12)
V (r) = 0  r > a  (13)

where V_{0} is a negative constant, so |V_{0}| is the depth of the potential well, while a is its width. As is the case for any solution of the Schrödinger equation, we must demand that our solution of (6) be normalizable, i.e. that

∫ ψψ^{∗}dV = 1  (14)

so that there is unit probability that our deuteron is found somewhere in the universe. In the current spherically symmetric case this means that we must have

4π ∫_{0}^{∞} r^{2}ψ(r)^{2}dr = 4π ∫_{0}^{∞} u(r)^{2}dr = 1  (15)

Clearly, a necessary condition for normalizability is that

lim_{r→∞} u(r) = 0  (16)

and this, in fact, is one of the boundary conditions for our ODE. Furthermore, we assert that for fixed values of the parameters of the model (a, M and V_{0}), a normalizable solution will only exist for certain discrete values of E—the eigenvalues of our Schrödinger equation. Before we consider the solution of (11) using shooting, we rewrite the equation in an equivalent form in which the minimum number of free parameters (which, if desired, can be made explicitly dimensionless) becomes evident. To this end, we define a rescaled radial coordinate, x

x ≡ \sqrt{2M}r  (17)

so that, among other things, we have

\frac{d^{2}u}{dr^{2}} → 2M \frac{d^{2}u}{dx^{2}}  (18)

Further, we choose units such that h̄ = 1 and V_{0} = −1 (you should establish that this is always possible if it isn’t immediately obvious to you).

With these choices, we are left with one free parameter—the width, a, of the potential well. Given that we have adopted the rescaled radial coordinate x, it is more convenient to use x_{0}, defined by

x_{0} ≡ \sqrt{2M}a  (19)

as the free parameter. Thus, our Schrödinger equation (11) becomes

\frac{d^{2}u(x)}{dx^{2}} + (E − V (x)) u(x) = 0  (20)

where

V (x) = −1  0 ≤ x ≤ x_{0}  (21)
V (x) = 0  x > x_{0}  (22)

and where we again note that E = E(x_{0}) is an eigenvalue of (20); i.e., for a given value of x_{0}, only for discrete values of E will we have a normalizable wave function.

The boundary conditions for (20) are derived from the demands that 1. ψ(r) be regular (analytic) at r = 0.

2. lim_{r→∞} u(r) = 0. The regularity condition at r = 0 means that ψ(r) admits an expansion

lim_{r→0} ψ(r) = ψ_{0} + r^{2}ψ_{2} + O(r^{4})  (23)

where the ψ_{i} are constants (the power series expansion cannot have terms which are odd in r, since ψ(r) would not have a well-defined derivative at r = 0 in that case). From this, it follows that u(r) has an expansion

lim_{r→0} u(r) = rψ(r) = rψ_{0} + r^{3}ψ_{2} + O(r^{5})  (24)

Thus, we have

u(0) = 0  (25)
\frac{du}{dr}(0) = ψ_{0}  (26)

We now observe that (6) (like all Schrödinger equations) is linear; given any solution, ψ(r) we have that cψ(r), where c is an arbitrary positive constant, is also a solution. The particular solution we seek is fixed by the normalization condition:

∫ ψψ^{∗}dV = 1  (27)

Operationally, this means that we can choose ψ(0) arbitrarily (say ψ(0) = 1 for convenience), then, for specified x_{0}, vary E until we find a solution which satisfies

lim_{r→∞} u(r) = lim_{r→∞} rψ(r) = 0  (28)

(In other words, the eigenvalue is independent of the normalization of the eigenfunction).

In preparation for a solution of our problem using LSODA we rewrite (20) in canonical first order form by introducing

w(x) ≡ \frac{du(x)}{dx}  (29)

We then have

\frac{du(x)}{dx} = w(x)  (30)
\frac{dw(x)}{dx} = (V (x) − E)u(x)  (31)

subject to

u(0) = 0  (32)
w(0) = 1  (33)

and with E to be determined so that

lim_{x→∞} u(x) = 0  (34)

We will now assume that for any given value of x_{0}, we are able to determine values E_{−} and E_{+} (perhaps by trial-and-error) such that

lim_{x→∞} u(x; E_{−}) = −∞  (35)
lim_{x→∞} u(x; E_{+}) = +∞  (36)

where the notation u(x; E_{i}) means the trial solution u(x) computed using eigenvalue estimate E_{i}. We further assume that properties (35) and (36) hold for any values E_{HI} and E_{LO} that bracket the desired eigenvalue, namely:

lim_{x→∞} u(x; E_{LO}) = −∞  (37)
lim_{x→∞} u(x; E_{HI}) = +∞  (38)

with either

E_{LO} < E < E_{HI}  (39)

or

E_{HI} < E < E_{LO}  (40)

Given the initial bracket [E_{−}, E_{+}], then, we can compute the desired eigenvalue, E, accurate to some desired tolerance δ, using a bisection search (also known as a binary search). Here is a typical implementation of a bisection search written in pseudo-code:

E_{LO} := E_{−} E_{HI} := E_{+} while |E_{HI} − E_{LO}| > δ do

E_{MID} := (E_{HI} + E_{LO})/2

if u(x; E_{MID}) → +∞ as x → ∞ then

E_{HI} := E_{MID}

   else

E_{LO} := E_{MID}

   end if
end while

E := E_{MID} We note that the convergence rate of the bisection search is completely

pre-determined; the size of the bracketing interval after n bisections is

\frac{1}{2^{n}} (E_{+} − E_{−})  (41)

The solution of equations (30)-(33) using LSODA is implemented as

the program deut.f; the bisection search for the eigenvalues E(x_{0})

is implemented via the shell-script Shoot-deut, which itself uses the

following collection of scripts which can be used to implement shell-

level bisection searches:

  bsnew <lo> <hi>     # Initializes a new search
  bscurr              # Returns the current (mid) value
                      # for the search
  bslo                # Replaces the low bracket value
                      #   with `bscurr`
  bshi                # Replaces the high bracket value
                      #   with `bscurr`
  bsdone [<tol>]      # Returns completion code 0 if search
                      #   bracket has been narrowed to a
                      #   relative precision <tol> (which
                      #   defaults to 10(-14), returns
                      #   completion code 1 otherwise
  bsnotdone [<tol>]   # Logical negation of bsdone.

In practice, of course, we do not (can not!) integrate all the way to x = ∞, instead, we integrate to some finite x = x_{max} where x_{max} is chosen sufficiently large so that, for the specfied convergence criterion δ, we can determine for all possible E whether the solution u(x; E) is diverging to +∞ or to −∞ at x = x_{max}.

Finally, as with many of the problems we have discussed in this course, the toy deuteron problem can be solved “analytically”—however, as in the cases of those other exactly soluble problems we have discussed, the numerical technique which we use to approximately solve this BVP can be extended very easily to solve entire classes of problems which are not amenable to exact solution.