Numerical Integration

Why Numerical Integration Matters

A velocity sensor records how a vehicle's speed changes over time. Total distance is an integral, but the sensor provides samples rather than an antiderivative. A simulator has a similar problem when it can evaluate a response at a chosen input yet cannot integrate the response symbolically. Numerical integration turns those evaluations into an approximation of accumulated change.

This lesson computes a real, finite-interval integral $I=\int_a^b f(x)\,dx$. It covers the composite trapezoidal rule, composite Simpson's rule, adaptive Simpson integration, Gaussian quadrature, and Monte Carlo integration. You will trace their sampling and weighting decisions, implement all five in three languages, and distinguish an actual error, an estimated error, and statistical uncertainty.

The first four methods exploit regularity in a one-dimensional function. Monte Carlo recasts an integral as an expectation and is especially useful when extending to many dimensions. None of these implementations automatically handles infinite intervals or an undefined endpoint; such problems need a transformation or a specialized method.

Core Concepts and Terminology

For $a\lt b$, a quadrature rule has the form $Q=\sum_i w_i f(x_i)$. The nodes $x_i$ specify evaluations and the weights $w_i$ include interval scaling. This is a dot product of two vectors, not a linear system: no matrix elimination or pivoting is required. Reversing the endpoints changes the integral's sign, while an uncertainty magnitude stays nonnegative.

On $[0,1]$, Simpson's single-panel approximation to $f(x)=x^4$ uses nodes $0,1/2,1$ and weights $1/6,4/6,1/6$. Its result is $5/24\approx0.208333$, while the exact integral is $1/5=0.2$. The diagram shows both the sample vector and the weighted sum.

Integration is signed: regions below the axis contribute negatively. If the application asks for geometric area, integrate $\lvert f(x)\rvert$ or split at sign changes. Also distinguish errors in the measured or simulated function from errors introduced by quadrature; refining a grid cannot repair an inaccurate model.

How Numerical Integration Works, Step by Step

1. Establish the domain and numerical contract. Validate finite endpoints and a finite representable interval width. Require a positive subdivision count, an even $n\ge2$ for composite Simpson, a positive finite adaptive tolerance, or at least two Monte Carlo samples when estimating a variance. The programs normalize reversed endpoints internally and apply the orientation sign to the final estimate. A zero-length interval returns zero without evaluating the function after validating the method's parameters.

2. Choose how to sample and combine values. The five methods use different information:

• Composite trapezoidal: join adjacent samples by straight lines and integrate those segments. With $f_i=f(a+ih)$, use $T_n=h[f_0/2+\sum_{i=1}^{n-1}f_i+f_n/2]$. It works naturally with equally spaced measured data.
• Composite Simpson: fit a quadratic through each group of three nodes on an equally spaced grid. With even $n$, use $S_n=(h/3)[f_0+f_n+4\sum_{i\text{ odd}}f_i+2\sum_{i\text{ even},\,0\lt i\lt n}f_i]$. Although the interpolant is quadratic, symmetry makes the rule exact for cubics.
• Adaptive Simpson: compare a coarse Simpson panel with two half-panels, then refine only where their disagreement exceeds a local tolerance budget.
• Gaussian quadrature: choose nodes and weights jointly to maximize polynomial exactness. The implementations use a fixed three-node Gauss–Legendre rule, optionally repeated over equal panels.
• Monte Carlo: sample positions uniformly, average their function values, and multiply by interval length. The sample variance supplies an estimated standard error.

3. For adaptive Simpson, reuse three cached values and evaluate only two new ones. On a current interval $[\ell,r]$ with midpoint $m$, define

$$ S_1=\frac{r-\ell}{6}[f(\ell)+4f(m)+f(r)]. $$

Evaluate the quarter points $q_\ell=(\ell+m)/2$ and $q_r=(m+r)/2$. Compute one Simpson value on each half, and let their sum be $S_2$. Set $\Delta=S_2-S_1$, estimate $e=\lvert\Delta\rvert/15$, and use the corrected value $Q=S_2+\Delta/15$ when accepting that interval.

Why divide by $15$? In the smooth asymptotic regime, halving Simpson's step divides its leading error by $2^4=16$. Thus $I-S_2\approx(S_2-S_1)/15$. This estimates the error of the uncorrected refined value; adding the correction often improves the approximation further. It is an estimate rather than a proof for an arbitrary integrand.

4. Accept or split without losing the interval's budget. If $e\le\tau_{\mathrm{local}}$, add the corrected value to the total. Otherwise split into two child intervals, give each half the parent's tolerance budget, and preserve cached endpoint and midpoint values. The code uses an explicit depth-first stack of ordinary arrays. It pushes the right child first so the left child is processed first.

The initial interval has remaining depth $D$. A remaining depth of zero still permits the two new evaluations needed for a comparison, but prevents further splitting. If its estimate exceeds its budget, return a best-effort contribution and set status depth. If rounding prevents distinct quarter points, keep the available coarse contribution, mark its error unavailable with infinity, and return status stalled. Underflow of a required child tolerance budget also produces stalled. An ok status means every terminal interval passed the estimator test; it does not certify the actual error.

