Homework 1; due on Tuesday, March 16
PHY336, Computational Physics, Spring 2021
Department of Physics, SUSTech
1. MONTE CARLO INTEGRATION
A sphere of radius r1 consists of two different materials, with densities ρ1 and ρ2. The material
with density ρ2 is located within a cylinder of radius r2, as illustrated in Fig. 1, and the material
of density ρ1 fills up the rest of the sphere. Write a program that calculates the two moments
of inertia of this sphere corresponding to rotation about the z and x axis. The inner cylinder is
centered around the z-axis, as also shown in the figure.
Figure 1: A solid sphere of radius r1 with an inner solid cylinder of radius r2. The cylinder consists
of a material with density ρ2; the rest of the sphere consists of a material with density ρ1.
Carry out the calculation using Monte Carlo sampling of the moment of inertia integral, (1)
where r⊥(x, y, z) is the perpendicular distance of the point (x, y, z) from the axis of rotation (here
the x or z axis, giving Ix and Iz). Enclose the sphere in a box with side L = 2r1 in order to
easily do the calculation using (x, y, z) points. Use the 64-bit linear congruential random-number
generator that you tested in homework assignment 1 (reading the seed integer from a file seed.in),
with the factor included to convert the integers to double-precision numbers in the range [0, 1) (the
generator with this factor included is available on the course web site).
The program should read the following input from a file read.in:
r1,r2,rho1,rho2,npt,nbi
1
where r1,r2 are the two radii (in m), rho1,rho2 are the densities (in kg/m3
), npt is the number of
random points generated per “bin” (for which bin averages are computed) and nbi is the number
of bins (on the basis of which the final average and statistical error are computed).
Bin averages should be computed for both the Iz and Ix moments of inertia and these should be
writen to a file bin.dat containing nbi lines, each with the bin number followed by the Iz and Ix
values (write these averages to the file after each bin is completed; it is not necessary to store the
data in the program). The final average and error bar (standard deviation of the mean) computed
using the bin averages should be written to a file res.dat.
As a specific case, do the calculation for a copper (8930 kg/m3
) sphere of radius 5 cm with an inner
gold (19320 kg/m3
) cylinder of radius 1 cm. Use 106 points per bin (npt) and do the calculation
for nbi=50,500,5000. For the final case, construct a histogram of the bin averages (with the width
of the histogram bins chosen in a reasonable way to get of the order tens of histogram bins with
significant weight). The report on this problem needs to contain only the final numerical results
(averages and standard deviations) for the three runs and the histogram for the last run. Comment
on the shape of the histogram.
2
2. Calculate the volume of a d (d=2,3,4) dimensional sphere of radius r=1 using Monte
Carlo (MC). Give an average value and an error estimate as a function of MC points.