Editing Simpson's rule (section)

== 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 and is the [[#Composite Simpson's 1/3 rule|'''composite Simpson's 1/3 rule''']] evaluated for <math>n = 2</math>. Simpson's 1/3 rule is as follows:
<math display="block">
\begin{align}
  \int_a^b f(x)\, dx
  &\approx \frac{b - a}{6} \left[f(a) + 4f\left(\frac{a + b}{2}\right) + f(b)\right]\\
  &= \frac{1}{3} h\left[f(a) + 4f\left(a + h\right) + f(b)\right],
\end{align}
</math>
where <math>h = (b - a)/n</math> is the step size for <math>n = 2</math>.

The error in approximating an integral by Simpson's rule for <math>n = 2</math> is
<math display="block">-\frac{1}{90} h^5f^{(4)}(\xi) = -\frac{(b - a)^5}{2880} f^{(4)}(\xi),</math>
where <math>\xi</math> (the [[Xi (letter)#Lowercase|Greek letter xi]]) is some number between <math>a</math> and <math>b</math>.{{sfn|Atkinson|1989|loc=equation (5.1.15)}}{{sfn|Süli|Mayers|2003|loc=§7.2}}

The error is asymptotically proportional to <math>(b - a)^5</math>. However, the above derivations suggest an error proportional to <math>(b - a)^4</math>. Simpson's rule gains an extra order because the points at which the integrand is evaluated are distributed symmetrically in the interval <math>[a,\ b]</math>.

Since the error term is proportional to the fourth derivative of <math>f</math> at <math>\xi</math>, this shows that Simpson's rule provides exact results for any polynomial <math>f</math> 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 function|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 <math>f''</math> exists and is [[convex function|convex]] in the interval <math>(a,\ b)</math>, then
<math display="block">(b - a) f\left(\frac{a + b}{2}\right) + \frac{1}{3} \left(\frac{b - a}{2}\right)^3 f''\left(\frac{a + b}{2}\right)
\leq \int_a^b f(x)\, dx
\leq \frac{b - a}{6} \left[f(a) + 4f\left(\frac{a + b}{2}\right) + f(b)\right].
</math>

=== Derivations ===
==== Quadratic interpolation ====
One derivation replaces the integrand <math>f(x)</math> by the [[quadratic polynomial]] (i.e. parabola) <math>P(x)</math> that takes the same values as <math>f(x)</math> at the end points <math>a</math> and <math>b</math> and the midpoint <math>a+h</math>, where <math>h = (b - a)/2</math>. One can use [[Lagrange polynomial|Lagrange polynomial interpolation]] to find an expression for this polynomial,
<math display="block">P(x) = f(a) \frac{(x - a - h)(x - b)}{-h(a - b)} + f(a+h) \frac{(x - a)(x - b)}{h(a+h-b)} + f(b) \frac{(x - a)(x - a - h)}{(b - a)(b - a - h)}.</math>

Define:
<math display="block">
A_{1}(x) = f(a) \frac{(x - a - h)(x - b)}{-h(a - b)},
</math>
<math display="block">
A_{2}(x) = f(a+h) \frac{(x - a)(x - b)}{h(a+h-b)},
</math>
and
<math display="block">
A_{3}(x) = f(b) \frac{(x - a)(x - a - h)}{(b - a)(b - a - h)}.
</math>

By making the substitution <math>u = x - a - h</math>, one can calculate:

<math display="block">
\begin{align}
\int_a^b A_{1}(x)\, dx = &\int_{-h}^{h}A_{1}(u)\, du\\
 &= f(a)\left(\int_{-h}^{h}\frac{u(u-h)}{2h^2}\, du\right) \\
 &= \frac{f(a)}{2h^2}\left(\int_{-h}^{h}u^2 - hu\, du\right) \\
 &= \frac{f(a)}{2h^2}\left[\frac{1}{3}u^3 - \frac{1}{2}hu^2 \right]_{u=-h}^{u=h} \\
 &= \frac{f(a)}{2h^2}\left[\left(\frac{1}{3}h^3 - \frac{1}{2}h^3 \right) - \left(-\frac{1}{3}h^3 - \frac{1}{2}h^3 \right) \right] \\
 &= \frac{f(a)}{2h^2}\left(\frac{1}{3}h^3 - \frac{1}{2}h^3  + \frac{1}{3}h^3 + \frac{1}{2}h^3 \right)  \\
 &= \frac{f(a)}{2h^2}\left(\frac{2}{3}h^3\right) \\
 &= \frac{h}{3}f(a).
\end{align}
</math>

<math display="block">
\begin{align}
\int_a^b A_{2}(x)\, dx = &\int_{-h}^{h}A_{2}(u)\, du\\
 &= f(a+h)\left(\int_{-h}^{h}\frac{(u+h)(u+a+h-b)}{h(a+h-b)}\, du\right)\\
 &= f(a+h)\left(\int_{-h}^{h}\frac{(u+h)(u-h)}{-h^2}\, du\right) \\
 &= - \frac{f(a+h)}{h^2}\left(\int_{-h}^{h}u^2 - h^2\, du\right) \\
 &= - \frac{f(a+h)}{h^2}\left[\frac{1}{3}u^3 - h^2u \right]_{u=-h}^{u=h} \\
 &= - \frac{f(a+h)}{h^2}\left[\frac{1}{3}h^3 - h^3 - \left(- \frac{1}{3}h^3 + h^3 \right) \right]  \\
 &= - \frac{f(a+h)}{h^2}\left(-\frac{2}{3}h^3 + \frac{1}{3}h^3 - h^3 \right)  \\
 &= - \frac{f(a+h)}{h^2}\left(\frac{-4}{3}h^3\right) \\
 &= \frac{4h}{3}f(a + h).
\end{align}
</math>

<math display="block">
\begin{align}
\int_a^b A_{3}(x)\, dx = &\int_{-h}^{h}A_{3}(u)\, du\\
 &= f(b)\left(\int_{-h}^{h}\frac{u(u+h)}{(b-a)(b-a-h)}\, du\right)\\
 &= f(b)\left(\int_{-h}^{h}\frac{u(u+h)}{2h^2}\, du\right)\\
 &= \frac{f(b)}{2h^2}\left(\int_{-h}^{h}u^2 + hu\, du\right)\\
 &= \frac{f(b)}{2h^2}\left[\frac{1}{3}u^3 + \frac{1}{2}hu^2 \right]_{u=-h}^{u=h}\\
 &= \frac{f(b)}{2h^2}\left[\frac{1}{3}h^3 + \frac{1}{2}h^3 - \left(-\frac{1}{3}h^3 + \frac{1}{2}h^3 \right) \right] \\
 &= \frac{f(b)}{2h^2}\left(\frac{5}{6}h^3 - \frac{1}{6}h^3 \right) \\
 &= \frac{f(b)}{2h^2}\left(\frac{2}{3}h^3\right)\\
 &= \frac{h}{3}f(b).
\end{align}
</math>

Therefore,

<math display="block">
\begin{align}
\int_a^b f(x)\, dx \approx&\int_a^b P(x)\, dx\\
 &= \int_a^b A_{1}(x)\, dx + \int_a^b A_{2}(x)\, dx + \int_a^b A_{3}(x)\, dx\\
 &= \frac{h}{3}f(a) + \frac{4h}{3}f(a + h) + \frac{h}{3}f(b)\\
 &= \frac{1}{3}h\left[f(a) + 4f(a + h ) + f(b)\right].
\end{align}
</math>
Because of the <math>1/3</math> 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 [[rectangle rule|midpoint rule]]
<math display="block">M = (b - a)f\left(\frac{a + b}{2}\right)</math>
and the [[trapezoidal rule]]
<math display="block">T = \frac{1}{2} (b - a)\big(f(a) + f(b)\big).</math>

The errors in these approximations are
<math display="block">\frac{1}{24} (b - a)^3 f''(a) + O\big((b - a)^4\big)</math>
and
<math display="block">-\frac{1}{12} (b - a)^3 f''(a) + O\big((b - a)^4\big)</math>
respectively, where <math>O\big((b - a)^4\big)</math> denotes a term asymptotically proportional to <math>(b - a)^4</math>. The two <math>O\big((b - a)^4\big)</math> 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]]
<math>\frac{2M + T}{3}.</math>
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]]''
<math display="block">\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).</math>

