NATIONAL BUREAU OF STANDARDS REPORT

A STATES TECAL METHOD FOR DETERMINING THE LOWEST

EIGENVALUE OF SCHRODINGER 1 S EQUATION

by

Mark Kac

Cornell University

and

National Bureau of Standards

and

Michael Cohen

Cornell University

and

National Bureau of Standards

U. S. DEPARTMENT OF COMMERCE

NATIONAL BUREAU OF STANDARDS

1#3
U. S. DEPARTMENT OF COMMERCE

Charles Sawyer, Secretary

NATIONAL BUREAU OF STANDARDS

A. V. Astin, Director THE NATIONAL BUREAU OF STANDARDS

The scope of activities of the National Bureau of Standards is suggested in the following listing of

the divisions and sections engaged in technical work. In general, each section is engaged in specialized research, development, and engineering in the field indicated by its title. A brief description

of the activities, and of the resultant reports and publications, appears on the inside of the back

cover of this report.

Electricity. Resistance Measurements. Inductance and Capacitance. Electrical Instruments.

Magnetic Measurements. Applied Electricity. Electrochemistry.

Optics and Metrology. Photometry and Colorimetry. Optical Instruments. Photographic

Technology. Length. Gage.

Heat and Power. Temperature Measurements. Thermodynamics. Cryogenics. Engines and

Lubrication. Engine Fuels. Cryogenic Engineering.

Atomic and Radiation Physics. Spectroscopy. Radiometry. Mass Spectrometry. Solid State

Physics. Electron Physics. Atomic Physics. Neutron Measurements. Infrared Spectroscopy.

Nuclear Physics. Radioactivity. X-Rays. Betatron. Nucleonic Instrumentation. Radiological Equipment. Atomic Energy Commission Instruments Branch.

Chemistry. Organic Coatings. Surface Chemistry. Organic Chemistry. Analytical Chemistry.

Inorganic Chemistry. Electrodeposition. Gas Chemistry. Physical Chemistry. Thermochemistry. Spectrochemistry. Pure Substances.

Mechanics. Sound. Mechanical Instruments. Aerodynamics. Engineering Mechanics. Hydraulics. Mass. Capacity, Density, and Fluid Meters.

Organic and Fibrous Materials. Rubber. Textiles. Paper. Leather. Testing and Specifications. Polymer Structure. Organic Plastics. Dental Research.

Metallurgy. Thermal Metallurgy. Chemical Metallurgy. Mechanical Metallurgy. Corrosion.

Mineral Products. Porcelain and Pottery. Glass. Refractories. Enameled Metals. Concreting Materials. Constitution and Microstructure. Chemistry of Mineral Products.

Building Technology. Structural Engineering. Fire Protection. Heating and Air Conditioning. Floor, Roof, and Wall Coverings. Codes and Specifications.

Applied Mathematics. Numerical Analysis. Computation. Statistical Engineering. Machine

Development.

Electronics. Engineering Electronics. Electron Tubes. Electronic Computers. Electronic

Instrumentation.

Radio Propagation. Upper Atmosphere Research. Ionospheric Research. Regular Propagation Services. Frequency Utilization Research. Tropospheric Propagation Research. High

Frequency Standards. Microwave Standards.

Ordnance Development. These three divisions are engaged in a broad program of research

Electromechanical Ordnance, and development in advanced ordnance. Activities include

Electronic Ordnance. basic and applied research, engineering, pilot production, field testing, and evaluation of a wide variety of ordnance materiel. Special skills and facilities of other

NBS divisions also contribute to this program. The activity is sponsored by the Department of

Defense.

Missile Development. Missile research and development: engineering, dynamics, intelligence,

instrumentation, evaluation. Combustion in jet engines. These activities are sponsored by the

Department of Defense. • Office of Basic Instrumentation • Office of Weights and Measures.
NATIONAL BUREAU OF STANDARDS REPORT

NBS PROJECT NBS REPORT

1101-11-5100 March 17, 1952 1553

A STATISTICAL METHOD FOR DETERMINING THE LOWEST

EIGENVALUE OF SCHRODINGER ' S EQUATION*

tflT Mark Kac

Cornell University

• and

National Bureau of Standards

and

Michael Cohen

Cornell University

and

National Bureau of Standards

PREPRINT

*The preparation of this paper was sponsored (in part) by the

Office of Naval Research

The publication, reprlr

unless permission is ob 25, D, C, Such permls

cally prepared If that

Approved for public release by the

Director of the National Institute of

Standards and Technology (NIST) on October 9, 2015

part, Is prohibited

ndards, Washington

ort has been speclfi-

>rt for Its own use.

Table of Contents

Page

I . Summary 1

II. Heuristic outline of the theory 1

III. Estimation of the error involved in the method 9

IV. Extension to the hydrogen atom * • 16

V. The slope method; results for hydrogen, helium, and

n double hydrogen” , . 22

VI. Techniques for speeding the sampling process; high speed

computing and importance sampling 3k

Appendices ....... ••

As Computation of discretization error for the harmonic

oscillator U6

Bs Further statistical properties of the random walks * * * 52

C« Statistical method for finding the lowest eigenfunction . 58

Bibliography 61
.

- ... ...

£ • ; I

-

OQ !
NBS 1553

3-17-52

PREPRINT COPI

A Statistical Method for Detennining the Lowest

•• ^

Eigenvalue of Schrodinger ! s Equation

by

Mark Kac

Cornell University

and

National Bureau of Standards

and

Michael Cohen

Cornell University

and

National Bureau of Standards

I. SUMMARY

The experiments discussed here are a continuation of those reported by

Donsker and Kac in [2], Atomic potentials, mainly that of hydrogen, are

dealt with here. Improved techniques for accumulating and treating the data

are presented; coupled with a high speed digital computer, these techniques

should make possible the rapid determination, to within a few percent, of

the lowest eigenvalue of many potentials# Considerable attention is given

to the errors in the method, and an elementary account of the underlying

theory is included.

II. HEURISTIC OUTLINE OF THE THEORT

With the introduction of an imaginary time variable, the time-dependent Schrodinger equation for a single particle in a one-dimensional

The preparation of this paper was sponsored (in part) by the Office

of Naval Research.

Now at the California Institute of Technology
r

: >

.

.

.

' :

• ' ' '• "

' < "

'

.

*

' ' '

' '

'

' :

' ’

'

'

'

r .
- 2 - potential field V(x) can be brought into the form

(1) — = 1 2 - v(x)P(x,t) . 3

1

^

Dtl

We shall shortly see that this form of the equation admits of a simple

probabilistic interpretation. There is some instructive value in the probability approach, inasmuch as it adds to our feeling for the meaning of

the wave-function and the equation. Furthermore, the probability model for

equation (l) yields a new computational method for finding the lowest eigen*

»• value of the time-independent Schrodinger equation. We shall concentrate

our attention mainly on this computational method.

In order to find the probabilistic meaning of (l) , let us consider the

problem of a single particle performing a random walk in one dimension subject to the following conditions:

A) At equispaced time intervals TJ the particle takes a step of

length h either to the left or the right, with equal probabilities „

B) If the particle stays at a given point x during a time interval dt, it is subject to destruction there with probability

V(x)dt + o(dt) . (We now restrict ourselves to the case

V(x) * 0, but shall later drop the restriction.)

If we let hP(x,t) represent the probability that the particle is at

the point x at time t (P(x,t) is thus a probability "density"), clearly

?(x,t) must satisfy the difference equation
•

•

.

•

'

.

.

.
- 3 - P(x,t) - |P(x - h, t -X)(l - 7(x - h)X + o(C))

+ fP(x + h, t - X)(l - v(x + h)t + o(-C))

(We have merely decomposed the event of being at the point x at time t

into the two alternative ways in which it could have occurred*) Subtracting P(x,t - X) from both sides and dividing throughout by X} we obtain

P(xJ t)-P(x,t-'C) _ 1 P (

x

-h . t-X) -2P (x . t~X)+P ( h . t-ta

X 2 % —

- |

P(x - h, t - X)v(x - h) - i P(x + h, t - r)7(x + h) + SEEL .

<- 1?

The left side, of course, is a finite difference approximation to If

we now require that the space and time increments in the random walk be related by the equation h

2

t;

f then the first term on the right is a second

3

2

P difference approximation to * . Then, in the limit as h 0 (and hence ^ yi ^ °)

9

the difference equation approaches the differential equation (l)

.

Thus, equation (l) is the limiting form of the difference equation

describing a random walk with destruction, where P(x,t)dx is the probability

that the particle is between x and x + dx at time t. Equivalently, the

equation describes a diffusion process with destruction, where P(x,t) is

the particle density at point x at time t, and V(x)P(x,t)dx is the time

rate of destruction of particles in the interval (x x + dx)

.

Under the former interpretation, if we start the particle from the

origin, we must require that our solution of equation (l) satisfy the subsidiary condition
. I . • • I

*' I I I I I I 11 1 . X J . ' • • • , , • •••. - 2 - r r; ' ' • ' • ' ' H I I H III . .
- h

-

(la)

00

for all t a o

, j

P(x,t)dx * 1

-oo

and the initial condition

As t 0

+

, P(x,t) -4 5(x) . (lb)

00

(£>(x) = 0 if x / 0 . / 8 (x)dx = 1) -oo

If V(x) oo as x —> - oo (we shall later relax this restriction),

the solution of (l) can be expanded in a series

oo “ X.t

P(x,t) = .2 c.e 3 ^-Cx)

0=1 1

1 0

where the X (>o) and are the eigenvalues and normalized eigenfunctions of

2

(2) - Av|/(x) = \ - v(x) yyu)

dx

subject to the conditions

oo / y (x)dxl < 05 , / -oooo

> -00

Y (*) 'dx < oo