flowchart TD
    A["Validate domain, tolerance, and depth"] --> B{"Zero interval?"}
    B -->|Yes| Z["Return zero without evaluations"]
    B -->|No| C["Cache endpoints and midpoint; push interval"]
    C --> D["Pop one interval"]
    D --> E{"Distinct quarter points?"}
    E -->|No| F["Keep coarse value; flag stalled"]
    E -->|Yes| G["Evaluate two points; compute refined value and error"]
    G --> H{"Error within local budget?"}
    H -->|Yes| I["Accumulate corrected value and error estimate"]
    H -->|No| J{"Depth and child budget available?"}
    J -->|No| K["Accumulate best effort; flag depth or stalled"]
    J -->|Yes| L["Push children with half budgets and cached values"]
    F --> M{"Stack empty?"}
    I --> M
    K --> M
    L --> M
    M -->|No| D
    M -->|Yes| N["Apply orientation; return estimate, error, count, status"]

The diagram describes the adaptive method. If even the initial midpoint is not representable strictly inside the interval, the code uses the two endpoint values for a best-effort trapezoid and reports stalled. A nonfinite function value or checked arithmetic result raises an exception in every implementation.

5. Map Gaussian nodes correctly. On the reference interval $[-1,1]$, three-node Gauss–Legendre has nodes $t=-\sqrt{3/5},0,\sqrt{3/5}$ and weights $5/9,8/9,5/9$. For a physical panel $[\ell,r]$, let $c=(\ell+r)/2$ and $d=(r-\ell)/2$. Then

$$ \int_\ell^r f(x)\,dx\approx d\sum_{j=1}^{3} w_j f(c+d t_j). $$

The factor $d$ is the change-of-variable Jacobian. Forgetting it scales every answer incorrectly. The code evaluates the equivalent expression $\ell+d(1+t_j)$ to avoid summing two large endpoints. General $q$-node Gauss–Legendre uses roots of the degree-$q$ Legendre polynomial and is exact through degree $2q-1$; the supplied programs tabulate only $q=3$.

6. Track Monte Carlo uncertainty alongside its average. For independent uniform $U_i$ on $[0,1]$, set $X_i=a+(b-a)U_i$ and $Y_i=f(X_i)$. Then

$$ \widehat I_N=(b-a)\overline Y,\qquad s^2=\frac{1}{N-1}\sum_{i=1}^N(Y_i-\overline Y)^2,\qquad \widehat{\operatorname{SE}}=\lvert b-a\rvert\sqrt{\frac{s^2}{N}}. $$

The code uses Welford's online mean and centered-sum update, avoiding storage of all samples. Its small deterministic generator exists to make the three languages reproduce the same sample stream; the statistical theory assumes independent sampling with finite variance. An approximate normal interval such as $\widehat I_N\pm1.96\widehat{\operatorname{SE}}$ needs those assumptions and an adequate sample size, and is not a guaranteed enclosure.

Worked Example with a Full Trace

Compute $I=\int_0^1x^4\,dx=1/5$. For a small fixed grid with $n=2$, $h=1/2$ and samples $0,1/16,1$, the trapezoidal rule gives $T_2=9/32=0.28125$. Simpson gives $S_2=5/24\approx0.2083333333$. For the adaptive trace, call that coarse single-panel value $S_1$; here the subscripts distinguish coarse and refined estimates rather than the composite grid count.

Run adaptive Simpson with total tolerance $\tau=10^{-4}$ and maximum depth $D=10$. The table lists every distinct evaluated node, sorted by position. The actual evaluation order in the programs is $0,1,1/2,1/4,3/4,1/8,3/8,5/8,7/8$.

The complete decision trace has one rejected root interval and two accepted children:

The final sum is $1/160+31/160=1/5$. The reported error estimate is $1/30720\approx0.0000325520833$, and the evaluation count is $9$. Both leaves pass their $0.00005$ budgets, so the status is ok. For this quartic, the correction cancels the entire Simpson truncation error in exact arithmetic; the positive reported diagnostic is not the actual error of the corrected result.

The three-node Gaussian rule also integrates this quartic exactly in exact arithmetic, using only three evaluations. Its physical nodes are $[1-\sqrt{3/5}]/2$, $1/2$, and $[1+\sqrt{3/5}]/2$, with physical weights $5/18,4/9,5/18$.

For a tiny Monte Carlo arithmetic example, suppose the two sampled uniform values happen to be $1/4$ and $3/4$. Their function values are $1/256$ and $81/256$, so the integral estimate is $41/256=0.16015625$. The sample variance is $25/512$, and its estimated standard error is $5/32=0.15625$. This deliberately small sample illustrates the calculation; it is separate from the seeded generator used by the programs.

Implementation in Python, C++, and Java

Each panel is a complete program with the same arithmetic and sample configuration. trapezoidal and simpson take the number of uniform subintervals; gauss3 takes the number of Gaussian panels, with three evaluations per panel. Adaptive Simpson takes an absolute tolerance and a maximum depth in $[0,30]$. Monte Carlo takes a sample count and an integer seed in $[1,2147483646]$.

Adaptive output contains the estimate, the sum of terminal error diagnostics, the number of function evaluations, and status ok, depth, or stalled. A stalled interval whose quarter points cannot be formed contributes an infinite error diagnostic. Fixed rules return a scalar without pretending to supply a checked error bound. Monte Carlo returns the estimate and a nonnegative standard error. Java uses the names adaptiveSimpson and monteCarlo; their logic matches the corresponding snake-case functions.

