Applications of Double Integrals

Learning Objectives

By the end of this section, you will be able to:

1. Compute areas of planar regions using double integrals by integrating the constant function 1
2. Calculate volumes under surfaces and between surfaces using double integration
3. Find mass, center of mass, and moments of laminas (thin flat plates) with variable density
4. Compute moments of inertia and understand their physical significance in rotational dynamics
5. Determine surface area of surfaces defined by $z = f(x,y)$
6. Apply double integrals to probability, computing probabilities from joint density functions
7. Connect these concepts to machine learning, including expected values, feature averaging, and integration over parameter spaces

The Big Picture: From Abstraction to Reality

"Double integrals transform abstract mathematical regions into measurable physical quantities—area, volume, mass, and probability."

In the previous sections, we learned how to compute double integrals over rectangular and general regions. But why do we care about these calculations? The power of double integrals lies in their applications—they provide the mathematical machinery to compute quantities that matter in physics, engineering, probability, and machine learning.

Each application follows the same pattern: we identify a quantity that can be broken into infinitesimal pieces, express each piece as a function times $dA$, and then integrate to get the total. Whether we are computing the mass of a thin plate, the volume of a solid, or the probability of an event, the underlying principle is the same.

Historical Context

The applications of double integrals emerged alongside the development of calculus itself. Isaac Newton and Gottfried Leibniz both understood that integration could compute areas and volumes, but the systematic treatment of multiple integrals came later.

Leonhard Euler (1707-1783) made extensive use of double integrals in mechanics and probability. His work on the rotation of rigid bodies required computing moments of inertia—quantities that depend on the distribution of mass relative to an axis.

Joseph-Louis Lagrange (1736-1813) developed the mathematical foundations of mechanics, where double and triple integrals appear constantly. His "Analytical Mechanics" treats the motion of extended bodies using integral calculus throughout.

In probability theory, Pierre-Simon Laplace (1749-1827) and Carl Friedrich Gauss (1777-1855) used double integrals to study two-dimensional probability distributions. The bivariate normal distribution, so important in statistics and machine learning, requires double integration for almost every calculation.

Area of a Region

The simplest application of double integrals: computing the area of a planar region. The idea is straightforward—if we integrate the constant function 1 over a region $D$, we count up all the infinitesimal area elements:

This might seem trivial, but it gives us a powerful tool: we can compute areas of regions that would be difficult to find by other methods. The same double integral machinery we use for more complex functions works for the constant function.

For a Type I region (vertical slices) with $a \le x \le b$ and $g_1(x) \le y \le g_2(x)$:

This recovers the familiar formula for area between curves!

Interactive: Area Visualizer

Explore how double integration computes the area of various regions. Watch how we sum up the widths of horizontal strips (Type II) or vertical strips (Type I):

Volume Under a Surface

When $f(x, y) \ge 0$, the double integral gives the volume of the solid region between the surface $z = f(x, y)$ and the xy-plane, bounded laterally by the region $D$:

This is the natural generalization of the single-variable integral formula for area under a curve. Where a single integral sums infinitesimal rectangles of height $f(x)$ and width $dx$, a double integral sums infinitesimal rectangular prisms of height $f(x,y)$ and base area $dA$.

Interactive: Volume Visualizer

Explore the volume under various surfaces. Toggle the Riemann columns to see how we approximate the volume with rectangular prisms:

Mass and Density

Consider a thin flat plate (a lamina) occupying a region $D$ in the xy-plane. If the density at each point is $\rho(x, y)$ (mass per unit area), then the total mass is:

For a uniform lamina where $\rho(x,y) = \rho$ is constant, the mass simplifies to $M = \rho \cdot \text{Area}(D)$. But for non-uniform density, we must integrate.

Center of Mass

The center of mass (or centroid for uniform density) is the point $(\bar{x}, \bar{y})$ where we could balance the lamina on a pin. It is computed using the first moments:

Interactive: Mass and Center of Mass

Explore how density distribution affects mass and center of mass. See how non-uniform density shifts the center of mass away from the geometric center:

Moments of Inertia

The moment of inertia measures how difficult it is to rotate an object about a given axis. It depends not just on total mass, but on how that mass is distributed relative to the axis. Mass far from the axis contributes more because it has a larger lever arm.

The perpendicular axis theorem for laminas states:

This makes sense geometrically: the squared distance from the origin is the sum of squared distances from the two coordinate axes.

Interactive: Moments of Inertia

Explore how distance from the axis affects the moment of inertia. Notice how mass far from the axis (brighter regions) contributes more:

Surface Area

To find the surface area of a surface $z = f(x, y)$ over a region $D$, we need to account for how the surface "stretches" relative to its projection onto the xy-plane.

Consider a small patch of the surface. When projected onto the xy-plane, it covers area $dA$. But the actual surface area is larger because the surface is tilted. The ratio depends on the partial derivatives $f_x$ and $f_y$:

The derivation comes from computing the cross product of the tangent vectors:

Interactive: Surface Area

Visualize how the surface area element varies across the surface. Toggle tangent patches to see the local tangent planes and normal vectors:

Probability Applications

In probability theory, double integrals compute probabilities for continuous bivariate distributions. If $f(x, y)$ is a joint probability density function (PDF), then:

Key properties of joint PDFs:

• Non-negativity: $f(x, y) \ge 0$ everywhere
• Normalization: $\iint_{-\infty}^{\infty} f(x, y) \, dA = 1$
• Marginal PDF: $f_X(x) = \int_{-\infty}^{\infty} f(x, y) \, dy$
• Expected value: $E[g(X, Y)] = \iint g(x, y) f(x, y) \, dA$

Interactive: Joint Probability

Explore different joint probability distributions and see how integrating over different regions changes the probability:

Machine Learning Connections

Double integral applications appear throughout machine learning:

1. Expected Values and Loss Functions

The expected loss over a data distribution is a double (or higher) integral:

This is why we use empirical risk minimization: we approximate the integral with a sum over training samples.

2. Bayesian Integration

In Bayesian inference, we often need to integrate over parameter space:

This predictive distribution integrates out the parameters, averaging predictions over all possible parameter values weighted by their posterior probability.

3. Kernel Methods and Inner Products

In kernel methods, we work with inner products in function space:

Python Implementation

Let's implement the key applications in Python:

Mass, Center of Mass, and Moments

1import numpy as np
2from scipy import integrate
3
4def double_integral_area(f_lower, f_upper, y_min, y_max, n=1000):
5    """
6    Compute area of region bounded by x = f_lower(y), x = f_upper(y)
7    from y = y_min to y = y_max.
8
9    Area = integral from y_min to y_max of [f_upper(y) - f_lower(y)] dy
10    """
11    y = np.linspace(y_min, y_max, n)
12    width = f_upper(y) - f_lower(y)
13    area = np.trapz(width, y)
14    return area
15
16def volume_under_surface(f, x_bounds, y_bounds, nx=100, ny=100):
17    """
18    Compute volume under z = f(x,y) over rectangular region.
19
20    Volume = double integral of f(x,y) dA
21    """
22    x = np.linspace(x_bounds[0], x_bounds[1], nx)
23    y = np.linspace(y_bounds[0], y_bounds[1], ny)
24    dx = x[1] - x[0]
25    dy = y[1] - y[0]
26
27    X, Y = np.meshgrid(x, y)
28    Z = f(X, Y)
29
30    volume = np.sum(Z) * dx * dy
31    return volume
32
33def center_of_mass_2d(rho, x_bounds, y_bounds, nx=100, ny=100):
34    """
35    Compute center of mass of a lamina with density rho(x,y).
36
37    x_bar = (1/M) * integral of x * rho(x,y) dA
38    y_bar = (1/M) * integral of y * rho(x,y) dA
39    """
40    x = np.linspace(x_bounds[0], x_bounds[1], nx)
41    y = np.linspace(y_bounds[0], y_bounds[1], ny)
42    dx = x[1] - x[0]
43    dy = y[1] - y[0]
44
45    X, Y = np.meshgrid(x, y)
46    R = rho(X, Y)
47
48    mass = np.sum(R) * dx * dy
49    moment_x = np.sum(X * R) * dx * dy
50    moment_y = np.sum(Y * R) * dx * dy
51
52    return {
53        'mass': mass,
54        'x_bar': moment_x / mass,
55        'y_bar': moment_y / mass
56    }
57
58def moment_of_inertia(rho, x_bounds, y_bounds, axis='origin', nx=100, ny=100):
59    """
60    Compute moment of inertia of a lamina.
61
62    axis='x': I_x = integral of y^2 * rho dA
63    axis='y': I_y = integral of x^2 * rho dA
64    axis='origin': I_0 = integral of (x^2 + y^2) * rho dA = I_x + I_y
65    """
66    x = np.linspace(x_bounds[0], x_bounds[1], nx)
67    y = np.linspace(y_bounds[0], y_bounds[1], ny)
68    dx, dy = x[1] - x[0], y[1] - y[0]
69
70    X, Y = np.meshgrid(x, y)
71    R = rho(X, Y)
72
73    if axis == 'x':
74        return np.sum(Y**2 * R) * dx * dy
75    elif axis == 'y':
76        return np.sum(X**2 * R) * dx * dy
77    else:  # origin
78        return np.sum((X**2 + Y**2) * R) * dx * dy
79
80# Example: Compute properties of uniform unit disk
81def uniform_disk(x, y):
82    """Density = 1 inside unit disk, 0 outside"""
83    return np.where(x**2 + y**2 <= 1, 1.0, 0.0)
84
85# Center of mass of uniform disk (should be at origin)
86result = center_of_mass_2d(uniform_disk, [-1, 1], [-1, 1])
87print(f"Mass: {result['mass']:.4f} (expected: π ≈ 3.1416)")
88print(f"Center of mass: ({result['x_bar']:.4f}, {result['y_bar']:.4f})")
89
90# Moment of inertia of unit disk
91I_0 = moment_of_inertia(uniform_disk, [-1, 1], [-1, 1], axis='origin')
92print(f"Polar moment I_0: {I_0:.4f} (expected: π/2 ≈ 1.5708)")