and the c. are so chosen as to satisfy the initial condition lb. (If we

0

look for a solution of (l) in the fom of a superposition of functions of

-At x the type ><p

r (x)T(t) ? we find that T(t) * e and that the space-dependent part vj/" satisfies equation (2).) The eigenfunctions

’vp

r ;(x) form

J

an orthonormal set, and the initial condition (lb) yields c. = ^.(0).

, 0 0
" <

-

. .

.

.

.

-

•

'

’
- 5 -

•We note here that (2) is Schrodinger 1 s equation for stationary states,

where A represents E, the energy of the state.

Returning to the limiting form of the random walk problem, we now have

the result that the probability of survival of the particle till time t is

oo -A .t oo oo - A .t

(3) /

P(x,t)dx - ? e

0 Yi(°) f Aj/ (x)dx s 2 b.e 0 4o 0 4x> J a*i 3

If we arrange the A f s in order of increasing magnitude, the first

term of the series will dominate as t —> oo and hence

oo

(a A = - lira

1 t->oo

log f P(x,t)dx

-00

We shall employ the original discrete random walk to obtain an estioo mate of J* P(x,t)dx for large finite t. Then equation (U) enables us to

-oo estimate . Returning to the discrete random walk, we inquire: what is the probability p(x,9) that a particle survives the time interval 9(9 ^ X) , given

that the particle spends the entire interval at the point x? It follows

directly from our description of the random walk that this probability

must satisfy the differential equation

3p(*,9)

Te V(x)p(x,9) .

Hence

p(x,e)

•v(x)e

Now, suppose a particle performs the discrete random walk, starting

from the origin and taking steps of length h * -7as at time intervals

vn
. .

'

-

'

'

" •

*
- 6 - X = 1. Given that the particle follows a particular path x(kt), the probability that it survives till time t(=nt^) on this path is

nt-1

(Ua)

I I" p(x(k*),X) = exp < - k|Q

V(x(kZ))Z

k=o

nt-1

exp

nt-I -

k*o%(*(!£) )ll . n n

| From the definition of the random walk,9 it follows that x(-) v * S, v - , n' k n

, where is the sum of k random variables X-^Xg, ••• each X.. taking

the value 1 or -1 with equal probabilities.

We have computed the probability of survival till time t on a particular path. The probability of survival till time t, irrespective of path*

is the average of (Ua) taken over all 2 possible paths x(kT? . Since

all paths are equally likely, the average we are interested in is the simple

arithmetic average . We symbolize this average by

(lib)

Then, for large n, we

(lie) <^exp

Thus, by equations (U) and (l±c), we can estimate by

( 5) X ~

-log (ex? {-I At7^))

It is not necessary actually to compute the average (Ub) over all

possible paths. We can get a statistical estimate of the average by
»

, V.!;v- ;

. Wl M jfi **•'( T*® fl »vi

.

'
- 7 -

considering the average over the paths m independent random walks.

In equation (?) we have now derived the main result which we shall

use for computational purposes. The next few paragraphs are intended to

render precise some of our remarks in later sections,^

It follows from the central limit theorem of probability that in the

g limit as n 4 oo and k -4 oo, the distribution of D

k/\/n becomes Gaussian

k with mean 0 and variance - ,

i.e. Prob

{

x <^ <x + dx|

V2lf k/n

exp

2 —

x

2k/n f dx

Then, as n 4 oo , the distribution function of the random variable

i s r

n k<nt °k/v^i) approaches that of the random variable J V(x(‘C))d‘£* here

the path x(T) is chosen from a space in which the measure of a set of paths

is the probability that a particle performing a Gaussian random walk will

follow a path in that set, (By a Gaussian random walk we mean a continuous

random walk, with independent increments, with the displacement x(£) at

time % governed by the law

(5a)

Prob y < x(T) < y + dy j- - -/ - y

2/2"Cdy

J.^>

.

V(

Sk//n) <exp - -J V(x(C))dV

Thus, as n -4 oo ,

where the latter average represents the average value of

exp ^

- vfk V(x(C“))dE^, taken over the space of continuous paths with the

^Tor the rigorous development of the theory, see [1],
. , * . *
- 8 -

measure described above. This average is just the Wiener integral of the

functional exp

{

- I

v(x(t))drj. , and the left side of (5a) msy be thought

of as the analogue of a Riemann sum approximating the Wiener integral. By

(Uc) and (5a)

t

(6)

|

- J* V(x('C))dt Y>'C P(x,t)dx

For the purpose of computing A^, it would clearly be better to get a sampling estimate of the right side of equation (5a) than of the left side.

However, in our actual sampling process we cannot deal with the space of

continuous paths xtC) . We can, nevertheless, speed the convergence of the

distribution of - k<nt

^(°k/\/n) to that of V(x(tO)(TE by taking ^k/\/n

directly as a Gaussian variable with mean 0 and variance k/n, rather than

waiting for k and n to become large enough for the central limit theorem to

apply. In practice, we achieve this by letting 5^ * X^ + * • • + X^, where

each X^ is Gaussian (rather than discrete) with mean 0 and variance 1.

» Returning to equation ( 5), we can use an IBM computer to find

<\

exp f"n klnt v( Sk/\/n)}^ and thus X^, as follows:

1) Pick values for n and t (errors involved in choosing finite

n and t will be discussed bel ow)

}

2) Choose a random number from a card containing Gaussian

deviates with mean 0 and variance 1;

3) Compute V(-=)s

Ml

U) Choose a random number Xr>|
. IM

-

•..us

•

*

-
?) Compute V( \

- 9 -

Choose a random number X,

X-+...+X,

Compute V( p -)j

,/5

Continue as long as k < nt;

Compute

as long

i{v& n

-l -

) . vffiS)

/

Compute e We have now found one term in the average ^

exp

namely the term corresponding to the path determined by our particular

, Vn J

{- 5 At «7n’ } > '

choice of the X^'s. Each choice of a sequence X^,Xg, •••

, corresponds

-A

to a random walk$ and if we take many random walks and average e for all

the walks , we should converge to an accurate value of

« kint v^}> • As an initial test of the method, the potentials V(x) = |x| and

2

V(x) * x (harmonic oscillator) were studied by means of 100 random walks

with n = 100, t^ = 3.7?, tg = ?.^ The experimental values for were 0.81

and 0.7?, compared with theoretical values 0.81 and 0.71 respectively.

Ill . ESTIMATION OF THE ERROR INVOLVED IN THE METHOD

The feasibility of the method for the calculation of eigenvalues of

complicated potentials depends, of course, on the accuracy which can be

*h?o minimize a numerical error discussed in Section III it proves

desirable to consider two different va.ass of t. For a detailed account

of the experiments on these potentials, see [2].
'

.

!

. " * ; '

.

.. • *c 1

•

1

\

•

.
- 10 -

obtained in a reasonable number of random walks. There seem to be three

major possible sources of error:

A)

B)

0

Error involved in discretizing the problem and sampling the

1 °k\ values of the finite sum - , 2 . V(-“) instead of J V(x(r))dXr # n n*t> o

i.e. ^expj-i

<^exp{-f V(x (D)dT

.

Error due to using equation (U) for finite t.

Sampling error, due to the fact that we sample only m paths

out of the space of all possible paths

.

Errors B and C, which we shall consider first, are related. It will

be seen in the consideration of statistical errors that it is impractical

to make t very large. Then the numerical error involved in using equation

(U) for finite t is twofold:

log b

1

l) We neglect the term —^ (see eqn.(3)). This error, which may

be serious for moderate t , is easily removed by considering two

different times t^ and Then it follows trivially from equations (3) and (6) that

^

exp ^

- I

' 1 V(x(Z))dZ

(7) A - - i to-t log

1

exp j- |

2

v(x(t-))dtr ^

and that the influence of b^ in the ratio becomes negligible even

for moderate t.

In the expansion of J P(x3 t)dx we omit the higher exponentials

in t (see eqn.(3)). Ignoring the ration l/b9 , we see that the

2)
q n . . " J - 1 \ • ! " . ' ?

•\ ,
- 11 -

-(

A

2

-A )t

fractional error caused by this omission is approximately e ±

. For the harmonic oscillator, A.. = , and hence the fractional error

in estimating the series by its first term is about e For t = b}

this is less than 0.5°/6 . The sampling error C arises from the fact that we try to approximate

the exact average sT-u- £V(x(£) )d£ (henceforth called > ) 1 ® ^

*s by the average - exp -j - ^

V(x^(tO )dtjtaken over only m independent

paths x^CC), ••• y x

m

(tr) . The expected square of the deviation from