The Gaussian rule and the adaptive stack use ordinary nested arrays. Quadrature has no matrix pivoting step. Every function evaluation is checked for a finite value; checked nonfinite arithmetic and invalid parameters raise an exception. The examples assume magnitudes where the remaining floating-point sums are meaningful. These implementations do not turn arbitrary singular functions into valid finite-interval inputs.

All five integrate $x^4$ on $[0,1]$. Use $1000$ trapezoidal subintervals, $10$ Simpson subintervals, adaptive tolerance $10^{-8}$ with depth $20$, one three-node Gaussian panel, and $100000$ Monte Carlo samples with seed $1$. The adaptive sample uses a tighter tolerance than the hand trace. Its reported diagnostic is about $7.9473\times10^{-9}$ with $65$ function evaluations.

import math
from dataclasses import dataclass


def finite(x):
    if not math.isfinite(x):
        raise ValueError("nonfinite function value or arithmetic result")
    return x


def interval(a, b):
    if not math.isfinite(a) or not math.isfinite(b) or not math.isfinite(b-a):
        raise ValueError("finite endpoints and finite interval width required")
    return [a, b, 1.0] if a <= b else [b, a, -1.0]


def count(n, minimum):
    if not isinstance(n, int) or n < minimum:
        raise ValueError("invalid subdivision or sample count")


def trapezoidal(f, a, b, n):
    count(n, 1)
    a, b, sign = interval(a, b)
    if a == b:
        return 0.0
    h = (b-a)/n
    total = 0.5*finite(f(a)) + 0.5*finite(f(b))
    for i in range(1, n):
        total += finite(f(a+i*h))
    return finite(sign*h*total)


def simpson(f, a, b, n):
    count(n, 2)
    if n % 2:
        raise ValueError("Simpson requires an even subdivision count")
    a, b, sign = interval(a, b)
    if a == b:
        return 0.0
    h = (b-a)/n
    total = finite(f(a)) + finite(f(b))
    for i in range(1, n):
        total += (4.0 if i % 2 else 2.0)*finite(f(a+i*h))
    return finite(sign*(h/3.0)*total)


def gauss3(f, a, b, panels=1):
    count(panels, 1)
    a, b, sign = interval(a, b)
    if a == b:
        return 0.0
    r = math.sqrt(3.0/5.0)
    # Rows are [reference node, weight]; no matrix pivoting is involved.
    rule = [[-r, 5.0/9.0], [0.0, 8.0/9.0], [r, 5.0/9.0]]
    total = 0.0
    for i in range(panels):
        left = a+(b-a)*(i/panels)
        right = b if i == panels-1 else a+(b-a)*((i+1)/panels)
        half = (right-left)/2.0
        subtotal = 0.0
        for node, weight in rule:
            subtotal += weight*finite(f(left+half*(1.0+node)))
        total += half*subtotal
    return finite(sign*total)


def panel(a, b, fa, fm, fb):
    return finite((b-a)*(fa/6.0+(2.0/3.0)*fm+fb/6.0))


@dataclass
class Adaptive:
    value: float
    error: float
    evaluations: int
    status: str


def adaptive_simpson(f, a, b, tol=1e-8, max_depth=20):
    if (not math.isfinite(tol) or tol <= 0.0 or not isinstance(max_depth, int)
            or not 0 <= max_depth <= 30):
        raise ValueError("positive finite tolerance and depth in [0, 30] required")
    a, b, sign = interval(a, b)
    if a == b:
        return Adaptive(0.0, 0.0, 0, "ok")
    fa, fb = finite(f(a)), finite(f(b))
    mid = a+(b-a)/2.0
    if not a < mid < b:
        return Adaptive(finite(sign*(b-a)*(fa/2.0+fb/2.0)), math.inf, 2, "stalled")
    fm = finite(f(mid))
    evaluations = 3
    # Stack rows: [a, b, fa, fm, fb, coarse Simpson, tolerance, depth].
    stack = [[a, b, fa, fm, fb, panel(a,b,fa,fm,fb), tol, max_depth]]
    total, error = 0.0, 0.0
    depth_hit, stalled = False, False
    while stack:
        lo, hi, fl, fm, fr, coarse, budget, depth = stack.pop()
        mid = lo+(hi-lo)/2.0
        ql, qr = lo+(mid-lo)/2.0, mid+(hi-mid)/2.0
        if not lo < ql < mid < qr < hi:
            total += coarse
            error = math.inf
            stalled = True
            continue
        fql, fqr = finite(f(ql)), finite(f(qr))
        evaluations += 2
        left = panel(lo,mid,fl,fql,fm)
        right = panel(mid,hi,fm,fqr,fr)
        delta = finite(finite(left+right)-coarse)
        estimate = abs(delta)/15.0
        corrected = finite(finite(left+right)+delta/15.0)
        if estimate <= budget or depth == 0 or budget/2.0 == 0.0:
            total += corrected
            error += estimate
            if estimate > budget:
                depth_hit = depth_hit or depth == 0
                stalled = stalled or budget/2.0 == 0.0
        else:
            # Push right first so left is processed first.
            stack.append([mid,hi,fm,fqr,fr,right,budget/2.0,depth-1])
            stack.append([lo,mid,fl,fql,fm,left,budget/2.0,depth-1])
    status = "stalled" if stalled else "depth" if depth_hit else "ok"
    return Adaptive(finite(sign*total), error, evaluations, status)


