Monte-Carlo Integration

Although Matlab’s trapz and quad are useful functions,physics problems often demand surface or volume integrals. Some physicists even claim that the structure of matter at the most fundamental level is 11 dimensions and perform integrals over dimensions that are hard to even imagine. It’s safe to say that at some point you’ll want an easy way to integrate over multiple dimensions. Matlab provides several 2D integrators, but it’s sometimes better to learn one method that extends to integrals in any number of dimensions. Monte-carlo methods extend naturally to any number of dimensions, so let’s start with one dimension and then go to more dimensions.

MC Integrals in 1D:

In one dimension, an integral can be written. .

To get some practice, let’s choose a=0, b=20, and f(x) =x2. Instead of picking the x values carefully with even spacing, dx, this time we’ll just randomly choose some number of x values using the rand function.

clear all; close all;

N = 1000000

x = rand(1,N);

plot(x(1:200),'.')

The x values of the first 200 indices is plotted above. Each xi has a random value between 0 and 1 which is displayed on the y-axis in the left-hand plot. This use of random samples and its relation to the statistics of gambling was the original inspiration for the name Monte Carlo.

The right-hand plot shows the original values scaled between our integration limits, a = 0 and b = 20. Here’s how we scale what comes out of the rand function to our interval, a to b.

a = 0; b = 20;

x = a + (b-a)*rand(1,N);

plot(x(1:200),'.');

Because the x values are randomly distributed, we can define the average value of the function:

The integral is the average value of the function multiplied by the total interval, (b-a).

.

The figure on the right shows red markers representing the x2 function at a series of randomly chosen x values. The blue line shows the average of the function, <f>. The area below the red markers is approximately the same as the area below the blue line. The area below the blue line is (b-a)*<f>. In Matlab the integral is simply:

Integral_MC = (b-a)/N*sum(x.^2)

How does this compare to Matlab’s fancy built-in routines? It takes a lot more samples, but eventually gets reliable results. Generally the error in the integral will go as the statistical error in the sampling, i.e. error ~. To get 3 significant figures, we’ll need N ~ 103x2 samples. The good news is that this error stays the same no matter how many dimensions are involved in the integral.

More Dimensions:

In more dimensions, the integration takes the form

What is important here is that we can still find the average of the function and then multiply it by the total volume: . In 3D, we can integrate over a box:.

N = 1e6;

L = 1; W = 1; H = 1;

x = L*rand(1,N); y = W*rand(1,N); z=H*rand(1,N);

Integral_MC = (L*W*H/N)*sum(sqrt(x.^2+y.^2+z.^2))

If our volume was spherical, we would use spherical coordinates:.

This time . Notice that we had to average the function multiplied by r2sin. To demonstrate, let’s integrate the function,.

N = 10000;

phiMax = 2*pi; thetaMax = pi; rMax = 1;

phi = phiMax*rand(1,N);

theta = thetaMax*rand(1,N);

r=rMax*rand(1,N);

Integral_MC = (phiMax*thetaMax*rMax/N)* ...

sum( (cos(phi)./r.^3).*r.^2.*sin(theta))

Electric potential of a plane with constant charge density:

Let’s consider an example, the potential from a plane of charge as shown in the figure.

The integral comes from summing all the potentials from the chunks of “point” charges. Just as we summed two charges for a dipole, now we’ll sum many charges. The important part is to understand how dQ, modeled as a point charge is related to dA where  is the charge density = charge/area (). Once the physics is understood, then we can consider doing the integral numerically.

sigma = 1e-6; % charge density is 1.0 microC/m^2

N = 1e6;

x0 = 2; % P(x0,y0)

y0 = 3;

xi = -1 + 2*rand(1,N); % -1<x<1 integration limits

yi = 0 + 2*rand(1,N); % 0<y<2 integration limits

% electric potential at P from a chunk of charge at xi, yi

Vi = pointPotential(x0-xi, y0-yi, sigma);

% integrate Vi for chunks of charge in x and y directions. Vtot = (2*2/N)*sum(Vi) % range is 2 in x and y.

Assignment: M8_MC_Integration

1. Use Monte-Carlo integration to verify that in spherical coordinates.

2.Compute the electric potential, V, from a disk, centered on the origin with radius, 1.2 cm, and charge density,  = 2 nC/m2.

The integral is:

3. Integrals of mass density lead to three interesting quantities:

where M is the total mass, is the center of mass position, and Icm is the moment of inertia about the center of mass, and the mass density, , may be a function of the special variables.

Consider a rectangular box: length=0.2 m, width = 0.2 m, and height = 1.0 m centered on the origin, (0,0,0) and with a mass density,for y and z in meters. Use numerical integration to find M, , and Icm. Think of a way to test the algorithm using a constant density. Verify that this gives reasonable results. Determine the number of significant figures in the computation. Report M, , and Icm with a reasonable number of significant digits.

(Be careful with Icm. The radial distance is computed in a plane perpendicular to the rotation axis. You’ll need to find Ixy around the z-axis, Iyz around the x-axis, and Ixz around the y-axis.)

4. Consider half a sphere, radius = 0.2 m, sitting below the x-y plane, as shown. Let the density be a constant,. Use numerical integration to find M,, and Icm.

Think of a way to test the algorithm using a complete sphere. After this is working, change the integration limit to get half a sphere.

5. The earth can be modeled as a sphere with a variable density profile. The internal densities of the earth are primarily known from measurements of earthquakes.

layer / radial range
(km) / density
(g/cm3)
Inner core / 0 -1,221 / 13.1 – 12.8
Outer core / 1222 - 3480 / 12.2 - 9.9
Lower mantle / 3481 - 5621 / 5.6 – 4.4
Upper mantle / 5622 - 6341 / 4.4 – 3.4
crust / 6342 - 6371 / 2.9 – 2.2

Determine the mass and moment of inertia about an axis through the center of the earth. (radius = 6371 km)

Compare your result with the “mean density” calculation.  = 5.52 g/cm3)

Page 1