<5**> of a single value exp {- l

7( x

k(

t')) dr

}

is given by < c*- 2

> - <PC >^‘

,

the usual expression for the mean square deviation of a single observation

from the mean® Then the mean square deviation from <oC > of the average of 12 2

m observations is - ( > - Thus, the relative fluctuation of

m

9 the average of m observations is given by

i_ y«^>- < >

Vfi >

1_

<c*>2 - 1

Fortunately, we are able to estimate < <*->* for

<oc

2

> - <^expj^-2 J

?(x(t))dtj^) .

Then, applying the results of equations (3) and (6) to the potential 2V,

*L A.t this point we are ignoring the discretization error A. Othern

t

wise, we should properly replace J V(x^(V))dV by a Riemann sum.

o
. ' • . ' T:

Vl *\ ... -r f li, J&tlLtfJ . / • 13 I !

«. . . • - • . ' r - I ' • sxrr «

9V, ' . -
- 12 -

we have ^ .2. . "Pf

<.&** > cL

e

where is the lowest eigenvalue of the potential 2V; ^-i(x) is the lowest

_ ATM

eigenfunction of the potential 2V; « y^(°) J

<^(x) dx.

-oo The relative fluctuation in the mean is thus __

T

( 8 )

, 2 -2 A't

b

l

e 1

- 1

1 (2A

1

-

1i1

)t/2

which, for large t, is of the order -7=, e . For any positive

vra 1 potential V, 2 ^ 0, and hence t must be kept relatively small,

even at the price of an increase in error B.

In the case V(x) * x , we know a simple change of scale

yields » 1. Then the relative fluctuation in <^J> when m - 100, t = 5

is approximately

i e (^-l) 5/2 „ 3o<£ . .

This result is not so disastrous as it looks, since the method of calculation involves the logarithm of <©c> and is not so sensitive to errors in

^“This fact is easily proved by the variational method. On the other

hand, one can convert our reasoning into a probabilistic "proof” that

2 A, - \L ^ 0 as follows : if u, > 2 , then for sufficiently large t we

-fLt -2 2 2 should have e x < < e and thus <oC ^ — <<*> <03 but this would

imply that for sufficiently large t, the mean square deviation of a single

observation of c* from the mean is negative. This is impossible.
! a t

I

. ''

‘

V 4 -: 5 i T I

-

• * I

^

fl

i n i $© »o*' *•" '

'i

,

•

1 f

'

. ! ; a

:i,

.

' ' '' '•

i • ,

.

;

: '

• • - \

-

. . $ :

'

; v

• JV-
- 13 - <©^> * Equation (5) says A^ ** r log <£p<> . Then

(9)

dX /V X — d (log <«*•>) /V 1 d<*>

tx <*<•>

For t = 5, the expected percentage error in our determination of <oc >

for the harmonic oscillator is thus

dX

A iOI7

30C

£ ~ 9% •

The observed relative error in for this case was about nThe error A is not easy to estimate. What we should like to find is

the relative error

<^exp

j

V(x('C))d'c|’^> - <^exp

|;

5 klnt

v(Sk/^]i>

<^exp

-jf-

|

V(xCC))dt }>

W

5

A method for estimating KTr is proposed here, but not rigorously justified*

In the case of the harmonic oscillator, however, we are able to calculate

K

2 exactly (Appendix A) . The answer obtained for K

2 agrees with that

x^ x^

predicted by our more general method. We shall discuss the method briefly]

it is worth noting that in practice the error seems unimportant.

With a given continuous path x(^) x^e can associate another path

x^Ct), which is a step-function approximation to x("C) . We define

Then ^exp {-*

X (t) - X (-)

n n

k<nt ^

when - ^ T* <

n n

* <^exp

|

7(x

n

(t.))dt 'y

.
. : : : ' ; ‘ - ' -

•• . . !!
- m -

,1 ,1

By definition x^O = x(- [ntj, and thus both x(t) and x (z) are Gaussian

random variables. If we plot the time X versus the difference (

X

- -[nt^)

n

of the time arguments of the two random variables, we get a sawtooth function which grows linearly from 0 to l./n with period l/n. On the average,

the step function lags behind the continuous function in time by an amount

l/2n. Thus, with respect to functionals which depend on the behavior of

x (r) over periods of time large compared to l/n, we may replace x

n

(£) hy

x(t*- . (This might be invalid if V is very rapidly varying, but seems

permissible for V(x) * x . In fact, it is easily verified that this scheme

gives x

n

k

('C)dr^> correctly, to the first order in l/n, for any k.)

We therefore conjecture that, to the first order in l/n,

f ^ • f

"t "fc <^exp < - f

V(x

n

(t))dt - <^exp < - ^

V(x(tBut the latter quantity is the same as By equations (3) and (6), for large t

<r { - ?

r" V(x(t))dX

J^>

.

t- i.

- f

2n V(x(t))dt

O

e- A

i

(t- k>

Thus

K

vV - x

i

(t-

— e=7£E

^1- - 2n

to the first order in l/n« For n =* 100, this error is negligible.

[7] * the greatest integer less than or equal to y.
' , . ? . ' tj i .
- 15 - Our estimate of the error K^. is considerably less than that of

Fortet [3], who proves rigorously for the potential V(x) = x , that

K

2 t^j:

2

4

5 V — X' - ^

and regards the right side as a good estimate of the left, Fortet ! s estimate follows directly from his preliminary proof that

x

2 (X)dt^> — <4:exp - f x

n

2 (c)dr} >1 ^ |

£/2

Vn

This result is true, but the chain of the inequalities used in proving it

involves some very wasteful inequalities, which operate to make the right

side a very high bound on the error. The highness of the bound is most

easily seen by considering the case of fixed n and large t. Then the right

side approaches + oo as t —> oo, while clearly the left side approaches 0

(since each of the terms on the left goes to 0)

.

We prove in Appendix A that K —, correct to the first order

1

‘ x ^n

in (This result agrees with our conjecture). For the interesting case

of n - 100, t = 5, we thus have

K '

If we take Fortet* s bound on K

2 as an estimate, we obtain for t

n = 100

- 5,

K

2 - 16 - n.600^5

Thus the latter estimate appears too high by a factor of ?000 and, if

correct, would indicate that the method is totally useless.
... •

.

‘

- '

' •

..

•• .

1 .

• ' :,i

'

. n

| o

. . T ‘

- ' -

:

t

. •

. ;j C i

. /

" ' ' * !

iWI
- 16 -

IV. EXTENSION TO THE HYDROGEN ATOM

The harmonic oscillator , of course, represents one of the simplest

potentials! the results obtained in the study of the hydrogen potential

are of more significance in determining the real value of the method. In

Cartesian coordinates and in the proper units , Schrodinger ! s equation for

the hydrogen atom is

i2-ur -2.,. ^2

1 f9~Y ry ry

V X 'of 't z Y - - > Y • x +y +z

We encounter here the following complicating features, which were not

present in the case of the harmonic oscillator:

1 1

X. V(r) = - depends on three coordinates. Hence

x +y +z

the random walks must be performed in three dimensions rather

than one.

2 e V is negative, has a singularity at the origin, and does not approach oo as r oo (If we want to continue thinking in terms of

random walks, we should now think of V as representing a probability of multiplication, rather than destruction)

3. The spectrum of eigenvalues is partly discrete and partly continuous, the continuous spectrum arising from the fact that V

does not approach oo as r *4 oo 9 The eigenvalues are now

\ * - rj (n * 1, 2, 3, ••• ) and the entire positive real axis,

n

Results analogous to those of the one-dimensional case can still be

obtained. Corresponding to equation (3) we have
.

'

l&M

d-oo

• 4 - 1 '

‘
- 17 -

p

dtr

J W) \oo-X.t oo oo \ , e

3 ^(0,0,0) fff'jf3

(x,y,z)dx dy dz + f e-A’tC(A)dA

where the integral is the contribution of the continuous spectrum. For convenience, we redefine X as i i and can now write

n <L d

where

Cl - Y^OjO.o) jff y (x,7,z) dx dy

-00

dz

The replacement of

J

¥(r(t?))d'C' by a finite sum proceeds as before,

72“’72T2”

s, +s, s. s, with J — —— substituted for — . Here S. is the sum of k inden Vn

^

pendent and identically distributed random variables X^, •••

, X^., each with

mean zero and variance one (similarly and S^) • Thus the three-dimensional path is determined by three independent one-dimensional random walks.

-t A dt . Then -

J*

V(r{X))&V = J r(t7 by

1 2

n k<nt

2

“

T" mT~

sf +sf +sf kx ky kz

n

~ 2

/

2 2 2 sr+sf +sf kx ky kz

As before, we approximate X^ by the relation
- . b /i | , , ' ' ’• . F ' • ; ' ' ' ' " v - . J '
- 18 - </^p>

<t^>

Again we use two different values of t to eliminate the error B,

In practice, some sampling time is saved by the following scheme. We

note that x(ZT) , being Gaussian with mean 0 and variance IT, has the same

distribution as \/t x(£) , Therefore r(t») has the same distribution as t

t

Vt r(^) . It follows that

J — d'T has the same distribution as 1 i ir J vt

J' rfgy

^ • Then, what we actually sample is the random variable

c<- n 1

-J ?5tT

dtr

100

/ —

2

—T

sf +sf +sf kx ky kz

We then compute \ by the relation

t2“tl

log </^>

The values of t used must be large enough to make error B (discussed

in section HI) small* The magnitude of this error, as we have already

seen, is about e (the sign of the exponent differs from that in

the similar expression for the harmonic oscillator because the solutions

now increase, rather than decrease, exponentially). For the hydrogen

atom, - A

2 * 0,2?. The corresponding quantity for the harmonic

oscillator is l.U. When t^ * 9, the fractional error in
1 '

"i - , „
- 19 - \e ) due to B is thus about 10% , and the consequent fractional error in about 2°/0 . For t

2

23 16, the accuracy is better.

What happens to the sampling error C when we use such large values of

t? If we sample m paths, we recall that the relative fluctuation of the

X

m

o rk

(*

)dt' 1 ( f

I1

-2 *1

) t/2

average - 2 e is approximately -7

= e $ u-,. the lowest m k=± vm 1 2 eigenvalue of the potential - - , is easily calculated as 2. Thus

Pi

- 2 ^i - 1, The corresponding term in the statistical error computation

for the harmonic oscillator was .I|lU, Hence the statistical fluctuation

for a given number of samples m and given t is considerably larger for the

1 2 potential - - than for the potential x „ Intuitively, one can explain this increase in the statistical error in

the following manner: The eigenvalue of the hydrogen potential depends

mainly on the behavior of the potential in the small region near the origin.

We see this fact immediately when we note that the main contribution to the

<•

J*

t 1

ritry

dr

average V e >> comes from the few paths r Cc) which stay close to

the origin, Plence, in order to calculate the eigenvalue accurately we must

get a detailed account of the behavior of the potential in the critical

small region about the origin. Naturally, many samples are required before

we get a sufficient number of paths staying close to the origin to sample

that region accurately. Equivalently, if we do not take a very large number of samples, there is a good chance that the few critical paths which

stay near the origin will not turn up In the proper frequency, thus

causing a large error in the eigenvalue , In the case of the harmonic
\

'

.

'

' '

' ' '

'

. i 9 •

* .: • •

•

.

;

.. a > m - § 'IWM

'

:

eH
- 20 - oscillator, however, there is no such small critical region which must be

carefully sampled* Remarks similar to these about hydrogen may also be

made about atomic potentials in general

.

Returning to our computations, we see that for 900 samples with t = 9

m

the expected relative fluctuation of - kli e

•t i

i

^dt

V 900

9/2 ~ 35 ~ 300°^, .

is about

If we compute more carefully, using the exact expression for the relative

fluctuation given in equation (8), we obtain 100°^ instead of 300^. It

follows that for m = 100 and t = 9,

d

X

x

1 d <£*->

t \

1

TZ5 100°7o 22®?, .

For t = 16, the error is of course much greater. For the initial group of

900 hydrogen random walks, with t^ = 9 and t^ = 16, computations according

to equation (10) yielded

.503 (true value *p00)

In view of the preceding remarks, this accuracy was quite surprising.

An examination of the data, however, seems to confirm our previous analysis. Of the 900 oaths sampled, I4. yielded extremely large values of

-1 x

J rTFT

^ * without these I4. paths, which clearly cannot be counted on to

recur regularly, the average e 4

the computed value for A about 20°/o smaller (note that the logarithmic

method of calculation is not too sensitive even to a large fluctuation in

6

Uoc

>).

would have been smaller, and
»

.

•

.

. .

:

:

'
- 21 -

Thus, there seem to be good grounds for attributing the high accuracy

in this determination of A^ to luck. Nevertheless, it seemed worthwhile

to run another group of samples. This was easily done in connection with

experiments on helium (described below), since each helium random walk

yields two hydrogen random walks as a by-product. Two more groups of 900

hydrogen walks were thus obtained.

The qualitative appearance of the data was similar to that of the first

group of walks. In the second group of 900, more than Q0°/o of the contri-

/ L°£

,V

bution to C e ycame from three walks yielding high values of ot , and

the calculated value of was 0.36. If a few more high values of had

been recorded, the calculated value of A^ would have been much closer to

0.50j but this again would have been purely a matter of luck. In the third

group, there were only two large values of o(, with the result that the low

value * 0,29 was obtained.

A more detailed analysis of some aspects of the hydrogen random walks

is given in section V and in Appendix 3.

The preceding discussion of hydrogen, which may be extended to other

atomic potentials, makes one fact rather clear; for t great enough to make

our numerical approximations good, the statistical fluctuation in

dt

/ £ Vpf \ (e K /

becomes prohibitively large unless we use a huge number of

samples. Equivalently, our computed value for is determined almost entirely by a very small group out of the 900 random walks. Hence, most of

I the data is wasted and our large group of samples acts statistically like

a small group.
.

.

: ~>d J

« » a* m .. Iml -Hu H M .

*’

.

U :#!

'

.

.
- 22 -

There now arises the ^question of whether there is any other method of

utilizing our data to find A_^, by means of some computation less subject

to statistical fluctuations. In the next section we shall deal with such

t a method, which we call the slope method , as distinguished from our previous

mean value method.

V. THE SLOPE METHOD 5 RESULTS FOR

HYDROGEN, HELIUM, AND n DOUBLE HYDROGEN’ 1

The method for finding A

^ which is to be presented here has the following advantages over the mean value method:

1) It makes use of the finer structure of the distribution of the

;

t d

IrTiET anc^ o ^ ' 2) It is quite insensitive to erratic individual observations of

the random variable* —

Consider the random variable

r^TfJ

* this random

variable assumes depends on the particular path r(r) which is chosen for

the path of integration. We have already observed (section IV) that

J

1 ^as the same density function as ]/i

J* «1

Thus, if we let /?(x^ be the density function of

j r i

-t 0

(i.e. /

> (x)dx = Prob -I x < J* ^

L o t J i

p dr

<J ,

,

Hr) < x + dx f ) we have

}

P

dr

\ft J vi'Z ) J' X

/(x) dx ~

o

A-,t

bn e t—^oo 1
' « / . c ti b n f ' .fSJi l ' ' • . '
- 23 -

or, writing t for t, we have

oo

o (n) J etX

/’(x)d

x

t 00

V

Alt'

One can verify, by a simple integration, that this asymptotic behavior of

oo tx

J

e y?(x)dx can be achieved by setting

( 12 ) /(x)

VhX1

X—00

2\TWT.

By an unproved Tauberian theorem which, however, is intuitively reasonable

and similar in form to the Hardy-Littlewood Tauberian Theorem, this is the

only asymptotic behavior of ^(x) which will make equation (ll) true. A

result similar to this can be obtained for any homogeneous potential function V.

Now we can determine A as follows:

1 1 dZ

l) Find the density function of

J* by replacing the integral

with a finite sum and performing random walks as before, to

get sample values of P • b ® 2

2) Letting x = J

1 \

, plot log /’(x) vs. x for large x. It is

o

^ to be hoped that a fairly straight line will result. The slope

of this line should be -

First we present uncritically an example of how this technique is applied to the hydrogen potential.

’The theorem can be proved in idle integral:|^form

//(Od f /
4

••

H . -

- & '••p£2SL

I I

’ ’ ^ ;f

' g

* Qm®T

.

-

.

.
- 21* -

Clearly, we cannot experimentally find a continuous density function

1 dZ

/?(x) . Our experimental data consists of 900 values of J

1 obtained

o rv / from 900 random walks. The closest we can come to plotting a density function is to plot a histogram, which shows the number of occurrences of

1 dZ

S rTcT

each of a sequence of intervals of length A; e.g., if

o

' ' we take A = 0.1, we obtain the histogram of Graph 1. The bar covering the

/

I j

'-T'T o

' ' lay between 1.5 and 1.6 (= 1.5 + A). Of course the choice of A does much to

determine the appearance of the histogram; a small A gives an irregular

curve, while a large A obscures the structure of the histogram by "smoothing” it too much. Since the purpose of our histogram is only to study the

behavior of the tail, it is clear that A = 0.1 makes the tail too irregular

for reliable use in our calculations. If we replot the histogram with

A = 0.2 (graph 2), the tail is considerably smoother but still retains some

fine structure.

2

Now, if we plot log f{x) vs . x for large x, we expect a straight

line with slope - • faking values of f from graph 2 and plotting

2 ^ them against x on semi-log paper, we obtain graph 3* Here the "tail” of

p is considered as beginning at x * 1.2, just to the right of the peak

of f . The question immediately arises: what is the best line to draw

through the points? We reason as follows (this argument and the general

character of the graph are typical of any potential):

1) The first few points at the left are not really in the tall

since x is not large enough. This is also clear from the
:• . i . .
90

80

70

60

5o

Uo

30

20

10

Graph 1

Histogram for 900 Samples

1

of J

dvZr('c)

A * 0.1

0.5 1.0 1.5 2.0 2.5 3.0 3.5 U .0

Values of J

dt/rfc)

Frequency

of

Occurreno

- 2‘Ub -

.ao

*

160

lho

Graph 2

Histogram for 900 Samples

1

Of / dVr(t)

120r\ !-

CO

100

80

60 -

Lo

A = 0,2

o.5 1*0 i c 5 2.0

n

X. Values of P dtr/r(f)

2.5' 3.0 3.5 U.O

Frequency

of

Occurrence

- 2 he -

200 r

- 25 -

graph, where it is seen that the points cannot be well fitted

by a line, i.e. they do not lie in the region where

—x2/4^2_

P(x) ~ const, e

2) The larger the x-coordinate " int, the more valid is the apThus, the data points becpme statistically less reliable as the numerical approximation improves. Consequently, we reason that the best line is

the one going through a set of points whose x-values are as large as possible, but which still represent enough data to be fairly reliable. Of

course, the test of the "reliability” of a point is that it lie on a line

with several other points. Near the mid-portion of the graph there generally appears a sequence of such points, followed by a number of "wild"

points which represent insufficient data. Drawing a line on graph 3 according to these principles and measuring its slope, we obtain A = .1*8

.

Obviously it is difficult to give a fair estimate of the error in such a method. A good upper bound, however, is obtained by looking at the slopes

of the most extreme lines one might attempt to draw through the data

points, and seeing how much those slopes differ from that of the original

line. In this case the answer is quite certainly less than 1096.

proximation ^(x) ^ const, e . But

3) the points which represent really large values of x represent very

little data, i.e. there are few occurrences of very large values

. Hence the points with the largest x-coordinates are

subject to large statistical errors.
y

< •

! :

• - I iv <

;

.

• , " -

' " * : 1

' \ ,,,J

-

,

.

;

•

... . . • :

• .* • •-

•

. r •. : w

\ft

‘ IS | m || I |

«

:

:

•

-

. • >

'

.

•

, ;

"

1 •
- 26 -

"The discussion gains "ibility who.. we consider the results for

the hydrogen potential obtained from 1800 mere random walks. The histogram

(graph U) has the same shape as graph 2 and has a peak at the same value

of x. When we plot log />(x) vs. x for large x, we get graph 5>, which shows

a marked improvement over graph 3. With the exception of the last two points

on the right, which represent insufficient data, the points quite certainly

lie near a line. Calculating the slope of a line fitted by eye, we find

A = .52 . Here the error seems less than sh°.

It is worth noting that about 30°/o of the random walks yielded values

-y

\ which are large enough to be considered in the "tail11 of the

o

^ ' distribution. Thus, in the case of the hydrogen atom, one can say roughly

that the slope method utilizes 30°7o of the data while the mean value method

t

makes use of only about 19& (since the average e° ^

is determined

largely by a handful of walks yielding large values of the integral) . This

difference in efficiency of the two methods seems critical in the case of

hydrogen, since the mean value method is infeasible for that potential

T<rhile the slope method is reasonably accurate. We can get an idea of the

limitations of the slope method by considering a more difficult potential.

In the proper units, the stationary SchrSdinger equation for helium

can be written as

| q

2

'k (5?, + \ V k

2

+ 1^1+ pjr’pi-“r

2

I /

( 13 )
mm !

'

'

'

•

. :

• *1 • < x ^

V

.

;

" • - - * -

- . ..3 .

.

' .

M-

" '

'

-

' -
Frequency

of

Occurrence

- 26a -

Graph U

Histogram for 1300 samples

of

J* d^/r(£)

o

2U0"

200 -

160 "

120 -

80-

Uo -

A » 0.2

1 1 1 < 0.5 1.0 1.5 - 2.0 2.5 3.0 3.5 u .0

Values of

j 6X/t{V)

o

- 26b -

- 27 -

In order to find X

^ by the slope method, we must determine, by sampling,

the distribution of the random variable

/ f-J— + —2 i—

where r^(£) and r^(t) are the paths of two independent Gaussian random

walks. As before, we discretize the walks and sample

.1

100

k*l

2

/

2 2 2

S. +sf +sf kx

2 ky^ kz

2

(1U) *

^ <VV*+( VV*+( VS)2

-

,Sky1,S

kZi

dependent three-dimensional random walks, as in Section IV.

Before we discuss the results for helium, let us consider a slightly

simpler potential, whose lowest eigenvalue can be calculated as a by-product of the helium computation. Consideration of this potential will

bring out the difficulties involved in the multi-electron potentials. The

potential we deal with is

where (S.kx1

) and (S^. ,Skz ) determine the paths of two in-

+

which corresponds to a helium atom without interaction between the electrons . This potential has been christened "Double Hydrogen". The
'

'

.

' '

•

• * k d ;

' "'.'-I

l

. *
28 -

Schrodinger equation for double hydrogen is

I «-» —v \ . ± n 2 . j ft (15') i

V ‘Xf (q,r

2

)H | C’fy,r2

) + (• lb

j£_)

>—* 1 ' ir^ i Ub ,v9

) « - +\rIt is easy to show that the eigenvalues of this equation are

'* = Pi + jv where p. and p^ are eigenvalues (not necessarily distinct) or

(16)

' £ The eigenvalues of (l6) are \i^_ * - (k - 1,2,3, ••• ) and the entire

k

positive real axis. The lowest eigenvalue of (l5) is thus = -Ii, and

the second lowest is Ap = -2.5. As before, we now redefine A so that

, ^

= U, “ 2.5. (It also follows simply from, our probabilistic approach that A

1

« 2ytj , since

t

o X1

, /

( nfw + ^m)ax

X ^ lQg-S e >

1

but r^(T) an(^ t

2

(^) are the paths of independent random walks, and hence

e e e .

The desired result follows.)

Of course, since we know that A_ - 2ii_ , one statistical method for

1 : JL finding would be to find by saiy rag the potential - , and then t< double the result. This scheme, however, breaks down as soon as the

-4 potential contains an interaction term between it. and r0 « Hence, since

_L
.

’

!

,
- 29 -

we want to test the general applicability of the method to multi-electron

potentials, we do not employ this short-cut, but sample the random

variable

Plotting the logarithm of the density function of this random

2

variable versus x , we obtain graph 6. The graph is based on 900 values

1/2 2

\

of S ( r Wj + r ~'(v) )

^ 9 °“ &bout 500 are represented on the graph

(the other values were too small to be in the "tail” of /’(x)). On the

basis of the solid line drawn, we find ^ = 5.6 (true value = U.0) . It

is evident that the data points do not clearly determine a line. If we

ignore the points on the left, arguing that they are not yet in the tail,

and draw the dotted line based on the right points, we get A^ = U.6. This

result, however, is hardly reliable. If we increase A from 0.6 to 0.8, we

obtain graph 7, which is considerably smoother. The slope of the line

through the first four points on the left yields = 8.5. Clearly

these points are not yet in the tail. The next three points represent

only 29* 8, and 6 walks respectively and are subject to too great statistical fluctuations . The expected statistical fluctuation of a point representing n walks is about \/n. Hence the limits of inaccuracy cf the

data points are represented approximately by the vertical bars we have

drawn. We see that the points are not inconsistent with the true value

A = U, but they hardly give us any further information.
.

-

' r

'

-

«q. 'i. - j V

* T

«) ^1

.

. r. -.

• - • \ r . ’

•

•

' ‘ * •

'jHVa

;

' "

,

•

• .

.

•
8.5

Graph

7

Logarithmic

Histogram

29a

OOU0JCJTLOOO jo Aauenbain;

30 -

The redeeming feature of the slope method, in the case of the hydrogen potential (and thus also for the potential -) was that the exponential

” tail" of the distribution started early enough so that the points in the

tail represented enough walks to be statistically reliable. This feature

does not seem to hold for double hydrogen. We now proceed to a more precise discussion of this point.

A simple extension of the reasoning leading to equation (12 ) yields

(17) r(x) ~ —

-WCD 2/WX

-x

2

/Ut

2\ZrirX a

”v " ^ 2

Thus, an estimate of the goodness of the approximation (12) is given by

% the ratio of the first two terms on the right side of (l7)j this ratio is

approximately

2/1 1 \ x ( ur up

When x is large enough to make this ratio large, we can say we are in the

tail® Thus, the criterion for being in the tail is

(18) tj ( t2 ~ a~

}

58 M

where M is some constant, whose value depends on how good we want the

approximation (12 ) to be.

For the hydrogen potential -,1we have -r- 11“* “ 8-2 = 6, whereas for 1111r Ar) A-i

<- ±

double hydrogen " j; “ . Thus, if the tail for hydrogen

starts at x

q

, the tail for double hydrogen starts at \/IJ0 Now we ask
. . •* ' ‘I . • y : l
- 31 -

the question: how mar$r random walks for double hydrogen must we perform

in order to obtain an accuracy comparable to that obtained for hydrogen?

If we perform random walks for hydrogen, the number of walks

yielding values of f

in the tail of /H(x) is about J

/^(x)dx.

n

x

o To achieve a comparable accuracy for double hydrogen, we want to get the

same number of points in the tail. Hence we require

00 pCD

/ fn{x) & = "dh J^r ^dh n (x)

X

where = number of double hydrogen random walks

/

?I)H

(x) = density function of

J* + wfrj dV •

Ignoring the constant factor in equation (12), we thus have the ordercf-magnitude equation

^DH

x

s

oo -x / 2

dx

N,H

u f

-x /l6 , e ' dx

x

Using the asymptotic series for the two integrals, we have, if x

q

is large.

2/2

rr -X 1 ^DH /p; e 0 _A/ Vh-O N

H -UOx

1 2

~X

L do

5713 ~ 6 e

Looking at graph 3, we see that the tail seems already to have started at

2

x =2. The above calculation is inaccurate for such a small value of

o (we could calculate ^DH/N^ more exactly, but it hardly seems worthwhile);
! s t . ' .
- 32 -

it yields the qualitatively correct result that, for a given accuracy in

the lowest eigenvalue, we require many more samples for double hydrogen

than we do for hydrogen.

Formally , the reason for the increase in the number of samples required is the creeping together of the reciprocals of the first and second

gen, we sample a single random variable X* for double hydrogen we sample

the sum X + Y, where X and Y are independent and identically distributed.

Roughly speaking, if p is the probability that a value of X lie in the

tail of the distribution of X, then the probability that a value of X + Y

2 lie m the tail of the distribution of X + Y is only about p . This argument can be made quantitative and leads to our previous result.

We conclude that 900 random walks constitute far too small a sample

for the double hydrogen potential or, more generally, for ary multi-electron

potential. The 900 helium random walks confirm this conclusion. The logarithmic histograms with A = O.U and A = 0.8 respectively are given by

graphs 8 and 9. The graphs are consistent with the true value X - 2.90,

but not enough points are in the tail to allow an accurate determination.

In the next section we discuss some schemes for making our method

feasible for the study of more complicated potentials. One more point is

worth mentioning, however, in connection with double hydrogen. We plotted

a logarithmic histogram (graph 7) for this potential, based on 900 values

eigenvalues (see eqn. (18)). n Physically” , the reason is this: for hydroLet us denote these

900 values by An +- A

2 (1=1, o «*,900) where

“i fc i 1

etc
: .

.

.

.

v

,

;
- 32a -

©

o

0u

u

3

oo

<Ho

r*5 O

£

0

&

©

£

20 30 HO So 60

-33 -

However, since /o

1

are independent, A^ + A,

i

2 (i ^ j) or

a

A-, + A-. (i / j) or A0 + A* (i / 3 ) also constitute legitimate sample

JL 13> 1 » tma 13+ €a. % values of the random variable X. Thus, out of our 1800 values A^ , , i i

we can manufacture 1800 x 1799/2 ** 160,000 values of the random variable X,

instead of merely 900 values* Of course, these 160,000 values are not all

a maximum of only 900 independent values of X in the set of 160,000. However, one feels intuitively that the dependence is not too strong and

that something might be gained by using the additional data. (The whole

situation can be treated quantitatively, but we omit the treatment here.

An analogous problem is that of sampling the distribution of the sums

rolled on a pair of dice. We roll 900 pairs, but then combine the 1800

single rolls into all possible pairs.)

Actually, we did not combine the data into all 160,000 possible pairs,

but considered only 18,000 pairs. Graph 10 is the logarithmic histogram

based on these pairs. The improvement over graphs 6 and 7 is marked, and

we claim about 6°/o accuracy for the result X * 3.99. One might tend to

challenge this result on the ground that we are getting n something for

nothing”, and thus violating what some consider to be a fundamental law of

nature. However, we should remember that the data already contained

enough information to find the lowest eigenvalue of ” single” hydrogen to

about 5% (graph 5), and thus (by our remarks on p, 28) the lowest eigenrvalue of double hydrogen to the same accuracy. Graphs 6 and 7 merely

represented inefficient use of the data.

There are
*

.

*
- 33a -

i

I aouajjnooo jo JCousri

- 3h -

This device of "multiplying" the data fails, of course, when an interaction terra is present (as in the case of helium) . Our final section deals

with procedures feasible for this case,

VI. TECHNIQUES FOR SPEEDING THE SAMPLING PROCESS

;

HIGH-SPEED COMPUTING AND IMPORTANCE SAMPLING

The experiments described in the previous sections were perfomed

mainly out of curiosity, to see whether equation (5) is verifiable within

a reasonable number of samples. At a later stage, one is naturally led to

ask whether this sampling method is of practical use, e.g. as a possible

competitor with the variational method for finding the lowest eigenvalue.

The authors are unaware of just what potentials are of real interest but

unamenable to variational treatment. However, there are a number of techniques which should make the present scheme feasible for the treatment,

without excessive sampling time, of any system whose dimensionality is

small. These techniques, which will be dealt with briefly here, are

1) The use of a high-speed digital computer, eliminating

entirely the "slow" punched-card operations;

2) The employment of sophisticated sampling schemes, which

we call "importance sampling"

.

Both techniques can be used together, effecting a further reduction of

sampling time

,

/
••

.

'

.1

i' i , . ml if m

• -> C

, ' ‘ •

1 * * ^

1 s> *

. :

, - ri ' ;

‘

•

1 . ,

• "

; 7 n ,1

1 $<t

1

.

...

.

'

.

v *

i ‘flu it r*

•

1
- 35 -

For an indication of the times involved in the use of IBM methods, we

note that 900 hydrogen random walks of 100 steps each require about 90 IBM

hours. To accumulate the same data on the National Bureau of Standards

Western Automatic Computer (SWAC), about one hour is required. No hydrogen

walks have yet been performed on the SWAC, but harmonic oscillator results

will soon be available. One can say, in general, that for our purposes the

SWAC is faster than IBM by a factor of at least 5>0 . The programming of

these random walk routines is quite elementary, involving mainly the iteration of simple operations . The problem of supplying random numbers , to

govern the random walks is of some interest. External sources, such as prepared tables, are unfeasible, since the machine must generate its own

random numbers if it is to operate at electronic speeds. A number of

arithmetic schemes for the internal generation of "random” numbers have

been proposed (the numbers so generated are only pseudo -random, since they

are completely determined by the numbers inserted into the arithmetic

scheme), and a particularly simple method of Dr. J, B, Rosser has been

found to produce digits which satisfy the usual tests of randomness.

Regardless of whether SWAC or IBM equipment is being used, one can

gain another factor of 10 or more by the use of importance sampling. This

topic has been discussed in recent literature (see [5]) and is of enough

simplicity and interest to warrant a brief review.

Suppose wTe have a random variable X which is continuously distributed

b

interval (a,b) with the probability

a over some
-

.

C . <fu K- or

' I( is

' ^yffH

, . - .•

: f

--

' : .

'

.

.

.
- 36 -

' and suppose we are interested in finding the expected value f of some

b

function f (X) • Then f = J f (x)/? (x)dx. (Alternatively, and for greater

a clarity, we sometimes write f = E^?(f(X)), i.e. f is the expectation of

f (X) over a population in which X has the density function /?(x) .) The integral can, of course, be estimated by sampling: we merely pick many values

of X from a population governed by the density function /?(x), compute f(X)

for, each value of X, and average the results. Suppose we also compute

(f - f) d ~ J* (f(x) “ f) 2

/

? (x)6x = (Ty p . a

Then it is well known that the expected error in our estimate of f, if we

take as our estimate the average of n sample values, is -i- * • s/rT s >f

Importance sampling is based merely on the recognition of the fact that

it is possible to bias the sampling procedure in such a way as to leave f

unchanged but to make smaller. Specifically, suppose y(x) is another

probability density function defined on the interval (a,b) with

b

^(x)dx * 1. We certainly have

a

(19) e

p

(f(x) ) = y

3 f( x)/(x)dx=y (x)dx ,

Thus, in our notation,

£U)P(

w (x)

00 E

p

(f (X) ) = E

^

i.e. the expectation of f(X) over a population governed, by the density

function ^

is the same as the expectation of f

f/(f over a population
'

..

•

*

<

'

•

• %

•

1 rl^sa

*
- 37 -

governed by the density function The interpretation of this fact is

f ’ simple: we bias the density function of the population by the factor

P

and hence must ’’weight" our values of f by the factor if we are to get

the correct answer. Let us compute the expected error in the average of n _ ]_ samples, if we estimate f by this scheme. The error is — „ . where

oy,

i,

dx

By the proper choice of the function y(x), ^

may be made much

smaller than - . In particular, if we knew T (which is of course what

we are trying to find), we could let

( 20 ) (^(x) = p{-x) f(x)

and thus have no error in our estimate. In practice, we do not know f, and

it is not easy to concoct a population governed by an arbitrary density

function fU). We can, however, use whatever we know about f to improve

our sampling procedure, for (20) tells us that the error will be decreased

if we distort p into a new distribution with more mass centered in the

regions where f(x) is large. In practice, it is not hard to construct a population with this general characteristic.

A simple illustration may clarify this

discussion. Suppose f(x) is a function of the

JL type pictured, and we wish to find f f(x)dx

X)

by a sampling procedure. This is the same as

"l 0 * X * 1

finding Eo(f(x)) where f(x) = ' $

0 otherwise
.

' 1C <

••

.

i

.

.

-

...

.
- 38 -

i.e. we pick values of X from a rectangular distribution on (0,1), and

find the average value of f(x). It is clear that most of the data are

wasted in this procedure, since the major contribution to the expectation

comes from the few values of X which we happen to choose from the neighborhood of 0.5, One easily sees that a better procedure would be to bias the

sampling scheme so that most of the points sampled lie in the vicinity of

0.

5. It is then necessary to weight the results somehow, so as to compensate for the bias introduced. The previous paragraph shows, quantitatively,

how to do this. The term "importance sampling" arises from the fact that

the amount of sampling time we spend in studying various parts of the population is proportional to the importance of those parts in determining the

average

.

The applicability of this technique to our problem is immediate. In

f |

dS/r(t)

the case of the hydrogen atom, we are interested in E \ e

1.

e. we sample from a space whose elements are paths rft)