The coefficients {{math|α}}, {{math|β}} and {{math|γ}} 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.)

=== Composite Simpson's 1/3 rule ===
If the interval of integration <math>[a, b]</math> is in some sense "small", then Simpson's rule with <math>n = 2</math> 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 <math>[a, b]</math>. 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 <math>[a, b]</math> into <math>n > 2</math> 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 <math>[a, b]</math> is split up into <math>n</math> subintervals, with <math>n</math> an even number. Then, the composite Simpson's rule is given by

Dividing the interval <math>[a, b]</math> into <math>n</math> subintervals of length <math>h = (b - a)/n</math> and introducing the points <math>x_i = a + ih</math> for <math>0 \leq i \leq n</math> (in particular, <math>x_0 = a</math> and <math>x_n = b</math>), we have
<math display="block">
\begin{align}
  \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{align}
</math>
This composite rule with <math>n = 2</math> corresponds with the regular Simpson's rule of the preceding section.

The error committed by the composite Simpson's rule is
<math display="block">-\frac{1}{180} h^4(b - a)f^{(4)}(\xi),</math>
where <math>\xi</math> is some number between <math>a</math> and <math>b</math>, and <math>h = (b - a)/n</math> is the "step length".{{sfn|Atkinson|1989|pp=257{{ndash}}258}}{{sfn|Süli|Mayers|2003|loc=§7.5}} The error is bounded (in absolute value) by
<math display="block">\frac{1}{180} h^4(b - a) \max_{\xi \in [a, b]} \left|f^{(4)}(\xi)\right|.</math>

This formulation splits the interval <math>[a, b]</math> 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]].