Simpson's Rule

Replace f by a parabola

The trapezoidal rule replaces the integrand f(x) by the straight line through (a, f(a)) and (b, f(b)) and integrates the resulting trapezoid. Simpson's rule replaces f instead by the parabola p(x) through three equally spaced points — the two endpoints and the midpoint c = (a + b)/2 — and integrates that parabola exactly.

Let h = (b − a)/2 be the half-width. The Lagrange polynomial of degree 2 through (a, f(a)), (c, f(c)), (b, f(b)) is:

p(x) = f(a) · (x − c)(x − b) / (2h²)
     − f(c) · (x − a)(x − b) / h²
     + f(b) · (x − a)(x − c) / (2h²)

Integrate p from a to b. The algebra collapses with the substitution u = x − c:

∫_a^b p(x) dx = (h / 3) · (f(a) + 4 f(c) + f(b))

That is Simpson's rule. The coefficients (1, 4, 1) arise from integrating the three Lagrange basis polynomials over [a, b]; the prefactor h/3 absorbs the interval length.

Why Simpson is exact for cubics (not just parabolas)

Simpson fits a degree-2 polynomial — you'd expect exactness for polynomials of degree ≤ 2. Surprisingly, the rule is exact for cubics as well, giving a degree of precision of 3. The reason is symmetry. Any cubic f(x) on [a, b] decomposes as p(x) + α (x − c)³ where p is the parabola through the three Simpson points and (x − c)³ is the odd-cubic correction. The correction integrates to exactly 0 on the symmetric interval [a, b] (odd about c), and the Simpson sum hits (x − c)³ at x = a, c, b, where the cubic part is (−h)³, 0, h³ respectively. The 4f(c) weight kills the 0, and the symmetric (1, 1) weights on f(a), f(b) cancel −h³ and h³. So Simpson sees the cubic correction as 0 — exactly as the true integral does.

The first thing Simpson misses is the quartic term. Expand f around c using Taylor:

f(x) = f(c) + f'(c)(x−c) + (f''(c)/2)(x−c)² + (f'''(c)/6)(x−c)³ + (f⁽⁴⁾(c)/24)(x−c)⁴ + ...

The 0th, 1st, 2nd, 3rd terms are integrated exactly by Simpson. The 4th term contributes the leading error:

∫_a^b (x−c)⁴ dx − Simpson_estimate[(x−c)⁴]
   = (2h⁵/5) − (h/3)(h⁴ + 0 + h⁴) = (2h⁵/5 − 2h⁵/3) = −4h⁵/15

So the per-panel error is −4h⁵/15 · f⁽⁴⁾(c)/24 = −(2h)⁵/2880 · f⁽⁴⁾(c). Equivalently in terms of the panel width b − a = 2h:

E = −(b − a)⁵ / 2880 · f⁽⁴⁾(ξ)    for some ξ ∈ (a, b)

This is the famous Simpson's-rule error formula. It demands f has a continuous fourth derivative — and the error grows with the fourth derivative's magnitude.

Worked example — ∫_0^1 e^x dx

Exact value e − 1 ≈ 1.7182818. Apply Simpson with a = 0, c = 0.5, b = 1, so h = 0.5:

f(0) = 1
f(0.5) = e^0.5 ≈ 1.6487213
f(1) = e ≈ 2.7182818

Simpson = (0.5 / 3) · (1 + 4 · 1.6487213 + 2.7182818)
        = (0.5 / 3) · (1 + 6.5948851 + 2.7182818)
        = (0.5 / 3) · 10.3131669
        = 1.7188611

Error ≈ 1.7188611 − 1.7182818 = +0.0005793, about 3 × 10⁻⁴. Predicted bound: |E| ≤ (1⁵/2880) · max|f⁽⁴⁾| = (1/2880) · e ≈ 0.000944. Our actual error fits well inside that bound.

Compare the trapezoidal rule with the same 3 evaluations chained as two panels of width 0.5:

