Change of variables and applications of multiple integrals

Lesson 5 evaluated double and triple integrals in rectangular, polar, cylindrical and spherical coordinates, and used them for mass, centroids and the second moment of area. This lesson completes the picture. It shows that polar, cylindrical and spherical coordinates are special cases of a general change of variables governed by the Jacobian determinant, which lets any suitable region be mapped to a simple one. It then applies multiple integrals to three more engineering questions: the force and centre of pressure of water on a dam gate, the mass moment of inertia of rotating parts such as flywheels, and probabilities for two random quantities at once, such as positioning errors in two directions. It ends with Monte Carlo integration, the numerical method that handles integrals in awkward shapes and many dimensions.

The Jacobian and change of variables

In one variable, substitution \(x = g(u)\) changes \(dx\) into \(g'(u)\,du\): the derivative measures how much a small interval is stretched. In two variables, a transformation \(x = x(u, v)\), \(y = y(u, v)\) maps a small rectangle \(du \times dv\) in the uv-plane to a small, nearly parallelogram-shaped patch in the xy-plane. Its area is \(|J|\,du\,dv\), where J is the Jacobian determinant:

This is the determinant of the matrix of the linear approximation, and MATH-125 shows that a determinant measures how a linear map scales area. The change of variables formula follows:

where S is the region in the uv-plane that maps onto R. For polar coordinates, \(x = r\cos\theta\), \(y = r\sin\theta\), the Jacobian is \(\cos\theta \cdot r\cos\theta - (-r\sin\theta)\sin\theta = r\), which is where the familiar \(r\,dr\,d\theta\) comes from. In three dimensions the Jacobian is a 3 × 3 determinant; for spherical coordinates it is \(\rho^2\sin\phi\).

Hydrostatic force on a submerged surface

Water pressure increases with depth, \(p = \rho g h\), so the force on a submerged gate or wall is an integral of pressure over its area, and the force acts lower than the centroid. For a plane vertical surface with depth h measured downward from the free surface:

Here \(\bar{h}\) is the depth of the centroid and \(I_c\) the second moment of area about a horizontal axis through the centroid (Lesson 5). The second form uses the parallel axis theorem. The depth \(h_{cp}\) of the centre of pressure is where a single equivalent force would act, which sets the moment that hinges and actuators must resist.

Mass moments of inertia

The resistance of a body to angular acceleration about an axis is its mass moment of inertia, the integral of mass times squared distance from the axis:

It plays the role of mass in rotation: torque equals \(I\alpha\), and kinetic energy is \(\tfrac{1}{2}I\omega^2\). It differs from the second moment of area of Lesson 5, which is a property of a cross-section used in bending, although the integrals look alike.

Composite parts and the parallel axis theorem

Machine parts are usually built from simple shapes, and their moments of inertia add, as long as all are taken about the same axis. For a piece whose own centroidal moment is \(I_c\) but which sits a distance d from the axis of rotation, the parallel axis theorem gives \(I = I_c + Md^2\), the mass version of the theorem used for areas in Lesson 5. It follows from expanding \(r_\perp^2\) about the centroid: the cross term integrates to zero because the first moment about the centroid vanishes.

Probability in two variables

If two random quantities X and Y have a joint probability density \(f(x, y)\), the probability that the pair falls in a region R is the double integral of f over R, and the density integrates to 1 over the whole plane. Expected values are integrals too: \(E[g(X, Y)] = \iint g\,f\,dA\). STAT-235 develops these ideas; here the point is that the methods of Lessons 5 and 6 apply directly.

The most famous example needs polar coordinates. The integral \(I = \int_{-\infty}^{\infty} e^{-x^2}dx\) has no elementary antiderivative, but its square is a double integral over the plane:

So \(I = \sqrt{\pi}\), which is the normalizing constant of the normal distribution, the most widely used model of measurement error.

Monte Carlo integration

For awkward shapes or integrals in many dimensions, nested integration becomes impractical. Monte Carlo integration estimates an integral by random sampling: scatter N random points uniformly through a simple box of volume \(V_{box}\) that contains the region, and count the fraction that fall inside. The volume is approximately that fraction times \(V_{box}\), and an integral of f is approximately \(V_{box}\) times the average of f over the points that land inside (counting f as zero outside). The error falls only as \(1/\sqrt{N}\), slowly, but independently of the number of dimensions, which is why the method dominates in radiation transport, financial risk and reliability analysis.

import numpy as np

rng = np.random.default_rng(1)
exact = 16 / 3   # volume common to two perpendicular cylinders of radius 1

for n in [1_000, 100_000, 10_000_000]:
    p = rng.uniform(-1, 1, size=(n, 3))          # box from -1 to 1 in each direction
    inside = (p[:, 0]**2 + p[:, 2]**2 <= 1) & (p[:, 1]**2 + p[:, 2]**2 <= 1)
    estimate = 8 * inside.mean()                 # box volume times fraction inside
    print(f"N = {n:>10,}: volume = {estimate:.4f}, error = {abs(estimate - exact):.4f}")
print(f"Exact value: {exact:.4f}")

N =      1,000: volume = 5.2800, error = 0.0533
N =    100,000: volume = 5.3373, error = 0.0039
N = 10,000,000: volume = 5.3326, error = 0.0007
Exact value: 5.3333

The exact value, 16/3, can be found with a careful triple integral, but the Monte Carlo estimate needed no thought about limits of integration. Multiplying the number of samples by 100 cuts the typical error by a factor of about 10, as the \(1/\sqrt{N}\) rule predicts; individual runs vary, because the error is itself random.