Derive Simpson's 1/3 rule

Simpson's 1/3 rule replaces f(x) on two sub-intervals by a parabola (2nd-degree interpolating polynomial) through three equally spaced points and integrates it exactly.

Derivation

Take points x₀, x₁ = x₀ + h, x₂ = x₀ + 2h with values y₀, y₁, y₂. Put p = (x − x₀)/h, so x = x₀ + ph and dx = h dp. Newton's forward interpolation up to second differences gives

f(x) ≈ y₀ + p·Δy₀ + p(p−1)/2 · Δ²y₀

Integrate from x₀ to x₂, i.e. p = 0 to 2:

∫ f dx = h ∫₀² [ y₀ + pΔy₀ + (p²−p)/2 · Δ²y₀ ] dp
       = h [ y₀·p + Δy₀·p²/2 + Δ²y₀·(p³/6 − p²/4) ]₀²
       = h [ 2y₀ + 2Δy₀ + Δ²y₀·(8/6 − 1) ]
       = h [ 2y₀ + 2Δy₀ + (1/3)Δ²y₀ ]

Use Δy₀ = y₁ − y₀ and Δ²y₀ = y₂ − 2y₁ + y₀:

∫ f dx = h [ 2y₀ + 2(y₁ − y₀) + (1/3)(y₂ − 2y₁ + y₀) ]
       = (h/3) [ 6y₀ + 6y₁ − 6y₀ + y₂ − 2y₁ + y₀ ]
       = (h/3) [ y₀ + 4y₁ + y₂ ]

This is Simpson's 1/3 rule for one pair of sub-intervals. Its error is −(h⁵/90) f⁗(ξ).

Composite rule

For n (even) sub-intervals, h = (b − a)/n, adding the n/2 pairs: the interior odd-indexed points appear once with weight 4 and the interior even-indexed points appear twice (end of one pair, start of next) with weight 2:

∫ₐᵇ f dx ≈ (h/3) [ y₀ + yₙ + 4(y₁ + y₃ + … + yₙ₋₁) + 2(y₂ + y₄ + … + yₙ₋₂) ]

Total error ≈ −(b−a)h⁴/180 · f⁗(ξ), so the rule is exact for polynomials up to degree 3.