Trap = 0.5 · (f(0)/2 + f(0.5) + f(1)/2)
     = 0.5 · (0.5 + 1.6487213 + 1.3591409)
     = 0.5 · 3.5078622 = 1.7539311
Error = +0.0356, about 60× larger

Same 3 function evaluations, Simpson is 60× more accurate on this integral. That's the practical payoff of jumping from second- to fourth-order accuracy.

Composite Simpson — chain rule for many panels

For higher accuracy, divide [a, b] into n equally spaced subintervals (n must be even) and apply Simpson to each pair of subintervals:

h = (b − a) / n
x_k = a + k h,   k = 0, 1, ..., n

Composite Simpson:
∫_a^b f ≈ (h/3) · [f₀ + 4(f₁ + f₃ + f₅ + ...) + 2(f₂ + f₄ + f₆ + ...) + f_n]

Global error: E = −(b − a) · h⁴ / 180 · f⁽⁴⁾(ξ)

The h⁴ scaling is the fourth-order accuracy. To shrink the error by 16, halve h (i.e., double the number of panels). For most smooth integrands, 100 panels gives ~10⁻⁸ accuracy.

JavaScript — Simpson's rule, single and composite

// Single Simpson panel on [a, b]
function simpsonPanel(f, a, b) {
  const c = (a + b) / 2;
  const h = (b - a) / 2;
  return (h / 3) * (f(a) + 4 * f(c) + f(b));
}

// Composite Simpson with n (even) panels
function simpsonComposite(f, a, b, n = 100) {
  if (n & 1) n++;                       // require even
  const h = (b - a) / n;
  let sum = f(a) + f(b);
  for (let k = 1; k < n; k++) {
    sum += (k & 1 ? 4 : 2) * f(a + k * h);
  }
  return sum * h / 3;
}

// Test: ∫_0^1 e^x dx = e - 1 ≈ 1.7182818
const f = x => Math.exp(x);
console.log(simpsonPanel(f, 0, 1));               // 1.7188611 (error 5.8e-4)
console.log(simpsonComposite(f, 0, 1, 4));        // 1.71832 (error 5e-5)
console.log(simpsonComposite(f, 0, 1, 16));       // 1.71828183 (error 4e-8)

// Adaptive Simpson — refine where error is large
function adaptiveSimpson(f, a, b, tol = 1e-8) {
  function go(a, b, fa, fc, fb, S, depth) {
    const c  = (a + b) / 2;
    const lc = (a + c) / 2;
    const rc = (c + b) / 2;
    const flc = f(lc), frc = f(rc);
    const Sl = ((c - a) / 6) * (fa + 4 * flc + fc);
    const Sr = ((b - c) / 6) * (fc + 4 * frc + fb);
    const S2 = Sl + Sr;
    if (depth <= 0 || Math.abs(S2 - S) < 15 * tol) {
      return S2 + (S2 - S) / 15;        // Richardson extrapolation
    }
    return go(a, c, fa, flc, fc, Sl, depth - 1)
         + go(c, b, fc, frc, fb, Sr, depth - 1);
  }
  const c = (a + b) / 2;
  const fa = f(a), fc = f(c), fb = f(b);
  const S = ((b - a) / 6) * (fa + 4 * fc + fb);
  return go(a, b, fa, fc, fb, S, 30);
}

console.log(adaptiveSimpson(f, 0, 1));            // ~1.7182818

Where Simpson sits in the Newton–Cotes family

Rule	Points	Per-panel formula	Error per panel	Degree of precision
Rectangle (midpoint)	1	(b−a) · f((a+b)/2)	(b−a)³/24 · f''(ξ)	1
Trapezoidal	2	(b−a)/2 · (f(a) + f(b))	−(b−a)³/12 · f''(ξ)	1
Simpson's 1/3	3	(h/3)(f(a) + 4f(c) + f(b))	−(b−a)⁵/2880 · f⁽⁴⁾(ξ)	3
Simpson's 3/8	4	(3h/8)(f₀ + 3f₁ + 3f₂ + f₃)	−(b−a)⁵/6480 · f⁽⁴⁾(ξ)	3
Boole's rule	5	(2h/45)(7f₀ + 32f₁ + 12f₂ + 32f₃ + 7f₄)	−(b−a)⁷/1935360 · f⁽⁶⁾(ξ)	5
Higher Newton–Cotes	6+	tabulated weights	weights become negative for n ≥ 8	n − 1 (even n) or n (odd n)

