Part 2: Quadrature - From Photon Counts to Dark Matter Halos

Learning Outcomes¶

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

• Derive rectangular, trapezoidal, and Simpson’s rules from Taylor series

• Understand optimal point placement through Gaussian quadrature basics

• Apply Monte Carlo integration for high-dimensional problems

• Choose appropriate methods based on smoothness, dimension, and noise

• Analyze error scaling and predict convergence behavior

The Fundamental Challenge¶

While differentiation is a local operation (needing only nearby points), integration is global—we must consider the entire domain. Quadrature (numerical integration) transforms the continuous integral:

into a discrete sum:

The art lies in choosing:

• The points $x_i$ (where to sample)

• The weights $w_i$ (how much each sample contributes)

• The number $n$ (balancing accuracy vs computation)

Physical Motivation: Why Astronomers Integrate¶

Integration is everywhere in astronomy because we observe integrated quantities:

1. Spectroscopy: We don’t see individual photons but integrated flux

   $$F = \int_{\lambda_1}^{\lambda_2} F_\lambda d\lambda$$

   (3)

2. Photometry: CCD pixels integrate photons over exposure time

   $$N_\text{photons} = \int_0^T \Phi(t) dt$$

   (4)

3. Structure: Mass distributions require volume integrals

   $$M(<r) = \int_0^r 4\pi r'^2 \rho(r') dr'$$

   (5)

4. Cosmology: Distances involve integrating through expanding space

   $$d_L = (1+z) \int_0^z \frac{c\,dz'}{H(z')}$$

   (6)

Each context demands different accuracy, and the function properties (smooth vs oscillatory, bounded vs singular) determine the best method.

Building Integration Methods: From Rectangles to Optimality¶

The Riemann Sum Foundation¶

The definition of the Riemann integral suggests the simplest approximation:

where $x_i^* \in [x_i, x_{i+1}]$ and $\Delta x_i = x_{i+1} - x_i$.

Different choices of $x_i^*$ give different methods:

• Left endpoint: Rectangle rule (left)

• Right endpoint: Rectangle rule (right)

• Midpoint: Midpoint rule

• Average of endpoints: Trapezoidal rule

Method 1: Rectangle Rule - The Starting Point¶

Why we care: The simplest possible approximation — helps understand error analysis that applies to all methods.

Mathematical Formulation¶

For uniform spacing $h = (b-a)/n$:

Error Analysis via Taylor Series¶

On each subinterval $[x_i, x_{i+1}]$:

Evaluating:

The rectangle rule uses only $f(x_i)h$, so the local error is:

Over $n$ intervals:

This is a first-order method—halving $h$ halves the error.

Implementation¶

def rectangle_rule(f, a, b, n):
    """Left rectangle rule"""
    h = (b - a) / n
    total = 0
    for i in range(n):
        x = a + i * h
        total += f(x)
    return h * total

Why Not to Use Rectangle Rule¶

• Only first-order accurate

• Systematically over/underestimates for monotonic functions

• Poor for oscillatory functions

• Mainly pedagogical value

Method 2: Trapezoidal Rule - The Workhorse¶

Why we care: The workhorse of experimental data integration—robust and simple.

Geometric Intuition¶

Instead of rectangles, connect consecutive points with straight lines, creating trapezoids. The area of each trapezoid averages the function values at its endpoints.

Mathematical Formulation¶

The area of a trapezoid with parallel sides $f(x_i)$ and $f(x_{i+1})$ and height $h$ is:

Summing over all intervals:

Error Analysis¶

Using Taylor expansion around the midpoint $x_i + h/2$:

The trapezoidal approximation on $[x_i, x_{i+1}]$:

Since the exact integral is $hf(x_i + h/2) + O(h^3)$:

Local error: $E_{local} = -\frac{h^3}{12}f''(\xi)$

Global error: $\boxed{E_{total} = -\frac{(b-a)h^2}{12}f''(\xi)}$ for some $\xi \in [a,b]$

This is a second-order method!

Common Failure Modes and Fixes¶

Method 3: Simpson’s Rule - Parabolic Perfection¶

Why we care: When your function is smooth and you want high accuracy with moderate computational cost.

The Key Insight¶

