Physics 6730 Crank-Nicholson-Crout Algorithm for the Time-Dependent Schrödinger Equation

The time-dependent Schrödinger equation in one dimension

$$\frac{\partial \psi(x,t)}{\partial t} = \frac{i\hbar}{2m}... ...l^2 \psi(x,t)}{\partial x^2} - \frac{iV(x)}{\hbar} \psi(x,t)$$ (1)

is a parabolic partial differential equation. We usually seek a solution on an interval $x \in [a,b]$ and $t > 0$. The solution is uniquely determined from boundary conditions:

$$\psi(a,t) = \psi(b,t) = 0 \mbox{and} \ \psi(x,0) = g(x).$$ (2)

One method for numerical solution solves for the values of the wavefunction on a regular grid of dimension $h = (b-a)/N_x$ in $x$ and $\tau$ in $t$:

$$\psi_j^{(k)} = \psi(a + jh,k\tau).$$ (3)

The derivatives are replaced by simple finite differences. The right side of the equation at the grid point $(i,j)$ is then

$$\frac{i\hbar}{2m h^2} (\psi_{j+1}^{(k)} + \psi_{j-1}^{(k)} ... ...h)/hbar \psi_j^{(k)} = \sum_{m = 0}^N i H_{j,m} \psi_m^{(k)}$$ (4)

where $H$ is a real symmetric tridiagonal matrix (provided $V(x)$ is real). The left side of the equation can be replaced either by a forward difference

$$\frac{\psi_j^{(k+1)} - \psi_j^{(k)}}{\tau},$$ (5)

which, when combined with the rhs gives the explicit algorithm

$${\bf\psi}^{(k+1)} = (1 + i H \tau) {\bf\psi}^{(k)},$$ (6)

or by a backward difference

$$\frac{\psi_j^{(k)} - \psi_j^{(k-1)}}{\tau},$$ (7)

leading to

$${\bf\psi}^{(k)} = {\bf\psi}^{(k-1)} + i H \tau {\bf\psi}^{(k)},$$ (8)

or (with the replacement $k \rightarrow k+1$), the implicit algorithm,

$$(1 - i H \tau){\bf\psi}^{(k+1)} = {\bf\psi}^{(k)}.$$ (9)

The Crank Nicholson Algorithm averages Eqs (6) and (9):

$$(1 - i H \tau/2){\bf\psi}^{(k+1)} = (1 + i H \tau/2) {\bf\psi}^{(k)},$$ (10)

This method is a second order algorithm in $t$, i.e. the discretization error decreases as $\tau^2$. The finite difference representation of the second derivative $d^2/dx^2$ is also good to second order in $h^2$. The Crank-Nicholson Algorithm also gives a unitary evolution in time. That is especially useful for quantum mechanics where unitarity assures that the normalization of the wavefunction is unchanged over time.

The algorithm steps the solution forward in time by one time unit, starting from the initial wave function at $t = 0$. According to the Crank-Nicholson scheme, the time stepping process is half explicit and half implicit. The implicit part involves solving a tridiagonal system. That solution is accomplished by Crout reduction, a direct method related to Gaussian elimination and LU decomposition.

To simplify the algorithm we have chosen units in which the Planck constant $\hbar = 1$, time step $\tau = 1$ and the spatial separation $h = 1$. This can always be arranged by an appropriate redefinition of mass and potential:

$$m = \mbox{(SI mass)} h^2/(\tau \hbar) \ and V = \mbox{(SI V)} \tau/\hbar.$$ (11)

We now present the algorithm in pseudocode. Note, the algorithm defines psi[0] and psi[nx] to be zero to make the coding more economical at the expense of a tiny bit of extra calculation at the ends.

  
Let nx = the number of data points + 1
Let psi[x] be the complex wave function for x = 0,...,nx
    with psi[0] = psi[nx] = 0.
Let m be the mass
Let V(x) be the potential 

Step 1 	INPUT psi[x] = input data for x = 1,...,nx-1
	Set psi[0] = psi[nx] = 0
 	Set lambda = i/(2*m)
	Set u[0] = 0

Step 2 	(Initialization of C-N tridiagonal matrix)
	For x = 1 to nx-1
		Set l[x] = 1 + lambda + iV(x)/2 + lambda/2 * u[x-1] 
		Set u[x] = -lambda/(2*l[x]) 

Step 3	(Iteration of C-N algorithm for nt time steps)
	For t from 1 to nt
		Set z[0] = 0
		(Solution of tridiagonal system by Crout reduction)
		For x from 1 to nx-1
			Set z[x] = ((1-lambda) psi[x] - iV(x)/2 * psi[x] + 
	     		  lambda/2 * (psi[x+1] + psi[x-1] + z[x-1]))/l[x] 
		(Back substitution)
		For x from nx-1 to 1
			Set psi[x] = z[x] - u[x]*psi[x+1] 

Step 4	OUTPUT psi[1] through psi[nx-1]