$ \newcommand\vu{\mathbf{u}} \newcommand\vv{\mathbf{v}} \newcommand\vw{\mathbf{w}} \newcommand\vi{\mathbf{i}} \newcommand\vj{\mathbf{j}} \newcommand\vk{\mathbf{k}} \newcommand\vo{\mathbf{o}} \newcommand\vr{\mathbf{r}} \newcommand\vt{\mathbf{t}} \newcommand\vx{\mathbf{x}} \newcommand\vy{\mathbf{y}} \newcommand\vz{\mathbf{z}} \newcommand\va{\mathbf{a}} \newcommand\vb{\mathbf{b}} \newcommand\vc{\mathbf{c}} \newcommand\ve{\mathbf{e}} \newcommand{\vE}{\mathbf{E}} \newcommand{\vS}{\mathbf{S}} \newcommand{\vk}{\mathbf{k}} \newcommand{\vq}{\mathbf{q}} \newcommand{\vzero}{\mathbf{0}} \newcommand{\vomega}{\mathbf{ω}} \newcommand{\mI}{\mathbf{I}} \newcommand{\mM}{\mathbf{M}} \newcommand{\mR}{\mathbf{R}} \newcommand{\mP}{\mathbf{P}} \newcommand{\mK}{\mathbf{K}} \newcommand{\mW}{\mathbf{W}} \DeclareMathOperator{\Tr}{Tr} \newcommand{\abs}[1]{\left\lvert {#1} \right\rvert} \newcommand{\norm}[1]{\left\lVert {#1} \right\rVert} \newcommand{\ensemble}[1]{\left\langle {#1} \right\rangle} \newcommand{\der}[2]{\frac{\partial {#1}}{\partial {#2}}} \newcommand{\Der}[2]{\frac{\delta {#1}}{\delta {#2}}} \newcommand{\lap}{\nabla^2} $

Simpson's Rule and Other Fourth-Order Quadrature Formulas

Research Numerical Algorithm Integration Formula

In this post we will list some popular closed-, open-, and semi-open form of numerical integration formulas for Simpson’s rule and other 4th order algorithms.

1. Closed formulas

By dividing the closed interval \([a, b]\) into \(N-1\) uniform subintervals with length of \(h=(b-a)/(N-1)\), the extended Simpson’s rule is given by (1)

\[\begin{aligned} \label{simpson} \int_a^b f(x)dx = & \frac{h}{3} \left( f_0 + 4f_1 + 2f_2 + 4f_3 + \cdots + 2f_{N-3} + 4f_{N-2} + f_{N-1} \right) \\ & -\frac{b-a}{180}f^{(4)}(\xi)h^4 \end{aligned}\]

where \(f_i = f(x_i)\) for \(x_i = ih\), and \(\xi\) is in \((a, b)\). The number of points, \(N\), must be odd, meaning that there are even number of subintervals. This requirement is a consequence of the fact that the coefficients alternate in a specific pattern, namely there are \(n\) pairs of \((4, 2)\) coefficient terms plus an extra \(4\)-coefficient term and two end-point terms. Therefore, there are \(2n+3\) terms in total, while \(2n+3\) is always odd.

To relax the constraint of odd number of points, which is important for situations when it is not convenient to freely choose the number of points such as in experiments or numerical simultaions, we can use

\[\begin{aligned} \label{three-term} \int_a^b f(x)dx = h \bigl( &\frac{9}{24}f_0 + \frac{28}{24}f_1 + \frac{23}{24}f_2 + \\ & f_3 + f_4 + \cdots + f_{N-5} + f_{N-4} + \\ & \frac{23}{24}f_{N-3} + \frac{28}{24}f_{N-2} + \frac{9}{24}f_{N-1} \bigr) \\ & +O(h^4) \end{aligned}\]

Note that the coefficients are no longer alternate except the first and last three terms. See Arif Zaman’s Note (2) for a derivation of the exact form of the error. Other non-alternate forms of fourth-order quadrature formulas are possible. An example from Wolfram Mathworld (3)

\[\begin{aligned} \label{four-term} \int_a^b f(x)dx = h \bigl( &\frac{17}{48}f_0 + \frac{59}{48}f_1 + \frac{43}{48}f_2 + \frac{49}{48}f_3 + \\ & f_4 + f_5 + \cdots + f_{N-6} + f_{N-5} + \\ & \frac{49}{48}f_{N-4} + \frac{43}{48}f_{N-3} + \frac{59}{48}f_{N-2} + \frac{17}{48}f_{N-1} \bigr) \\ & +O(h^4) \end{aligned}\]

2. Open and semi-open formulas

Sometimes open and semi-open formulas are required. Below is a fourth-order open formula for \((0, N-1)\)

\[\begin{aligned} \label{open-three-term} \int_a^b f(x)dx = h \bigl( &\frac{55}{24}f_1 - \frac{4}{24}f_2 + \frac{33}{24}f_3 + \\ & f_4 + f_5 + \cdots + f_{N-6} + f_{N-5} + \\ & \frac{33}{24}f_{N-4} - \frac{4}{24}f_{N-3} + \frac{55}{24}f_{N-2} \bigr) \\ & +O(h^4) \end{aligned}\]

And here is a fourth-order semi-open formula for \((0, N-1]\)

\[\begin{aligned} \label{semiopen-three-term} \int_a^b f(x)dx = h \bigl( &\frac{55}{24}f_1 - \frac{4}{24}f_2 + \frac{33}{24}f_3 + \\ & f_4 + f_5 + \cdots + f_{N-5} + f_{N-4} + \\ & \frac{23}{24}f_{N-3} + \frac{28}{24}f_{N-2} + \frac{9}{24}f_{N-1} \bigr) \\ & +O(h^4) \end{aligned}\]

Note that above formulas share the same discrete scheme as the closed formulas.

Reference

  1. Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numercial Recipes: The Art of Scientific Computing; 3rd ed.; Cambridge University Press: Cambridge, 2007.
  2. Arif Zaman’s Note
  3. An example from Wolfram Mathworld