Monte Carlo integration on a sphere

I meant to embed a Mathematica notebook into last week’s post about generating random points on the surface of a sphere, but I was having trouble with the Wolfram Cloud. Those problems have been fixed, so here’s a Mathematica notebook that does the various calculations necessary to generate the points and to produce some of the images included in that post:

If you scroll down to the end of the notebook, you’ll notice that in addition to plotting the random points to see if they look like they’re distributed uniformly over the surface, I also use the points to estimate the moments of inertia of the sphere about each of the three axes and compare those estimates with the analytical solution. The estimate is done through a technique called Monte Carlo integration.

The moment of inertia of a thin shell about the x-axis is

Ixx=Aρ(y2+z2)dA

where A is the surface of the sphere and ρ is the areal density (mass per unit area) of the shell material. We’ll use a unit density, so

Ixx=A(y2+z2)dA

There are similar equations for the moments of inertia about the other axes, Iyy and Izz.

The idea behind Monte Carlo integration is to generate a set of random points over the integration domain and add up the weighted sum of the integrand at all those points. The weights are the proportion of the domain represented by each point. In our case, the domain of integration is the surface of the sphere, we are generating n points on that surface, and each point represents one nth of the total surface area. So

Ixxi=1n4πR2n(yi2+zi2)=4πR2ni=1n(yi2+zi2)

Similarly,

Iyy4πR2ni=1n(xi2+zi2)

and

Izz4πR2ni=1n(xi2+yi2)

These are the equations you see in the notebook for ixxEstimate, iyyEstimate, and izzEstimate. The Total function does the summation.

We get the exact moments of inertia by rewriting the equation in terms of θ and φ,

Ixx=π/2π/2ππ[(Rcosϕsinθ)2+(Rsinϕ)2]R2cosϕdθdϕ

which leads to

Ixx=83πR4

Because of the symmetry of the sphere, the moments of inertia about the other two axes, Iy and Izz, have the same formula, which is called iAnalytical in the notebook.

The Monte Carlo estimates are all within one percent of the exact value, which is another indication that the points are indeed uniformly distributed over the surface. Each time the notebook is reevaluated, the estimates will change (because the random points will change), but the estimates will always be close to the exact value.

You might think these estimates are pretty poor, considering we used 5000 points to get them. And you’d be right. Monte Carlo integration tends to be much less efficient than other methods of numerical integration when the dimension of the problem is low. But the number of integrand evaluations needed by those other methods typically rises exponentially with the number of dimensions in the domain, something that doesn’t happen with Monte Carlo integration. So the Monte Carlo method is a good tool to have when you need to integrate over a high-dimensional domain. And it’s a decent way to check your work when efficiency isn’t a concern.