Numerical Integration

Numerical quadrature using trapezoidal rule, Simpson's formulas, and adaptive integration methods.

What it is

Numerical integration, historically termed numerical quadrature, is the computational framework that evaluates definite integrals: $I = \int_a^b f(x) dx$ when the integrand $f(x)$ lacks an elementary closed-form anti-derivative (such as $e^{-x^2}$, $\frac{\sin x}{x}$, or $\sqrt{1 - k^2 \sin^2\theta}$), or when $f(x)$ is supplied as discrete experimental data points from sensors.

Quadrature algorithms split into two primary architectures:

1. Newton-Cotes composite rules (equispaced grid spacing $h = \frac{b-a}{N}$):
   • Composite Trapezoidal Rule: Connects adjacent samples via linear chords. Global truncation error decays at second order: $I_{\text{trap}} = \frac{h}{2} \left[ f(a) + 2 \sum_{i=1}^{N-1} f(x_i) + f(b) \right], \qquad E_{\text{trunc}} = -\frac{(b-a) h^2}{12} f''(\xi) = \mathcal{O}(h^2)$
   • Composite Simpson's 1/3 Rule (even $N$): Fits quadratic parabolic arcs across successive point triplets. Global truncation error decays at fourth order: $I_{\text{Simp}} = \frac{h}{3} \left[ f(x_0) + 4 \sum_{i=1,3,\dots}^{N-1} f(x_i) + 2 \sum_{i=2,4,\dots}^{N-2} f(x_i) + f(x_N) \right], \qquad E_{\text{trunc}} = -\frac{(b-a) h^4}{180} f^{(4)}(\xi) = \mathcal{O}(h^4)$
   • Romberg Integration: Applies Richardson extrapolation across trapezoidal estimates with step sizes $h, h/2, h/4, \dots$, systematically cancelling $h^2, h^4, h^6$ error terms.
2. Gaussian Quadrature (optimal non-uniform collocation nodes):
   • Gauss-Legendre Quadrature: Chooses evaluation nodes $x_i$ and weights $w_i$ based on the roots of orthogonal Legendre polynomials $P_n(x)$. With just $n$ function evaluations, it integrates any polynomial of degree up to $2n - 1$ exactly: $\int_{-1}^1 f(t) dt \approx \sum_{i=1}^n w_i f(t_i)$

Why it appears in physics

Integration is fundamentally intertwined with physical conservation and accumulation laws:

• Classical mechanics: Mechanical work performed by variable forces $W = \int_{x_1}^{x_2} F(x) dx$; center-of-mass coordinates and rotational moments of inertia require integrating continuous density fields $\rho(\mathbf{r})$.
• Electromagnetism: Electrostatic potentials and magnetic vector potentials around realistic conductor geometries require integrating Coulomb and Biot-Savart kernels that lack elementary primitives.
• Thermodynamics and statistical physics: The Stefan-Boltzmann law ($u = a T^4$) stems from integrating Planck's spectral density across all frequencies $\int_0^\infty \frac{x^3}{e^x - 1} dx = \frac{\pi^4}{15}$; quantum partition functions and Debye lattice heat capacities require definite integrals.
• Nonlinear dynamics: The exact period of a large-amplitude simple pendulum is a complete elliptic integral of the first kind $T = 4 \sqrt{\frac{L}{g}} \int_0^{\pi/2} \frac{d\phi}{\sqrt{1 - \sin^2(\theta_0/2) \sin^2\phi}}$, impossible to express in basic trigonometric terms.

Visual or geometric idea

Geometrically, the definite integral defines the signed area enclosed between $f(x)$ and the horizontal coordinate axis over $[a, b]$.

• Trapezoidal quadrature tiles the area with adjoining vertical trapezoids capped by straight chords.
• Simpson's 1/3 rule tiles the area with parabolic segments that mold closely to local curvature.
• Gauss-Legendre quadrature clusters evaluation nodes toward interval boundaries where polynomial approximation error typically flares, maximizing accuracy without globally refining uniform grid density.