,

and compute the

«f ar/rCc)

average value of the functional F(r) = e .In practice, the

space of paths from which we sample contains only a finite number of elements r^ (i ® 1, •••

, N), since the paths are determined by the choice of

a finite number (300) of Gaussian deviates, which themselves have a basic

—3 th increment of 10 . If p(r^) is the probability associated with the i

path, we can write

- 39 -

N

ill F(r

±

) P( r

i) = E

p

(F(r)) .

We have observed in Section IY that the major contribution to this expectation comes from a few paths r^ which stay near the origin, i.e. for which

F(r.) is very large. It follows from the foregoing discussion that we can

decrease the error in our estimate by associating a new set of probabilities

q(r^) with the paths. The q(r^) will be an improvement over the p(r^) if

q(r^) * p(r^)f(r^) where f(r^) is, roughly, proportional to F(r^) (it is

sufficient that f be large when F is large, and small when F is small).

Then, as before, we have

(21) Ep (F) = -l F(r

i

)p(r.) = & [

F (r.)

^

]

q(r.) = Sq(F^-) ,

but the variance of the q-scheme is less than that of the p-scheme.

There arises the practical question of how to find a distribution q

having the desired property, without hopelessly complicating the calculations. The scheme presented here is, perhaps, the simplest possible one.