Past n = 7, Newton–Cotes weights start to include negative values, which creates numerical instability for high-amplitude integrands. Simpson's 1/3 hits the sweet spot of high accuracy with positive, intuitive weights — which is why it became the standard textbook rule.

Simpson vs Gauss-Legendre with the same evaluation count

	Simpson's 1/3 (3 evaluations)	Gauss-Legendre n = 3 (3 evaluations)
Node locations	Equally spaced (a, c, b)	Roots of P_3 (0, ±√(3/5))
Weights	(h/3, 4h/3, h/3)	(5/9, 8/9, 5/9) × (b−a)/2
Degree of precision	3	5
Convergence (h refinement)	O(h⁴) globally	O(h⁶) globally if h-refinement applied
Works with tabular data	Yes (equally spaced)	No (irregular nodes)
Best for	Equally spaced data; simple production code	Smooth integrand, free node choice

For a smooth integrand with full control over evaluation points, Gauss is more accurate per evaluation. For tabular data on a fixed grid, Simpson is the natural rule.

Where Simpson's rule actually gets used

• Engineering production code. Composite Simpson is the default in many CAD/CAM and finite element preprocessors when integrating geometric properties (areas, centroids, moments).
• Probability and statistics. Computing tail probabilities of distributions with no closed-form CDF — normal, gamma, etc.
• Scientific instrument software. Numerical integration of measured curves — spectra, NMR signals, calorimetry — on equally spaced sample points.
• Computer graphics. Filtering operations, motion blur integrals, lens-camera convolutions — often evaluated by Simpson over the sample grid.
• Adaptive quadrature. The recursive adaptive Simpson algorithm subdivides only where the error estimate (Richardson-extrapolated) exceeds tolerance — the basis of SciPy's quad in early days and many textbook codes.
• Pricing path-dependent options. The "Asian option" payoff integrals are smooth and benefit from composite Simpson over many time steps.
• Spacecraft trajectory analysis. Integrating perturbation forces along an orbit — equally spaced time samples → Simpson.

Common mistakes

• Using composite Simpson with odd n. Composite Simpson pairs subintervals two at a time, so the total panel count must be even. With odd n, either skip a sub-panel (biased) or use Simpson's 3/8 on the last 3 sub-panels and 1/3 on the rest (hybrid).
• Applying Simpson when f⁽⁴⁾ is unbounded. The error formula relies on a finite f⁽⁴⁾. For integrands with singularities (1/√x near 0, sqrt-singular at endpoints), Simpson's stated convergence rate disappears and the actual error decays slowly.
• Believing higher Newton–Cotes is always better. Newton–Cotes of order ≥ 8 has negative weights and can amplify roundoff dramatically. Use composite low-order rules (Simpson) or switch to Gauss for very high accuracy.
• Confusing Simpson's 1/3 with Simpson's 3/8. Both due to Thomas Simpson. The 1/3 rule uses 3 points; the 3/8 rule uses 4 points and is exact through degree 3 (same degree of precision, slightly worse error constant). 1/3 is the default; 3/8 is used to handle remainders when n is not a multiple of 2.
• Trying Simpson on oscillatory or rapidly varying integrands. ∫₀^N cos(100x) dx — Simpson's polynomial fit is a terrible model. Use Filon-type rules or specialised oscillatory quadrature instead.
• Floating-point error accumulation. Summing many large terms in composite Simpson can lose precision. Use Kahan summation or Neumaier compensation when n is huge.