def monte_carlo(f, a, b, n, seed=1):
    count(n, 2)
    if not isinstance(seed, int) or not 1 <= seed < 2147483647:
        raise ValueError("seed must be an integer in [1, 2147483646]")
    a, b, sign = interval(a, b)
    if a == b:
        return [0.0, 0.0]
    state, mean, m2 = seed, 0.0, 0.0
    for i in range(1, n+1):
        state = (16807*state) % 2147483647
        u = state/2147483647.0
        y = finite(f(a+(b-a)*u))
        delta = y-mean
        mean += delta/i
        m2 += delta*(y-mean)
        finite(mean)
        finite(m2)
    variance = max(0.0, m2/(n-1))
    return [finite(sign*(b-a)*mean), finite((b-a)*math.sqrt(variance/n))]


def main():
    def f(x):
        return (x*x)*(x*x)

    exact = 0.2
    t = trapezoidal(f,0.0,1.0,1000)
    s = simpson(f,0.0,1.0,10)
    a = adaptive_simpson(f,0.0,1.0,1e-8,20)
    g = gauss3(f,0.0,1.0)
    mc, se = monte_carlo(f,0.0,1.0,100000,1)
    checks = [abs(t-exact)<=4e-7, abs(s-exact)<=1.4e-5,
              a.status=="ok" and abs(a.value-exact)<=1e-10,
              abs(g-exact)<=1e-12, abs(mc-exact)<=5.0*se]
    for name, result, good in zip(["trapezoidal","simpson","adaptive","gauss3","monte-carlo"],
                                  [t,s,a.value,g,mc],checks):
        print(f"{name}: integral = {result:.6f}, check = {'PASS' if good else 'FAIL'}")
    print(f"adaptive status = {a.status}")
    print(f"monte-carlo standard error = {se:.6f}")
    print(f"all methods = {'PASS' if all(checks) else 'FAIL'}")


if __name__ == "__main__":
    main()

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <functional>
#include <iomanip>
#include <iostream>
#include <limits>
#include <stdexcept>
#include <string>
#include <vector>
using Fn=std::function<double(double)>;