Instead of linear interpolation (trapezoids), use quadratic interpolation (parabolas) through consecutive triplets of points.

Derivation via Interpolation¶

For three equally-spaced points $x_0, x_1, x_2$ with $x_1 = x_0 + h$ and $x_2 = x_0 + 2h$, the unique parabola through $(x_i, f(x_i))$ is:

Integrating this parabola from $x_0$ to $x_2$:

These are the famous Simpson weights: 1, 4, 1.

Composite Simpson’s Rule¶

For $n$ intervals (must be even):

The pattern of weights is:

• First point: 1

• Odd indices (1,3,5,...,n-1): 4

• Even indices (2,4,6,...,n-2): 2

• Last point: 1

Error Analysis and the Fourth-Order Miracle¶

Through careful Taylor series analysis, the local error for Simpson’s rule on interval $[x_i, x_{i+2}]$ is:

Global error: $\boxed{E_{total} = -\frac{(b-a)h^4}{180}f^{(4)}(\xi)}$

Simpson’s rule is fourth-order accurate despite using only parabolas!

Why Simpson’s Rule is Fourth-Order¶

The “miracle” of Simpson’s fourth-order accuracy comes from symmetry. When we integrate the interpolating parabola, the $f'''(\xi)$ term (which would give $O(h^3)$ error) has odd symmetry about the midpoint and integrates to zero:

This cancellation is automatic—we get fourth-order accuracy “for free” by using symmetric intervals! Since the error term involves $f^{(4)}$, and the fourth derivative of any cubic polynomial is zero, Simpson’s rule integrates cubic polynomials exactly.

When Simpson’s Rule Shines¶

• Smooth functions with continuous fourth derivative

• Periodic functions

• When high accuracy is needed with moderate $n$

• Functions from physical models (usually smooth)

Common Failure Modes and Fixes¶

Method 4: Gaussian Quadrature - Optimal Point Placement¶

Why we care: Sometimes being clever about WHERE you sample beats taking MORE samples.

The Revolutionary Idea¶

All previous methods used equally-spaced points. But what if we could choose both the points AND weights optimally?

Gauss’s insight: For $n$ points, we can make the method exact for all polynomials up to degree $2n-1$ by choosing specific locations. These points are the roots of Legendre polynomials, which have special orthogonality properties that maximize accuracy.

Example: 2-Point Gaussian Quadrature¶

Instead of evaluating at equally-spaced points, Gaussian quadrature uses special points. For 2 points on $[-1,1]$:

• Points: $x_1 = -\frac{1}{\sqrt{3}} \approx -0.577$, $x_2 = \frac{1}{\sqrt{3}} \approx 0.577$

• Weights: $w_1 = w_2 = 1$

This seemingly arbitrary choice makes the method exact for all cubic polynomials!

In Practice¶

You won’t derive Gauss points yourself. Instead:

from scipy.special import roots_legendre

# Get n Gauss points and weights for [-1,1]
points, weights = roots_legendre(n)

# Transform to general interval [a,b]
x_gauss = (b-a)/2 * points + (a+b)/2
w_gauss = (b-a)/2 * weights

# Compute integral
integral = sum(w * f(x) for x, w in zip(x_gauss, w_gauss))

When Gaussian Quadrature Matters¶

• When you control where to sample (not experimental data)

• For very smooth functions

• When function evaluations are expensive (fewer points needed)

• As the foundation for spectral methods in advanced simulations

Key Takeaway: This principle—optimal placement beats more samples—appears throughout computational physics, from choosing timesteps in orbit integration to placing grid points in PDE solvers.

Method 5: Monte Carlo Integration - When Dimensions Explode¶

Why we care: The only practical method for high-dimensional integrals that appear in statistical mechanics, quantum mechanics, and Bayesian inference.

Understanding Dimensions in Integration¶

When we talk about “dimensions” in Monte Carlo integration, we mean the number of integration variables, not spatial dimensions:

• 1D integral: $\int_a^b f(x) dx$ (one integration variable)

• 2D integral: $\int \int_R f(x,y) dxdy$ (two integration variables)

• 6N-D integral: Phase space integral over N particles with 3 position and 3 momentum components each

For example, computing the partition function for 100 atoms requires a 300-dimensional integral (3 position coordinates per atom)!

The Curse of Dimensionality¶

For a $d$-dimensional integral with $n$ points per dimension, deterministic methods need $n^d$ evaluations:

• 1D: 100 points

• 2D: 100² = 10,000 points

• 3D: 100³ = 1,000,000 points

• 10D: 100¹⁰ = 10²⁰ points (impossible!)

This exponential growth makes grid-based methods impossible for $d > 5$.

The Monte Carlo Solution¶

Randomly sample $N$ points $\vec{x}_i$ in the domain $\Omega$:

where $V(\Omega)$ is the volume of the domain.

Error Analysis¶

For the Monte Carlo estimator $\hat{I} = \frac{V(\Omega)}{N}\sum_{i=1}^N f(\vec{x}_i)$:

Step 1: The expected value of our estimator is the true integral:

Step 2: The variance of the estimator is:

where $\sigma_f^2 = \text{Var}[f] = \mathbb{E}[f^2] - \mathbb{E}[f]^2$ is the variance of $f$ over the domain.

Step 3: The standard error (standard deviation of the estimator) is:

Key insight: Error $\propto N^{-1/2}$ independent of dimension $d$!

This assumes:

• Independent random samples

• Finite variance of the integrand ($\sigma_f < \infty$)

• Sufficient samples for CLT to apply (typically N > 30)

When Monte Carlo Wins¶

Above ~4 dimensions, Monte Carlo becomes superior!

Common Failure Modes and Fixes¶

Adaptive Methods and Richardson Extrapolation¶

Adaptive Quadrature¶

When functions vary rapidly in some regions but are smooth elsewhere:

def adaptive_simpson(f, a, b, tol):
    """Recursively refine where needed"""
    # Compute with n and 2n points
    I1 = simpson(f, a, b, n=10)
    I2 = simpson(f, a, b, n=20)
    
    if abs(I2 - I1) < tol:
        return I2
    else:
        mid = (a + b) / 2
        left = adaptive_simpson(f, a, mid, tol/2)
        right = adaptive_simpson(f, mid, b, tol/2)
        return left + right

Richardson Extrapolation¶

If your method has error $I(h) = I_{exact} + Ch^p + O(h^{p+1})$, compute at two resolutions:

This cancels the leading error term! For trapezoidal rule (p=2):

This gives fourth-order accuracy from a second-order method!

Real Application: Galaxy Luminosity from Spectrum¶

Astronomers measure galaxy spectra as discrete samples at wavelengths $\lambda_i$ with fluxes $F_i$. The total luminosity requires integration:

Challenges¶

1. Irregular wavelength spacing from different instruments

2. Noise in flux measurements

3. Missing data from atmospheric absorption

4. Emission lines need high resolution

Solution Strategy¶

For irregular spacing, use trapezoidal rule with variable h:

def integrate_spectrum(wavelengths, fluxes):
    """Integrate spectrum with irregular spacing"""
    total = 0
    for i in range(len(wavelengths) - 1):
        dλ = wavelengths[i+1] - wavelengths[i]
        f_avg = (fluxes[i+1] + fluxes[i]) / 2
        total += dλ * f_avg
    return total

For emission lines: ensure ~10 points across line width or fit Gaussian profile and integrate analytically.

Choosing the Right Integration Method¶

Method Comparison Summary¶

Practical Guidelines¶

1. Start simple: Try trapezoidal first—it often suffices

2. Check smoothness: Plot your function before choosing

3. Test convergence: Double n and see if result changes

4. Consider cost: Is f(x) expensive to evaluate?

5. Remember libraries: Use scipy.integrate.quad for production

Debugging Integration¶

Common integration bugs and their diagnosis:

Bridge to Part 3: The Unity of Numerical Methods¶

You’ve now mastered two fundamental classes of problems: finding where functions equal zero and measuring areas under curves. These aren’t isolated techniques—they’re connected by deep mathematical principles.

In Part 3, we’ll explore these connections. You’ll see how root finding and integration are inverse operations, how the same error analysis principles govern both, and how to combine methods for complex problems. Most importantly, you’ll develop the intuition to choose the right method for any computational challenge.

The synthesis ahead will prepare you for the dynamic problems of Module 3, where these static methods become building blocks for simulating evolution through time.

Next: Part 3 - Synthesis