Probability Applications

1import numpy as np
2from scipy import stats
3from scipy.integrate import dblquad
4
5def joint_probability(pdf, x_range, y_range):
6    """
7    Compute P((X,Y) in region) by integrating joint PDF.
8
9    pdf: joint probability density function f(x,y)
10    x_range, y_range: tuples (min, max) for integration
11    """
12    result, error = dblquad(
13        pdf,
14        y_range[0], y_range[1],  # y limits
15        lambda y: x_range[0],    # x lower limit
16        lambda y: x_range[1]     # x upper limit
17    )
18    return result, error
19
20def marginal_pdf_x(joint_pdf, y_support, x_value, n=1000):
21    """
22    Compute marginal PDF of X: f_X(x) = integral of f(x,y) dy
23    """
24    y = np.linspace(y_support[0], y_support[1], n)
25    integrand = np.array([joint_pdf(x_value, yi) for yi in y])
26    return np.trapz(integrand, y)
27
28def expected_value_xy(joint_pdf, x_support, y_support, func='x'):
29    """
30    Compute E[g(X,Y)] = integral of g(x,y) * f(x,y) dA
31
32    func: 'x' for E[X], 'y' for E[Y], 'xy' for E[XY], 'x2' for E[X²]
33    """
34    if func == 'x':
35        g = lambda x, y: x
36    elif func == 'y':
37        g = lambda x, y: y
38    elif func == 'xy':
39        g = lambda x, y: x * y
40    elif func == 'x2':
41        g = lambda x, y: x**2
42    else:
43        raise ValueError(f"Unknown function: {func}")
44
45    integrand = lambda y, x: g(x, y) * joint_pdf(x, y)
46    result, _ = dblquad(
47        integrand,
48        x_support[0], x_support[1],
49        lambda x: y_support[0],
50        lambda x: y_support[1]
51    )
52    return result
53
54# Example: Bivariate normal distribution
55mu_x, mu_y = 0, 0
56sigma_x, sigma_y = 1, 1
57rho = 0.5  # correlation coefficient
58
59def bivariate_normal(x, y):
60    """Joint PDF of bivariate normal with correlation rho"""
61    z = ((x - mu_x)/sigma_x)**2 - 2*rho*((x - mu_x)/sigma_x)*((y - mu_y)/sigma_y) + ((y - mu_y)/sigma_y)**2
62    norm = 2 * np.pi * sigma_x * sigma_y * np.sqrt(1 - rho**2)
63    return np.exp(-z / (2*(1 - rho**2))) / norm
64
65# Probability that X > 0 AND Y > 0 (first quadrant)
66prob, err = joint_probability(bivariate_normal, (0, 5), (0, 5))
67print(f"P(X > 0, Y > 0) = {prob:.4f} (expected: {0.25 + np.arcsin(rho)/(2*np.pi):.4f})")
68
69# Expected values
70E_X = expected_value_xy(bivariate_normal, (-5, 5), (-5, 5), 'x')
71E_XY = expected_value_xy(bivariate_normal, (-5, 5), (-5, 5), 'xy')
72print(f"E[X] = {E_X:.4f} (expected: 0)")
73print(f"E[XY] = {E_XY:.4f} (expected: ρ = {rho})")

