2 Documentation

 The Schrödinger applet provides a numerical solution of the time-independent one-dimensional Schrödinger equation,

  

$${\frac{\hbar^2}{2m}} {\frac{d^2\,}{dx^2}} \psi \psi \psi$$ (1)

The wave function is set to zero at the ends of the range to define the boundary conditions.

2.1 Panel Documentation

 The problem to be solved is specified by the contents of the seven panel windows on the left side of the applet. Each can be changed by the user.

  The input fields (from top to bottom) are:

• the input value of ${\frac{\hbar^2}{2m}}$ .
• the minimum and maximum values of x which give the range of the wave function; the wave function is zero at X_min and X_max .
• the number of grid points used by the differential equation integrator, uniformly distributed over the range.
• the number of nodes (zeros) in the solution, i.e. the excited state number.
• the scaling factor used to plot the potential and the energy.
• the analytic potential.

The analytic potential can be quite complicated or can be a simple expression such as 1/2*x**2. The rules for parsing (translating to computer code) the potential are given in Sec. 2.2.  

The floating point numbers ${\frac{\hbar^2}{2m}}$ , X_min , X_max , the energy scaling factor, and any numbers in the potential can be specified either in fixed point (e.g. 1.0 or 7.678, etc.) or by using scientific notation as in Fortran, (e.g. 1.0e-5, or 1.0e4, etc., where the integer after the e is the exponent).  

While the user may specify any input field, the applet also provides an ``Example Inputs'' bar. Clicking the bar, located at the top of the right side, will produce a menu of sample inputs to be selected. Choosing an item from the list will write a set of values in the input boxes that will allow you to calculate that system by clicking the Solve button. This list merely provides a survey of the types of problems which the user can examine and some examples for the potential parser. Go to Sec. 3 to explore the examples in detail.

Next to the Solve button is a button marked Stop   normally this button will be disabled and will be ``greyed out.'' When you click Solve, the Stop button will be enabled and by clicking it, you can stop the calculation. This can be useful if you have inadvertently typed in input data that will cause the integrator to run for a long time, e.g. by using an enormous number of grid points.

Below the plotting area are three buttons labeled Print,   Read,   and Write.   The Print button will bring up a dialog that will allow you     to print the current applet appearance. The preferred method of printing is to use the Browser's print function for those Browsers that support printing applets (e.g. Netscape 4.05 preview release 1 (AWT1.1.5)). This applet Print button will allow printing on some other browsers, and when the applet is being run as a Java application or in a separate Window. The quality of the printing is platform dependent. Printing the applet in Landscape orientation normally provides the best quality.

The Read   and Write buttons   bring up file dialogs. Selecting Write and giving a file name will write a file with the values of the 7 input fields of the applet, the calculated energy, followed by the values of x , the wave function, and the potential. The file is a standard ASCII text file.   Here is the start of a file generated for a double harmonic oscillator:

# 0.5  hbar^2/(2m)
# -10.0 xmin
# 10.0 xmax
# 1000 ngrid
# 0 nnode 
# 0.1 vscale
# 0.49977564535774377  Energy
#a = 3 
#0.5*min((x-a)^2,(x+a)^2)
-10.0 0.0 0.0
-9.97997997997998 6.367057436351622E-12 24.36006026046066
-9.95995995995996 1.2855968537246738E-11 24.220521322122924
-9.93993993993994 1.9589595863550545E-11 24.081383184986787
   .
   .
   .

The first 6 lines contain the first 6 input fields from the applet, the next line is the computed energy, followed by the potential definition. All of these lines begin with a # character. The remaining lines are x , $\psi$(x) , and v(x) , at the grid points. The plotting package gnuplot views the # character as a comment, so these files can be directly plotted using gnuplot. These files allow users to save cases of interest, with numerical tables of wave functions that can be used to calculate matrix elements etc. in other programs.  

Clicking the Read button,       and then selecting a previously saved file, will read the 7 input fields from the file. You need to click Solve to recalculate and display the results of this input.

These three buttons are the only parts of the applet that use Java 1.1 features and will require you to use one of the latest browsers. Early browsers such as Netscape 2.0 or 3.0 will run the rest of the applet, but will not allow printing or file manipulation. For more information, see section 2.7

2.2 Parser Function

   The applet parser allows users to define their own potential functions for which the applet will solve the Schrödinger equation. By clicking in the panel window labeled V(x) = , the user can type one or more lines of assignment statements in the style of Fortran or C. For example, the parser can recognize expressions like the set

L=0
D0=3145
alpha=2.352709
r0=2.413683
ex=exp(-alpha*(x-r0))
D0*ex*(ex-2.)+.96415129*L*(L+1)/x^2

which defines a Morse potential for a pair of Cl atoms.

The operators that are recognized are + , - , * , / , ^ , and ** . Both ^ and ** can be used as the exponent operator. In addition, the functions sin, cos, exp, log, max, min, atan, sqrt, abs, H, and pulse are recognized.

The max, min and pulse functions have 2 arguments, and the rest take just one. H is the Heaviside function; H(y) is zero when y is negative and 1 when y is positive. The pulse(y,z) function is 1 when |y| < |z| and zero otherwise.

