W e conclude this chapter with a few remarks concerning the M o n t e Carlo method:
1. D o not use the M o n t e Carlo method if any other method is available.
It is a method of last resort.
2. Use variance reduction if and only if it is really necessary. Statistical estimates of the- variance using computed results should be considered only as suggestive of the true variance. Great caution should be exercised in attributing precision to a result merely from low variance of past histories, since the variance estimate is an asymptotic function of the number of histories. U s e known expected values wherever possible.
3. Particles must be followed for many nuclear lifetimes to remove any hereditary influences of the initial population. T h e fluctuations within one lifetime will be much smaller than those over many lifetimes, and these small fluctuations can give a very false impression of reliability and accuracy.
276 V I . T H E M O N T E C A R L O M E T H O D
4. Check the coded program by all possible methods, and check it again.
T h e checking of M o n t e Carlo coded programs is especially difficult (1) because of the intensely complex logic, (2) because this logic is not deterministic, (3) because events wrongly calculated may be sufficiently rare to hide the error and yet large enough to influence the results, and (4) because the statistical nature of the calculations masks such errors.
Several suggestions are offered to make any errors obvious and to isolate them.
1. Check each individual subroutine before assembling it into the larger routine. Check as many aggregates of subroutines as possible.
2. Run a series of simple, extreme problems so designed as to bound the class of all possible problems for which the M o n t e Carlo will be used.
Artificial substances in which neutrons experience only one type of nuclear event are particularly useful to check the routines with respect to that event. Everything should be checked, such as tallies, the velocities and locations of particles, the conservation of neutrons. Some-times the injection of particles at one point with one velocity will facilitate precise checks; at other times statistical checks on velocity distributions, for example, on a reasonable sample of particles must be run. Input data should be as troublesome to the computer as possible.
Run problems involving simple and complex geometries in which the matter consists of a vacuum to check the crossing of boundaries. I f possible, observe the trajectories of the particles on an on-line oscilloscope to check these routines. Boundary crossings are particularly troublesome because of round-off error within the machine. W h e n a particle hits a boundary, design the boundary crossing routine such that a small amount is added to its coordinates to force it across regardless of this round-off.
References
The literature on the Monte Carlo method is surprisingly sparse. There are many articles and chapters devoted to the topic, but no comprehensive review exists to the authors' knowledge. For reasonably broad coverage of the method, References 1-5 are recommended. Certain special topics such as random walk, random number generation, etc. are discussed in References 6-11. References 12 and 13 are two very readable texts on statistics, including the central limit theorem.
/. Cashwell, E. D., and Everett, C. J., " T h e Monte Carlo Method for Random Walk Problems." Pergamon, N e w York, 1959.
2. Kahn, H., in "Symposium on Monte Carlo Method" ( H . A. Meyer, ed.).
Wiley, N e w York, 1956.
3. Kahn, H., Applications of Monte Carlo. Rand Report A E C U - 3 2 5 9 (1954).
4. Goertzel, G., and Kalos, Μ . H., Monte Carlo method in transport problems, in "Progress in Nuclear Physics," Vol. 2, pp. 315-369. Pergamon Press, 1958.
5. Mayne, A. J., "Symposium on Monte Carlo Methods," pp. 103, 123, 176, 249.
Wiley, N e w York (1953).
6. Bauer, W . F., The Monte Carlo method. J. Soc. Ind. Appl. Math. 6, 438-451 (1954).
7. Curtis, J. H., Sampling methods applied to differential and difference equa
tions, in "Seminar on Scientific Computation." International Business Machines Corp., Inc., Nov. 1949.
8. Kahn, H., Random sampling (Monte Carlo) techniques in neutron attenua
tion problems. Nucleonics 6 (1950).
9. Kahn, H., Modification of the Monte Carlo method, in "Seminar on Scientific Computation." International Business Machines Corp., Inc., Nov. 1949.
10. Troost, M . , "Study of Neutron Transport by Monte Carlo Methods."
Sc.D. Thesis, Massachusetts Inst. Technol., Cambridge, Massachusetts, June, 1958.
11. Moshman, J., The generation of pseudorandom numbers on a decimal calculator. J. Assoc. Comp. Mach. 1; 88.
12. Hoel, P. C , "Introduction to Mathematical Statistics." Wiley, N e w York, 1954.
13. Cramer, H., "Mathematical Methods of Statistics." Princeton Univ. Press, Princeton, N e w Jersey, 1946.
1. A linear interpolation between tabular points F(x}) and F{Xj^x) assumes the random variable χ is uniformly distributed over the interval to Xj. If F(x) is known to be concave upward or downward more accurate parabolic interpolation may be used. Show that for F(x) concave upward we have
where κ is a random number. What are the advantages and disadvantages of the parabolic interpolation methods ?
2. Prove Eq. (6.3.16).
3. Show that the direction cosines of a random direction may be selected by choosing a point randomly within a semicircle and are given by
Problems
and for F(x) concave downward
2
2 7 8 V I . T H E M O N T E C A R L O M E T H O D
where 171 and η2 are acceptable random numbers, scaled between 0 and 1 and between — 1 and + 1, respectively, if and only if
vl + vt
< 1and where ηζ is any random number scaled between —1 and + 1 . Hint:
Consider the trigonometric relations sin 2<p and cos 2<p, and sin ψ and cos φ.
Show that the efficiency of selection of the random numbers is 78.5%.
4. By selecting a point at random that lies within a unit sphere, show that the cosines ocd , βα , yd of a random direction are given by
Show that the efficiency of selecting the cosines of this random direction is 52.4%.
8. Suppose the probability that y has a value in the range Ay centered at y is f{y)Ay) where/(y) = 3y2 and suppose that 0 ^ y ^ 1. Show that y = η/3,
where 0 ^ 17 ^ 1.
9. The differential probability distribution of y is f(y) = 4/π(\ + y2) , where 0 ^ y ^ 1. Show that y = tan 77-17/4, where 0 ^ η ^ 1.
10. The differential probability distribution of y is f(y) = 2/rr V y( l + y).
Devise a rejection technique for choosing y from this distribution. Find the efficiency.
11. Suppose y = η\ where 0 < 77 < 1. Show that the probability that y has a value in the interval Ay centered at y is
1 . . .
12. Prove Eq. (6.4.24).
13. Prove Eqs. (6.6.11) and (6.6.12).
14. Flow chart a Monte Carlo routine in which isotropic elastic scattering, anisotropic elastic scattering, capture, inelastic scattering, (w, 2n) events, and fissions are each considered individually.
15. Devise a boundary crossing routine for zones consisting of concentric spheres only, taking round-off error into account.
16. Devise a boundary crossing routine for zones consisting of only coaxial cylinders infinitely long, taking round-off error into account.
17. W e have coded a program that has no provision for recording the coordinates of the particles that leak out of a reactor but that has a provision for recording the coordinates of each neutron in each zone within the reactor. Devise a simple method for finding the spectrum and angular distribution of neutrons that leak out of the reactor.
18. Program a subroutine to calculate the logarithm of a random number according to the rejection technique discussed.
19. Program a Monte Carlo routine in which isotropic elastic scattering and anisotropic elastic scattering are considered as one nuclear event and capture,
fission, (w, 2n) events, inelastic scattering are considered as another event.
20. Devise a generalized boundary crossing routine for zones consisting of spheres, cylinders, planes, cylindrical cones, elliptical cones, elliptical cylinders, and ellipsoids. Hint: Consider that such boundaries can be described by a formula of the type
ax2 + by2 + c(z - z{))2 - Κ = 0. ( 1 )
Show that the distance to any boundary is given by
— e - f j Λ/ e2 — hr
and j 1 if the expression (1) changes from negative to positive as the neutron crosses the boundary and j = — 1 if the expression (1) changes from positive to negative. A neutron is inside the a zone if jr > 0 for any of the boundaries.
21. Program a Monte Carlo routine that takes delayed neutron emitters into account. Hint: Devise a routine to consider representative times, then extra
polate past results. Calculate what the sources will be and use them to inject neutrons.
22. Suppose we have an infinite plane slab used as a shield. The attenuation of neutrons through the slab is very great.
280 V I . T H E M O N T E C A R L O M E T H O D
(a) Program a routine that will automatically locate particle splitting planes stably.
(b) Program a routine that will calculate accurately the number of neutrons that leave the slab through a small spot on the surface.
23. Discuss the merits and faults of the following criteria for locating splitting planes:
(a) The planes are so located that the number of deaths in the zone on the left of the prospective position, i.e., the space between itself and the most recently placed boundary, is one-half the original number.
(b) The prospective boundary is to be so placed that the number of deaths to its right will be half the number of particles that emerge from the preceding boundary.
(c) The prospective boundary is so placed that as many more deaths above the initial number of particles injected occur in the region to its right as occurred short of the required number to the right of the preceding boundary.
All of these schemes for locating splitting planes are inferior to the one mentioned in the text.