JOHN
edlted by W MOORE
-
courses in quantnm chemistry
Numerical Solution of the One-Dimensional Schrodinger Equation Uldis Blukis and J. M. Howell College, Clty Unlverslty of New York, Brooklyn, NY 11210
Brooklyn
The one-dimensional Schrodinger equation may be written as
where the independent variable is a spatial coordinate X, p is the reduced mass, V(X) is the potential energy of a model, and the other symbols have their usual meanings. To solve numerically a differential equation of the form y" = F(x,y) one may use the simple but accurate Noumerov formula (1,2). When the function F(x,y) = f(x)y, as in eqn. (I), the Nonmerovformula can be written
understandine of the difference between an eiaenfunction and other solutions of the Schrodinger eqnation, as well as the effect of the nature of fix) (hence E, and V(x)) on the s h a ~ e and amplitude of eigenfunctions. Reduced Equations
For numerical solutions, reduced equations are a practical necessity since they contain a minimum of parameters characteristic of a molecule or atom and of univer&l constants and, therefore, minimize the number of numerical solutions that one needs to calculate. However, they have another advantage: they separate the universal properties of a given model from those of a particular molecule to which the model is being applied. The separation, purchased for the price of a slight increase in abstraction, can simplify learning about a given model. The reduction of an equation may be viewed as a problem in dimensional analysis: given a set of molecular parameters and of universal constants a, b , c, d , what values must the powers k , I , m , n have if E , i n E , = a k b l c m d d E(see Table 1, last column) be dimensionless? All one-dimensional Schrodinger equations may he transformed to the particular dimensionsless form y', = f(x)y = [-2E,
where x, - x,-1 = h. The integer n labels the successive points a t which y is computed, counted from the initial point a t n = 0. The residue R,+1 is approximately given by Compared to most formulas for numerically solving differential equations the Noumerov formula (eqn. (2)) has an advantage in application (3):f(x,+i) can be computed exactly, unlike F(x,+,,y,+~) which can be only estimated. Thus, iterative methods, as the commonly used predictor-corrector methods ( 4 ) ,are not necessary for solving eqn. (1). A simple program based on the Noumerov formula for any computer or for a pocket calculator can determine accurate eigenfunctions of eqn. (1)with a great variety of potentials V(X).We DruDose . . that this ca~ahilitvshould be verv useful for teaching introductory quantnm mechanics. In the remainder of this paper we discuss some advantages ( 5 )of the numerical approach to introductory quantum mechanics, present some specifics of using Noumerov's formula (dimensionless parameters and coordinates, initiation, errors, general properties of f(X)),and present some results. Advantages
The analytical approach to the one-dimensional Schrodinger eqnation requires that one learn various rather sophisticated methods for finding exact truncated series solutions, and for approximating solutions. The Noumerov formula is the basis of a sinele and s i m ~ l method e for finding accurate numerical approximations to both exact and ap-' oroximate analvtic solutions and eieenvalues. I t ~ e r m i t the s beg~nningstudent to generate hlsbwn eiaenvaiues, eigenfunctions, and other solutions rather than merely look at somebody else's results. Therefore, he should develop a keener
+ 2V(x)]y
(4)
where all variables are dimensionless or reduced. The unreduced and reduced potential energies, the reduced energy E,, and the reduced s ~ a t i acoordinate l x for unrednced equations of the form (1)are given in Table 1. Some one-dimensional equations in their usual form contain a first derivative. Such equations are the radial Schrodinger equations of models for the hydrogen atom or oneelectron atomic ions, for the diatomic vihrotor, and for the total angular momentum equation. A transformation of the dependent variable that has the general form y = g(x)$ (6) plus reduction (if necessary) will produce an eqnation of form (4). In Table 2, we list nntransformed and transformed equations in this group. Solutions of the Schrodinger Eauation
by the constraints imposed upon the solutions. Using the harmonic oscillator as our example we consider the effect of f(x) = -2E, + x2 on the solutions of the harmonic oscillator Schrodinger equation. We note that f(x) is negative when lx 1 < (2E,)1/2 and f(x) is positive when Ix 1 > (2E,)lI2. At x = f(2E,)'/2 the solution has its outermost inflection points. In an interval on x where f(x) is positive y"(-Ay'/Ax) and y have the same sign. This means that for positive y the slooe v it . "v' increases with increasing x while for negative decreases with increasing x. In sum, when f(x) is positive the solution y always curves away from the x-axis with increasing magnitude of x. On the other hand, when f(x) is negative the solution Y alwavs curves toward (or oscillates about) the xaxis. ~ h maghitude k of f(x).y(x) determines the rate of curvinr.. a t point x. Therefore, in the negative f(x) region near .
x = 11 ut.expet.1 i , . v i l l . ~ l ~ t ,(11 ~ ~111f - . S ~ ~ ~ I Iil0t)ol I I ~ ) U 1Iw I - i i \ ~ i that will III. I ht 111,~t.e ,.I illt.l\,ll~.~aed ~lle 1hinht.r I ht. \oiu? 01 L,,
but that will space out as-f(4) decreasestoward zero at the Volume 60
Number 3
March 1983
207
outermost inflection points. Within this region the solution has further inflection points whenever y (x) = 0, i.e., on the x-axis. In the region of positive f(x) the solution will curve away from the x-axis (Fig. 1). We now consider the constraints that must be imposed upon the solutions of the harmonic oscillator equation to generate eigenfunctions. The constraints must uniauelv . " determine the that characterize a particular eigenfunction. There are three such parameters: the two arbitrarv constants in the general solution of a second order differ-. ential equation that appear during the first step (n = 1) of the numerical integration as yo and y, in the Noumerov eqn. (2) and the reduced energy E,. In the prevailing, the Born interpretation, the parameterdetermining constraints arise through the association of the eigenfunctions with probability concepts. Namely, the probable distribution of harmonic oscillators along x is given by y*(x)y(x). Therefore, the total probability of finding the harmonic oscillators anywhere along x, given by integrating over the probable distribution, must be finite:
from the quadratic integrability constraint hut from the symmetry of the eigenfunction. The procedure for determining the eigenfunctions and eigenenergies of the harmonic oscillator is then as follows. An antisymmetric function must have y(0) = 0. If we start the numerical integration a t x = 0,then yo = y(0) = 0. For a symmetric function we can obtain a relation between y, and yo = y(0) via Taylor's expansion. Since the derivatives y(")(O) = 0 when n is odd and all even derivatives can he obtained by successive differentiations of the Schrodinger equation, the Taylor expansion becomes
Equation (5) requires that the eigenfunction y(x) be quadratically integrable. If we further require that the total probabilitv be eaual to one, then the eieenfunction is both parameters, on> of which is E,, whie normalization yields the remaining unknown parameter. With symmetric potentials, such as the harmonic oscillator potential, it is convenient to determine one parameter not Table 1.
v(X)
Model
-1
Reduced Oscillators X
-
Potentlai Constants
Equations With First Derivatives Transformed Equation
w.4
Total Angular Momentum
*"+
ctn
m2 m 'sin2x - *= J ( J +
I)'#
siniJ2x
-%
208
Journal of Chemical Education
Er -
E
X
Untransformed Equation
I
Figure 1. Numerical approximations (h = 0.1)to normalized harmonic oscillator eigenfunctions wilh low ( n = 0) and high ( n = 10)energy.
V(d
Table 2.
C
ctn 2~
I
+sin2x "' y
Reduced Parameters and Variables
Three numerical approximations (h tisymmetric harmonic oscillator solution. Figure 2.
= 0.1) to
the firsf
-1C ( n = 1) an-
Figure 4. Numerical approximations (h = O.l)to the ground state eigenfunctions for the quartic (E, = 0.66600. V(d = fl)and double-well quarlic (E, = 0.86952 0.00001, V(X) = 1 + f l - 2x2).
*
The eigenfunctions and eigenenergies of all symmetric potentials that do not possess singularities a t finite values of x, e.g., the quartic potential V = x4, or the quartic potential with a quadratic hump, V = x4 - ax2, can be determined by the ahnve nrocednre (Fies. 3 and 41. Potenriala uthcr than yrnmvtrwone-. s,ilhout .;~ngolnriti(.; fall into t w o ~lcl4st.iwit" rt.wect I < , I hc milnner i n n Ihich t hr ir numerical eigenfunctionsmay he determined: (1) unsymmetric ootentials without sineularities at finite x; (2) potentials with " singularities a t finite x. An examule of a potential in the first class is the Morse
~.
-1
t
Figure 3.Numerical approximations (h = 0.1) to the ground slate elgenfunctians for the harmonic (E, = 0.50000 i0.00005. V(x) = xZ/Z)andthe quartic (E, = 0.66800 f 0.00004, V(x) = fl) oscillators. of eigenfunctions for two potentials in this class are shown in Figures 6 and 7. T-. o the second classinotentials with sinrmlarities at finite ~--x-belong the H-like and the total angular momentum equations in Table 2. We shall discuss the determination of the eigenfunctions and eigenenergies of both equations. For the H-like eauation we have chosen to start the numerical integration arbitrarily close to its regular singular point at x = 0'. Therefore, for the first use of eqn. (2), we must , = f(O)y(O) = ~ " ( 0 1y~ , = know the values y o = ~ ( 0 1f(xO)yo y(x1). For any stationary state solution of the H-like Schrodinger equation general proof that ~~~~
[
=y,l
h4 h8 hzb.,+-(1t2E:)--(2E:t7Er)+ 12 180
.]
(61
etc. The first line of eqn. where fo = f(O),6 = [dzf(x)/d~2],=~, (6) is valid for any symmetric potential. The second line is the Taylor expansion for the harmonic oscillator. The requirement that symmetric or antisymmetric harmonic oscillator solutions he quadratically integrable will determine the energy parameter E,. Numerical solutions, due to errors inherent in numerical procedures, can do no more than approximate a quadratically integrable solution. The best numerical solutions are selected on grounds of their hehavior as they approach the limits of their range on x. The criterion, applied through visual inspection of best behavior, simply selects the harmonic oscillator solution that stays nearest to the x-axis for the longest distance as the limits are approached (Fig. 2). The third parameter, y l for antisymmetric solutions and, most conveniently, yo for symmetric solutions of the harmonic oscillator, determines the y-scale of the solution. During the determination of the other two parameters it may he set at any convenient value. When the other two parameters have been determined, the eigenfunction may he normalized by requiring the numerical integration of eqn. ( 5 ) to yield the value "one."
~
~
-
~
~
(71
~(0= 10
is possible. For a stationary state the expectation value of any linear momentum must equal zero. Thus, for the expectation value of the radial momentum we have, assuming real eigenfunctions,
'
in the subsequent proof, we follow the traditions of quantum mechanics texts in not distinguishing the point r = 0 from points arbitrarily close tor = 0.Since the solution of a differential equation is undefined at r = 0,a rigorous treatment would be unable to say anything about y(0). However, since the point at r = 0 is of measure zero, it will not contribute anything to the integrals on which our proof depends, nor will it affect the subsequent numerical integration.
Volume 60
Number 3
March 1983
209
Figure 6. A numerical approximation ( h = 0.1) to the normalized ground state eigenfunction (E, = 0.489625 f 0.000005) of the cubic anharmonic oscillator with V ( x ) = x2/2 - $112 Since Me potentiai diverges toward negative infinity 85 x goes toward positive infinity, the oscillatory part of the solution at large xthat diverges toward infinity (not shown in the figure) is not part of the eigenfunction.
Result E, = 0 4896252 0 000005 y - y o = (8 1715 +O 0 0 0 5 ) ~
0 4895 Figure 5. The determination of the numerical approximation. shown in Figure 6, to a cubic anharmonic osciiiatw eigenfunction: (a) First ten guesses or "shots:" (b) Last seven shots with scales expanded. The arrows connect successive guesses. numbered sequentially, of the best values for the eigenenergy, E,. and of thedifference between the first two values of the unnormaiized (yo = I ) eioenfunction. .vl - .vn. - me best vaiues are found at the intersection of the dashed ~ reg on 5 me IUO ~ncs,cil len I 1 p .nes, tn.2' 'iepilralc me four te9 o n w.&n 0, .,piw of erers me Irsl referrmg l o ih? l e n w o t dn llle n q a : r e r-ilk 5 and !he s ! ~ w alo 'he "ens$ or un inr poz ! ,e x - 3 r a oeflnru ny m e f0.r ~ J S S I DC, combinations of solution behavior on the minus x and plus xaxes. On the minus x-axis. the behavior of the interest is found where the solution begins to diverge as xapproaches negative infinity. Here the solution may either cross the xaxis after the last extremum (indicated by the first letter C in the region label) or not cross it (first letter N). Figure 4 shows an example of each behavior. On the plus x-axis. the crossing or non-crossing of interest is located at the boundary of the physically meaningful potential energy region, roughly wherethe potential has its finite maximum (see Fig. 6).
The upper limit of the integral equals zero because of quadratic integrability. Therefore lim ~ ~ (=00) r-0
and, hence, y(0) = 0 (cf. ref. (8))Q.E.D. To find y"(0) and yix), we use exact analytic solutions of simplified Schrijdinger equations that are valid near x = -, and x = 0. Near x = m, the H-like Schrodinger equation hecomes 210
Journal of Chemical Education
Figure 7. Numericalapproximations(h = 0.1) tothe narmaiiredfirst excited state eigenfunctions of the harmonic oscillator (E,= 1.50000 i 0.00001) and the Morse osciiiator ( E = 1.4550 0.001, y = 5).
*
Y"= -2E,y which, for anegative E, has the solutiony Near x = 0 for 1 > 0, we get
(9) = e~p[-(-2E,)'/~x].
which has the obvious solutions, y = x l + l and y = x - ! The second solution must be rejected if 1 > % because it is not quadratically integrable. The two limiting solutions suggest a solution that might he valid a t any x, namely, Y = zl+'exp(--4
(11)
Substitution of this solution into the H-like equation (Table 2) shows that eqn. (11)is an eigenfunction even for 1 = 0 if E, = -1/[2(1 + and thatyr'(0) is -2 if 1 = 0,2 if 1 = 1,a n d 0 if 1 > 1.In summary: We have found analytic eigenfunctions and the corresponding energies for one set (the ls,2p,3d, . . . set) of states; yo = 0, y l from eqn. (11) and the above values of yr'(0) = f(x0)yo can he used to find numerically the re-
Figure 8. Numerical approximations ( h = 0.1) to the normalized H-like eigenfunctions of the 2s (E, = 0.12496 f 0.00001: exact E, = -0.125) and the 3 p (E, = -0.3125) states.
maining eigenfunctions and the corresponding eigenenergies (Fie. ,- ~m 8). ~, The total angular momentum has an effective potential that is symmetric about x = K/2 and that results in a differential equation with regular singular points at the limits of the range x = 0,n. Therefore, the numerical solutions may be initiated from xo = ~ 1 in2the same manner as for the harmonic oscillator. Some of the resulting approximate solutions are shown in Figure 9. T o choose the best approximation to an eigenfunction. we note that eieenfunctions must ohev viO) = v i r i = 0. This follows from thkfact thaty*(O)y(O)aniy*(n);(nj are the ~rohahilitiesof findine the electron on the zero volume positive and negative z-axes, respectively (8).
1 Figure 9.Numerical approximations (h = 0.01)to the normalized total angular momentum eigenfunctions: (a) lml = 2 (assumed), J = 3.000 -+ 0.001: (b) rn = 0 (assumed). J = -0.2, 0, +0.2. + 0 5 Figure b shows that due totheir negative f(4.eigenfunctions with rn = 0 are hard todistinguish from othersolutions.As (mi increases toward 'ts maximum value, m( = J, k eigenvalue of J becomes easier to distinguish (as shown by Figure a).
1
Errors
Programs
The value of v , + ~ from the nth steD is affected hv the round-off error characteristic of the calc;lational device heing used. For many computers, this error may he diminished by going to double precision calculations. A second error that affects the results for step n is the truncation error-the neglect of the residue R in eqn. (2). A third error, the initiation error, occurs in YT - vn if the initiation uses an an~roximate
A program to perform the calculations described in this paper was written in the PASCAL language (UCSD version) for the Apple I1 microcomputer. Computations can be performed on polynomial, Morse, hydrogen-like, and effective angular momentum potentials. The program could be easily modified to accept other potentials. The user of the program must repeatedly input the unknown parameters which constrain the eieenfunction and search for their best values. The solutions touthe differential equations are graphically displayed either on a television monitor or on a Tektonix 4662 hard copy plotter (for which an assembly language interface routine was written). Conies of the Droeram on 5% in. f l o o ~ v diskette and documentation will heprovided if a blank di'kette is sent to the authors. We have also written for the Hewlett-Packard programmahle pocket calculator 19C (it has printout) programs (available from the authors) that produce eigenfunctions and eigenenergies of symmetric, H-like, and total angular momentum models. The speed of the programs is about 8 points on the solution curve per minute. Therefore, a program for use with the two-dimensional shooting method that is necessary for unsymmetric potentials would he impractically slow.
two errors can be decreased to values below the round-off error simply by decreasing h. A final error, unlike the previous ones, grows with the number of steps. This is the cumulative error-due to the round-off, initiation, and the truncation errors of individual steps the numerical solution gradually departs from the particular solution that was selected hy the initial conditions and by the guessed eigenvalue. The Schrodinger equation has solutions that. at small x values. are near the eieenfunctions and that contain exponential teirns vith p s i t i & exponents. The cumulative error will cause all numerical solutions to be dominated at large magnitudes of x by the positive exponent terms and, hence, to diverge. However, our examples show that the divergence is easily made to occur past the point where the exact solution has an apureciable amvlitude and. n is sufficiently large, the increase is exponential with step number n and with step size h (10). A final source of error may become important near singularities in f(x) if the denominator in eqn. (2) diverges from values near one. With our numerical method, this happens near the singular points of the total angular momentum equation. The H-like equation, however, is stable in our numerical procedure even near x = 0.
Literature Cited (1) Noumerov, B.,&bl. Cenlrol Artrophys. Obs. Ru~sio,voi.7,1923, Moscow;Royol as^ t i o n Soc., M.N., 84,592 (1924): Scheid, F., "Theory and Prohlems of Numerical Analysis.). MeCraw-Hill,New York,1968, pp. 224,212-233. (2) Boleman, J. S..Amer J. Phys.40.1511 (1872). User Noomernvformula tosolve the
1970, pp. 130rI. A delightful and practical int;~duefionto n;merical methods. (51 Sherwin, C. S., "lntrodoction tn Quantum Mechmcis," Holt, Rinehart and Winston. New York, 1969, p. 20. Sherwin uses numerical methods to teach the quantum af the hydrogen atom. " ~ e d ~ ~or~ the i < harmonicoscillatdrsnd s
Volume 60
Number 3
March 1983
211
(6) Margenau, H., and Murphy, G. M.."The Mathematicsof Chemistryand Physics: D. Van Nostrand.Princeton. NJ.1956.2nd Ed. ( I ) Daniel, d. W., and Msure. R. E., "Computation and Theory in Ordinary Differential Equations." W. H. heeman. San Frsnciscu. 1970. (8) Messiah, A,, "Quantum Mechanics: John Wiley, New York, 1961, vol. I, p. 346. : ~Or~w~Hlll, (9) Pading, L., and Wibon. E.B., ~'Intttdddtiiit o Quantum M ~ ~ h h i i i M NewYork. 1935.p~.59-63. (10) Hildehiand, F. B., "introduction to Numericai Analysis: McGiaw~Hill,New York. 1956. Formulae6.12.18and6.12.19.SseslsoCuuley, J. W., Math. Cornput. 15,361 (1961).
Illustration of Quantizationand Perturbation Theory Using Microcomputers Christian Kubach Laboretoire des Collisions Atomiques et Moleculaires, Bit. 351, UniversitC de Paris Sud, 91405 Orsay, France Microcomputers allow students to become more familiar with the concepts of quantum mechanics and also to ro heyond prohlemshavin&alytic solutions. The purposesbf the present contribution are to exhibit the notion of quantization and to reinforce students' knowledge of the behavior of wavefunctions and the results of perturbation theory. Determination of vibrational energy levels of diatomic molecules, along with the associated wavefunctions, is of particular interest to chemists. This involves solution of a one-dimensional Schrodinger equation, which can he integrated either by linear variation or by numerical methods (see e.g., (1) and (2)). The Noumerov algorithm, which pertains to numerical methods (see e.g., ( 3 ) )has been selected here for its simplicity. Outline of the Noumerov Algorithm The one-dimensional Schrodinger equation may be reduced in suitable units to:
In the Noumerov algorithm the solution of eqn. (1)is calcnlated for equally-spaced values of x defined by r,+l-x,=x,-x,-l=...=h
(2)
Using a Taylor expansion off through order hS one obtains: where fk = f(xh ). Differentiating eqn. (3) twice gives fl, since: h2(f,+,+ fA-l) = 2h2< + h4%+ 0(h6)
(4)
Finally using eqn. (1) one obtains:
Equation (6) allows the determination of f,+l if the solution is known for the two previous points: x, and x , _ ~ . Application to the Potential Well For hound states of a potential well there are two classically forbidden regions defined by V(x) > E. The physical solutions of the Schrodinger equation vanish at points far into these regions. This property is used to begin the construction off a t a point sufficiently far into the first classically forbidden
' The solution of the Schrodinger equation is not normalized, and so may have any value compatible with the precision of the microcomputer used. f
212
Journal of Chemical Education
region that f may be chosen to have zero value. More preciselv.
bidden region. Selection of Parameters The harmonic potential (V(x) = x2),along with the known eigenvalues E, = 2(u + 1/2), u = 0 , l . . ., is used both to determine the integration parameters (xo,h . . . ) and to be a first exercise for the students. The following energy-dependent parameters have been selected to ensure sufficient precision for the construction of the lowest energy functions of the harmonic oscillator (v < 6) within a small computational time.
Eauation (7) shows that a leneth of 2 units of x is chosen as a starting point in the first c l a k a l l y forbidden region, (the quantity 2 a represents the length of the classicdy allowed region). Equation (8) gives the point X M where the construction of the fnnction f is achieved and is checked for convergence (the points xo a n d x are ~ symmetric with respect to the location of the potential minimum). Finally, eqn. (9) gives, for a fnnction having N nodes, 10 (N+ 1) integration points for the determination of the solution in the classical allowed region (since E, = 2 (u 112) for a wavefunction with u nodes). Other scaling parameters have also been determined in order to ensure the visualization of both the function and the potential on the screen of a Radio Shacka TRS-80 for a large range of input energy 0.1 < E < 15. The parameters of eqns. (7)-(9) can also he used for nonharmonic potentials, provided the ootentials chosen are nearlv harmonic.
+
" .
-"
E defining the integration parameters eqns. (7)-(9). The results aw~earinr .. - on the screen show. for each calculated value of x , the solution f(x) of eqn. (11, r h potential ~ V(x) and the location of the classicallv allowed reeion. In eeneral the mathematical solution f isdivergirig in the seconl classically forbidden region (x 5 XM). . .. excewt for articular values E,, (eigenvalues) where wavefunctions converge in this region. The eigenvalues may be bracketed using the relation between the signs of the divergence and of E - E, (see Fig. 1). Problems Given to the Students The harmonic wotential is of articular interest since it is a first approximation to the lowest vibrational levels of a diatomic molecule, and since its eieenvalues and behavior of eigenfunctions are usually knownby the students. With this potential the user observes a clear difference between physical (non-diverging) and only mathematical (diverging) solutions of eqn. (I), illustrating the quantization concept (see Fig. 1). One also visualizes the lowest energy wavefnnctions and verifies the law E, = 2(u 112) with a precision better than n 79" The addition of a perturbation term to the harmonic oscillator offers two benefits as an instructional exercise: i t allows the student to review perturbation theory, and it provides a closer approximation to the potential energy curves of diatomic molecules. A comparison between the results obtained for the two following potentials:
+