Since the potential panel is editable, as are the other input panels,   the arrow, insert, delete, home, end, etc. keys, as well as text selection with the mouse, provide simple editing features.

 

2.3 Plotter

   The middle right of the applet panel is the plotting area which, upon successful completion of a solution, displays the wave function in red and the potential in blue with a green horizontal line at the energy value. (Note that the wave function is scaled so that its maximum magnitude is 1.0 , and is not normalized so that its square integral is unity. Recall that the energy values shown need to be multiplied by the V_scale value that was typed in.)

The plotting routines are from     ptplot by Edward A. Lee and Christopher Hylands, copyright University of California.

When the results have been displayed, the user can zoom,   by using the mouse to click and pull down and to the right to produce a small square of the plot, which will be expanded into the full plotting region. You can do this repeatedly. This feature is extremely useful for viewing fine details of a wave function. You can restore the plot to its original form by clicking Solve again.  

2.4 Natural Units

  Natural units are exploited in several of the example inputs. One way to extend the range of numerical solutions is to use dimensional analysis. In the harmonic oscillator Hamiltonian,

$${\frac{\hbar^2}{2m}} {\frac{d^2}{dx^2}} {\textstyle\frac{1}{2}} \omega^{2}_{}$$ (2)

we have three parameters, $\hbar$ , m , and $\omega$ . Nature tells us the value of $\hbar$ in Joule-seconds. However, we aren't required to use SI as our measuring system. A more ``natural'' set of units can be found by looking at the dimensions of our three quantities: $\hbar$ has units of mass length ² time ^{- 1}, m has units of mass, and $\omega$ has units of time ^{- 1}. For any particular harmonic oscillator, we can choose our natural units to be
 

$$\sqrt{ \frac{\hbar}{m \omega} }$$ Length Unit =

Mass Unit = m

$${\textstyle\frac{1}{\omega}}$$ Time Unit = (3)

The numerical values of $\hbar$ , m and $\omega$ in natural units are all equal to 1 ( m is 1 natural mass unit. $\omega$ is 1 natural inverse time unit, etc.). We can therefore do the calculations in these units, and at the end convert to physical units. For example, we will find the ground state energy of the harmonic oscillator is 0.5 in our natural units. A natural energy unit is

 

$$\hbar \omega$$ (4)

so the value of the ground state energy is 0.5$\hbar$$\omega$ as usual. In natural units, the Harmonic Oscillator Schrödinger Equation is

 

$$\left [ -\frac{1}{2}\frac{d^2}{dx^2}+ \frac{1}{2} x^2 \right]\psi \psi_{n}^{}$$ (5)

which has no other free parameters. Therefore all harmonic oscillators are described by this equation.  

2.5 Influence of Boundary Conditions

  If you solve the above harmonic oscillator you will find that as you solve for higher and higher excited states, the wave function has a larger extent in x. For the first 20 or so states, you will find that you can enforce boundary conditions that $\psi_{n}^{}$(x) = 0 at x = $\pm$ 10 without changing these solutions. Note that this boundary condition introduces a new parameter, the distance where we set the wave function to zero. You will get different answers for different values if you increase n enough.  

2.6 Numerical Errors

    As described in the solution algorithm section, the differential equation integration algorithm uses a finite difference formula for computing the wave function. If the interval is divided into N steps, the error in integration is of order N ^{- 5} if the potential is well behaved, i.e., if the potential can be expanded in a Taylor series. Obviously discontinuous potentials such as finite square wells or steps will produce larger errors. The recommended way to estimate the step size error is to repeat the calculation with, for example, double and half the number of integration points to make sure that the values have converged to an acceptable error.

Similar tests can be done to make sure that the boundaries are not influencing the results for cases where the original problem's solution extends to infinity.  

2.7 Error Messages and Failure Modes

   The applet is relatively stable, but there are ways to break it.

• Parser Errors  

  The parser, Sec. 2.2, accepts simple algebraic expressions and requires a fairly rigid syntax. If there is an input or parsing error, a message will be written to the energy panel.

• Failure to Find a Solution

    The numerical integration may fail because its searching procedure (integrating inward from the two end points) cannot match the the solutions at the midpoint. If the integration fails, an ``Integration failed'' message will be written to the energy panel. Generally this is caused by X_min and X_max chosen too far apart so that the wave function overflows or underflows during the calculation. It can also be caused by too few grid points. Adjusting the grid points or range specified for the solution usually fixes the problem unless a pathological potential is being solved.  

• Java Errors

    Java is an evolving language. We have tested the applet using various versions of sun's appletviewer, Netscape and Internet Explorer. In addition the applet can be run as an application if you have a Java runtime environment installed. For full functionality, the applet requires Java version 1.1. The Print, Read, and Write buttons will not work under version 1.0, but the rest of the applet should. Browsers that support the Java 1.1 features we use include Internet Explorer 4.0 and Netscape 4.05 preview release 1 (AWT1.1.5). The standard Netscape 4.05 will not work properly. It must support Java 1.1, so you must use one of the preview releases that have JDK1.1 installed. Currently, these are freely available at http://developer.netscape.com/software/jdk/download.html To reiterate, Browsers that don't or only partially support Java 1.1 may not recognize the print, read and write functions, but everything else will likely work fine.  