Minimal symbolic language

1. Composite Trapezoidal Rule: $T(h) = h \left[ \frac{f(a) + f(b)}{2} + \sum_{i=1}^{N-1} f(a + i h) \right]$
2. Composite Simpson's 1/3 Rule: $S(h) = \frac{h}{3} \left[ f(x_0) + 4 \sum_{i=1,3,\dots}^{N-1} f(x_i) + 2 \sum_{i=2,4,\dots}^{N-2} f(x_i) + f(x_N) \right]$
3. Romberg recurrence: $R_{k,j} = \frac{4^j R_{k, j-1} - R_{k-1, j-1}}{4^j - 1}$
4. Coordinate mapping to $[-1, 1]$ for Gauss-Legendre: $x = \frac{b - a}{2} t + \frac{a + b}{2}, \qquad dx = \frac{b - a}{2} dt$ $\int_a^b f(x) dx = \frac{b - a}{2} \int_{-1}^1 f\left( \frac{b-a}{2}t + \frac{a+b}{2} \right) dt$

What the result means

The computed integral provides the accumulated physical quantity (joules of work, volts of electrostatic potential, seconds of period, coulombs of charge). Asymptotic convergence guarantees that doubling grid intervals $N$:

• Under Trapezoidal integration, error divides by $2^2 = 4$.
• Under Simpson's rule, error divides by $2^4 = 16$.
• With 5-point Gauss-Legendre quadrature, smooth integrands frequently achieve 12 to 14 accurate decimal digits using just 5 function calls.

When not to use it

1. Integrands with endpoint singularities: If $f(x) \to \infty$ at boundary points (e.g. $\int_0^1 \frac{1}{\sqrt{x}} dx$), closed Newton-Cotes schemes crash from division by zero at $x=0$. Algebraic regularization or open Gauss-Jacobi quadrature must be used.
2. Highly oscillatory integrands: If $f(x) = g(x) \sin(\omega x)$ with rapid phase $\omega \gg 1$, standard formulas require millions of grid nodes to sample each cycle ($h \ll 2\pi/\omega$). Filon or asymptotic methods are mandatory.
3. High-dimensional integrals ($D > 4$): In statistical mechanics or quantum Feynman path integrals, grid quadrature suffers the curse of dimensionality ($N^D$ operations). In high dimensions, Monte Carlo integration is the only viable path.

Operational summary

• Inspect integrand $f(x)$: is it smooth and continuous? Are endpoint singularities present?
• For tabulated laboratory data on uniform intervals, apply composite Simpson's rule (even $N$) or the trapezoidal rule (odd $N$ or irregular spacing).
• For continuous analytic expressions, deploy Gauss-Legendre or adaptive Simpson quadrature for maximum accuracy per function evaluation.
• Perform grid doubling tests ($N$ vs $2N$) to verify convergence and extract a posteriori error estimates via Richardson extrapolation.
• Ensure the final value incorporates composite physical dimensions $[f] \cdot [x]$.

Procedure for Numerical Quadrature in Physics

Before calculating

Before evaluating numerical integrals:

• Inspect integration domain $[a, b]$ for vertical poles ($f(x) \to \pm \infty$) or step discontinuities.
• For experimental data, check whether spatial or temporal sampling $\Delta x$ is strictly uniform or irregular.
• For analytic functions with parity symmetry (even or odd), exploit domain reductions prior to computing.

Step 1: identify the mathematical object

Select the quadrature scheme matching input data characteristics:

• Equispaced tabulated data with even subinterval count $N$: Composite Simpson's 1/3 rule.
• Tabulated data with non-uniform spacing: Non-uniform composite Trapezoidal rule.
• Smooth continuous analytic equations: Gauss-Legendre quadrature.

Step 2: choose the operation

For composite Simpson's 1/3 rule with $N$ subintervals ($N$ even, $N+1$ nodes):

