A numerical solution of the one- dimensional, time-independent Schroedinger equation.

Reference: S. E. Koonin, Computational Physics, Addison-Wesley Publishing Company Inc, (1986).

The program n3.txt is written in QuickBasic.  It demonstrates a method for solving the time-independent Schroedinger equation in one dimension.  This is a boundary value problem. We solve for the eigenvalues E of the Hamiltonian for which the wave function and its derivative are continuous and for which the wave function takes on given values at the boundaries of the integration volume.  The potential in the program can easily be modified.

Explanation:

We want to solve the eigenvalue equation

or

,

with

.

We want to solve it in a finite region between x=0 and x=L.

We assume that V(x)=0 at x=0 and x=L.  We also assume that V(x)<0 somewhere in the region between x=0 and x=L.  Because we cannot extend our numerical solution to infinity, we assume that V(x)=¥ for x<0 and x>L.

We want to find the energies E<0 for which the wave function y(x) and its derivative are continuous and the boundary conditions y(0)=y(L)=0 are satisfied.  These are the energies of the bound states of the system.

We start by expanding y(x) in a Taylor series expansion.

.

.

Combining the two equations above yields

      [equation (a)]

Here O((Dx)⁴) denotes terms of higher than fourth order in Dx.

From the Schroedinger equation

we obtain for y''''

or

.

Inserting this expression into equation (a) and neglecting terms of higher than fourth order in Dx we obtain

.

For our numerical work, let {x_n} denote points on a grid defined in the region from x=0 to x=L, n=0, 1, …, N.

Define y_n=y(x_n), k_n=k(x_n), dx=x_n+1-x_n. Then

.

Integrating the wave function using this algorithm, i.e. finding its value on the next grid point from its value at the previous two grid points, is called integrating using the Numerov method.

To solve for the bound states of the system, we pick y₀=0, y₁=c₁, where c₁ is some small number.  This number c₁ determines the overall normalization of the wave function.  We now calculate y_2, y₃, etc.  We start integrating in the classically forbidden region, where the magnitude of the wave function increases approximately exponentially.  As we reach the classically allowed region, the wave function becomes oscillatory at the classical turning point. As we pass through the second turning point and enter again into the classically forbidden region, the integration becomes numerically unstable, because it can contain an admixture of the exponentially growing solution.  Integration into a classically forbidden region is likely to be inaccurate.

We therefore generate a second solution, picking y_N=0, y_N-1=c₂, and calculating y_N-2, y_N-3, etc.

For both integrations we pick E=V_min+DE, i.e. we pick E just slightly larger than V_min.  To determine whether the energy is an eigenvalue, we compare the results of our integrations at a matching point x_m in the classically allowed region.

We first re-normalize the wave functions y¹ and y² obtained from our two integrations such that y¹_m=y²_m, i.e. such that the two functions have the same value at the matching point.  We then check if the derivatives of the two functions are equal at x_m, i.e. if .

We use

.

We already have made y¹_m=y²_m.  We therefore only have to check if y¹_m-1=y²_m-1. We define

.

If f¹0 within a chosen limit of tolerance, then we repeat our integrations after incrementing E by DE.

If f¹0 for two successive iterations, but f has changed sign, then the eigenvalue has been overshot.  We decrease the magnitude and change the sign of DE and start decrementing E.  We repeat until f=0 within the chosen limit of tolerance.

How do we adjust DE to find the value of E for which f=0?

We want to find the root of f(E).

Take an initial guess E_i.  Expand f(E) about E_i.  Let

.

E_i+1 is the next better value for the root of f(E).

.

.

Therefore

.