SYLLABUS Previous: 2.5.1 Advectiondiffusion equation
Up: 2.5 Implicit CrankNicholson
Next: 2.6 Computer quiz
2.5.2 Schrödinger equation
In quantum mechanics, the evolution of a particle is generally modelled
with a complex wave function
. If the particle does't
decay with time, the total probability of finding this particle somewhere
in space of course remains exactly unity at all times
.
Solving the Schrödinger equation normalized to have Planck's constant
, the particle mass
and the static potential
it is clearly an advantage to formulate a numerical scheme that conserves
the probability by keeping the evolution operator unitary even as it is
discretized. This can be achieved by writing a formal solution
in the socalled Cayley form with
This approximation is unitary, accurate to secondorder in time
. It suggests an evolution of the form
Replacing the second order derivative in the Hamiltonian operator
with finite differences centered in space
, one
obtains a scheme that is stable, unitary and is in fact the same
CrankNicholson method applied here for the Schrödinger equation:
which is finally cast into the linear system

(2.5.2#eq.17) 
Exploiting the tridiagonal structure of the matrix, the scheme has been
implemented in JBONE
using complex arithmetic
BandMatrixC a = new BandMatrixC(3, h.length); //Complex objects
BandMatrixC b = new BandMatrixC(3, h.length);
Complex[] c = new Complex[h.length];
Complex z = new Complex();
double[] V = physData.getPotential(); //Heavyside(xL/2) times
double scale = 10.*velocity; // an arb. scaling factor
double dtodx2 = timeStep/(dx[0]*dx[0]);
Complex pih = new Complex(0., 0.5*dtodx2);
Complex mih = new Complex(0.,0.5*dtodx2);
Complex pip1 = new Complex(1., dtodx2);
Complex mip1 = new Complex(1., dtodx2);
for (int j=0; j<=n; j++) {
z = new Complex(0.,0.5*scale*timeStep*V[j]);
a.setL(j,mih); //Matrix elements
a.setD(j,pip1.add(z));
a.setR(j,mih);
b.setL(j,pih); //Right hand side
b.setD(j,mip1.sub(z));
b.setR(j,pih);
}
c=b.dot(h); //Right hand side with
c[0]=c[0].add(b.getL(0).mul(h[n])); // with periodicity
c[n]=c[n].add(b.getR(n).mul(h[0]));
hp=a.solve3(c); //Solve linear problem
for (int j=0; j<=n; j++){ //Plot norm & real part
fp[j]=hp[j].norm();
gp[j]=hp[j].re();
}
The code again relies on the complex
BandMatrixC.solve3()
method to solve the linear system efficiently with
complex operations using the standard LUfactorization
that will be studied in 3.
The applet below illustrates the scheme
with a calculation of the scattering of a wave packet on a square
potential barrier.
JBONE applet: press Start/Stop
to compute the scattering of a wave packet on a square potential
barrier with a height specified with PhysDataValue and a
width extending throughout the right side of the periodic domain.
The real part of the wave function (blue line) starts to oscillates
and the evolution shows how the probability density (black line)
splits into a reflected and a transmitted part.

The perfect conservation of the lowest order moment shows that the
total probability is indeed conserved for all times.
SYLLABUS Previous: 2.5.1 Advectiondiffusion equation
Up: 2.5 Implicit CrankNicholson
Next: 2.6 Computer quiz