1. Compute step size: $h = \frac{b - a}{N}$.
2. Generate grid points: $x_i = a + i h$ for $i = 0, 1, \dots, N$.
3. Group terms weighted by Newton-Cotes pattern ($1, 4, 2, 4, 2, \dots, 4, 1$).

Step 3: preserve units and signs

• Integrand values $f(x_i)$ carry dimensions $[f]$.
• Step size $h$ carries dimensions of coordinate variable $[x]$.
• The integral carries dimensions $[f] \times [x]$ (e.g. force in newtons times position in meters produces joules).
• Preserve signs under reversed limits: $\int_b^a f(x) dx = -\int_a^b f(x) dx$.

Step 4: calculate

1. Simpson 1/3 evaluation: $I(h) = \frac{h}{3} \left[ f(a) + 4 \sum_{k=1}^{N/2} f(x_{2k-1}) + 2 \sum_{k=1}^{N/2 - 1} f(x_{2k}) + f(b) \right]$
2. Refined evaluation with doubled grid ($2N$ subintervals): Compute $I(h/2)$ using twice as many nodes.
3. A posteriori error estimation (Richardson extrapolation): $E_{\text{approx}} \approx \frac{I(h/2) - I(h)}{15}$ High-order refined estimate: $I_{\text{refined}} = I(h/2) + \frac{I(h/2) - I(h)}{15}$

Step 5: interpret the result

• Contrast error bound $E_{\text{approx}}$ against tolerance targets ($10^{-6}$ for engineering, $10^{-12}$ for precision science).
• If error exceeds tolerance, double interval density again ($4N$) until stable convergence is confirmed.
• Format the final value with appropriate physical units.

Self-check checklist

• Is subinterval count $N$ strictly even for Simpson's 1/3 rule?
• Did I scale by $h/3$ rather than $h/2$?
• Do dimensions of the answer strictly match $[f] \cdot [x]$?
• Did I confirm the absence of hidden zero-denominators inside $[a, b]$?
• Does error divide by approximately 16 when doubling nodes in Simpson's rule?

Guided Example: Oscillation Period of a Large-Angle Nonlinear Pendulum

Physical Situation

A simple pendulum of length $L = 1.000 \text{ m}$ in gravitational field $g = 9.80665 \text{ m/s}^2$ is released from rest at large amplitude $\theta_0 = 60.0^\circ = \frac{\pi}{3} \text{ rad}$. The small-angle approximation ($\sin\theta \approx \theta$) predicts: $T_0 = 2\pi \sqrt{\frac{L}{g}} = 2.006409 \text{ s}$ For $\theta_0 = 60^\circ$, nonlinearity delays oscillation. Exact energy conservation leads to the complete elliptic integral: $T = 4 \sqrt{\frac{L}{g}} \int_0^{\pi/2} \frac{d\phi}{\sqrt{1 - k^2 \sin^2\phi}}, \qquad k = \sin(30^\circ) = 0.500$ We compute exact period $T$ using composite Simpson's 1/3 rule with $N = 4$ intervals.

Mathematical Translation

Front constant: $C = 4 \sqrt{\frac{1.000}{9.80665}} \approx 1.277319 \text{ s}$ Integrand on $[0, \pi/2]$: $f(\phi) = \frac{1}{\sqrt{1 - 0.25 \sin^2\phi}}$ With $N = 4$, interval width is $h = \frac{\pi}{8} \approx 0.3926991 \text{ rad}$. Nodes: $\phi = [0, \pi/8, \pi/4, 3\pi/8, \pi/2]$.

Step-by-Step Calculation

1. Nodal Evaluations:
   • $\phi_0 = 0$: $f(\phi_0) = 1.0000000$
   • $\phi_1 = \pi/8$: $\sin^2(\pi/8) \approx 0.1464466 \implies f(\phi_1) \approx 1.0188048$
   • $\phi_2 = \pi/4$: $\sin^2(\pi/4) = 0.5000000 \implies f(\phi_2) = \frac{1}{\sqrt{0.875}} \approx 1.0690450$
   • $\phi_3 = 3\pi/8$: $\sin^2(3\pi/8) \approx 0.8535534 \implies f(\phi_3) \approx 1.1297629$
   • $\phi_4 = \pi/2$: $\sin^2(\pi/2) = 1.0000000 \implies f(\phi_4) = \frac{1}{\sqrt{0.75}} \approx 1.1547005$
