Numerical Integration

Approximate ${\int }_a^{b}f(x)dx$ when no closed form is available.

Simpson’s rule

Fit a parabola through the endpoints and the midpoint:

${\int }_a^{b}f\approx \frac{b-a}{6}(f(a)+4f\left(\frac{a+b}{2}\right)+f(b)).$

Composite Simpson with $n$ (even) subintervals:

double simpson(function<double(double)> f, double a, double b, int n) {
    double h = (b - a) / n, s = f(a) + f(b);
    for (int i = 1; i < n; i++)
        s += f(a + i * h) * (i % 2 ? 4 : 2);
    return s * h / 3;
}

Error $O(h^4)$ — exact for cubics, and excellent for smooth functions.

Adaptive Simpson — the one to use

Recursively subdivide only where the estimate is inaccurate:

double simp(function<double(double)> f, double a, double b) {
    double c = (a + b) / 2;
    return (b - a) / 6 * (f(a) + 4 * f(c) + f(b));
}
 
double adaptive(function<double(double)> f, double a, double b, double eps, double whole, int depth) {
    double c = (a + b) / 2;
    double left = simp(f, a, c), right = simp(f, c, b);
    if (depth <= 0 || fabs(left + right - whole) <= 15 * eps)
        return left + right + (left + right - whole) / 15;      // Richardson correction
    return adaptive(f, a, c, eps/2, left, depth-1)
         + adaptive(f, c, b, eps/2, right, depth-1);
}
 
double integrate(function<double(double)> f, double a, double b, double eps = 1e-9) {
    return adaptive(f, a, b, eps, simp(f, a, b), 50);
}

It spends effort where the function is difficult and almost none where it is flat. The /15 term is Richardson extrapolation, which upgrades the accuracy for free.

The methods

Method	Error	Exact for	Note
Rectangle (midpoint)	$O(h^2)$	linear	simplest
Trapezoid	$O(h^2)$	linear
Simpson	$O(h^4)$	cubic	the default
Simpson 3/8	$O(h^4)$	cubic
Romberg	$O(h^{2k})$	high degree	trapezoid + repeated Richardson
Gauss-Legendre ($n$ points)	exact to degree $2n-1$		best per evaluation, needs tabulated nodes
Monte Carlo	$O(1/\sqrt{N})$	—	the only option in high dimensions

When integration appears in contests

Problem	Integrand
Area under a curve	$f(x)$
Area of a region bounded by curves	the difference
Arc length	$\sqrt{1+f^{\prime }(x{)}^2}$
Volume of revolution	$\pi f(x{)}^2$
Volume by cross-sections	area of the slice at $x$
Expected value of a continuous distribution	$xp(x)$
Probability of a region	the density
Union of circles / shapes	slice area as a function of $x$

The volume by cross-sections pattern is the most common: compute the area of the 2D slice at each $x$ (often itself a geometric computation) and integrate over $x$ with adaptive Simpson.

Pitfalls

Simpson assumes smoothness

Adaptive Simpson silently gives wrong answers on functions with discontinuities or kinks, because the error estimate compares two smooth approximations that happen to agree.

Fix: split the interval at every known breakpoint (where shapes start or stop overlapping, where a maximum changes which curve is on top) and integrate each smooth piece separately. This is usually the difference between a correct solution and a mysteriously wrong one.

Other traps: singularities at the endpoints (substitute to remove them), and oscillatory integrands (Simpson can alias badly — increase the base subdivision).

Monte Carlo

For high dimensions, deterministic quadrature is hopeless ($n^d$ points). Monte Carlo converges at $O(1/\sqrt{N})$ independently of dimension:

${\int }_{\Omega }f\approx \frac{\mathrm{vol}(\Omega )}{N}{\sum }_{i=1}^{N}f(x_i).$

Slow (4× the samples for one more digit) but dimension-free. Variance reduction (importance sampling, stratification, quasi-random sequences) helps considerably.

Exact alternatives — check first

Many contest “integration” problems have exact answers:

• polygon area → shoelace, exact;
• circle-polygon intersection → closed-form sectors and triangles;
• piecewise-linear functions → sum of trapezoids, exact;
• polynomial integrands → integrate symbolically.

Reach for numerical integration only when the integrand genuinely has no closed form.

See also: Floating Point · Polygon Area · Probability