Schrödinger's Equation
(draft version)

Terry Robb
<tdr@vaxc.cc.monash.edu.au>

Introduction

This Notebook is in response to a message that I posted to the MathGroup on 26-Feb-1993. The animation is in the Evolution section below.

Date: 26 Feb 1993 15:49:52 +1100
From: TDR@vaxc.cc.monash.edu.au
Subject: Schroedinger Equation and InterCall

Dear Mathgroup,

>Date: 20 Feb 1993 21:45:04 -0600
>From: gmk@pegasos.ccsr.uiuc.edu (Gottfried Mayer-Kress)
>Subject: MathLink and Integration of Schroedinger Equation
>
>We want to use mathlink to integrate the one-dimensional
>time-dependent Schroedinger equation with arbitrary potential.
>We think about using the numerical recipe algorithm (pg 851 of the
>most recent edition). Has anybody implemented that or can recommend
>some other package?

I've dusted off one of my old Crank-Nicholson solvers (written in fortran)
and adjusted it so that it can solve Schroedinger-type equations of the form:
        I*D[psi,t] == -D[psi,{x,2}] + V[x,t]*psi
where psi[x,t], a<x<b, t>0 is the complex wave-function, and
the equation is solved subject to some initial condition:
        psi[x,0]==psi0[x] for a<x<b,
and some linear time-dependent boundary conditions, ie.
        a0[t]*psi[a,t]+a1[t]*psi'[a,t] == ac[t],
        b0[t]*psi[b,t]+b1[t]*psi'[b,t] == bc[t].

By importing the fortran subroutine (named SCHROED) into Mathematica using
InterCall I was able to do a few experiments involving shooting a wave-packet
at a potential barrier, and watching an animation of the packet being
mostly reflected but also partially tunnelling through the wall.

I could supply the potential function V[x,t] (required by the fortran
SCHROED routine) by writing it as a Mathematica function which gave me
a lot of flexibility. Typically  V=Function[{x,t}, If[x<4 || x>6,0,40]]
for a potential wall, but I could set it to any Mathematica expression.

Also the initial condition for a wave-packet was computed in Mathematica
too. The fortran subroutine SCHROED just advanced the solution an amount
dt of time per call, and so I simply called the routine inside a Mathematica
loop and did a ListPlot of the evolution.

I used 400 points to represent the solution, and the rendering of each plot
was the slow part. Calling SCHROED from within Mathematica using InterCall
was very fast.

If anyone is interested in the details, I'll put a Notebook at MathSource
about how easy it was to do using the InterCall package.

Terry.

  ###################################################################
  Dr Terry Robb                              Email  : tdr@vaxc.cc.monash.edu.au
  Department of Mathematics           Phone: +61 (3) 565-5666
  Monash University                        Fax    : +61 (3) 565-4403
  Clayton 3168, Melbourne, AUSTRALIA.
  ###################################################################

Discussion

About this Notebook

This Notebook uses Mathematica and InterCall to communicate with an external fortran subroutine that inplements the Crank-Nicholson scheme for numerically solving an equation such as Schrödinger's.

Steps involved in doing this are:

• load the InterCall package into Mathematica

• import the desired external routine named SCHROED

• call SCHROED, sending initial condition

• call SCHROED in a loop, to evolve the wave function

The rest of the Notebook is just calling Mathematica graphics and stuff not directly related to calling or solving the equations.

You may be interested in timing information. For the calculation in this Notebook I have used 2001 points and 1000 time-steps to produce 26 frames. On my machine (a DECstation 5000) the cpu-time spent inside the SCHROED routine is about 0.15 seconds per time-step but involves also about 0.3 seconds of communication. Each frame of the animation thus consumes about 10 seconds of time to call SCHROED and a further 10 seconds to render the plot on a Macintosh IIfx. Overall it therefore takes about 10 minutes of real-time for the 26 frames, with most of the time being used by Mathematica graphics.

The good old days

For an interesting comparison of how an animation of Schrödinger's equation was done 25 years ago, see the paper by Goldberg et al. on Computer Movies.

