In numerical integration, Simpson's rules are several approximations for definite integrals, named after Thomas Simpson (1710–1761).
The most basic of these rules, called Simpson's 1/3 rule, or just Simpson's rule, reads
In German and some other languages, it is named after Johannes Kepler, who derived it in 1615 after seeing it used for wine barrels (barrel rule, ). The approximate equality in the rule becomes exact if is a polynomial up to and including 3rd degree.
If the 1/3 rule is applied to n equal subdivisions of the integration range [a, b], one obtains the composite Simpson's 1/3 rule. Points inside the integration range are given alternating weights 4/3 and 2/3.
Simpson's 3/8 rule, also called Simpson's second rule, requires one more function evaluation inside the integration range and gives lower error bounds, but does not improve on order of the error.
If the 3/8 rule is applied to n equal subdivisions of the integration range [a, b], one obtains the composite Simpson's 3/8 rule.
Simpson's 1/3 and 3/8 rules are two special cases of closed Newton–Cotes formulas.
In naval architecture and ship stability estimation, there also exists Simpson's third rule, which has no special importance in general numerical analysis, see Simpson's rules (ship stability).
Simpson's 1/3 rule
Simpson's 1/3 rule, also simply called Simpson's rule, is a method for numerical integration proposed by Thomas Simpson. It is based upon a quadratic interpolation. Simpson's 1/3 rule is as follows:
where is the step size.
The error in approximating an integral by Simpson's rule for is
where (the Greek letter xi) is some number between and .
The error is asymptotically proportional to . However, the above derivations suggest an error proportional to . Simpson's rule gains an extra order because the points at which the integrand is evaluated are distributed symmetrically in the interval .
Since the error term is proportional to the fourth derivative of at , this shows that Simpson's rule provides exact results for any polynomial of degree three or less, since the fourth derivative of such a polynomial is zero at all points. Another way to see this result is to note that any interpolating cubic polynomial can be expressed as the sum of the unique interpolating quadratic polynomial plus an arbitrarily scaled cubic polynomial that vanishes at all three points in the interval, and the integral of this second term vanishes because it is odd within the interval.
If the second derivative exists and is convex in the interval , then
Derivations
Quadratic interpolation
One derivation replaces the integrand by the quadratic polynomial (i.e. parabola) that takes the same values as at the end points and and the midpoint . One can use Lagrange polynomial interpolation to find an expression for this polynomial,
Using integration by substitution, one can show that
Introducing the step size , this is also commonly written as
Because of the factor, Simpson's rule is also referred to as "Simpson's 1/3 rule" (see below for generalization).
Averaging the midpoint and the trapezoidal rules
Another derivation constructs Simpson's rule from two simpler approximations: the midpoint rule
and the trapezoidal rule
The errors in these approximations are
and
respectively, where denotes a term asymptotically proportional to . The two terms are not equal; see Big O notation for more details. It follows from the above formulas for the errors of the midpoint and trapezoidal rule that the leading error term vanishes if we take the weighted average
This weighted average is exactly Simpson's rule.
Using another approximation (for example, the trapezoidal rule with twice as many points), it is possible to take a suitable weighted average and eliminate another error term. This is Romberg's method.
Undetermined coefficients
The third derivation starts from the ansatz
The coefficients , and can be fixed by requiring that this approximation be exact for all quadratic polynomials. This yields Simpson's rule. (This derivation is essentially a less rigorous version of the quadratic interpolation derivation, where one saves significant calculation effort by guessing the correct functional form.)