Simpson's rule

One derivation replaces the integrand ${\displaystyle f(x)}$ by the quadratic polynomial (i.e. parabola) ${\displaystyle P(x)}$ that takes the same values as ${\displaystyle f(x)}$ at the end points ${\displaystyle a}$ and ${\displaystyle b}$ and the midpoint ${\displaystyle m=(a+b)/2}$. One can use Lagrange polynomial interpolation to find an expression for this polynomial,

$${\displaystyle P(x)=f(a){\frac {(x-m)(x-b)}{(a-m)(a-b)}}+f(m){\frac {(x-a)(x-b)}{(m-a)(m-b)}}+f(b){\frac {(x-a)(x-m)}{(b-a)(b-m)}}.}$$

Using integration by substitution, one can show that[3][2]

$${\displaystyle \int _{a}^{b}P(x)\,dx={\frac {b-a}{6}}\left[f(a)+4f\left({\frac {a+b}{2}}\right)+f(b)\right].}$$

Introducing the step size ${\displaystyle h=(b-a)/2}$, this is also commonly written as

$${\displaystyle \int _{a}^{b}P(x)\,dx={\frac {1}{3}}h\left[f(a)+4f\left({\frac {a+b}{2}}\right)+f(b)\right].}$$

Because of the ${\displaystyle 1/3}$ factor, Simpson's rule is also referred to as "Simpson's 1/3 rule" (see below for generalization).

Another derivation constructs Simpson's rule from two simpler approximations: the midpoint rule

$${\displaystyle M=(b-a)f\left({\frac {a+b}{2}}\right)}$$

and the trapezoidal rule

$${\displaystyle T={\frac {1}{2}}(b-a){\big (}f(a)+f(b){\big )}.}$$

The errors in these approximations are

$${\displaystyle {\frac {1}{24}}(b-a)^{3}f''(a)+O{\big (}(b-a)^{4}{\big )}}$$

and

$${\displaystyle -{\frac {1}{12}}(b-a)^{3}f''(a)+O{\big (}(b-a)^{4}{\big )}}$$

respectively, where ${\displaystyle O{\big (}(b-a)^{4}{\big )}}$ denotes a term asymptotically proportional to ${\displaystyle (b-a)^{4}}$. The two ${\displaystyle O{\big (}(b-a)^{4}{\big )}}$ 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 ${\displaystyle {\frac {2M+T}{3}}.}$ 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.

The third derivation starts from the ansatz

$${\displaystyle {\frac {1}{b-a}}\int _{a}^{b}f(x)\,dx\approx \alpha f(a)+\beta f\left({\frac {a+b}{2}}\right)+\gamma f(b).}$$

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.)

If the interval of integration ${\displaystyle [a,b]}$ is in some sense "small", then Simpson's rule with ${\displaystyle n=2}$ subintervals will provide an adequate approximation to the exact integral. By "small" we mean that the function being integrated is relatively smooth over the interval ${\displaystyle [a,b]}$. For such a function, a smooth quadratic interpolant like the one used in Simpson's rule will give good results.

However, it is often the case that the function we are trying to integrate is not smooth over the interval. Typically, this means that either the function is highly oscillatory or lacks derivatives at certain points. In these cases, Simpson's rule may give very poor results. One common way of handling this problem is by breaking up the interval ${\displaystyle [a,b]}$ into ${\displaystyle n>2}$ small subintervals. Simpson's rule is then applied to each subinterval, with the results being summed to produce an approximation for the integral over the entire interval. This sort of approach is termed the composite Simpson's 1/3 rule, or just composite Simpson's rule.

Suppose that the interval ${\displaystyle [a,b]}$ is split up into ${\displaystyle n}$ subintervals, with ${\displaystyle n}$ an even number. Then, the composite Simpson's rule is given by

Dividing the interval ${\displaystyle [a,b]}$ into ${\displaystyle n}$ subintervals of length ${\displaystyle h=(b-a)/n}$ and introducing the points ${\displaystyle x_{i}=a+ih}$ for ${\displaystyle 0\leq i\leq n}$ (in particular, ${\displaystyle x_{0}=a}$ and ${\displaystyle x_{n}=b}$), we have

$${\displaystyle {\begin{aligned}\int _{a}^{b}f(x)\,dx&\approx {\frac {1}{3}}h\sum _{i=1}^{n/2}{\big [}f(x_{2i-2})+4f(x_{2i-1})+f(x_{2i}){\big ]}\\&={\frac {1}{3}}h{\big [}f(x_{0})+4f(x_{1})+2f(x_{2})+4f(x_{3})+2f(x_{4})+\dots +2f(x_{n-2})+4f(x_{n-1})+f(x_{n}){\big ]}\\&={\frac {1}{3}}h\left[f(x_{0})+4\sum _{i=1}^{n/2}f(x_{2i-1})+2\sum _{i=1}^{n/2-1}f(x_{2i})+f(x_{n})\right].\end{aligned}}}$$

