Background Theory:Using the Shooting Method to Solve the Time-Independent Schrödinger Equation for a 1D Quantum Well¶
Go back to the interactive notebook
Shooting method¶
In numerical analysis, the shooting method is a method for solving a boundary value problem by reformulating it as an initial value problem. Roughly speaking, we 'shoot' out trajectories in different directions from an initial trial value until we find a trajectory that has the desired boundary value. You can check out the following link for additional information on the method https://en.wikipedia.org/wiki/Shooting_method.
For the specific example of the one-dimensional time-independent Schrodinger equation with a quantum well potential, we know that the wavefunction will converge to zero at both the far left and right boundaries in order for the wavefunction to be normalizable (i.e. $\psi(x_{\pm \infty})=0$). By keeping the boundary value at the left hand side equal to zero, one can try different eigenvalues of the Schrödinger equation and obtain the corresponding eigenfunctions (by means of a numerical integrator such as the Numerov algorithm discussed below). Only the true eigenvalue will result in the solution wavefunction converging to zero at the right hand side. By scanning over the possible trial energies and monitoring the solution wavefunction at the right hand boundary, we can find all allowed eigenvalues and their corresponding wavefunctions. This numerical method is referred to as the shooting method. Through its use can obtain the eigenvalues and eigenfunctions of the Schrödinger equation for this 1D quantum well.Numerical integration and the Numerov algorithm¶
The time-independent Schrödinger equation is an ordinary differential equation (ODE) of second order, where the 1st-order term does not appear in the equation, i.e. it assumes the following structure: $\large \dfrac{d^2 y}{d x^{2}} = -g(x)y(x) + s(x)$
In the particular case of the time-independent Schrödinger equation, we have:
$\large \left[ -\dfrac{\hslash^2}{2m} \, \dfrac{\partial^2}{\partial x^{2}} + V(x)\right] \psi(x) = E\psi(x)$ (1)
and so $g(x)= \frac{2m}{\hslash^2}(E-V(x))$ and $s(x)=0$.
For a one dimensional system, the second-derivative can be evaluated numerically via the following formula
$\large \psi ''(x_{i})= \dfrac{1}{\delta x^2}\left[ \psi(x_{i+1})-2\psi(x_i)+\psi(x_{i-1}) \right]$ (2)
where $x_i$ gives the position at the i-th point on a discretized grid of $i=1,...,N$ points representing space in the x-dimension and $\delta x = x_{i+1}-x_{i}$ is the grid spacing.
Substituting equation 2 into equation 1, we can create an iterative procedure for generating the wavefunction $\psi(x)$:
$\large \psi(x_{i+1}) =\delta x^2 \psi ''(x_{i}) +2\psi(x_i)-\psi(x_{i-1}) = -\dfrac{2m \delta x^2}{\hslash^2} \left[E-V(x_i)\right]\psi(x_i) +2\psi(x_i)-\psi(x_{i-1})$
I.e. if we know the value of $\psi$ at two preceding points $x_i$ and $x_{i-1}$, then we can obtain the value of $\psi$ at the next point $x_{i+1}$. Carrying this out for all values of $i$, we obtain our solution wavefunction.
However, the values of the first two starting points are unknown. For the square well potential shown in the interactive notebook, we can assume $\psi(x_0)$ is zero and $\psi(x_1)$ is a very small positive (or negative) number. See task 4 in the interactive notebook for further discussion of the issue of initial conditions.
There are occasions where the above approximation to the derivative is simply not accurate enough for the problem at hand. In this case, higher-order approximations must be employed.
The Numerov method is one such higher-order method. It is used to specifically solve the kind of
ODE which has a form like that of the time-independent Schrödinger equation, i.e., one having the form $\dfrac{d^2 y}{d x^{2}} = -g(x)y(x) + s(x)$. The method capitalizes on this particular form to approximate the solution to order $O((\delta x)^6)$, where $\delta x$ is the step size for the integration. The method works by allowing one to relate the value of the solution at a given point on a discretized grid representing space, $y_{n+1}$, to the two previous points, $y_n$ and $y_{n-1}$, through the relationship:
$\large y_{n+1}\left(1+{\frac {(\delta x)^{2}}{12}}g_{n+1}\right)=2y_{n}\left(1-{\frac {5(\delta x)^{2}}{12}}g_{n}\right)-y_{n-1}\left(1+{\frac {(\delta x)^{2}}{12}}g_{n-1}\right)+{\frac {(\delta x)^{2}}{12}}(s_{n+1}+10s_{n}+s_{n-1})+O((\delta x)^{6})$
where $s_n = s(x_n)$ and $g_n = g(x_n)$.
See https://en.wikipedia.org/wiki/Numerov's_method for a detailed derivation of the method.