Essentially, we bias the random numbers chosen so that at any step the

’‘particle” is more likely to step tox^rard the origin than away from it.

Most of the paths sampled will, therefore, be the interesting ones x-jhich

stay near the origin. Furthermore, the biasing is conducted In such a fashion that the factor p(i\)/q(r^) is easily computed. ¥e recall from

Section IY that in the unbiased x^alks one step (say the k ) consisted in

choosing three numbers from the Gaussian deviates with mean 0,

N P(r±

)
' • ' . . : .-J: < . • :
- Uo -

variance 1, and basic increment .001. These three numbers represent the

components of the step on each of the Cartesian axes . Then the probability

p(r^) to be associated with a 100-step path r^ determined by the 300

choices X

1± , Y

1±9 Z1±, ••• , x

iQ0i* Y

100i*

Z100i ls

As fan as we are concerned, the only important property of p(r.

)

is that it

has. the form

100

= Vi * 10-9

where g is a probability density, and an even function. Suppose now that

we bias the walk in the following fashion: let oC>/3> o with = 1;

"fell at the k step, pick as before from a Gaussian distribution, but ignore

the sign of X, ; then we play some game of chance with two possible outcomes

K " k-1

A and B, with nrobabilities <=^ and /3 ; if .ZL X. ^ o and the outcome is A,