Their figure 5 is the equivalent of the animation given in this Notebook.

[Goldberg A, Schey H, Schwartz J. 1967, American Journal of Physics vol. 35, p177-186.]

Setup

<<InterCall.m

We wish to import the fortran subroutine named SCHROED. To do this it is necessary to know the detailed datatype structure of the arguments to SCHROED and these are given in the header for the fortran routine, ie

subroutine SCHROED(PSI,dt, t,dx,a,b,n, wrk, Vfunction,LinearBC)
        complex*16 PSI(n),wrk(n)
        real*8 dt,t,dx,a,b,Vfunction
        integer n
        external Vfunction,LinearBC
        ...

The SCHROED subroutine advances the array PSI (of length n) from time t to t+dt subject to some user-supplied potential function and boundary conditions.

The potential V[x,t] is supplied as a function called Vfunction with header:

real*8 function Vfunction(x,t)
        real*8 x,t
        ...

And the coefficients involved in the boundary conditions are supplied as the subroutine called LinearBC with header:

subroutine LinearBC(a0,a1,ac,b0,b1,bc,t)
        complex*16 a0,a1,ac,b0,b1,bc
        real*8 t
        ...

With knowledge of the above we can compose an appropriate default setting that mentions the datatype structure of SCHROED to InterCall. Thus

AddDefault[
   SCHROED[$PSI->$PSI,$DT,$T->$T,
           $DX->($B-$A)/($N-1),$A,$B,$N->ROWS[$PSI],
           $WRK:>Null,$VFUNC->$VFUNC,
           $LINEARBC->Function[{a0,a1,ac,b0,b1,bc,t},
             a0=1; a1=0; ac=0; b0=1; b1=0; bc=0;] ],
   S[Z[$N],R,R,R,R,R,I,Z[$N],RF[R,R],S[Z,Z,Z,Z,Z,Z,R]],
   "schroed.o"];

The AddDefault commands tells InterCall three things. Firstly it specifies some sensible default settings to save the user typing any unnecessary arguments. Secondly in declares the datatype structure of the routine, and finally it mentions the object file where the external code is to found. Having told InterCall these things about SCHROED we can proceed to import it into Mathematica.

Import[{SCHROED}];

Once imported a routine named schroed will be available from within Mathematica.

Initial condition

We want to setup an initial wave-packet, ideally over an infinite region. Numerically an infinite region is not an easy thing to represent, so we take infinity to be somewhat smaller (20) and break the resulting region [-20,20] into 2001 equally spaced steps. This gives a resolution of 0.02 units along the x-axis. (Thoughout this Notebook units are chosen so that h/(2π)=1 and mass=1/2.)

inf = 20.;
dx=0.02; dt=0.001;

Now we compute an initial wave-packet psi0 of width 2a, group velocity 2k, and centred at point x0 along the x-axis. For consistency we make sure it satisfies the boundary conditions used later, and also we normalize psi0 so that the total probability is 1 exactly.

k = 15; a = 1; x0 = -inf/2.;
psi0=Table[Exp[I*k*(x-x0)-(x-x0)^2/(2*a)],
           {x,-inf,+inf,dx}];

psi0[[1]]=0; psi0[[Length[psi0]]]=0;     (* apply BC  *)
psi0=psi0/Sqrt[psi0.Conjugate[psi0]*dx]; (* normalize *)

Here is a plot of the real-part of psi0. Note there are roughly k*Sqrt[a] waves inside the packet.

