Computational Physics

A7 Differential equations 579

with respect toxon thejth position as follows:

∂^2 ψ(xj,t)

∂x^2

=

1

x^2

[ψj− 1 (t)+ψj+ 1 (t)− 2 ψj(t)]+O(x^2 ). (A.66)

We denote the discretised Hamiltonian byHD, so that the Schrödinger equation can
be written as

i

dψj(t)

dt

=HDψj(t). (A.67)

As a next step, we discretise time using intervalst. Indicating thenth time step
with an upper indexnto the wave functionψ, the time derivative ofψjnat this time

step can be approximated by(ψjn+^1 −ψjn)/t(with an error of ordert), so that
the discretised form ofEq. (A.64)is given by

i

1

t

(ψjn+^1 −ψjn)=

1

x^2

[−ψjn− 1 (t)−ψjn+ 1 (t)+ 2 ψjn(t)]+Vjψjn (A.68)

or, in shorthand notation:

ψjn+^1 =( 1 −itHD)ψjn. (A.69)

The stability of this method can be investigated using the so-called Von Neumann
analysis. In this analysis, one considers the analytical solution to Eq. (A.69) forVj
constant:

ψjn=ξnexp(ikjx). (A.70)

The wave vectorkdepends on the boundary conditions, andξandkare related –
the relation can be found by substituting this solution inEq. (A.68). If there exists
such a solution with|ξ|>1, it follows that a small perturbation in the solution
tends to grow exponentially, because an expansion of a generic small perturbation
in terms of the solutions(A.70)will almost certainly contain the|ξ|>1 modes.
For our algorithm, we find

ξ= 1 −it

[

4

^2 x

sin^2 (kx/ 2 )+Vj

]

(A.71)

which means that in all cases|ξ|>1, indicating that this method is always unstable.
Therefore, we shall consider a modification of the method.
First, the wave function occurring in the right hand side of(A.69)is calculated
at timet=(n+ 1 )t. We then arrive at a form which reads in shorthand:

ψjn+^1 =ψjn−itHDψjn+^1. (A.72)

This algorithm seems impractical, however, because to determineψn+^1 on the left
hand side of(A.72)it is needed on the right hand side: it is animplicitform. We can