Surface Area Computation

1import numpy as np
2from scipy.integrate import dblquad
3
4def surface_area(f, f_x, f_y, x_bounds, y_bounds, region=None):
5    """
6    Compute surface area of z = f(x,y) over a region.
7
8    S = integral of sqrt(1 + f_x^2 + f_y^2) dA
9
10    f: function z = f(x,y)
11    f_x, f_y: partial derivatives
12    region: optional function (x,y) -> bool for non-rectangular regions
13    """
14    def integrand(y, x):
15        if region is not None and not region(x, y):
16            return 0.0
17        dzdx = f_x(x, y)
18        dzdy = f_y(x, y)
19        return np.sqrt(1 + dzdx**2 + dzdy**2)
20
21    result, error = dblquad(
22        integrand,
23        x_bounds[0], x_bounds[1],
24        lambda x: y_bounds[0],
25        lambda x: y_bounds[1]
26    )
27    return result, error
28
29# Example 1: Paraboloid z = x² + y² over unit disk
30f = lambda x, y: x**2 + y**2
31f_x = lambda x, y: 2*x
32f_y = lambda x, y: 2*y
33disk = lambda x, y: x**2 + y**2 <= 1
34
35area, err = surface_area(f, f_x, f_y, [-1, 1], [-1, 1], disk)
36exact = np.pi * (5*np.sqrt(5) - 1) / 6
37print(f"Paraboloid surface area: {area:.4f} (exact: {exact:.4f})")
38
39# Example 2: Hemisphere z = sqrt(1 - x² - y²)
40def hemisphere_area():
41    f = lambda x, y: np.sqrt(np.maximum(0, 1 - x**2 - y**2))
42
43    # Partial derivatives (handle boundary carefully)
44    def f_x(x, y):
45        val = 1 - x**2 - y**2
46        return -x / np.sqrt(val) if val > 0.01 else 0
47
48    def f_y(x, y):
49        val = 1 - x**2 - y**2
50        return -y / np.sqrt(val) if val > 0.01 else 0
51
52    # Use smaller disk to avoid singularity at boundary
53    disk = lambda x, y: x**2 + y**2 <= 0.99
54
55    area, err = surface_area(f, f_x, f_y, [-1, 1], [-1, 1], disk)
56    return area
57
58hemisphere = hemisphere_area()
59print(f"Hemisphere surface area: {hemisphere:.4f} (exact: 2π ≈ {2*np.pi:.4f})")
60
61# Example 3: Saddle z = x² - y² over unit square
62f = lambda x, y: x**2 - y**2
63f_x = lambda x, y: 2*x
64f_y = lambda x, y: -2*y
65
66area, err = surface_area(f, f_x, f_y, [-1, 1], [-1, 1])
67print(f"Saddle surface area over [-1,1]²: {area:.4f}")

Test Your Understanding

Summary

Double integrals transform abstract mathematical computation into concrete physical and probabilistic quantities. The key applications covered in this section are:

Key Takeaways

1. Unifying pattern: All applications follow the same structure—integrate an appropriate function over a region
2. Physical interpretation: The integrand represents what we are accumulating (height, mass, probability, etc.)
3. Moments measure distribution: First moments give center of mass; second moments give moment of inertia
4. Surface area uses derivatives: The factor $\sqrt{1 + f_x^2 + f_y^2}$ accounts for surface tilt
5. ML ubiquity: Expected values, Bayesian marginalization, and many other ML concepts are fundamentally double integrals

Coming Up Next: In the next section, we extend to three dimensions with Triple Integrals. You'll learn to compute volumes and masses of solid regions, where the setup mirrors what we've done here but with an additional dimension of integration.