2. Simpson Composite Weighting: $I = \frac{h}{3} [f_0 + 4f_1 + 2f_2 + 4f_3 + f_4]$ $\text{Sum} = 1.0000000 + 4(1.0188048) + 2(1.0690450) + 4(1.1297629) + 1.1547005$ $\text{Sum} = 1.0000000 + 4.0752192 + 2.1380900 + 4.5190516 + 1.1547005 = 12.8870613$ $I = \frac{0.3926991}{3} \times 12.8870613 = 0.1308997 \times 12.8870613 \approx 1.6869125$
3. Period Computation: $T = 1.277319 \times 1.6869125 \approx 2.154725 \text{ s}$

Reading the Result

The large-angle period is $T = 2.155 \text{ s}$, which is $7.4\%$ longer than the harmonic estimate ($2.006 \text{ s}$). The exact analytical value is $2.15473 \text{ s}$, demonstrating 5-figure precision with only 4 Simpson intervals.

Common Error Avoided

Using the trapezoidal rule without checking periodic boundary derivatives or using a coarse Riemann sum that drastically underpredicts peak pendulum slowdown near release.

Inverse Example: Deducing Total Injected Beam Charge from Current Pulses

Physical Situation

An accelerator beam Faraday cup records current $I(t)$ during a rapid pulsed laser extraction: $I(t) = I_0 \left(\frac{t}{\tau}\right) e^{-t/\tau}, \quad \tau = 2.0 \text{ ns}, \quad I_0 = 5.0 \text{ A}$ Total extracted charge is $Q = \int_0^\infty I(t) dt$.

Given Result

Integrated charge measured: $Q = 20.0 \text{ nC}$.

Analytical Reconstruction

Integrating by parts: $\int_0^\infty t e^{-t/\tau} dt = \tau^2$ $Q = \frac{I_0}{\tau} \tau^2 = I_0 \tau = (5.0 \text{ A}) \times (2.0 \times 10^{-9} \text{ s}) = 1.0 \times 10^{-8} \text{ C} = 10.0 \text{ nC}$ If measured charge is $20.0 \text{ nC}$, actual peak current or pulse duration was double nominal.

What Physical Question It Can Answer

Does pulse integration provide total particle count? Yes: $N_e = Q / e = 10.0 \times 10^{-9} / 1.602 \times 10^{-19} \approx 6.24 \times 10^{10}$ electrons per bunch.

Boundary Verification

As $\tau \to 0$, $Q \to 0$; as $I_0 \to 0$, $Q \to 0$.

Common Mistakes in Numerical Quadrature

1. Running Simpson's 1/3 rule on an odd number of subintervals

