Mathematical foundations of Simpson Quadrature Rules.
When Thomas Simpson published his mathematical rules in 1743, he revolutionized numerical integration by replacing linear secants with parabolic arcs. By passing a quadratic polynomial through three consecutive equidistant points, Simpson's 1/3 rule achieves a remarkable mathematical bonus: because the odd cubic term in the Taylor expansion is symmetric about the midpoint, its integral vanishes identically, granting exactness for all cubic polynomials without requiring any additional function evaluations. The resulting composite rule achieves fourth-order O(h⁴) error convergence with a remarkably small error constant (1/180). In this laboratory, we implement both the 1/3 and 3/8 rules, analyze the Peano error kernel, and verify the O(h⁴) convergence rate under systematic mesh refinement.
Viva voce questions & analytical theory
How is Simpson's 1/3 rule derived?
Simpson's 1/3 rule approximates f(x) on each double sub-panel [x_{2i}, x_{2i+2}] by a quadratic Lagrange polynomial P_2(x) passing through (x_{2i}, f_{2i}), (x_{2i+1}, f_{2i+1}), and (x_{2i+2}, f_{2i+2}). Integrating this parabola analytically yields (h/3) [f_{2i} + 4 f_{2i+1} + f_{2i+2}].
Why is Simpson's 1/3 rule exact for cubic polynomials even though it fits a quadratic parabola?
By symmetry, the cubic term (x - x_{2i+1})³ is an odd function with respect to the midpoint x_{2i+1}. Its integral over the symmetric interval [-h, h] vanishes identically: ∫_{-h}^h t³ dt = 0. Thus, Simpson's rule unexpectedly gains an extra degree of algebraic precision (d = 3).
What is the leading error term for composite Simpson's 1/3 rule?
The composite error is given by E_S = - [(b - a) / 180] · h⁴ · f^{(4)}(ξ), where ξ ∈ [a, b]. The factor 1/180 is remarkably small compared to the trapezoidal error factor 1/12.
When is Simpson's 3/8 rule used instead of Simpson's 1/3 rule?
Simpson's 1/3 rule strictly requires an even number of panels N. If an integration grid has an odd number of panels, one typically applies Simpson's 1/3 rule over the first N - 3 panels and uses Simpson's 3/8 rule (which fits a cubic polynomial across 3 panels) over the remaining 3 panels.
How do the error coefficients of the 1/3 and 3/8 rules compare?
For Simpson's 1/3 rule on double panel 2h, error is -(h⁵/90) f^{(4)}. For Simpson's 3/8 rule on triple panel 3h, error is -(3h⁵/80) f^{(4)}. Per unit interval, the 1/3 rule is slightly more accurate and requires fewer function evaluations.
What is a Newton-Cotes formula?
A Newton-Cotes formula is a numerical quadrature scheme constructed by integrating a polynomial interpolant on equally spaced nodes. Closed formulas include the endpoints (trapezoidal, Simpson), while open formulas omit endpoints (midpoint rule).
Why are higher-order Newton-Cotes formulas (e.g. n ≥ 8 panels) rarely used in practice?
Due to Runge's phenomenon, polynomial interpolants on equally spaced nodes develop negative quadrature weights for degree n ≥ 8. Negative weights cause catastrophic subtractive cancellation of floating point numbers, destroying numerical stability.
Algorithm precautions and stability rules
- Verify that panel count N is strictly even before invoking Simpson's 1/3 rule; use the 3/8 rule on the final triplet if N is odd.
- Avoid recalculating function values at internal nodes; evaluate odd and even index nodes in separate accumulator loops.
- Ensure floating point precision remains 64-bit IEEE 754 to prevent loss of significance when h⁴ becomes small (e.g. h⁴ ≈ 10⁻⁸).
- Scale by h/3 only once at the end of the accumulation loop to prevent repeated multiplication roundoff errors.
- Verify that the integrand is at least four times continuously differentiable (C⁴) on [a, b] to guarantee fourth-order convergence.
- Do not use high-order Newton-Cotes rules (n ≥ 8) with negative weights; compound Simpson's rule over small panels instead.
Theoretical limitations
- Simpson's rule requires equidistant node spacing; if data points are non-uniformly distributed, spline quadrature or Gaussian quadrature must be used.
- For integrands with discontinuous derivatives (such as |x - 0.5|), convergence degrades from O(h⁴) to O(h) or O(h²).
- Simpson's rule is restricted to 1D integrals; extending to multi-dimensional domains suffers from the curse of dimensionality, motivating Monte Carlo methods.