ListPlot[psi0//Re, PlotRange->All,PlotJoined->True];1;

Wave-packets of initial form Exp[I*k*x-x^2/(2*a)] naturally spread out over time. Initially the spread or half-width of the packet is a, but after a time t its spread becomes Sqrt[a^2+(2*t)^2].

The wave-packet will be aimed at a potential wall of strength u0, placed between the points xp0 and xp1 on the x-axis. The energy of the incoming wave-packet is k^2 (in the units mentioned above) and so if the strength u0 of the potential is made large enough it will make it hard for the wave-packet to penetrate.

u0=250; xp0=2.; xp1=3.33;
potential = Function[{x,t}, If[x<xp0 || x>xp1, 0, u0]];

Some routines

Now we setup a better default for SCHROED so that we need not always specify the region. In addition parameters for a clamped boundary condition are pre-compiled (with the InterCall compiler) so that $LINEARBC is not unnecessarily re-compiled on each call to schroed.

AddDefault[CLAMP->Function[{a0,a1,ac,b0,b1,bc,t},
             a0=1; a1=0; ac=0; b0=1; b1=0; bc=0;],
           S[Z,Z,Z,Z,Z,Z,R] ];
SetDefault[schroed, $A-> -inf, $B->inf, $LINEARBC->CLAMP]

With this setting we can define a routine called evolve[deltatime] that will advance the current solution by an amount deltatime of time. Internally schroed advances the solution only by a time dt so we must call schroed a number of times (m) and then do a final advance to make a total of exactly deltatime.

evolve[deltatime_] := Block[{m=Floor[deltatime/dt]},
    Do[schroed[dt],{m}]; schroed[deltatime-dt*m]; ]

The routine drawit[] simply displays the solution $PSI at the current iteration. It displays a sample of only 500 points (roughly) because otherwise ListPlot would take a long time to plot if the full list were used. The InterCall commands Peek[$PSI] and Peek[$T] pull back the current solution and time from the external schroed routine.

drawit[] := Block[{psi,t,pot,p1,p2,n,dn,nps=500},
   psi = Peek[$PSI]; n=Length[psi]; dn=(n-1)/nps;
   psi = Map[First,Partition[psi,dn]]; t=Peek[$T];
   pot = Table[potential[x,t]/k^2*.7,{x,-inf,inf,dn*dx}];
   p0=ListPlot[pot,      PlotJoined->True,
               PlotStyle->RGBColor[0,0,1],
               DisplayFunction->Identity];
   p1=ListPlot[psi//Re,  PlotJoined->True,
               DisplayFunction->Identity];
   p2=ListPlot[psi//Abs, PlotJoined->True,
               PlotStyle->RGBColor[1,0,0],
               DisplayFunction->Identity];
   Show[p0,p1,p2, Graphics[Line[{{nps,-1},{nps,1}}]],
               PlotRange->{-0.8,1}, Ticks->None,
               PlotLabel->PaddedForm[t,{3,2}],
               DisplayFunction->$DisplayFunction]; ]

Evolution

First we communicate the initial condition and potential function to the external routine named schroed.

schroed[0, $PSI->psi0,$T->0,$VFUNC->potential];

And now we draw what schroed has just been told:

drawit[];

To animate this, simply select the outer most cell and type y

The frames that are being animated are the one above, plus a further 25 computed next.

Do[evolve[0.04]; drawit[], {25}]

• The rapid oscillations when the packet hits the potential wall are accurate – they are NOT numerical error or instability.

• In the animation, the potential wall is high enough so as to reflect most of the incoming wave-packet – although part of it does quantum tunnel through a bit.

• It's interesting to see a standing wave form inside the potential wall too, due to a resonance effect slowing its decay.

• Also wave-packets naturally spread out over time, as mentioned earlier. This too can be seen in the animation.

Miscellaneous

The Crank-Nicholson scheme guarantees that the solution is unitary throughout the evolution. That is the integral of all the probabilities,

Integrate[psi[x,t]*Conjugate[psi[x,t]], {x,-Infinity,Infinity}]

is 1 for all time no matter how poor the numerical resolution. Numerically psi[x,t] is represented by the array $PSI and the above integral just corresponds to a dot product multiplied by the grid spacing dx.

Peek[$PSI].Conjugate[Peek[$PSI]]*dx

References

InterCall

InterCall was developed by Terry Robb. For more information on InterCall send email to intercall_forward@wri.com. An order form is available via MathSource – email the line send 0202-587-0011 to mathsource@wri.com for details. Information can also be posted or faxed to Analytica below.

Analytica                     Fax: +61 (9) 386 5666
P.O. Box 522
Nedlands, WA 6009
AUSTRALIA.