2.8 Solution Algorithm

   The applet solves the Schrödinger equation numerically on a discrete grid. The interval between X_max and X_min is divided into N - 1 equal intervals of width $\Delta$ = (X_max - X_min)/(N - 1) and, as explained below, the Numerov method is used to convert the differential equation into a difference equation for the values of the wave function on each grid point.

The method begins by rewriting the Schrödinger equation as an expression for the second derivative in terms of the value of the wave function,

 

$$\psi {\frac{2m}{\hbar^2}} \left [ V(x) - E \right]\psi \equiv \tilde{v} \psi$$ (6)

Here $\psi$'' is the second derivative of $\psi$ , and the function $\tilde{v}$(x) is defined to simplify notation. The Numerov method uses this form to integrate numerically the second derivative starting from the known value of $\Psi$ at a boundary and using a guessed value of E .

The Numerov method is derived by writing the Taylor series expansions

 

$$\psi \pm \Delta \psi \pm \Delta \psi {\frac{\Delta ^2}{2}} \psi \pm {\frac{\Delta ^3}{6}} \psi {\frac{\Delta ^4}{24}} \psi \pm \Delta^{5}_{}$$ (7)

and

 

$$\psi \pm \Delta \psi \pm \Delta \psi {\frac{\Delta ^2}{2}} \psi \pm \Delta^{3}_{}$$ (8)

The two equations defined by taking opposite signs in the second equation above can be added and solved for $\psi$'''' to give the usual second difference formula,

 

$$\psi {\frac{\psi''(x+\Delta )-2\psi''(x)+\psi''(x-\Delta )}{\Delta ^2}} \Delta^{2}_{}$$ (9)

Notice that all the odd order terms in $\Delta$ cancel so the error is O($\Delta^{2}_{}$) . Adding the two equations defined by taking opposite signs in the first equation and substituting for the fourth derivative gives the result

 

$$\psi \Delta \psi \Delta \psi \Delta^{2}_{} \psi {\frac{\Delta ^2}{12}} \left [ \psi''(x+\Delta ) - 2 \psi''(x) + \psi''(x-\Delta ) \right].$$ (10)

We haven't used the differential equation yet. Substituting Eq. 6 for $\psi$'' gives the Numerov three term recursion relation

 

$$\psi \Delta {\frac{\Delta ^2}{12}} \tilde{v} \Delta \psi \Delta {\frac{\Delta ^2}{12}} \tilde{v} \Delta \psi {\frac{10 \Delta ^2}{12}} \tilde{v}$$ (11)

The error in this expression is O($\Delta^{6}_{}$) . By starting with the values of $\psi$ at two adjacent grid points, with a guessed value of E , the Numerov recursion relation produces the wave function at all the grid points.  

The code enforces the boundary conditions that $\psi$ = 0 at X_min and X_max .   You can think of the system as being placed in an infinite square well potential that forces the wave function to zero at these points. The values are chosen by you. Other input that you need to provide are the value of ${\frac{\hbar^2}{2m}}$ , the number of grid points which sets the step size $\Delta$ , the number of nodes (zero crossings) that you want in the wave function, which determines the corresponding eigenstate, and the potential.

Since $\psi$ = 0 at the end points, the value at the neighbor grid point is arbitrary since it varies with the normalization. The code uses zero at the end points and an arbitrary small value at the neighboring points to start the Numerov recursion relation to generate the value of the wave function everywhere. The code integrates from X_min to a point near the center, and from X_max to just beyond the same point. The logarithmic derivative from the left and from the right at the match point are approximately calculated. These use the same Taylor series expressions above, but you subtract rather than add them to get a formula for the derivative. If the E value initially guessed is exactly right the two approximate derivatives from the left and right will be equal, and there will be exactly the correct number of nodes. If the energy is too high, there will either be too many nodes found, or the correct number of nodes and the difference of the approximate logarithmic derivatives will be positive. If the energy is too low, there will be either too few nodes, or the correct number of nodes and the difference in the approximate logarithmic derivatives will be negative. Therefore, we always know if the guessed energy is too high or too low. The integrator starts out with an energy of approximately the energy of an infinite square well from X_min to X_max with the same number of nodes as the solution we want, and doubles this energy until it is too high. It then uses an energy that is the negative of the infinite square well energy and doubles it until it is too low. The actual starting values are unimportant as long as the order of magnitude is sensible. Once the correct energy is bracketed, a binary chop is used to home in on the correct energy. The average of the upper and lower bounds is used, and after integrating this value replaces either the old upper or lower bound. This method by itself gives 1 binary place per iteration, so typically around 50 iterations would give machine accuracy. To speed this up, when both the upper and lower bounds have the correct number of nodes, Newton's method is used to predict the new energy. Typically the energy converges within 10 or 20 iterations.