double finite_value(double x) {
    if (!std::isfinite(x)) throw std::invalid_argument("nonfinite function value or arithmetic result");
    return x;
}
std::vector<double> interval(double a,double b) {
    if (!std::isfinite(a) || !std::isfinite(b) || !std::isfinite(b-a))
        throw std::invalid_argument("finite endpoints and finite interval width required");
    return a<=b ? std::vector<double>{a,b,1.0} : std::vector<double>{b,a,-1.0};
}
void count(int n,int minimum) {
    if (n<minimum) throw std::invalid_argument("invalid subdivision or sample count");
}
double trapezoidal(const Fn& f,double a,double b,int n) {
    count(n,1); auto v=interval(a,b); a=v[0]; b=v[1]; double sign=v[2];
    if (a==b) return 0.0;
    double h=(b-a)/n;
    double total=0.5*finite_value(f(a))+0.5*finite_value(f(b));
    for (int i=1;i<n;++i) total+=finite_value(f(a+i*h));
    return finite_value(sign*h*total);
}
double simpson(const Fn& f,double a,double b,int n) {
    count(n,2);
    if (n%2) throw std::invalid_argument("Simpson requires an even subdivision count");
    auto v=interval(a,b); a=v[0]; b=v[1]; double sign=v[2];
    if (a==b) return 0.0;
    double h=(b-a)/n, total=finite_value(f(a))+finite_value(f(b));
    for (int i=1;i<n;++i) total+=(i%2 ? 4.0 : 2.0)*finite_value(f(a+i*h));
    return finite_value(sign*(h/3.0)*total);
}
double gauss3(const Fn& f,double a,double b,int panels=1) {
    count(panels,1); auto v=interval(a,b); a=v[0]; b=v[1]; double sign=v[2];
    if (a==b) return 0.0;
    double r=std::sqrt(3.0/5.0);
    // Rows are [reference node, weight]; no matrix pivoting is involved.
    std::vector<std::vector<double>> rule={ {-r,5.0/9.0}, {0.0,8.0/9.0}, {r,5.0/9.0} };
    double total=0.0;
    for (int i=0;i<panels;++i) {
        double left=a+(b-a)*(static_cast<double>(i)/panels);
        double right=i==panels-1 ? b : a+(b-a)*(static_cast<double>(i+1)/panels);
        double half=(right-left)/2.0, subtotal=0.0;
        for (const auto& row:rule) subtotal+=row[1]*finite_value(f(left+half*(1.0+row[0])));
        total+=half*subtotal;
    }
    return finite_value(sign*total);
}
double panel(double a,double b,double fa,double fm,double fb) {
    return finite_value((b-a)*(fa/6.0+(2.0/3.0)*fm+fb/6.0));
}
struct Adaptive { double value,error; std::int64_t evaluations; std::string status; };
Adaptive adaptive_simpson(const Fn& f,double a,double b,double tol=1e-8,int max_depth=20) {
    if (!std::isfinite(tol) || tol<=0 || max_depth<0 || max_depth>30)
        throw std::invalid_argument("positive finite tolerance and depth in [0, 30] required");
    auto v=interval(a,b); a=v[0]; b=v[1]; double sign=v[2];
    if (a==b) return {0.0,0.0,0,"ok"};
    double fa=finite_value(f(a)), fb=finite_value(f(b)), mid=a+(b-a)/2.0;
    if (!(a<mid && mid<b))
        return {finite_value(sign*(b-a)*(fa/2.0+fb/2.0)),std::numeric_limits<double>::infinity(),2,"stalled"};
    double fm=finite_value(f(mid));
    std::int64_t evaluations=3;
    // Stack rows: [a, b, fa, fm, fb, coarse Simpson, tolerance, depth].
    std::vector<std::vector<double>> stack={ {a,b,fa,fm,fb,panel(a,b,fa,fm,fb),tol,static_cast<double>(max_depth)} };
    double total=0.0,error=0.0;
    bool depth_hit=false,stalled=false;
    while (!stack.empty()) {
        auto row=stack.back();stack.pop_back();
        double lo=row[0],hi=row[1],fl=row[2],fmid=row[3],fr=row[4],coarse=row[5],budget=row[6];
        int depth=static_cast<int>(row[7]);
        mid=lo+(hi-lo)/2.0;
        double ql=lo+(mid-lo)/2.0,qr=mid+(hi-mid)/2.0;
        if (!(lo<ql && ql<mid && mid<qr && qr<hi)) {
            total+=coarse;error=std::numeric_limits<double>::infinity();stalled=true;continue;
        }
        double fql=finite_value(f(ql)),fqr=finite_value(f(qr));evaluations+=2;
        double left=panel(lo,mid,fl,fql,fmid),right=panel(mid,hi,fmid,fqr,fr);
        double delta=finite_value(finite_value(left+right)-coarse);
        double estimate=std::abs(delta)/15.0;
        double corrected=finite_value(finite_value(left+right)+delta/15.0);
        if (estimate<=budget || depth==0 || budget/2.0==0.0) {
            total+=corrected;error+=estimate;
            if (estimate>budget) {
                depth_hit=depth_hit || depth==0;
                stalled=stalled || budget/2.0==0.0;
            }
        } else {
            // Push right first so left is processed first.
            stack.push_back({mid,hi,fmid,fqr,fr,right,budget/2.0,static_cast<double>(depth-1)});
            stack.push_back({lo,mid,fl,fql,fmid,left,budget/2.0,static_cast<double>(depth-1)});
        }
    }
    std::string status=stalled ? "stalled" : depth_hit ? "depth" : "ok";
    return {finite_value(sign*total),error,evaluations,status};
}
std::vector<double> monte_carlo(const Fn& f,double a,double b,int n,std::int64_t seed=1) {
    count(n,2);
    if (seed<1 || seed>=2147483647) throw std::invalid_argument("seed must be in [1, 2147483646]");
    auto v=interval(a,b);a=v[0];b=v[1];double sign=v[2];
    if (a==b) return {0.0,0.0};
    std::int64_t state=seed;double mean=0.0,m2=0.0;
    for (int i=0;i<n;++i) {
        state=(16807*state)%2147483647;
        double u=state/2147483647.0,y=finite_value(f(a+(b-a)*u));
        double delta=y-mean;mean+=delta/(i+1);m2+=delta*(y-mean);
        finite_value(mean);finite_value(m2);
    }
    double variance=std::max(0.0,m2/(n-1));
    return {finite_value(sign*(b-a)*mean),finite_value((b-a)*std::sqrt(variance/n))};
}
int main() {
    Fn f=[](double x){return (x*x)*(x*x);};
    double exact=0.2,t=trapezoidal(f,0,1,1000),s=simpson(f,0,1,10);
    Adaptive a=adaptive_simpson(f,0,1,1e-8,20);
    double g=gauss3(f,0,1);auto mc=monte_carlo(f,0,1,100000,1);
    std::vector<bool> checks={std::abs(t-exact)<=4e-7,std::abs(s-exact)<=1.4e-5,
        a.status=="ok" && std::abs(a.value-exact)<=1e-10,std::abs(g-exact)<=1e-12,
        std::abs(mc[0]-exact)<=5.0*mc[1]};
    std::vector<std::string> names={"trapezoidal","simpson","adaptive","gauss3","monte-carlo"};
    std::vector<double> results={t,s,a.value,g,mc[0]};bool ok=true;
    std::cout<<std::fixed<<std::setprecision(6);
    for (int i=0;i<5;++i) {
        ok=ok && checks[i];
        std::cout<<names[i]<<": integral = "<<results[i]<<", check = "<<(checks[i]?"PASS":"FAIL")<<'\n';
    }
    std::cout<<"adaptive status = "<<a.status<<'\n';
    std::cout<<"monte-carlo standard error = "<<mc[1]<<'\n';
    std::cout<<"all methods = "<<(ok?"PASS":"FAIL")<<'\n';
    return 0;
}

import java.util.ArrayList;
import java.util.Locale;
import java.util.function.DoubleUnaryOperator;

