\begin{center}
{\bf PHYS 410/555 Computational Physics}\\
{\em Solving Two-Point Boundary Value Problems Using ``Shooting''\\}
{\em A Toy Model for the Deuteron}
\end  {center}
%
Recall that, by definition, {\em 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)
{\em 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 {\em search}---in 
general we will {\em not} be able to construct an algorithm which in ``one go'' 
results  in the eigenvalue and the associated {\em 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 {\em 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 {\em 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:
\begin{equation}
\label{bvp1}
   u''(x) = f(u,u',x)  \quad\quad\quad 0 \le x \le 1 
\end{equation}
subject to the boundary conditions
\begin{eqnarray}
\label{bvs1a}
	u(0) & = & u_0 \\
\label{bvs1b}
   u(1) & = & u_1
\end{eqnarray}
where $u_0$ and $u_1$ are specified constants.
Assume that we have somehow determined a function $U(x)$ which 
satisfies~(\ref{bvp1}), (\ref{bvs1a}) and (\ref{bvs1b}).  We then note that 
if we now consider the {\em initial value problem} defined by~(\ref{bvp1})
subject to the {\em initial conditions}:
\begin{eqnarray}
\label{ics1a}
	u(0) & = & U(0) = u_0 \\
\label{ics1b}
   u'(0) & = & U'(0) 
\end{eqnarray}
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, $\theta$, of an artillery gun (i.e. $\theta$ is the angle 
the gun's barrel makes with the 
horizontal), so that the shell hits a target at some given range.
Assuming that $\theta < \pi/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 $\theta$ 
so that the next shot comes closer to the target.

Similarly, when integrating a BVP such as~(\ref{bvp1}) 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) \longrightarrow u(1) < u_1$ and 
$u'(0) < U'(0) \longrightarrow 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~(\ref{bvp1}) 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 {\em 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, $\lambda$, 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
$[\lambda_-,\lambda_+]$ 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 {\em 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 {\em highly} idealized; we assume that 
the deuteron wave function, $\psi({\vec r})$, is a {\em spherically symmetric}
solution of the time-independent Schr\"odinger equation, with the 
proton-neutron interaction described by a square wave potential.  

Thus we wish to solve
\begin{equation}
\label{sch1}
-\frac{\hbar^2}{2M} \nabla^2 \psi({\vec r}) + V({\vec r} \psi({\vec r})  = 
       E \psi({\vec r})
\end{equation}
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 
\begin{equation}
	\psi({\vec r}) \rightarrow \psi(r)
\end{equation}
and 
\begin{equation}
	\nabla^2 \psi(r) =
   \frac{1}{r^2}\frac{d}{dr}\left(r^2 \frac{d\psi}{dr}\right)
\end{equation}
Further, defining
\begin{equation}
	u(r) \equiv r \psi(r)
\end{equation}
we have (as you should verify)
\begin{equation}
	\nabla^2 \psi(r) \rightarrow \frac{1}{r} \frac{d^2 u(r)}{dr^2}
\end{equation}
Thus, (\ref{sch1}) can be rewritten:
\begin{equation}
\label{sch2}
	\frac{d^2u(r)}{dr^2} + \frac{2M}{\hbar^2}\left(E - V(r)\right)u(r) = 0
\end{equation}

As noted above, we will model the proton-neutron interaction as a 
finite-range square potential.  Thus, we take
\begin{eqnarray}
   V(r) & = &  V_0  \quad\quad\quad  0 \le r \le a \\
   V(r) & = &  0    \quad\quad\quad\,\,\,\,  r > a
\end{eqnarray}
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\"odinger equation, we must 
demand that our solution of~(\ref{sch1}) be normalizable, i.e. that 
\begin{equation}
\label{norm1}
	\int \psi \psi^* dV = 1
\end{equation}
so that there is unit probability that our deuteron is found {\em somewhere}
in the universe.  In the current spherically symmetric case this means that 
we must have 
\begin{equation}
\label{norm2}
	4 \pi \int_0^\infty r^2 \psi(r)^2 dr =
   4 \pi \int_0^\infty u(r)^2 dr = 1
\end{equation}
Clearly, a necessary condition for normalizability is that 
\begin{equation}
	\lim_{r\to\infty} u(r) = 0
\end{equation}
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\"odinger equation.

Before we consider the solution of~(\ref{sch2}) 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$
\begin{equation}
	x \equiv \sqrt{2M} r 
\end{equation}
so that, among other things, we have
\begin{equation}
	\frac{d^2u}{dr^2} \to 2M\frac{d^2u}{dx^2}
\end{equation}
Further, we choose units such that $\hbar=1$ and $V_0 = -1$ (you should 
establish that this {\em 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
\begin{equation}
	x_0 \equiv \sqrt{2M} a
\end{equation}
as the free parameter.

Thus, our Schr\"odinger equation~(\ref{sch2}) becomes
\begin{equation}
\label{sch3}
	\frac{d^2u(x)}{dx^2} + \left( E - V(x) \right) u(x) = 0
\end{equation}
where
\begin{eqnarray}
	V(x) & = & -1 \quad\quad\quad 0 \le x \le x_0 \\
	V(x) & = &  0 \quad\quad\quad\quad x > x_0
\end{eqnarray}
and where we again note that $E = E(x_0)$ is an {\em eigenvalue} 
of~(\ref{sch3}); 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~(\ref{sch3}) are derived from the demands that 
\begin{enumerate}
\item $\psi(r)$ be regular (analytic) at $r=0$.
\item $\lim_{r \to \infty} \psi(r) = 0$.
\end{enumerate}

The regularity condition at $r=0$ means that $\psi(r)$ admits an 
expansion
\begin{equation}
	\lim_{r\to0} \psi(r) = \psi_0 + r^2 \psi_2 + O(r^4)
\end{equation}
where the $\psi_i$ are constants (the power series expansion cannot have 
terms which are odd in $r$, since $\psi(r)$ would not have a well-defined 
derivative at $r=0$ in that case).  From this, it follows that $u(r)$ 
has an expansion
\begin{equation}
	\lim_{r\to0} u(r) = r\psi(r) = r\psi_0 + r^3\psi_2 + O(r^5)
\end{equation}
Thus, we have 
\begin{eqnarray}
	u(0) & = & 0 \\
	\frac{du}{dr}(0) & = & \psi_0 
\end{eqnarray}

We now observe that~(\ref{sch1}) (like {\em all} Schr\"odinger equations)
is {\em linear}; given any solution, $\psi(r)$ we have that $c \psi(r)$, 
where $c$ is an arbitrary positive constant, is also a solution. 
The {\em particular} solution we seek is fixed by the normalization 
condition:
\begin{equation}
	\int \psi \psi^* dV = 1
\end{equation}
Operationally, this means that we can choose $\psi(0)$ {\em arbitrarily}
(say $\psi(0) = 1$ for convenience), then, for specified $x_0$, 
vary $E$ until we find a solution which satisfies 
\begin{equation}
	\lim_{r \to \infty} \psi(r) = 0
\end{equation}
(In other words, the eigenvalue is independent of the normalization of 
the eigenfunction).

In preparation for a solution of our problem using {\tt LSODA} we 
rewrite~(\ref{sch3}) in canonical first order form by introducing 
\begin{equation}
	w(x) \equiv \frac{du(x)}{dx}
\end{equation}
We then have
\begin{eqnarray}
\label{ode1}
	\frac{du(x)}{dx} & = & w(x) \\
	\frac{dw(x)}{dx} & = & (V(x) - E) u(x) 
\end{eqnarray}
subject to 
\begin{eqnarray}
	u(0) & = & 0 \\
\label{iv2}
	w(0) & = & 1 
\end{eqnarray}
and with $E$ to be determined so that 
\begin{equation}
	\lim_{x \to \infty} u(x) = 0
\end{equation}

\def\ELO{E_{\rm LO}}
\def\EHI{E_{\rm HI}}
\def\EMID{E_{\rm MID}}
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 
\begin{eqnarray}
\label{eplus}
	\lim_{x \to \infty} u(x;E_-) = -\infty  \\
\label{eminus}
	\lim_{x \to \infty} u(x;E_+) = +\infty 
\end{eqnarray}
where the notation $u(x;E_i)$ means the trial solution $u(x)$ computed 
using eigenvalue estimate $E_i$.  We further assume that 
properties~(\ref{eplus}) and~(\ref{eminus}) hold for {\em any} values 
$\EHI$ and $\ELO$ satisfying $\EHI > E$ and 
$\ELO < E$ respectively, namely:
\begin{eqnarray}
	\lim_{x \to \infty} u(x;\ELO) = -\infty  \\
	\lim_{x \to \infty} u(x;\EHI) = +\infty 
\end{eqnarray}

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

\begin{obeylines} 
{\tt 
~~~~~~$\ELO := E_-$
~~~~~~$\EHI := E_+$
~~~~~~while $\EHI - \ELO > \delta$ do
~~~~~~~~~$\EMID := (\EHI + \ELO)/2$
~~~~~~~~~if $u(x;\EMID) \to +\infty$ as $x \to \infty$ then
~~~~~~~~~~~~$\EHI := \EMID$
~~~~~~~~~else
~~~~~~~~~~~~$\ELO := \EMID$
~~~~~~~~~end if
~~~~~~end while
~~~~~~$E := \EMID$
}
\end{obeylines}
~\\[0.05in]

We note that the convergence rate of the bisection search is completely 
pre-determined; the size of the bracketing interval after $n$ bisections
is 
\begin{equation}
  \frac{1}{2^n} \left( E_+ - E_- \right)
\end{equation}

The solution of equations~(\ref{ode1})-(\ref{iv2}) using {\tt LSODA}
is implemented as the program {\tt deut.f}; the bisection search 
for the eigenvalues $E(x_0)$ is implemented via the shell-script 
{\tt Shoot-deut}, which itself uses the following collection of scripts which
can be used to implement shell-level bisection searches:
\begin{verbatim}
  bsnew <lo> <hi>     # Initializes a new search
  bscurr              # Returns the current (mid) value for the search
  bslo                # Replaces the low bracket values with `bscurr`
  bshi                # Replaces the high bracket values 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.
\end{verbatim}

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

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 {\em not} amenable to exact solution.

