Editing Simpson's rule (section)

== Composite Simpson's rule for irregularly spaced data ==
For some applications, the integration interval <math>I = [a, b]</math> needs to be divided into uneven intervals &ndash; perhaps due to uneven sampling of data, or missing or corrupted data points. Suppose we divide the interval <math>I</math> into an '''even number <math>N</math> of subintervals''' of widths <math>h_k</math>. Then the composite Simpson's rule is given by{{sfn|Shklov|1960}}
<math display="block">\int_a^b f(x)\, dx \approx \sum_{i = 0}^{N/2 - 1} \frac{h_{2i} + h_{2i + 1}}{6} \left[\left(2 - \frac{h_{2i + 1}}{h_{2i}}\right) f_{2i} + \frac{(h_{2i} + h_{2i + 1})^2}{h_{2i} h_{2i + 1}} f_{2i + 1} + \left(2 - \frac{h_{2i}}{h_{2i + 1}}\right) f_{2i + 2}\right],</math>
where
<math display="block">f_k = f\left(a + \sum_{i = 0}^{k - 1} h_i\right)</math>
are the function values at the <math>k</math>th sampling point on the interval <math>I</math>.

In case of an '''odd number <math>N</math> of subintervals''', the above formula is used up to the second to last interval, 
and the last interval is handled separately by adding the following to the result:{{sfn|Cartwright|2017|loc=Equation 8. The equation in Cartwright is calculating the first interval whereas the equations in the Wikipedia article are adjusting for the '''last''' integral. If the proper algebraic substitutions are made, the equation results in the values shown}}
<math display="block">\alpha f_N + \beta f_{N - 1} - \eta f_{N - 2},</math>
where
<math display="block">
\begin{align}
  \alpha &= \frac{2h_{N - 1}^2 + 3h_{N - 1} h_{N - 2}}{6(h_{N - 2} + h_{N - 1})},\\[1ex]
  \beta &= \frac{h_{N - 1}^2 + 3h_{N - 1} h_{N - 2}}{6h_{N - 2}},\\[1ex]
  \eta &= \frac{h_{N - 1}^3}{6 h_{N - 2}(h_{N - 2} + h_{N - 1})}.
\end{align}
</math>

{| role="presentation" class="wikitable mw-collapsible mw-collapsed"
|'''Example implementation in [[Python (programming language)|Python]]'''
|-
|<syntaxhighlight lang="python">
from collections.abc import Sequence

def simpson_nonuniform(x: Sequence[float], f: Sequence[float]) -> float:
    """
    Simpson rule for irregularly spaced data.

    :param x: Sampling points for the function values
    :param f: Function values at the sampling points

    :return: approximation for the integral

    See ``scipy.integrate.simpson`` and the underlying ``_basic_simpson``
    for a more performant implementation utilizing numpy's broadcast.
    """
    N = len(x) - 1
    h = [x[i + 1] - x[i] for i in range(0, N)]
    assert N > 0

    result = 0.0
    for i in range(1, N, 2):
        h0, h1 = h[i - 1], h[i]
        hph, hdh, hmh = h1 + h0, h1 / h0, h1 * h0
        result += (hph / 6) * (
            (2 - hdh) * f[i - 1] + (hph**2 / hmh) * f[i] + (2 - 1 / hdh) * f[i + 1]
        )

    if N % 2 == 1:
        h0, h1 = h[N - 2], h[N - 1]
        result += f[N]     * (2 * h1 ** 2 + 3 * h0 * h1) / (6 * (h0 + h1))
        result += f[N - 1] * (h1 ** 2 + 3 * h1 * h0)     / (6 * h0)
        result -= f[N - 2] * h1 ** 3                     / (6 * h0 * (h0 + h1))
    return result
</syntaxhighlight>
|}

{| role="presentation" class="wikitable mw-collapsible mw-collapsed"
|'''Example implementation in [[R (programming language)|R]]'''
|-
|<syntaxhighlight lang="r">
SimpsonInt <- function(fx, dx) {
  n <- length(dx)
  h <- diff(dx)
  stopifnot(exprs = {
    length(fx) == n
    all(h >= 0)
  })
  res <- 0
  for (i in seq(1L, n - 2L, 2L)) {
    hph <- h[i] + h[i + 1L]
    hdh <- h[i + 1L] / h[i]
    res <- res + hph / 6 *
      ((2 - hdh) * fx[i] +
         hph ^ 2 / (h[i] * h[i + 1L]) * fx[i + 1L] +
         (2 - 1 / hdh) * fx[i + 2L])
  }

  if (n %% 2 == 0) {
    hph <- h[n - 1L] + h[n - 2L]
    threehth <- 3 * h[n - 1L] * h[n - 2L]
    sixh2 <- 6 * h[n - 2L]
    h1sq <- h[n - 1L] ^ 2
    res <- res +
      (2 * h1sq + threehth) / (6 * hph) * fx[n] +
      (h1sq + threehth) / sixh2 * fx[n - 1L] -
      (h1sq * h[n - 1L]) / (sixh2 * hph) * fx[n - 2L]
  }
  res
}
</syntaxhighlight>
|}