public class Main {
    static double finite(double x) {
        if (!Double.isFinite(x)) throw new IllegalArgumentException("nonfinite function value or arithmetic result");
        return x;
    }
    static double[] interval(double a,double b) {
        if (!Double.isFinite(a) || !Double.isFinite(b) || !Double.isFinite(b-a))
            throw new IllegalArgumentException("finite endpoints and finite interval width required");
        return a<=b ? new double[]{a,b,1.0} : new double[]{b,a,-1.0};
    }
    static void count(int n,int minimum) {
        if (n<minimum) throw new IllegalArgumentException("invalid subdivision or sample count");
    }
    static double trapezoidal(DoubleUnaryOperator f,double a,double b,int n) {
        count(n,1);double[] v=interval(a,b);a=v[0];b=v[1];double sign=v[2];
        if (a==b) return 0.0;
        double h=(b-a)/n,total=0.5*finite(f.applyAsDouble(a))+0.5*finite(f.applyAsDouble(b));
        for (int i=1;i<n;++i) total+=finite(f.applyAsDouble(a+i*h));
        return finite(sign*h*total);
    }
    static double simpson(DoubleUnaryOperator f,double a,double b,int n) {
        count(n,2);
        if (n%2!=0) throw new IllegalArgumentException("Simpson requires an even subdivision count");
        double[] v=interval(a,b);a=v[0];b=v[1];double sign=v[2];
        if (a==b) return 0.0;
        double h=(b-a)/n,total=finite(f.applyAsDouble(a))+finite(f.applyAsDouble(b));
        for (int i=1;i<n;++i) total+=(i%2!=0 ? 4.0 : 2.0)*finite(f.applyAsDouble(a+i*h));
        return finite(sign*(h/3.0)*total);
    }
    static double gauss3(DoubleUnaryOperator f,double a,double b,int panels) {
        count(panels,1);double[] v=interval(a,b);a=v[0];b=v[1];double sign=v[2];
        if (a==b) return 0.0;
        double r=Math.sqrt(3.0/5.0);
        // Rows are [reference node, weight]; no matrix pivoting is involved.
        double[][] rule={ {-r,5.0/9.0}, {0.0,8.0/9.0}, {r,5.0/9.0} };
        double total=0.0;
        for (int i=0;i<panels;++i) {
            double left=a+(b-a)*((double)i/panels);
            double right=i==panels-1 ? b : a+(b-a)*((double)(i+1)/panels);
            double half=(right-left)/2.0,subtotal=0.0;
            for (double[] row:rule) subtotal+=row[1]*finite(f.applyAsDouble(left+half*(1.0+row[0])));
            total+=half*subtotal;
        }
        return finite(sign*total);
    }
    static double panel(double a,double b,double fa,double fm,double fb) {
        return finite((b-a)*(fa/6.0+(2.0/3.0)*fm+fb/6.0));
    }
    static class Adaptive {
        double value,error;long evaluations;String status;
        Adaptive(double value,double error,long evaluations,String status) {
            this.value=value;this.error=error;this.evaluations=evaluations;this.status=status;
        }
    }
    static Adaptive adaptiveSimpson(DoubleUnaryOperator f,double a,double b,double tol,int maxDepth) {
        if (!Double.isFinite(tol) || tol<=0 || maxDepth<0 || maxDepth>30)
            throw new IllegalArgumentException("positive finite tolerance and depth in [0, 30] required");
        double[] v=interval(a,b);a=v[0];b=v[1];double sign=v[2];
        if (a==b) return new Adaptive(0.0,0.0,0,"ok");
        double fa=finite(f.applyAsDouble(a)),fb=finite(f.applyAsDouble(b)),mid=a+(b-a)/2.0;
        if (!(a<mid && mid<b))
            return new Adaptive(finite(sign*(b-a)*(fa/2.0+fb/2.0)),Double.POSITIVE_INFINITY,2,"stalled");
        double fm=finite(f.applyAsDouble(mid));long evaluations=3;
        // Stack rows: [a, b, fa, fm, fb, coarse Simpson, tolerance, depth].
        ArrayList<double[]> stack=new ArrayList<>();
        stack.add(new double[]{a,b,fa,fm,fb,panel(a,b,fa,fm,fb),tol,maxDepth});
        double total=0.0,error=0.0;boolean depthHit=false,stalled=false;
        while (!stack.isEmpty()) {
            double[] row=stack.remove(stack.size()-1);
            double lo=row[0],hi=row[1],fl=row[2],fmid=row[3],fr=row[4],coarse=row[5],budget=row[6];
            int depth=(int)row[7];mid=lo+(hi-lo)/2.0;
            double ql=lo+(mid-lo)/2.0,qr=mid+(hi-mid)/2.0;
            if (!(lo<ql && ql<mid && mid<qr && qr<hi)) {
                total+=coarse;error=Double.POSITIVE_INFINITY;stalled=true;continue;
            }
            double fql=finite(f.applyAsDouble(ql)),fqr=finite(f.applyAsDouble(qr));evaluations+=2;
            double left=panel(lo,mid,fl,fql,fmid),right=panel(mid,hi,fmid,fqr,fr);
            double delta=finite(finite(left+right)-coarse),estimate=Math.abs(delta)/15.0;
            double corrected=finite(finite(left+right)+delta/15.0);
            if (estimate<=budget || depth==0 || budget/2.0==0.0) {
                total+=corrected;error+=estimate;
                if (estimate>budget) {
                    depthHit=depthHit || depth==0;
                    stalled=stalled || budget/2.0==0.0;
                }
            } else {
                // Push right first so left is processed first.
                stack.add(new double[]{mid,hi,fmid,fqr,fr,right,budget/2.0,depth-1});
                stack.add(new double[]{lo,mid,fl,fql,fmid,left,budget/2.0,depth-1});
            }
        }
        String status=stalled ? "stalled" : depthHit ? "depth" : "ok";
        return new Adaptive(finite(sign*total),error,evaluations,status);
    }
    static double[] monteCarlo(DoubleUnaryOperator f,double a,double b,int n,long seed) {
        count(n,2);
        if (seed<1 || seed>=2147483647L) throw new IllegalArgumentException("seed must be in [1, 2147483646]");
        double[] v=interval(a,b);a=v[0];b=v[1];double sign=v[2];
        if (a==b) return new double[]{0.0,0.0};
        long state=seed;double mean=0.0,m2=0.0;
        for (int i=0;i<n;++i) {
            state=(16807L*state)%2147483647L;
            double u=state/2147483647.0,y=finite(f.applyAsDouble(a+(b-a)*u));
            double delta=y-mean;mean+=delta/(i+1);m2+=delta*(y-mean);
            finite(mean);finite(m2);
        }
        double variance=Math.max(0.0,m2/(n-1));
        return new double[]{finite(sign*(b-a)*mean),finite((b-a)*Math.sqrt(variance/n))};
    }
    public static void main(String[] args) {
        Locale.setDefault(Locale.ROOT);
        DoubleUnaryOperator f=x -> (x*x)*(x*x);
        double exact=0.2,t=trapezoidal(f,0,1,1000),s=simpson(f,0,1,10);
        Adaptive a=adaptiveSimpson(f,0,1,1e-8,20);
        double g=gauss3(f,0,1,1);double[] mc=monteCarlo(f,0,1,100000,1);
        boolean[] checks={Math.abs(t-exact)<=4e-7,Math.abs(s-exact)<=1.4e-5,
            a.status.equals("ok") && Math.abs(a.value-exact)<=1e-10,Math.abs(g-exact)<=1e-12,
            Math.abs(mc[0]-exact)<=5.0*mc[1]};
        String[] names={"trapezoidal","simpson","adaptive","gauss3","monte-carlo"};
        double[] results={t,s,a.value,g,mc[0]};boolean ok=true;
        for (int i=0;i<5;++i) {
            ok=ok && checks[i];
            System.out.printf("%s: integral = %.6f, check = %s%n",names[i],results[i],checks[i]?"PASS":"FAIL");
        }
        System.out.println("adaptive status = "+a.status);
        System.out.printf("monte-carlo standard error = %.6f%n",mc[1]);
        System.out.println("all methods = "+(ok?"PASS":"FAIL"));
    }
}

