580 Appendix A
make this explicit by an inversion:
ψjn+^1 =∑
j′( 1 +itHD)−jj′^1 ψjn′. (A.73)As the matrix in brackets is tridiagonal, the inversion can be done efficiently – it
requires onlyO(L)steps, i.e. of the same order as the other steps of the algorithm.
The inversion method will be discussed in Appendix A8.1.
Performing the Von Neumann analysis for this method, we find
ξ=1
1 +it[( 4 /^2 x)sin^2 (kx/ 2 )+Vj]. (A.74)
This looks more promising than the first method: foreverychoice oftandx
one gets|ξ|<1, so that the method is stable. It is important to taketsmall as
the time-derivative ofψihas a first order accuracy intonly. The accuracy of the
discretised second derivative with respect toxis still of orderx^2.
We have neglected an important point: we would likeψnto be normalised at
every time step. This is unfortunately not the case for the methods discussed up to
now. Writing the exact solution to the continuous differential equation (A.64) as
ψ(x,t)=exp(−itH/)ψ(x,0), (A.75)we notice that for all timest, the norm ofψis equal to that att = 0 since
exp(−itH/)is a unitary operator:
〈ψ(t)|ψ(t)〉=〈exp(−itH/)ψ( 0 )|exp(−itH/)ψ( 0 )〉=〈ψ( 0 )|ψ( 0 )〉. (A.76)The operator 1±itHDis in general not unitary. If we are able to find a unitary
operator which gives us at the same time a stable method, we have the best method
we can think of (apart from finite accuracy inherent to discretisation). It turns out
that the operator
1 −itHD/ 2
1 +itHD/ 2
(A.77)
does the job, and surpasses expectations: not only is it unitary and stable, it is even
correct tosecondorder int. This can be seen by realising that using this matrix
is equivalent to taking the average ofHacting onψnandψn+^1 , or by noticing that
the matrix in (A.77) is a second order approximation to exp(−itH/)int.
This method is recommended for solving the time-dependent Schrödinger equa-
tion using finite difference methods. It is known as theCrank–Nicholsonmethod.
Note that the presence of the denominator in the operator (A.77) makes this an impli-
cit method: a matrix inversion must be carried out at every time step, requiringO(L)
calculations, which is the overall time scaling of the algorithm.