k-X 0 3

interpret X, as negative; if .21 X. ^ o and the outcome is B, interpret X, K J =

0 - - r a i as positive; if ,.2^ X.. < o and the outcome is

j ^ k interpret Xk

as

I positive \

| negativej

# Sirailarly, choose Gaussian deviates and play the auxiliary

game for the Y- and Z- coordinates. Clearly, the effect of this scheme is

that, at a given step, the absolute value of a particular coordinate decreases with probability oL and increases with probability . Hence, the

walks are systematically biased toward the origin. The two-outcome game

involves little extra effort; if we take = .6 and /3 = . 1;, we need
ij

'

'

:

.

•

-

. .

'

rr.

:

' I

'

.

.

•

'

; .. v _
- Ill -

merely pick a digit out of an equidistributed population of zeros

,

ones, •••

, nines. If the digit chosen is not greater than 5, we call the

outcome A, otherwise B.

Let us write down q(r^) , the probability of following the path r^ under

this biased system. Let

8 =

k-i

ft if has the same sign as ^2^

k-1

if has opposite sign to

j5l

x

3

(Similarly 9 (Yk

) , 9 (Zk

)) .

t It is easily seen then that

q(ri

} = [

2^X

ki> 9 t\i)] 9 (yid}

] [

2«^> 9 (Zki} ]

Then

/ \ 100 n OAA -N. „ N. -300

p(r. ) OAn t

| 1 o-300 > •- i /x i

iTJT7

= 2 31 ^ki59

^

9

^

= ^

where IL is the number of times the auxiliary game yielded the outcome A.

We can rewrite this last result as

p(r.)vr \W_. i

J /I Vi /_i\ ^7 W v

2 ft)

X 300-I'k

The computation of this factor adds little complexity to the IBM work;

we merely keep a tally on N^, the number of outcomes A, and read the

weight factor p(r^)/q(r^) from a table. The extra IBM time involved in
.

•

1

.
-lousing this importance sampling scheme is only about ±Q°/o or 15% .

As a preliminary experiment, mainly to test the qualitative features

of the method, 65 hydrogen random walks with importance sampling were performed. ^ was computed by means of equations (10) and (21), with

t^ = 9 and t^ = 16. The results seem quite promising.

One point concerning our treatment of the data, however, should be mentioned at the outset. The possibility exists that we may sample a path

whose weight factor p(r^)/q(r^) is exceptionally large because the auxiliary

game yielded the outcome B improbably often. If we sample enough paths, the

effect of such aberrations will be negligible (unless, of course, the importance scheme was incorrectly chosen and virtually neglects significant

regions of the sample space, which make themselves felt only through such

aberrations) . However, if we sample only a few paths and happen to include