All three programs print:

• trapezoidal: integral = 0.200000, check = PASS
• simpson: integral = 0.200013, check = PASS
• adaptive: integral = 0.200000, check = PASS
• gauss3: integral = 0.200000, check = PASS
• monte-carlo: integral = 0.199877, check = PASS
• adaptive status = ok
• monte-carlo standard error = 0.000842
• all methods = PASS

Each deterministic check compares against the independently known $I=0.2$ with an explicit threshold in the driver. In particular, the trapezoidal actual error is about $3.333333\times10^{-7}$, so rounding its answer to six decimal places hides its remaining error. The Monte Carlo check compares this particular seeded run's actual error with five estimated standard errors; it is a reproducible demonstration check, not an accuracy guarantee for future samples.

The reproducible sample stream uses $s_{k+1}=16807s_k\bmod2147483647$ and $u_k=s_k/2147483647$. The first four states from seed $1$ are $16807$, $282475249$, $1622650073$, and $984943658$. C++ and Java use 64-bit integer multiplication so the generator update does not overflow. This compact illustrative generator is kept explicit so a language's default random-number generator cannot silently change the sample input.

Complexity, Correctness, and Invariants

Let $C_f$ be the cost of one function evaluation, $L=\lvert b-a\rvert$, $h=L/n$, $P$ the number of Gaussian panels, $N$ the Monte Carlo sample count, and $D$ the maximum adaptive depth. The following errors concern exact-arithmetic quadrature under the stated smoothness assumptions; roundoff and input-model error are additional effects.

For fixed rules and a fixed sample count, best, average, and worst work are the same because the loop counts do not depend on function values. All methods return in $O(1)$ time on a validated zero-length interval. Adaptive work depends on where the integrand needs resolution; there is no distribution-free average tree size. The depth bound guarantees termination, including a reported failure to meet the diagnostic target. Unrepresentable quarter points can terminate some branches earlier than the full-tree count.