This composite rule with ${\displaystyle n=2}$ corresponds with the regular Simpson's rule of the preceding section.

The error committed by the composite Simpson's rule is

$${\displaystyle -{\frac {1}{180}}h^{4}(b-a)f^{(4)}(\xi ),}$$

where ${\displaystyle \xi }$ is some number between ${\displaystyle a}$ and ${\displaystyle b}$, and ${\displaystyle h=(b-a)/n}$ is the "step length".[4][5] The error is bounded (in absolute value) by

$${\displaystyle {\frac {1}{180}}h^{4}(b-a)\max _{\xi \in [a,b]}\left|f^{(4)}(\xi )\right|.}$$

This formulation splits the interval ${\displaystyle [a,b]}$ in subintervals of equal length. In practice, it is often advantageous to use subintervals of different lengths and concentrate the efforts on the places where the integrand is less well-behaved. This leads to the adaptive Simpson's method.

Dividing the interval ${\displaystyle [a,b]}$ into ${\displaystyle n}$ subintervals of length ${\displaystyle h=(b-a)/n}$ and introducing the points ${\displaystyle x_{i}=a+ih}$ for ${\displaystyle 0\leq i\leq n}$ (in particular, ${\displaystyle x_{0}=a}$ and ${\displaystyle x_{n}=b}$), we have

$${\displaystyle {\begin{aligned}\int _{a}^{b}f(x)\,dx&\approx {\frac {3}{8}}h\sum _{i=1}^{n/3}{\big [}f(x_{3i-3})+3f(x_{3i-2})+3f(x_{3i-1})+f(x_{3i}){\big ]}\\&={\frac {3}{8}}h{\big [}f(x_{0})+3f(x_{1})+3f(x_{2})+2f(x_{3})+3f(x_{4})+3f(x_{5})+2f(x_{6})+\dots +2f(x_{n-3})+3f(x_{n-2})+3f(x_{n-1})+f(x_{n}){\big ]}\\&={\frac {3}{8}}h\left[f(x_{0})+3\sum _{i=1,\ 3\nmid i}^{n-1}f(x_{i})+2\sum _{i=1}^{n/3-1}f(x_{3i})+f(x_{n})\right].\end{aligned}}}$$

While the remainder for the rule is shown as[6] ${\displaystyle -{\frac {1}{80}}h^{4}(b-a)f^{(4)}(\xi ),}$ we can only use this if ${\displaystyle n}$ is a multiple of three. The 1/3 rule can be used for the remaining subintervals without changing the order of the error term (conversely, the 3/8 rule can be used with a composite 1/3 rule for odd-numbered subintervals).

In the task of estimation of full area of narrow peak-like functions, Simpson's rules are much less efficient than trapezoidal rule. Namely, composite Simpson's 1/3 rule requires 1.8 times more points to achieve the same accuracy as trapezoidal rule.[8] Composite Simpson's 3/8 rule is even less accurate. Integration by Simpson's 1/3 rule can be represented as a weighted average with 2/3 of the value coming from integration by the trapezoidal rule with step h and 1/3 of the value coming from integration by the rectangle rule with step 2h. The accuracy is governed by the second (2h step) term. Averaging of Simpson's 1/3 rule composite sums with properly shifted frames produces the following rules:

$${\displaystyle \int _{a}^{b}f(x)\,dx\approx {\frac {1}{24}}h\left[-f(x_{-1})+12f(x_{0})+25f(x_{1})+24\sum _{i=2}^{n-2}f(x_{i})+25f(x_{n-1})+12f(x_{n})-f(x_{n+1})\right],}$$

where two points outside of the integrated region are exploited, and

$${\displaystyle \int _{a}^{b}f(x)\,dx\approx {\frac {1}{24}}h\left[9f(x_{0})+28f(x_{1})+23f(x_{2})+24\sum _{i=3}^{n-3}f(x_{i})+23f(x_{n-2})+28f(x_{n-1})+9f(x_{n})\right],}$$

where only points within integration region are used. Application of the second rule to the region of 3 points generates 1/3 Simpon's rule, 4 points - 3/8 rule.

These rules are very much similar to the alternative extended Simpson's rule. The coefficients within the major part of the region being integrated are one with non-unit coefficients only at the edges. These two rules can be associated with Euler–MacLaurin formula with the first derivative term and named First order Euler–MacLaurin integration rules.[8] The two rules presented above differ only in the way how the first derivative at the region end is calculated. The first derivative term in the Euler–MacLaurin integration rules accounts for integral of the second derivative, which equals the difference of the first derivatives at the edges of the integration region. It is possible to generate higher order Euler–Maclaurin rules by adding a difference of 3rd, 5th, and so on derivatives with coefficients, as defined by Euler–MacLaurin formula.