one with a veiy large weight, the average value E {ftis determined

largely by the single path with high weight, and the rest of the samples

are useless. It seems reasonable, in a small sample, to ignore such a path

(certainly so if we know, as in the present case, that it represents an

unimportant region of the sample space) . Consequently, we have neglected

a path which occurred with a weight factor 25 times greater than that of

ary other path, and 100 times greater than that of all but three other paths.

These other three paths could have been neglected, with almost no change in

our results. ("What we have done is statistically quite justifiable, and

corresponds to the rejection of measurements which fall several standard
. - • V L * ' ' .
- U3 -

deviations from the mean! the measurements may be correct, but the sample

is not large enough to give them meaning.) The following discussion refers

to 6b random walks , omitting the one with the large weight

.

On the basis of 6U hydrogen random walks with importance sampling,

« *5>l8 (true value ,5>00) . We recall from Section IV that three groups

of 900 random walks without importance sampling yielded the values ^ = «30, .36, and .29 respectively. The question immediately arises of

whether the accurate result yielded by importance sampling may not be mere

luck, like the result for the first 900 unbiased walks. To Investigate this

possibility, the 6I4. walks were divided into arbitrary subgroups, and

was recomputed on the basis only of the data in each subgroup. If our

result were due to luck, the computed values of A- should vary greatly from

one subgroup to another. Actually, a high degree of consistency is observed in the values of A^ so obtained. For example, a division of the

walks into four arbitrary groups of 16 yielded the values

\ = 0.1*9 , 0.39 , 0.58 , 0.53

respectively. The corresponding computation over six arbitrary groups

each consisting of 300 unbiased random walks, yields

X = 0.29 , 0.140 , 0.38 , 0.1*7 , 0.1*9 , 0.35

respectively. We can therefore safely state that groups of 16 randomwalks, biased according to our scheme, give at least as good an estimate

of (by the mean value method) as do groups of 300 unbiased walks.

JL
' H . ;« f

.

‘

'
- bh -

Importance sampling is also applicable to the M slope method11 of Section V. In constructing a histogram for unbiased walks, we plotted x

versus the number of observed paths r. for which x < F(r.

)

- x + A. If i x enough paths are observed, this latter quantity is proportional to

p(rj

2 p(r.) which is the same as 2>,. . . —7

-—r a(r.). But the

x<F(r^)-x+A x 1

' latter quantity can be estimated by the biased walks, and is in fact propor4 tional to the sum of the weights p(r^)/q(ru) of the paths yielding values

of F(r^) in the range (x,x + A)* Thus, we can plot a histogram for the unbiased walks on the basis of data obtained from biased walks ; for a fixed

number of walks, the biased scheme does not give a good picture of the entire histogram, but yields an accurate picture of the tail, where most of

the data are concentrated. Such a histogram was plotted on the basis of

the 6H biased walks; there was insufficient data, however, to define anything more than the order of magnitude of the slope of a logarithmic histogram, which yielded

With very little extra computation, the biased hydrogen walks furnish

us with data for biased haimonic oscillator random walks . We recall that

the lowest eigenvalue of the harmonic oscillator is computed from the cuan-

-J

,fc x

2

(t)dtr

tity E > e

°

*

, which we estimate by random walks • The walks

which spend much time near the origin contribute the larger values of

this exponential; hence, the biasing scheme for hydrogen, which stresses

these same paths, might be expected to be useful for the harmonic oscillator. On the basis of 100 biased harmonic oscillator random walks, we
. '{•.. •"

• : - • 5

.

’ ••

,* •• • >

. .

. V

•••

.

•:

.

.... .

• .

’

’ •

. ... r l .
obtain * 0.81 (true value 0.71). Division of the data into ten subgroups, each of ten walks, yields

0.66 0.85 0.86 0.78 0.66

0.8? 1.15 0.96 0.81 0*71

respectively, for each of the subgroups* In this case, our final result

for is less accurate than that obtained from unbiased walks, but our

subgroups are more consistent.

One can think of many other types of importance sampling, e.g, where

the strength of the biasing is a function of position, or where the biasing

scheme is successively improved by some sequential procedure. Even with

the present biasing scheme, lowest eigenvalue of helium should be obtainable

to within 2rfo in two hours 1 time on the SWAC. The authors will welcome

correspondence concerning unsolved potentials, to which the sampling method

may be applicable*

The authors wish to express their thanks to Dr. Wolfgang Wasow, of

the National Bureau of Standards, whose valuable comments brought to light

many obscurities in the original manuscript; and to Dr. Everett C. Yowell,

of the National Bureau of Standards, who not only supervised the IBM work

infallibly, but reduced card jams as fast as the junior author could

create them*

3/27/52
'

"
- 146 -

APFE1DIX A

COMPUTATION 0? DISCRETIZATION ERROR FOR THE HARMONIC OSCILLATOR

The discretization error K 0 was defined in Section III as

r nt

n

<

exp

j

- \ j£ S,J - <^exp

|

- f x

2 (z)dV

We define

<Zej®

|

-

J*

x

2 (E)dT

|

X(n,t) - <<xpj- \ S

k

2

}> .

We shall compute l(n,t) exactly. Since, by equation (5a)

lira I(n,t) = /exp <f ~ J

x

2 (r)dr V/> , n—>oo \ l o

J '

our calculation of l(n,t) will yield the exact value of the Wiener integral

as a byproduct. Then we can compute K 0 , x

We recall that S, - S, n is a Gaussian random variable with mean 0 and

k k-1

variance 1 . Therefore

Prob

{

x

l < Sx < ^

+ d^, x, < S

2 < X2 + a**, •••

, x

nt

< Snt < X

nt + dxntj

nt

2

(2i») exp{-![*1 + + (^-X

2

} + *" + 0 }

dx_ dx^, • a ’ dx , x 2 nt

Consequently,
' - > • . .. ; .
- hi -

I(n,t) - /exp {

- \ X S

fc2

j

f 1 2 ir 2 2

exn

- — kJi \ " 2 L^l

+ ( x2~x

i) + ••• + l n

dxi dx

nt .

(x , “X , i

)

nt nt-1

Integrals of the form

I = (TT)

m

2

m

Z a. . x.x.

, 1 , 0=1 10 1 0 dx^ dx

m

(where a. . is symmetric and positive definite) are easily evaluated by

10

transformation to principal axes. By a pure rotation of the form

y « B x (where - det B * 1) we can bring I into the form

The A’s are the elements of the diagonal matrix A - B A B~*^.

Hence

I = (det A

)

-1/2 (det (a. .)

V 10

Applying this result to l(n,t), we obtain

Let us consider the k x k determinant

Expanding according to the first row or column, we establish the recursion relation

D

k

(a,b) • a Dj

{_1

(a,b) - (a,b) .

If we define D * 1 and D- * b, this equation is true for k * 2*

o x We rewrite the recursion relation as Vl + k-1

a D, (k ^ 1)
.

'

’

-

: I

«

•

.

*
- h9 -

Utilizing the identities

cosh (k + 1) 9 + cosh (k - l) 0=2 cosh 9 cosh k 9

sinh (k + l) 0 + sinh (k - l) 9 = 2 cosh 9 sinh k 9

we can write

Dj (a,b) = A cosh k 9 + B sinh k 9

where

cosh 9 = a/2

To fit the starting conditions on D and Dn we set

o 1

A = 1 B =

Thus

I( n,« - [

cosh' (n t 9) +

2ba

/ a

2- It

1

sinh (n t 9) J

-1/2

where cosh 9=1+ 1/n * Solving for e , ve obtain

0

= 1 + l/n

2 + i \fz + l/'n

2

= l + 'i--+\+0 (n-3

)

1 n

As n —=> oo , en t 0

e"^^ « Therefore

lim I(n 5 t) = <^exp - J

x^CEjd'cJ’^ = £cosh t yfe
'

*

• • ,

• *
- 5o -

In estimating K

2

> we shall coniine ourselves to the case t > > 1,

x

t/n << 1. This is the case of interest for the harmonic oscillator, where

we took t = 5, n = 100. Therefore we neglect terras involving the factor

-t\/2

e

^ compared with terras involving e"^

v ^

. Thus

• /yp j

“ I

x

2 (T)dtjy> v'T

l(n,t)^ e

ntQ

2 rx* —i—

i

L

-1/2

-W2

We shall expand the difference e

^ in powers of l/n, retaining

the term of lowest order. It is convenient to deal with logarithms at

this point

The vanishing of terms of order n in the expression for 9 is significant, since it implies that the expression for log l(n,t) contains no

error term of order t/n. This is the computational reflection of our

qualitative conjecture that the error, to first order, involves only a shifting of the time origin, rather than a dilation of the time scale.
. . ' . 1 T.: . >• ' *•'
- 51 -

We obtain, then

log I(n,t) _ -VS _ JL. + ^

2

2\/Sn

-WS

I(n,.*) - s”^ Ti — - ^ 4. L 2\/5n

Hence

K

p

= -i— + 0(t/n'

x^ 2\/Sn

0(t/n2

) ,

0(t/n2

)l \/2

) .

- 52 -

APPENDIX B

FURTHER STATISTICAL PROPERTIES OF THE RANDOM WALKS

Since many properties of the hydrogen random walks are predictable

analytically, it seemed worthwhile to check the experimental results against

some of these predictions.

The most significant result obtained, and one which may represent a source of error, concerned the distribution of the end points r(l) of the

hydrogen walks. We recall

//Too T2 TToo TS 7Too \7

<« “

+ {v*i

z’->

(

100

Since \k»l XjJ , etc,, are Gaussian with mean 0 and variance 100, then

2 2 2 2 [r(l)r has the same distribution as X + Y + Z where X, Y, and Z are independent Gaussian random variables with mean 0 and variance 1.

Hence

2 2 2

n z p r ~ x +y * Prob / r(l) <^- 1 = (2lT) J J J e

l J 2 2 2 2

X +y +z

dx dy dz

= (27JT3^2 f e

2

U1T r

2

dr

o

2 oc - c ,

» = erf - —•— e = F(oc)

^tr

where the last line is obtained through an integration by parts.
;

.

• r "•••

•

. ~

: - ’

, ••

' ^

.

*

-

'
- 53 -

If we tabulate the end points r(l) of the first 900 random walks, the

expected number of walks with r(l) 4s* is 900 F(=*l) . Let us define

S(°0 * number of walks with r(l) (empirical)

It Is easy to show that

<[ ]> s(oc) - 900 F(=*-) I >

s (7- (ex.) = 900 F(o^)(l - F(cX-))

r

®<-) "I Furthermore, the probability that |s(®*) - 900F(C<-) | 2 c7_('>c)

J- is

3=K<^)J

r*681

is 4 .95 I

.

( .997

J

Hence, discrepancies of more than 2 cr(c*> ) between the empirical and theoretical distributions are to be regarded with suspicion. We compare the two

distributions in the table below*

DISTRIBUTION OF END POINTS OF 900 HYDROGEN WALKS

F(c=c) 900 F(<*) s(<*0 • 0• (pc) ,S(^)-900 F(«-) ,

1 {<=<

)

o*5 .030 27 28 5.1 0.2

1*0 .197 177 185 11.9 0.7

1.5 .U78 130 U27 i3.o 0.2

2.0 .739 665 6U5 13 .2 1.5

2.5 .900 810 753 9.0 6.3

3.0 .971 87U 836 5.0 7.6 •

3.5 ,99k 89U 867 2.3 11.7

U.o .999 899 880 1.0 19 .0

The tail of the empirical distribution S( evidently falls off

much less rapidly than that of a true Gaussian,, The magnitude of the
. • *'

. r .

• . , " ' -•

.

• V

-v . : 'i fi i: f, 'V i

' o »<flc

i

•

?

: v

'

'

,

%

•

, "

. •

.

•

' sjfc

• \

....

.

- •

• j

'

•

y
- 5U -

resultant error in and possible remedies, have not been investigated,

mainly because the SVJAC will not employ the Rand Corporation Gaussian

deviates. Qualitatively, the effect of the non-Gaussi an behavior is to increase the values of r(£), and thus to decrease the values of J

dt/rfC) *

- • o This may provide a partial explanation of another discrepancy between theory

and observation, discussed below.

Comparisons were also made between the theoretical and experimental

mean values of integrals arising in the hydrogen and helium computations.

Since the operations of averaging and summing (or integrating) commute, we

have

|

x at/r(-c)\> -

But

r(x)

- (2fTX)“3/2 f e-r2/2r

U7T

r

2

, /2 /l

r v IT V r

Hence

<f **«> - (

*

-1/2 dt

x/Trr

= 1.59

The empirical value for 1800 samples was 1.U3.

In order to appreciate the significance, if any, of this discrepancy,

2 p

1

one must know the variance O** of a single observation of I dT/rfC)

.

r, O

Knowing , we then know the variance of the average of 1800 observaf

° ^2 ^o2

tions , since o* =

.

t^qq
-

'

.

r

.

.

4
3y definition

^0

/,!

3

-' WoW

2 - ( P d't/r(r) - —i- v2 nr )

2

>- <(j>/W

2

>-^

a

If we think of f dX/r(?Z) as a Riemann sum of many non-independent random

±

variables , it is easily seen that

r?Eh;

M

\

d :/r\Z)

/ (2tr)3 J

0

Jc ?372 d\

dr

2 *

CD

f

6 r. ^ U rjz

*1 +7l+*l

2T

1

2 2

(xp-x,) +(y9-y

1

) +(s

2-z1

) 2 <Vv

2 2 >n 2 2

•00 /x, +y +Z-, /x0 +yn + z 9 V a. I. I V ^ £

dx, ••• ds0 J. u

= U log 2

1

Thus

<yf = I* log 2 - JL = .22

cr0

<r = -—2~~ - .on . \TTHoo

(In good agreement with the empirically

o observed <T^ " - .19*)

0

^Evaluation of th$ integral is left as an exercise for the reader.
«
- 56 -

The difference of *16 between the calculated ana .. crved value of

i

J

d'C/rff ) is therefore equal to l5 which is far t eyond the limits of

statistical inaccuracy. Dr. R* P. Feynman first called attention to a major systematic error, and gave the following method for estimating it e

r 1

,

, . What we actually sample is not the integral j d /r(''), but rather the

10( u

o Riemann sum *01

^ —

^loo5 Thus, the average to be considered is that of

the Riemann sum, which is somewhat smaller than that of the integral (the

error encountered here is the ’'discretization error*1 discussed in Section

It is geometrically evident that a fairly good estimate of the Riemann sum

in parentheses is

i

1.005’

oo5

ar

v'T

1.003 - .071 1.86

Thus <^.01 kJj_ ?

'(Y7iooy

^ = (compared with 1.5# for <^ | df/r(f^>)

and our observed value 1*U3 is in error by 5 cr, instead of 1 The remaining error may, perhaps, be accounted for by the non-Gaussian behavior

of the tail of the distribution of rfl) (discussed above), and the conse1

quent depression of values of f d'^/r ) .
’
- 57 -

One can also calculate the average value of the integral of the helium

integration term.

1

•\f(xg-xp 2

+ (^-yp ^(Zn-Zj ) 2

2 2 2

*1 *?1 * Z

1

2X

>y, +z,

2t dx-j • • • da,

The observed value was .93. The discrepancy is probably to be explained

along the same lines as above.

APPENDIX C

STATISTICAL METHOD FOR FINDING THE LOWEST EIGENFUNCTION

The reasoning in Section II leads easily to a statistical method for

finding (x) , the lowest eigenfunction. From Section II we have

~ A. ,t

P(x, t)dx - H C- .(x)e J dx . 0 J

Furthermore

P(x,t) cte «

Prob i particle survives till time t and is in (x,x + dx) at time t r

1

r i

J

Prob

'I particle picks a path such that x < x(t) < x + dx

j

x

Prob ^particle survives till t, given that x < x(t) < x + dx

J

But

Prob I particle picks a part such that

. x < x(t)

. < x + dx

It

l

v2Tf t

-x /2t . e ox $

and the probability of survival till t, given that x < x(t) < x + dx, is

just the average of the survival probability exp ^ /

V(x(£)dTj^ over

the set of paths such that x < x(t) < x + dx. Hence

P(x,t)dx = 1 -x

2

/2t x — .. A ' /3x <(eoq>\ -J V(x(t))dZ x < x( t) < x + dx

V2 Tf t >

The conditional, average in brackets can be found empirically by subdividing

the x-axis into short intervals, and computing the average of
'

.

• • .

'
P(x,t)dx for large t, we can thus determine b$ (up to a normalizing

factor) . By considering several different values of t and waiting until

the rate of time decay of P[x9 i) becomes the same at all points x, we can

ascertain when the space distribution has settled down to the correct eigenfunction. This settling-down process can be speeded by starting the walks

not from a delta-function distribution, but rather from a distribution which

approximates our best guess for (x) . Then the coefficients 0 o , •••

will be small •

No experiments have yet been performed with the above scheme, which

clearly requires very many random walks

.

If we start our walks from a point x^ instead of the origin, then

,t

P(x,t) = 2 ) Y.(x)e 3 o u J

and hence

N^(xo )e = j P(x,t) ^f/

1(x)dx = <Qaxp

j - V(x(x))dT

J* v|/^(x(t)^

If we have an empirical estimate of and A n from an earlier set of

walks, we may perform a new set of walks starting at x and use the last

xIf N paths x^(tr)(i =1, •••

, N) are sampled, an equivalent and computationally easier estimate of P is

P(x,t)dx = pT 2 exp* {

- I V(x.$))dt \ . N x<x,(tVx+dx L %

x J
.

,

.

jjpir
- 60 -

equation to improve our gueg^ for a.vx (x, ) . The osefnines s of th: pro' 5 - . cedars is doubtful, since it depends on the ,! approximate orthogonality” of

our final ‘^le true ''V2

> "¥y etc., and involves an inordinate

amount of effort to improve our guess for at a single point*
'
xbliography

K* Kao, On the distribution of certain Wiener functionals, Trans* Am#

Math. Soc. 65 (l9h9) PP* 1-13.

H. Kac and M. D# Donsker, A saiapling method for determining the lowest

eigenvalue and the principal eigenfunction of Schrodinger 1 s equation-,

Journal of Research of the National Bureau of Standards, UU (1950}

r-rs PP e I • H a Fortet, On the estimation of an eigenvalue by an additive functional

of a stochastic process, with special reference to the Kac-Donsker

method, to appear in Journal of Research of the National Bureau of

Standards

.

I. ¥. King, Stochastic Methods in Quantum Mechanics, Proceedings of

IBM Seminar on Scientific Computation, Nov# 19h? •

E. Kahn, Modification of the Monte Carlo Method, Proceedings of IBM

Seminar on Scientific Computation, Nov# 19U?

THE NATIONAL BUREAU OF STANDARDS

Functions and Activities

The functions of the National Bureau of Standards are set forth in the Act

of Congress, March 3, 1901, as amended by Congress in Public Law 619, 1950.

These include the development and maintenance of the national standards of

measurement and the provision of means and methods for making measurements

consistent with these standards; the determination of physical constants and

properties of materials; the development of methods and instruments for testing

materials, devices, and structures; advisory services to Government Agencies on

scientific and technical problems; invention and development of devices to serve

special needs of the Government; and the development of standard practices,

codes, and specifications. The work includes basic and applied research, development, engineering, instrumentation, testing, evaluation, calibration services

and various consultation and information services. A major portion of the

Bureau’s work is performed for other Government Agencies, particularly the

Department of Defense and the Atomic Energy Commission. The scope of activities is suggested by the listing of divisions and sections on the inside of the

front cover. Reports and Publications

The results of the Bureau’s work take the form of either actual equipment and

devices or published papers and reports. Reports are issued to the sponsoring

agency of a particular project or program. Published papers appear either in the

Bureau’s own series of publications or in the journals of professional and scientific

societies. The Bureau itself publishes three monthly periodicals, available from

the Government Printing Office: The Journal of Research, which presents complete papers reporting technical investigations; the Technical News Bulletin, which

presents summary and preliminary reports on work in progress; and Basic Radio

Propagation Predictions, which provides data for determining the best frequencies

to use for radio communications throughout the world. There are also five series

of nonperiodical publications: The Applied Mathematics Series, Circulars, Handbooks, Building Materials and Structures Reports, and Miscellaneous Publications.

Information on the Bureau’s publications can be found in NBS Circular 460,

Publications of the National Bureau of Standards ($1.00). Information on calibration services and fees can be found in NBS Circular 483, Testing by the National

Bureau of Standards (25 cents). Both are available from the Government Printing Office. Inquiries regarding the Bureau’s reports and publications should be

addressed to the Office of Scientific Publications, National Bureau of Standards,

Washington 25, D. C.
-