\documentclass[12pt]{article}
\usepackage{epsfig}
%\topmargin=+0.0in
%\oddsidemargin=0.25in
\textwidth=6.0in
%\textheight=8in
%\footskip=6ex
%\footheight=2ex
\begin{document}
\centerline{\bf Solving the Schoedinger Equation Numerically}
\bigskip
\begin{center}
{\bf Solving the Non-Relativistic Schroedinger Equation \\
for a spherically symmetric potential}
\end{center}
\bigskip
If the energy of a particle is non-relativistic, and its interaction
is described by a potential energy function, the "physics" is described
by solutions to the the time independent Schr\"{o}dinger equation:
\begin{equation}
-{{\hbar^2}\over{2m}} \nabla^2 \Psi + V(r) \Psi = E \Psi
\end{equation}
\noindent This equation is the starting point whether one is calculating cross sections
for a scattering experiment or calculating
bound state energies. The equation is the same for each of these types of
problems, but the boundary conditions for $\Psi (\vec{r})$ will be different..
\noindent For bound state solutions, the wavefunction $\Psi$ and the integral
$\int \Psi^* \Psi dV$ over all space must be finite.
Therefore, the "boundary conditions" at infinity are:
$r \rightarrow \infty$, $\Psi$ must approach zero
faster than $1/r$. Since the potential is spherically symmetric,
the angular dependence can be separated from the radial. Writing
$\Psi = R(r) Y_{lm}(\theta, \phi)$ as a product of a radial part times a
{\it spherical harmonic} ($Y_{lm}(\theta, \Phi)$), the above equation reduces to
\begin{equation}
-{{\hbar^2}\over{2m}} ({1 \over r^2} {{d} \over {d r}}
(r^2 {{d R} \over {d r}}) - {{l(l+1)} \over {r^2}} R(r))
+ V(r) R(r) = E R(r)
\end{equation}
A further simplification is obtained by writing $R(r)$ as
$R(r)=u(r)/r$. The radial part of the Schr\"{o}dinger equation finally
becomes
\begin{equation}
-{{\hbar^2}\over{2m}} ( {{d^2 u(r)} \over {d r^2}}
- {{l(l+1)} \over {r^2}} u(r) )
+ V(r) u(r) = E u(r)
\end{equation}
\noindent For $\Psi$ to be finite, $u(0)$ equals $0$, and for bound
states, $u(r)$ goes to zero as $r \rightarrow \infty$.
\noindent The integer $l$ is related to the particles {\it orbital
angular momentum}.
\newpage
In our second assignment, we will solve the Schroedinger equation
numerically. It will be necessary to solve it numerically, since
the potential V(r) will yield a simple analytic solution.
As derived previously, the radial part of the Schroedinger is:
\begin{equation}
-{{\hbar^2}\over{2m}} ( {{d^2 u(r)} \over {d r^2}}
- {{l(l+1)} \over {r^2}} u(r) )
+ V(r) u(r) = E u(r)
\end{equation}
\noindent where $l$ is the orbital angular momentum of the particle
about the center of the potential.
We will solve this differential equation numerically by using a simple
"finite difference" method. There are other algorithms, but here we can
obtain accurate results using the simplest of methods. The main idea is
to make the continuous variables and functions, $r$, $V(r)$, and $u(r)$
discrete. We do this by
making the radial coordinate $r$ discrete with a step size $\Delta$. That is,
$r \rightarrow r(i) = i*\Delta$ where $i$ is an integer. The variable
$r$ and the functions $V(r)$ and $u(r)$ become arrays: $r \rightarrow
r(i)$, $V(r) \rightarrow V(r(i)) \rightarrow V(i)$, and $u(r) \rightarrow
u(r(i)) \rightarrow u(i)$.
The first derivative of $u$ with respect to $r$ is approximately
\begin{equation}
{{du} \over {dr}} \approx {{u(i+1) - u(i)} \over \Delta}
\end{equation}
\noindent from the definition of the derivative. Actually the limit
as $\Delta \rightarrow 0$ gives the derivative. Taking $\Delta$
small enough will give a good approximation. We need the second derivative
in the Schroedinger equation. Using the finite difference approximation
once again with the first derivative yields:
\begin{eqnarray*}
{{d^2u} \over {dr^2}} & \approx & {{ {{u(i+1) - u(i)} \over \Delta} -
{{u(i) - u(i-1)} \over \Delta} } \over \Delta} \\
{{d^2 u} \over {dr^2}} & \approx & {{u(i+1)+u(i-1)-2u(i)} \over {\Delta^2}} \\
\end{eqnarray*}
\newpage
\noindent After substituting these expressions into the differential
equation and doing some algebra, one obtains a discrete version
of the Schr\"{o}dinger equation:
\begin{equation}
u(i+1) = 2u(i) - u(i-1) + \Delta^2 {{l(l+1)} \over {r^2}} u(i) +
{{2m \Delta^2} \over {\hbar^2}} (V(i) - E) u(i)
\end{equation}
\noindent This "discrete Schr\"{o}dinger equation" allows one to
determine all values of $u(i)$ if $u(0)$ and $u(1)$ are known. This
can be done by iteration, since $u(i)$ is completely determined
from $u(i-1)$ and $u(i-2)$.
\bigskip
\noindent {\bf Boundary Conditions}
\bigskip
The Schr\"{o}dinger equation will describe both bound state and
scattering situations. In both cases, the
"boundary condition" on $u(r)$ at $r=0$ is the same, namely $u(0)=0$. For a
bound state solution the other boundary condition is that
$u(i \rightarrow \infty) \rightarrow 0$. Only certain values of
$E$ will allow this limit to be satisfied. These values of $E$ are the
allowed "bound state" energies that the particle can have.\\
\noindent {\bf Bound States}\\
\begin{eqnarray*}
u(0) & = & 0 \\
u(r \rightarrow \infty) & \rightarrow & 0
\end{eqnarray*}
For a scattering situation, $u(r)$ will oscillate as $r \rightarrow \infty$.
We will consider scattering problems later on in the course. Now lets
see how we can solve for the energies of any bound states.\\
\newpage
\noindent {\bf Solving for the bound state energies}
\bigskip
We can determine the bound-state energies using a "bracket and half"
method. First a trial energy $E_t$ is chosen which lies
below the ground state energy. In the discrete Schr\"{o}dinger
equation, $u(i+1)$ is determined from the
values of $u(i)$ and $u(i-1)$. We assign $u(0)$ a
value of $0$, and $u(1)$ is assigned a non-zero value (e.g.
$u(1)=1.0$). The discrete Schroedinger equation can then be used to iterate $u(i)$ to a large value
of $i = imax$, well outside the range of the potential.
We assign $u(imax)$ the value $u(imax)\equiv test0$. Next, the
trial energy is increased by an amount $\delta E$ and the process is
repeated. The function $u(imax)$ will now have a different value which
we call $test1$. If $test0$ and $test1$ have the same sign, then
the trial energy is changed again by an amount $\delta E$, $test1
\rightarrow test0$, and the process is repeated. If $test0$ and $test1$ have opposite signs,
then the wave function at $r=imax*\Delta$ has changed sign and the
trial energy has passed over the ground state energy. The energy step is
reversed and halved, $\delta E \rightarrow -\delta E/2$,
$test1 \rightarrow test0$, and the process is repeated to the desired accuracy.
To determine the energy of the next excited state, one starts
with a trial energy just above the ground state energy. The
trial energy is stepped up in a similar manner until the energy
converges. The next higher allowed energy is found in a similar
manner. For our homework problem, you will solve for the
bound state energies of a valance neutron and a valance proton
in different nuclei.
\end{document}