Trapezoidal rule

The formula

One panel between two samples:

∫ₐᵇ f(x) dx ≈ (b−a) · [ f(a) + f(b) ] / 2

Reading: the area of a trapezoid — width times the average of the two end heights. That is the whole idea; everything else is bookkeeping.

The composite rule for n intervals of uniform width h = (b−a)/n:

∫ₐᵇ f(x) dx ≈ h · [ f₀/2 + f₁ + f₂ + … + f_{n−1} + fₙ/2 ]

Reading: every interior sample gets full weight, the two ends get half. Interior points sit on the edge of two trapezoids and get counted once from each side.

The version that matters for real DAQ records — no uniform-spacing assumption at all:

∫ f dt ≈ Σ (t_{i+1} − t_i) · (f_i + f_{i+1}) / 2

Reading: each pair of consecutive timestamps makes its own trapezoid. Jittered clocks, dropped packets, mixed sample rates — the rule doesn't care.

The error, for a twice-differentiable integrand on a uniform grid:

E = −(b−a)·h²/12 · f″(ξ)        for some ξ in [a, b]

Reading: second-order accurate — halve the spacing, error drops by 4. The sign is doing work too: where the curve is concave down (f″ < 0 makes E positive as written, meaning the rule reads low), the chords cut the corner and you undercount. Systematically.

The running-sum form, for velocity-from-acceleration and every other cumulative integral:

v_k = v_{k−1} + (Δt/2) · (a_k + a_{k−1})

Reading: new total equals old total plus one trapezoid. This recursion is also exactly what the bilinear transform (s ← (2/T)·(z−1)/(z+1)) does to an integrator — trapezoidal integration and Tustin's method are the same move wearing different badges ^[7].

Where you meet it

Velocity off an accelerometer. Sled test, shaker table, drop tower: the DAQ hands you acceleration at 10 kHz and the customer wants velocity change. scipy.integrate.cumulative_trapezoid(a, t) — or MATLAB's cumtrapz — is the whole analysis. Delta-V for a pyroshock event, terminal velocity for a drop, integrated once more for displacement. This is the single most common numerical integral run in any test group.

Charge from a current log. Amp-hours through a battery under cycling, charge delivered by a pyro firing circuit, dose accumulated on a shunt record — all of them are ∫ i dt, and all of them get computed as trapezoids on the logged samples. Coulomb counting for state-of-charge is this integral running continuously.

Total impulse from a thrust trace. Thrust versus time off the stand's load cell, numpy.trapezoid(F, t), done — total impulse in one line, and delivered Isp after dividing by propellant weight. Simpson's rule buys more digits on the same file, but the trapezoid version survives dropped samples and uneven timestamps without ceremony, which is why it's the default reach.

Inside your digital controller. Discretize the I-term of a PID with Tustin's method and the accumulator update that ships in the firmware is the trapezoid recursion above ^[7]. If a control review asks "forward Euler, backward Euler, or trapezoidal?" — this is what the third option means, and it's the one that maps the whole stable left half-plane inside the unit circle.

How it works

The rule replaces the curve with chords, so it pays wherever the curve bends. Integrating sin(x) over [0, π] (true value 2) with 10 intervals gives 1.9835 — low by 1.6·10⁻², low because a chord under a concave-down arc always misses the bulge. At 20 intervals the error is 4.1·10⁻³: a factor of 4.00, the h² contract holding exactly. All numbers here checked by direct computation.

That contract hands you a free convergence check. Compute the integral twice — once with every sample (I_h), once with every other sample (I_2h):

error in I_h ≈ (I_h − I_2h) / 3

Since halving the spacing cuts the error 4-fold, the gap between the two answers is about 3 times the fine answer's remaining error. On the sin example the estimate says 4.12·10⁻³ against a true error of 4.11·10⁻³. When a review board asks how good your integrated number is, this is the two-line answer.

The rule's honest advantage over Simpson isn't accuracy — it's robustness. Simpson's weights demand uniform spacing and an even interval count; feed them timestamped field data with a dropout and you get a confidently wrong number. The trapezoid rule takes each interval on its own terms. It also degrades gracefully at kinks: on √x over [0, 1], where every derivative blows up at zero, it still converges (6.2·10⁻³ at 10 intervals, 2.0·10⁻⁴ at 100, 6.5·10⁻⁶ at 1,000 — about n^1.5), just slower than advertised.

Then there's the case where the humble rule embarrasses everything: smooth periodic functions integrated over a full period. The endpoint errors cancel around the loop — the Euler–Maclaurin formula makes this precise — and convergence goes exponential for analytic integrands ^[1]^[6]. Integrating e^(cos x) over [0, 2π] with 8 intervals lands within 1.3·10⁻⁶; 16 intervals hits machine precision. One cycle of a power waveform, one full antenna-pattern cut, THD from an integer number of periods: plain trapezoid is the right tool, and "upgrading" to Simpson makes it thousands of times worse.

The mistake people make isn't in the rule — it's in what integration does to error. The h² term is truncation error and it's usually negligible at DAQ sample rates. What kills accelerometer-to-velocity work is bias: a constant 0.01 m/s² offset integrates to a 0.6 m/s velocity error after one minute, growing linearly forever. Zero-mean noise is gentler but still accumulates as a random walk, growing like √t (verified by simulation: integrating unit white noise at 1 kHz for 100 s gives a velocity spread of about 0.32 m/s, matching σ·√(t·Δt)). Remove the pre-event mean, high-pass, or re-zero at known quiet points — and when the integrated trace drifts, blame the sensor offset, not the quadrature.

Tooling note: NumPy 2.0 renamed np.trapz to numpy.trapezoid ^[8]; SciPy has cumulative_trapezoid for running integrals; MATLAB keeps trapz and cumtrapz. All of them accept a vector of sample locations, so uneven spacing costs nothing.

History

The oldest known use of this rule is also one of the oldest pieces of applied mathematics on record. In 2016, Mathieu Ossendrijver of Humboldt University published his analysis of Babylonian cuneiform tablets, dated between 350 and 50 BCE, showing astronomers in Babylon computing Jupiter's displacement along the ecliptic as the area of a trapezoid drawn in an abstract time-velocity space — the planet's apparent speed drops roughly linearly over 60 days, and they multiplied time by average velocity exactly as the one-panel formula says ^[2]^[3]. Historians had credited that idea — distance as area under a speed graph — to the scholars of 14th-century Oxford and Paris. The tablets moved it back some seventeen centuries ^[2]^[3].

The modern rule arrived as the simplest member of a family. Isaac Newton originated the approach of integrating a polynomial through equally spaced samples, and Roger Cotes — the Cambridge prodigy who edited the second edition of Newton's Principia between 1709 and 1713 — extended it into the tables now called the Newton–Cotes formulas; the trapezoid rule is the two-point case ^[4]^[5]. Cotes died in 1716 at thirty-three, and his quadrature work reached print only in the posthumous Harmonia mensurarum of 1722 ^[5].

The last chapter is recent. In 2014, Nick Trefethen and André Weideman published a 70-page SIAM Review paper on the "exponentially convergent trapezoidal rule," documenting how the supposedly crude method achieves geometric convergence for analytic functions on periodic intervals and the real line — and quietly powers algorithms for inverse Laplace transforms, special functions, and matrix functions ^[6]. Twenty-four centuries in, the simplest rule in numerical analysis is still publishing.

Related tools

/tools/specific-impulse
/tools/battery-life
/tools/rms-peak
/tools/capacitor-energy

Sources

Written by HE in our own words from the cited sources — engineering judgment included, your stamp still required. All entries →