It seems reasonable because users assume quadrature routines automatically average available intervals. It fails because Simpson's 1/3 rule groups intervals in pairs to fit parabolic arcs across triplets $(x_0, x_1, x_2)$, strictly requiring an even number of subintervals $N$ (odd number of nodes $N+1$). If $N$ is odd, the final panel is left unpaired. It is detected when software throws array indexing exceptions or applies invalid weight coefficients to endpoint panels. It is corrected by checking $N$ parity; if $N$ is odd, apply Simpson's 1/3 across the first $N-1$ panels and trapezoidal (or Simpson's 3/8) on the final subinterval. Mini-example: applying Simpson 1/3 to 5 subintervals ($N=5$, 6 nodes) leaves one subinterval without a parabolic partner, producing $10\%$ truncation errors.

2. Evaluating closed Newton-Cotes formulas on boundary singularities

It seems reasonable because improper integrals frequently converge to finite values (e.g. $\int_0^1 x^{-1/2} dx = 2$). It fails because closed rules (standard Trapezoidal and Simpson) evaluate the integrand at exact endpoints $a$ and $b$. At $x=0$, $1/\sqrt{0}$ triggers division by zero. It is detected via runtime exceptions for division by zero or infinite values. It is corrected by algebraic variable substitutions ($x = u^2$), adopting open quadrature formulas (midpoint rule), or using Gauss quadrature schemes that never sample boundary limits. Mini-example: for $\int_0^1 x^{-1/2} dx$, substitution $x = t^2$ yields smooth $\int_0^1 2 dt = 2$.

3. Discretizing rapid oscillatory integrands on coarse grids

It seems reasonable because if an integrand oscillates symmetrically between $-1$ and $+1$, intuition expects self-cancellation. It fails because if step size exceeds the oscillation period ($h > \lambda_{\text{wave}}$), spatial aliasing occurs: nodes systematically land on crests or troughs, generating massive false areas instead of zero cancellation. It is detected when doubling nodes swings the computed integral wildly in sign and magnitude. It is corrected by enforcing at least 10 to 20 nodes per full cycle ($h \leq \frac{2\pi}{10 \omega}$) or applying Filon quadrature. Mini-example: integrating $\int_0^{100\pi} \sin(10 x) dx$ with step $h = \pi/5$ evaluates $\sin(2\pi k) = 0$ everywhere, returning $I=0$ by pure arithmetic accident, but slight step shifts yield uncontrolled chaotic values.

4. Confusing local single-panel truncation with global composite order

It seems reasonable because elementary derivations prove single-panel Simpson error is $\mathcal{O}(h^5)$. It fails because spanning interval $[a, b]$ sums $N = (b-a)/h$ panels. Summing $N$ terms of order $\mathcal{O}(h^5)$ loses one power of $h$: $\sum_{i=1}^N \mathcal{O}(h^5) \approx \frac{b-a}{h} \mathcal{O}(h^5) = \mathcal{O}(h^4)$. It is detected when someone assumes doubling nodes slashes error by $2^5 = 32$ rather than the actual factor $2^4 = 16$. It is corrected by tracking global asymptotic orders: $\mathcal{O}(h^2)$ for composite Trapezoidal and $\mathcal{O}(h^4)$ for composite Simpson. Mini-example: doubling nodes in composite Simpson integration reduces error by 16 times, never 32.

5. Overlooking accumulated roundoff on ultra-dense meshes

It seems reasonable to assume $N = 10^8$ trapezoidal panels delivers near-infinite accuracy. It fails because summing $10^8$ floating-point terms in single precision (Float32) compounds roundoff noise that overwhelms truncation gains, degrading overall accuracy. It is detected when benchmark error stops dropping and begins rising with increasing $N$. It is corrected by adopting higher-order algorithms (Gauss-Legendre, Romberg) achieving target precision with dozens of nodes, or using Float64 with Kahan compensated summation. Mini-example: summing ten million trapezoidal slices in Float32 accumulates $\sim 10^{-4}$ roundoff drift, destroying single-precision accuracy.

Quick detection rule

• Is $N$ odd for Simpson's rule? Invalid: add one node or use 3/8 rule on the final panel.
• Is the function oscillating rapidly? Check if $h < \pi / \omega$. If not, your integral is aliased noise.

Physical Applications of Numerical Quadrature

1. Electromagnetism: electrostatic potential from continuous charge geometries

Appears in high-voltage engineering and accelerator physics when computing electrostatic potential $V(\mathbf{r})$ generated by metallic electrodes or continuous ion beams at arbitrary field coordinates via Coulomb's volume integral: $V(\mathbf{r}) = \frac{1}{4\pi\epsilon_0} \int_V \frac{\rho(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|} d^3\mathbf{r}'$ It is calculated by parameterizing boundary surfaces and applying two- or three-dimensional Gauss quadrature within the Boundary Element Method (BEM). The result maps spatial potential and evaluates electric field gradients $\mathbf{E} = -\nabla V$ to predict dielectric breakdown in air ($3 \text{ MV/m}$) and prevent corona discharge. The critical danger is evaluating potential directly on or near the charged surface; as observation coordinate $\mathbf{r}$ approaches source $\mathbf{r}'$, the kernel $|\mathbf{r} - \mathbf{r}'|^{-1}$ forms an integrable weak singularity that produces artificial infinite spikes on uniform grids. As a final sanity check, local coordinate transformations (Duffy's polar method) extract the kernel singularity analytically and Gauss's flux law is verified on an enclosing boundary.

2. Thermodynamics and astrophysics: total blackbody radiated power (Stefan-Boltzmann law)

Appears in stellar astrophysics when calculating total bolometric luminosity by integrating Planck's spectral radiance across all electromagnetic wavelengths from $\lambda = 0$ to $\lambda = \infty$: $B(T) = \int_0^\infty \frac{2 h c^2}{\lambda^5} \frac{1}{e^{\frac{h c}{\lambda k_B T}} - 1} d\lambda$ It is calculated by substituting dimensionless variable $x = \frac{h c}{\lambda k_B T}$, yielding: $B(T) = \frac{2 k_B^4 T^4}{h^3 c^2} \int_0^\infty \frac{x^3}{e^x - 1} dx$ The result evaluates the canonical integral $\int_0^\infty \frac{x^3}{e^x - 1} dx = \frac{\pi^4}{15}$, recovering the Stefan-Boltzmann law $u = \sigma T^4$ with $\sigma = \frac{2\pi^5 k_B^4}{15 c^2 h^3} \approx 5.67037 \times 10^{-8} \text{ W/(m}^2\cdot\text{K}^4)$. The critical danger is truncating the semi-infinite domain at low $x$ (e.g. $x=5$) when modeling infrared bands; truncating early discards over $15\%$ of thermal emission, biasing stellar effective temperatures. As a final sanity check, Gauss-Laguerre quadrature (tailored for semi-infinite integrals with weight $e^{-x}$) confirms convergence to the analytic $\pi^4/15 \approx 6.493939$.

3. Relativistic dynamics and particle physics: cross sections and phase-space integrals

Appears in high-energy physics at particle colliders (CERN) to predict particle production rates (such as Higgs boson cross sections or Møller scattering $e^- e^- \to e^- e^-$) by integrating Feynman invariant matrix elements $|\mathcal{M}|^2$ across relativistic Lorentz phase space: $\sigma = \frac{1}{2 s} \int |\mathcal{M}|^2 d\Phi_n$ It is calculated using adaptive multidimensional Monte Carlo algorithms like VEGAS with importance-sampled grid refinement. The result produces total cross section $\sigma$ in picobarns ($1 \text{ pb} = 10^{-40} \text{ m}^2$), establishing detection thresholds for beyond-Standard-Model physics. The critical danger is using fixed-step quadrature across narrow Breit-Wigner resonances ($|\mathcal{M}|^2 \propto \frac{1}{(s - M^2)^2 + M^2 \Gamma^2}$); if decay width $\Gamma$ is narrow (such as the $Z^0$ boson, $\Gamma_Z \approx 2.5 \text{ GeV}$ vs $M_Z \approx 91.2 \text{ GeV}$), fixed grids skip the resonance peak entirely, underestimating cross sections by orders of magnitude. As a final sanity check, multichannel variable flattening flattens resonance spikes prior to numerical quadrature.

4. Solid-state physics: density of states and Fermi level in conduction bands

Appears in condensed matter physics when evaluating conduction electron density $n$ and finite-temperature chemical potential $\mu(T)$ by integrating the product of density of states $g(E)$ and the Fermi-Dirac distribution: $n = \int_0^\infty g(E) f(E) dE = \int_0^\infty \frac{1}{2\pi^2} \left( \frac{2m^*}{\hbar^2} \right)^{3/2} \frac{\sqrt{E}}{e^{\frac{E - \mu}{k_B T}} + 1} dE$ It is calculated using Gauss-Fermi quadrature or Sommerfeld asymptotic low-temperature expansions ($k_B T \ll E_F$). The result determines chemical potential $\mu(T)$, Pauli spin susceptibility, and linear electronic specific heat $C_e = \gamma T$. The critical danger is evaluating integrals at cryogenic temperatures ($T \to 0$) without resolving the steep near-step cliff in the Fermi-Dirac function at $E = \mu$; derivatives reach $1/(4 k_B T)$, demanding sub-micro-eV adaptive grid refinement across the Fermi surface. As a final sanity check, the $T \to 0$ limit is checked against Sommerfeld's exact baseline $n = \frac{2}{3} g(E_F) E_F$.

5. Fluid mechanics and aerodynamics: airfoil lift via surface pressure integration

Appears in wind tunnel testing and aerospace engineering when evaluating aerodynamic lift force $L$ and pitching moments by integrating surface static pressure $P(s)$ and viscous skin friction $\tau(s)$ over closed wing profiles: $L = \oint_{\text{profile}} [ -P(s) \cos\theta(s) - \tau(s) \sin\theta(s) ] ds$ It is calculated by sampling static pressure ports distributed around suction and pressure surfaces, applying non-uniform composite trapezoidal integration along closed curvilinear contours. The result produces lift coefficient $C_L$, certifying aircraft takeoff, stall, and cruise margins. The critical danger is inadequate pressure tap density near the leading edge; at high angles of attack, flow accelerates sharply around the nose, forming a narrow suction spike. Missing the leading $5\%$ of chord underestimates total wing lift by over $30\%$, misidentifying aerodynamic stall. As a final sanity check, Kutta-Joukowski's theorem ($L' = \rho_\infty V_\infty \Gamma$) confirms that surface-integrated lift matches far-field circulation $\Gamma = \oint \mathbf{v} \cdot d\mathbf{r}$.

Guided Practice: Numerical Quadrature

Level 1: recognition

Select the best quadrature algorithm for each physical problem:

1. Evaluating magnetic flux $\Phi = \int_S \mathbf{B} \cdot d\mathbf{A}$ across 11 equispaced planar Hall sensor readings.
2. Evaluating the Gaussian error function $\operatorname{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2} dt$ to $10^{-14}$ precision in quantum simulations.
3. Integrating gravitational potential from non-uniform spatial core density samples across 8 coordinates.
4. Computing $\int_0^1 \frac{\cos(x)}{\sqrt{x}} dx$.

Solution:

1. Composite Simpson's 1/3 rule (11 points yield $N=10$ even intervals, optimal for $\mathcal{O}(h^4)$).
2. Gauss-Legendre quadrature with $n=8$ or $n=10$ nodes (exponential convergence on smooth entire functions like $e^{-t^2}$).
3. Non-uniform composite Trapezoidal rule (weights individual trapezoidal slices).
4. Regularizing variable substitution $x = u^2 \implies dx = 2u du$, transforming the integral into $2 \int_0^1 \cos(u^2) du$, which is completely smooth and pole-free, easily computed by Gauss-Legendre or Simpson.

Level 2: direct application

Apply composite Simpson's 1/3 rule with $N=2$ subintervals ($h = 0.5$) to approximate: $I = \int_0^1 \frac{1}{1 + x^2} dx$ Exact analytical solution is $\arctan(1) - \arctan(0) = \frac{\pi}{4} \approx 0.7853982$.

Solution:

1. Grid nodes: $a = 0, b = 1, N = 2 \implies h = 0.5$. $x_0 = 0.0, x_1 = 0.5, x_2 = 1.0$.
2. Function values $f(x) = \frac{1}{1 + x^2}$:
   • $f(0.0) = 1.0000000$.
   • $f(0.5) = \frac{1}{1 + 0.25} = 0.8000000$.
   • $f(1.0) = \frac{1}{1 + 1} = 0.5000000$.
3. Simpson 1/3 evaluation: $I_{\text{Simp}} = \frac{0.5}{3} [1.0000000 + 4(0.8000000) + 0.5000000] = \frac{0.5}{3} [4.7] \approx 0.7833333$
4. Contrast against exact value: $E_{\text{abs}} = |0.7853982 - 0.7833333| = 0.0020649 \quad (0.26\% \text{ relative error with only 3 nodes}).$

Level 3: tool selection

Compute crystal Debye heat capacity: $C_v(T) = 9 N k_B \left( \frac{T}{\Theta_D} \right)^3 \int_0^{\Theta_D/T} \frac{x^4 e^x}{(e^x - 1)^2} dx$ In the cryogenic limit ($T \to 0$), upper bound $\Theta_D / T \to \infty$. How should this be computed numerically without precision loss?

Solution: At cryogenic temperatures, the upper limit tends to infinity while the integrand decays as $x^4 e^{-x}$.

1. Analytically, the infinite integral converges to: $\int_0^\infty \frac{x^4 e^x}{(e^x - 1)^2} dx = \frac{4\pi^4}{15} \approx 25.97576$
2. Rather than forcing quadrature on an arbitrary large cutoff like 1000 (inviting floating-point overflow in $e^x$), replace the integral with its closed-form value $\frac{4\pi^4}{15}$, recovering Debye's $T^3$ law: $C_v(T) \approx \frac{12\pi^4}{5} N k_B \left( \frac{T}{\Theta_D} \right)^3$
3. For intermediate regimes ($0.1 < T/\Theta_D < 2$), apply 10-point Gauss-Legendre quadrature on $[0, \Theta_D/T]$.

Level 4: physical interpretation

Show why Simpson's 1/3 rule, derived by fitting 2nd-degree quadratics ($ax^2 + bx + c$), is exact for any cubic polynomial ($f(x) = a_3 x^3 + a_2 x^2 + a_1 x + a_0$). What physical symmetry explains this bonus accuracy order?

Solution: Consider integrating monomial $f(x) = x^3$ on symmetric interval $[-h, h]$:

1. Exact integral of odd functions over symmetric spans is zero: $\int_{-h}^h x^3 dx = 0$
2. Simpson's evaluation: $S(h) = \frac{h}{3} [(-h)^3 + 4(0)^3 + h^3] = \frac{h}{3} [-h^3 + 0 + h^3] = 0$ Simpson's output equals zero identically, matching the analytic integral. Physical interpretation: By parity symmetry, cubic truncation error is an odd function relative to the interval midpoint. Integrating across the symmetric bracket causes left-hand errors to cancel right-hand errors identically. The $h^3$ error term vanishes, promoting Simpson's rule to fourth-order $\mathcal{O}(h^4)$ accuracy.

Level 5: intentional mistake

A student calculates electrostatic potential at the origin from a uniform line of charge of length $L=1$: $V = \frac{1}{4\pi\epsilon_0} \int_0^1 \frac{\lambda}{x} dx$ The student runs composite trapezoidal integration with $x_0 = 0$. The code crashes with division by zero. Explain the physical and mathematical blunder and show how it is resolved.

Conceptual blunder: $f(x) = 1/x$ possesses a simple pole at $x=0$.

1. Mathematically, the improper integral diverges logarithmically: $\int_0^1 \frac{1}{x} dx = [\ln x]_0^1 = -(-\infty) = +\infty$
2. Quadrature cannot integrate non-integrable singular poles.
3. Physically, $x=0$ lies directly on an idealized zero-thickness line charge. Real wires carry finite radius $R > 0$. The physical axial potential is: $V = \frac{\lambda}{4\pi\epsilon_0} \int_0^1 \frac{dx}{\sqrt{x^2 + R^2}} = \frac{\lambda}{4\pi\epsilon_0} \ln\left( \frac{1 + \sqrt{1 + R^2}}{R} \right)$ Including physical conductor thickness $R$ renders the integrand finite and smooth across the entire interval, enabling clean numerical quadrature.