The trapezoidal and Simpson bounds imply $n=O(\varepsilon^{-1/2})$ and $n=O(\varepsilon^{-1/4})$, respectively, for fixed interval and derivative bounds when targeting absolute error $\varepsilon$. Composite three-node Gauss–Legendre has $P=O(\varepsilon^{-1/6})$ under its sixth-derivative bound. Plain Monte Carlo requires $N=O(\varepsilon^{-2})$ for a fixed root-mean-square target when its variance is finite. These rates compare different assumptions and constants; Monte Carlo's rate alone does not depend on dimension, but the variance and evaluation cost can depend strongly on dimension.

Why fixed rules are correct under their assumptions. Trapezoidal and Simpson rules integrate their local interpolation formulas, then sum over a partition. Their derivative remainder terms supply the stated error bounds. Gaussian nodes are roots of an orthogonal polynomial. Dividing any polynomial of degree at most $2q-1$ into a multiple of that degree-$q$ polynomial plus a lower-degree remainder explains the $2q-1$ exactness: the multiple integrates to zero by orthogonality and vanishes at the nodes, while the interpolatory rule integrates the remainder exactly.

Why adaptive bookkeeping is correct. Accepted leaves together with pending stack intervals form a disjoint partition of the original interval, apart from shared endpoints. Splitting replaces one interval by two children without a gap or overlap; halving the parent's tolerance preserves the total assigned budget. The accumulated integral is the sum of accepted leaf approximations, never the sum of both an ancestor and its descendants. If every leaf estimate passes, the sum of estimates is at most the original tolerance in exact arithmetic. That proves the budget bookkeeping, not that the estimator cannot miss a feature.

Why Monte Carlo is unbiased in its ideal sampling model. A uniform point $X$ on $[a,b]$ satisfies $\mathbb E[f(X)]=I/(b-a)$ for ordered endpoints. Linearity of expectation gives $\mathbb E[\widehat I_N]=I$. Independence gives the variance reduction by $N$; finite variance is needed for the stated root-mean-square result. A finite pseudorandom stream approximates this sampling model and does not make those probabilistic statements deterministic guarantees.

Variants, Pitfalls, and Optimizations

A small refinement difference can miss a narrow feature. Consider the continuous triangular bump $f(x)=\max(0,1-\lvert x-0.30\rvert/0.01)$ on $[0,1]$. It is nonzero only on $(0.29,0.31)$ and has exact integral $0.01$. Adaptive Simpson's first five nodes are $0,1/4,1/2,3/4,1$, all with value zero. Both estimates are zero, so the code returns ok with zero estimated error and an incorrect integral. Supplying breakpoints $0.29,0.30,0.31$ and integrating the linear pieces exposes the feature.

This is a limitation of finite sampling, not just of the triangular function: an unobserved smooth bump can produce the same issue. If known physics or data indicates a discontinuity, a corner, a narrow peak, or a rapid oscillation, split or transform the problem before relying on a local estimator. Refinement tests compare what has been sampled, not everything between samples.

Respect the input data. Standard composite Simpson requires equal spacing and an even subdivision count. It cannot be applied to arbitrary unevenly spaced measurements with the same $1,4,2,4,\ldots,1$ weights. A nonuniform trapezoid, $\sum_i (x_{i+1}-x_i)(y_i+y_{i+1})/2$, works directly with ordered irregular samples. Interpolation of noisy data adds a modeling decision; more quadrature evaluations of that interpolant do not create new measurements.

Use specialized variants when their assumptions fit. Adaptive Gauss–Kronrod pairs reuse samples to compare different-degree rules and are common in general-purpose integrators. Romberg extrapolates a sequence of trapezoidal grids on smooth problems. Periodic smooth integrands can make the trapezoidal rule exceptionally efficient. Endpoint singularities may benefit from a substitution, Gauss–Jacobi weights, or tanh–sinh quadrature. Highly oscillatory integrals often need oscillation-aware techniques rather than a generic tolerance decrease.

Improve Monte Carlo by reducing variance. Stratification, antithetic samples, control variates, and importance sampling can lower variance for suitable integrands. With an importance density $p(x)$, sample from $p$ and average $f(X)/p(X)$ over the domain where the integral receives contributions; the density must be positive there and the reweighted variance must be controlled. Quasi-Monte Carlo uses low-discrepancy points, and its error analysis differs from the independent-sample standard-error formula. Correlated samples also need a different uncertainty calculation.

Handle floating-point limits honestly. Add relative tolerance only with an appropriate absolute tolerance, especially near cancellation where the integral is close to zero. Compensated or pairwise summation can improve long sums. Cache shared nodes in adaptive rules; when doubling a trapezoidal grid, reuse all old samples. Never silently accept a maximum-depth exit as meeting the target. Treat an observed zero Monte Carlo sample variance cautiously: it can mean the samples missed a small region rather than that the integrand is constant.

Practice Problems and Real-World Use

Start with exact polynomial integrals so mistakes in weights, orientation, and scaling are visible. Then move to smooth transcendental functions, deliberately hidden features, and stochastic examples. Record evaluation counts as well as runtime: the cost of calling a simulation can dominate the arithmetic inside the integrator.

For further reading, NIST's quadrature chapter gives the classical rules, Gaussian quadrature, and remainder formulas. SciPy's Simpson documentation distinguishes sampled-data integration and its spacing-dependent exactness. Art Owen's Simple Monte Carlo chapter develops the sampling model and error estimation. Their assumptions should be checked against the data or callable function before interpreting a reported accuracy figure.